Deep learning time-series modeling for assessing land subsidence under reduced groundwater use

Land subsidence in Choshui delta

The Choshui delta, situated in central Taiwan, is one of the largest alluvial fans in the region, formed through the continuous deposition of sediments on the western side of a hilly terrain. Soil compression within the delta is a natural geological process, occurring due to the weight of accumulating sediments^4,5,6. However, human activities, particularly the extensive extraction of groundwater for agricultural and industrial use, have significantly exacerbated land subsidence in the area. The Choshui delta alone accounts for approximately 75% of Taiwan’s total land subsidence, covering an area of about 600 km² as of 2015. By the end of 2021, the region remained highly vulnerable to land subsidence, experiencing the highest rates in central Taiwan, with a maximum subsidence rate of 7.8 cm per year. This concerning trend highlights the urgent need for comprehensive monitoring and mitigation strategies to address the detrimental impacts on infrastructure, ecosystems, and local communities.

Figure 1 shows the location of the study area and the cumulative subsidence data from 2015 to 2023 (this figure was created with ArcGIS 10.8 software, obtained from https://www.arcgis.com/index.html), emphasizing the long-term effects of groundwater over-exploitation, which is a potential primary driver of subsidence in western Taiwan. Notably, the subsidence funnels have gradually shifted from coastal to inland areas, revealing the spatial and temporal evolution of subsidence patterns. This underscores the importance of effective groundwater management to mitigate further land instability. In addition, Fig. 1 is intended to provide context by showing the regions in Taiwan with severe land subsidence currently targeted for subsidence mitigation efforts and the locations of all MLCW stations within those regions as shown in the yellow diamonds in Fig. 1. The STES MLCW data were used in our study as a demonstration case. The purpose of including all MLCW locations in Fig. 1 is to highlight that STES is one of several stations in the critical subsidence areas, and to indicate that our methodology could be applied to other MLCW sites in these severely subsiding regions in the future.

Figure 2 presents thematic maps of the Choshui Delta, including: (a) accumulated subsidence from 2015 to 2023; (b) fine-grained soil content; (c) electric power consumption density distribution; (d) managed wells distribution. As illustrated in Fig. 2a, accumulated subsidence from 2015 to 2023 was assessed using leveling survey data collected (WRA, 2023). Accumulated subsidence was estimated using inverse distance weighting (IDW) interpolation within a GIS framework. Regarding spatial resolution, the IDW interpolation has a spatial resolution of 250 m × 250 m. Figure 2b presents the spatial distribution of fine-grained soil content, derived from comprehensive geological investigations that included 104 borehole logs. Each borehole was classified using the Unified Soil Classification System (USCS), where fine-grained soils are defined as materials with 50% or more passing through a No. 200 sieve, such as fine sand, silt, and clay. To generate a continuous surface from the point-based borehole data, the IDW interpolation method was employed within a GIS environment. This approach estimates values at unsampled locations by assigning greater influence to nearby data points, inversely related to their distance. The fine-grained soil percentage at each borehole was calculated by dividing the total thickness of fine-grained layers by the total drilling depth, resulting in a spatially continuous map of fine-grained soil distribution across the study area. Figure 2c and d show the spatial distribution of electric power consumption density and the locations of managed wells, respectively. Managed wells in this study refer to groundwater extraction wells that are officially registered under a government management system for simple agricultural wells. These wells differ from fully monitored or instrumented pumping wells, as they are manually operated by farmers and are not equipped with flow meters. Consequently, direct measurements of groundwater withdrawal volumes are unavailable. Instead, electricity consumption data, recorded via installed electricity meters, serve as a proxy for estimating groundwater use. The scale of agricultural groundwater extraction in the study area is substantial. Across the Choshui Delta, nearly 200,000 such managed wells have been registered. Specifically, within a 500-m radius of the Xiutan Elementary School MLCW (STES), a total of 141 officially registered wells are present, as illustrated in Fig. 4. All of these are simple agricultural wells primarily used for irrigation purposes. Electricity usage data from these wells are aggregated and analyzed to infer spatial and temporal patterns of groundwater abstraction in the region.

Since 2008, 31 MLCWs have been established, particularly in areas like Yunlin County’s Tuku, Yuanchang, and Huwei Townships, which experience the most significant subsidence. The MLCW is a specialized monitoring instrument established by Taiwan’s Water Resources Agency to assess land subsidence, especially in areas vulnerable to excessive groundwater extraction. MLCWs are equipped with multiple sensors installed at various depths, which measure changes in the vertical distance between layers over time. These depth-specific observations allow for precise monitoring of soil compaction and deformation at different strata, providing invaluable data for evaluating the mechanisms of land subsidence.

The MLCW at Xiutan Elementary School (STES) in Tuku Township, located in one of the most subsided areas, has been selected for machine learning-based time-series analysis to predict land subsidence caused by groundwater depletion. At STES, the subsurface profile extends to a depth of 300 m and consists of predominantly coarse to medium sand and fine-grained clay, with interbedded layers. Figure 3 illustrates that the soil consists of up to 19 layers, with varying thicknesses of coarse and fine-grained materials, including significant clay intervals between 140–163 m and 218–255 m. The complex geological profile poses challenges for establishing a conventional hydrogeological model based on physical consolidation mechanisms for the area.

According to Fig. 3, the analysis of the stratigraphic column and the compression observations relative to the bottom of the MLCW at STES indicates the occurrence of compression behavior across the depth range of 0 to 300 m. Particularly, the highest compression rate is observed within the depth range of 9 m to 63 m in recent years. Based on the monitoring results, the primary compression layer of the MLCW at STES is located in the shallow soil layer (between 0 and 63 m). The soil distribution at this depth primarily consist of alternating layers of medium sand, fine sand and silt. The compression at STES is predominantly influenced by fluctuations in the groundwater level, and it demonstrates a continuous increase as the water level gradually decreases over time. Figure 3 also shows that significant subsidence is observed within a depth of 63 m in terms of compression rate.

Land subsidence is a direct result of excessive groundwater extraction, necessitating a thorough examination of groundwater utilization patterns. While data on groundwater usage specific to well discharge is unavailable, the electricity consumption by wells can serve as a major indicator for assessing groundwater usage. The land subsidence observed at an MLCW is affected by groundwater usage within a specific buffer radius. To determine the appropriate buffer radius, the correlation coefficient was calculated between well electricity consumption and cumulative subsidence for different buffer radii ranging from 100 to 2000 m, considering subsidence depths from 0 to 60 m. A buffer radius of 500 m yielded the highest correlation coefficient of 0.91. Therefore, a buffer radius of 500 m was selected for further analyses. Figure 4 depicts the map of land use, managed wells, and monitoring wells within a 500 m radius of the STES. Within a radius of 500 m from the MLCW at STES, there are a total of 141 managed wells. According to Fig. 4, it appears that the land use types within a 500 m radius of the STES, primarily dominated by agricultural land use. The location of multi-layer groundwater level monitoring well (MLGLW) with a depth of 134 m below the surface, is situated adjacent to the MLCW at STES.

Figure 5 illustrates the relationship between drawdown and settlement at STES from 2017 to 2023. In this study, the drawdown and settlement values for 2017 are both set to zero to facilitate the comparison of changes in subsequent years. The results indicate that in 2021, following a severe drought, the drawdown reached its lowest recorded level in recent years. This substantial decline in groundwater levels led to an increase in settlement, resulting in significant and irreversible deformation. The results demonstrate the phenomenon of elastic–plastic deformation, which refers to the combined behavior of a material under stress.

The observed phenomenon—namely that settlement continues even after groundwater levels begin to recover—is primarily due to the time-lagged mechanical response of compressible soil layers³³. This means that settlement does not respond instantaneously to pore pressure changes, especially in fine-grained soil, which exhibit delayed consolidation. Moreover, the methodology used to construct.

Figure 5 draws on the concept of stress-compaction analysis, originally developed by³³ in their seminal work. In this framework, a correlation is established between cumulative drawdown and cumulative compaction (settlement). This figure helps to reveal whether the soil system undergoes elastic or inelastic (irreversible) deformation. In our case, the continued increase in settlement despite water level recovery suggests that irreversible compaction has occurred—an indication that permanent aquitard deformation has taken place. Therefore, even when groundwater levels rebound, the compacted soil layers do not recover their original thickness, and settlement continues due to prior stress history. Crucially, irreversible compaction occurs only when groundwater heads fall below the pre-consolidation level. Figure 5 is used to assess whether the subsurface system exhibits elastic or inelastic (permanent) deformation in response to groundwater level changes. The analysis indicates that groundwater heads have declined below the pre-consolidation threshold. Specifically, the data reveal that a drawdown exceeding approximately 6 m marks the onset of irreversible (inelastic) compression. This suggests that once groundwater levels fall below this threshold, permanent deformation of the soil strata is likely to occur.

Hydrogeologic characteristics of the STES

To evaluate the hydrogeologic characteristics of the monitoring well, we analyzed three borehole logs in the vicinity of the STES monitoring site, as shown in Fig. 6, and classified the subsurface stratigraphy based on grain size distribution. From this analysis, the aquifer system in the study area was delineated into four aquifers separated by three aquitards, as shown in Fig. 7. The groundwater level data used in the proposed machine learning model were collected from a monitoring well screened at a depth of 134 m, corresponding to Aquifer 2, which can be generally categorized as a leaky confined aquifer. However, based on site-specific in situ pumping tests and lithological analysis, we observed that the hydrogeological layering in this region is complex. The boundaries between aquifers and aquitards are not distinctly isolated, and fine-grained materials are interbedded within the aquifer zones.

To verify the hydraulic connectivity between aquifers, we conducted two controlled pumping tests. Under identical pumping conditions, the groundwater drawdown responses in Aquifers 1 and 2 yielded correlation coefficients of 0.86 and 0.94, respectively. These results, now presented in newly added Fig. 8 and Fig. 9, indicate strong vertical hydraulic connectivity between the upper and lower aquifers. This implies that although the water level data were obtained from a deeper well, they are influenced by near-surface pumping activity due to aquifer interconnection.

Datasets

In this study, multiple datasets were utilized, including accumulation of subsidence data, groundwater level fluctuations, electricity consumption of managed wells, and rainfall, as outlined in Table 1. The data availability, temporal resolution, and sources for each dataset are summarized in Table 1.

The data for the accumulation of subsidence, groundwater level, and electric power consumption are sourced from the Water Resources Agency (WRA), while precipitation data is obtained from the Central Weather Administration, Taiwan. It is worth mentioning that there is an absence of time series subsidence data for the period between 2012 and 2014 at STES.

Figure 10 illustrates the cumulative compression spanning from 2008 to 2023, utilizing data sourced from the WRA. Because subsidence data are missing from March 2012 to March 2014 at this site, the analyses in this study focus on the period with available data (2008–2023). We apply a data reconstruction technique to fill the 2012–2014 gap in the subsidence time series at STES. Figure 11 illustrates monthly variations in well electricity consumption and groundwater levels at STES. The well electricity use and groundwater level exhibit distinct seasonal patterns corresponding to the wet and dry seasons: higher electricity consumption, indicative of more pumping coincides with lower groundwater levels, and vice versa. Figure 12 shows monthly precipitation and groundwater level; groundwater levels tend to decline during the dry season and recover during the wet season, demonstrating the direct influence of rainfall on groundwater recharge. Detailed descriptions of these time series datasets are provided in the following sections.

Monthly compression data

Compression at different soil depths is monitored using MLCW, which measure depth variations up to 340 m with 1 mm precision⁹. In this study, we use monthly compression changes from MLCWs installed at STES—the area with the highest subsidence—as input for the LSTM model. Figure 10 shows monthly compression, rebound, and cumulative compression at STES. Since subsidence data from 2012 to 2014 are missing, we apply LSTM to reconstruct this gap.

Monthly groundwater level data

Groundwater use may be driving subsidence in the study area^7,8. To investigate this, groundwater level data were collected. At STES, the MLGLW reaches a depth of 134 m. Figure 11 shows the relationship between monthly electricity use of managed wells and groundwater levels, highlighting clear seasonal fluctuations. Between 2008 and 2023, groundwater levels varied by up to 15 m within the aquifer at that depth.

Monthly electric power consumption data

Monthly electricity consumption data were collected for 141 managed wells within a 500 m radius of the MLCW at STES. These wells range in depth from 8 to 150 m, with an average depth of 27.5 m. On average, they consume 78 kWh per month. As shown in Fig. 11, electricity use shows seasonal variation, mirroring fluctuations in groundwater levels—higher consumption aligns with lower groundwater levels, and vice versa.

Monthly rainfall data

Precipitation data were sourced from the Central Weather Bureau. As shown in Fig. 12, both total monthly precipitation and groundwater levels follow a clear seasonal pattern, with about 80% of annual rainfall occurring from June to September. Increases in rainfall correspond to rising groundwater levels, while decreases lead to declines, highlighting the direct influence of precipitation on groundwater variation.