Measuring and Interpreting Multilayer Aquifer‐System Compactions for a Sustainable Groundwater‐System Development

Ever decreasing water resources and climate change have driven the increasing use of groundwater causing land subsidence in many countries. Geodetic sensors such as InSAR, GPS and leveling can detect surface deformation but cannot measure subsurface deformation. A single‐well, single‐depth extensometer can be used to measure subsurface deformation, but it cannot delineate the depths of major compaction and provide insight about the deformation mechanism throughout a complex aquifer system, unless man extensometers at different depths are used. We present a multilayer compaction well (MLCW), installed in a borehole, that uses magnetic rings to detect stratum compaction at 25 depths as deep as 300 m below land surface. Our laboratory and field assessments indicate 1 mm precision and accuracy for one single‐depth magnetic reading. We tested the performance of MLCW by measuring aquifer‐system compaction over the proximal, middle, and distal fans of the Choushui River Alluvial Fan (CRAF) that has long experienced severe land subsidence. The MLCW measurements were used to create time‐depth diagrams of compaction, showing different compaction rates at different layers of aquifers and aquitards to identify the depths of major compactions. The elastic (reversible) and inelastic (irreversible) compactions from MLCW were used in stress‐strain analyses to estimate skeletal specific storages and the safe groundwater levels, below which groundwater extractions have caused irreversible compactions. The hydrogeological parameters derived from MLCW measurements can help governmental agencies to determine effective land‐use and water‐use policies, and ascertain the best strategy for utilizing artificial recharge to prevent land subsidence and achieve sustainable groundwater management.

rather than depth-specific deformation that can be indicative of the cause of the vertical deformation below the surface. In contrast, deformation measurements below the surface at specific depth intervals can provide clues for understanding the cause and mechanism of land subsidence, especially when combined with GWL and geological data.
Previous works on direct determination of subsurface deformation are based on measurements from extensometers (Metzger et al., 2002;Riley, 1969Riley, , 1984Ye et al., 2016). For example, Riley (1969) combined extensometer and hydrograph records in California to produce stress-strain diagrams for estimating hydrogeological parameters at specified depth intervals along a vertical profile. In the literature, extensometer measurements have been used to (1) describe the nature and the depth interval of aquifer-system compaction, (2) analyze the hydrogeological units responsible for the historical and current compaction rates, and (3) determine the recoverable and nonrecoverable storage properties and vertical hydraulic conductivity values of the compacting units (Burbey, 2003;Galloway et al., 1999;Liu & Helm, 2008).
To measure subsurface deformations at multiple depths, several extensometers must be installed. Installing multiple extensometers can be costly and requires a large land area. An alternative device for single-site, multiple-depth measurements is a device based on the fiber grating method . To our knowledge, the main challenge with optic fiber sensors is how an optical fiber or multiple optical fibers is (are) anchored into the soil in a stratum, so that vertical movements of the stratum are effectively detected by the fiber sensor. An optical fiber can be very thin, making it hard to lower it to deep layers and to stabilize it. Another alternative is multiple-layer compaction monitoring well(MLCW; Hung et al., 2010Hung et al., , 2012Hung et al., , 2018. A MLCW is a device that can collect readings at pre-anchored magnetic rings in a well to detect subsurface compaction between the magnetic rings. A MLCW is characterized by multilayered, flexible settings for monitoring depth-dependent compaction to a depth of 300 m, and provide data through time just like a snapshot of compaction across depth. In this study, we show an improved, semiautomatic version of the MLCW and a detailed method for installing it. The precision and accuracy of the improved MLCW are determined by laboratory and field tests. We use the MLCW measurements collected at several sites over the Choushui River Alluvial Fan (CRAF; see Section 4 for its location) in central Taiwan to demonstrate the usefulness of these measurements. The CRAF is the most groundwater-stressed region in Taiwan and has suffered from severe land subsidence since 1980. Specifically, we use MLCW measurements here to (1) identify the depths of major groundwater withdrawal, (2) determine stress-strain relationships associated with elastic and inelastic aquifer-system compactions, and (3) estimate the safe GWL above which the risk of irreversible land subsidence is low. We also discuss the potential application of MLCW measurements to monitoring the safety of Taiwan High Speed Rail (THSR). Figure 1 shows the workflow for installing a MLCW. A MLCW site is selected based on the following criteria:

Assessing a Candidate Site for Installation
1. The site is near the center of a subsidence-affected area 2. The site is near a groundwater monitoring well, continuous Global Navigation Satellite System (GNSS) station and a precision leveling benchmark 3. It is easy to obtain land access permits around the site 4. The site is permanent and will not be damaged by future constructions 5. The usable land around the site is sufficiently large to allow the construction of the MLCW Before installing a MLCW at the selected site, it is recommended to understand the aquifer formations near the well, especially the depths and the distribution of the layers with large compaction potential. This information is needed when determining the depths to deploy magnetic rings. At a MLCW site, the depths of the deployed magnetic rings are based on three models and datasets: (1) a hydrogeological model, (2) geotechnical well logging data, and (3) geophysical well logging data.
In the CRAF, the needed hydrogeological model has been constructed by the Central Geological Survey (CGS) of Taiwan (Chen & Jiang, 1999), based on 88 geological drillings made from 1992 through 1998. These drillings were all down to 300 m and have a total profile length of 17,051.5 m. Geotechnical well logging uses soils sampled at boreholes around a selected MLCW site to determine particle sizes, specific gravity, and the Atterberg limits.
Typically, a 100 m borehole is drilled within about 20 m of a selected MLCW site. The second borehole is the MLCW and is drilled to a depth of 300 m to analyze the aquifer formations along this borehole by geophysical well loggings. The 300 m borehole is new and is not the continuation of the 100 m borehole because the 100 m borehole often collapses after sampling the aquifer formations and is no longer suitable for further drilling. The two soil test datasets collected at the 100 m and 300 m boreholes are then used to determine the boundaries of aquifer formations at the 300 m borehole, along which the magnetic rings are then installed.
HUNG ET AL.  Geophysical well logging data are used to assess the hydrogeological properties of the soils and to distinguish between the different strata in the 300 m borehole. Specifically, for a geophysical well logging, we measure the following parameters: normal resistivity, single point (SP) resistivity, lateral resistivity, natural potential and spontaneous potential, natural gamma (NG), and temperature. These measurements can be used to detect the fine-grained sediments with the compaction potentials. This geophysical well logging method has been used by Reese and Wacker (2007). Figure 2 shows an example of magnetic ring installation at the MLCW site YWJS. The 25 rings (Figure 2b) are installed at the boundaries between the major aquifers defined by the CGS (layers B1-B3, Figure 2c) and at the boundaries between the fine-and coarse-grained soils defined by the geotechnical well logging data in the 100 m borehole and geophysical well logging data in the 300 m borehole ( Figure 2a, the borehole for the MLCW). Figure 2a shows different signal responses to different geological materials. For example, the fine-grained soil (clay and silt) results in low SP resistivity and high NG, whereas the coarse-grained soil (sand and gravel) results in high SP and low NG. In a coarse-grained soil, there is a marked difference between the short normal (SN) resistivity and the long normal (LN) resistivity due to high permeability, but the difference between SN and LN is small in a fine-grained soil. The variations in the differences also help to distinguish between the two kinds of soils. After selecting the depths for the magnetic rings, the rings are attached to PVC pipes which are then inserted into the borehole. After the pipe insertion, the magnetic rings are separated from the pipes by burning the materials between the pipes and the rings using electric short circuits. After this stage, the anchored magnetic rings can move with the sediment, rather than the pipes. A final step is to fill sands into the space between the pipes and the borehole.

Detecting the Position of a Magnetic Ring Center for Invar Tape Reading
The reading of the Invar tape when the probe is at the center of an anchored magnetic ring is the essential measurement of a MLCW. In our previous work (Hung et al., 2010), the detection of the ring center was made using analog signals produced by an electric current meter. In the improved MLCW, the detection of the ring center is made digitally (see Figure 3 for its components). Using the principle of electromagnetic interaction, we designed a measuring system that works as follows: when the probe's center (corresponding to the position of reading) is at the center of the magnetic ring, the electric current is zero and an Invar tape reading is made manually. As shown in Figure 3, a device is attached to the probe and broadcasts the current's signal at the rate 100 HZ to the receiver. The signals (current readings) are enhanced by the signal line in Figure 3. The signal line is 310 m long and is a specially made wire that can intercept the broadcast HUNG ET AL.   Chen and Jiang (1999). signal sent to the receiver. It is attached to the external side of the PVC pipes during the pipe-insertion stage to avoid later interruptions with the movement of the probe. The receiver receives the signal, which is displayed on the tablet PC for deciding whether the probe has reached the ring center (this happens when the current is zero). The broadcasting device is able to correct for the effect of temperature variation on magnetic fields. To ensure stable broadcasted signals at greater depths, we use three frequencies for the broadcasting device.

Assessing the Precision and Accuracy of MLCW Invar Tape Readings at a Laboratory and at Field
The precision and accuracy of depth readings of the Invar tape at the magnetic rings were assessed annually using laboratory and field tests described below. Currently, we have not detected any deterioration in any component of the MLCW. A depth reading of the Invar tape shows the depth of a magnetic ring ( Figure 3) and is affected by the precision and accuracy in determining the zero electric current position based on the signals from the broadcasting device (Figure 3), which in turn are affected by the ring's depth, groundwater temperature variation from the ring to the receiver, and systematic errors (if existing) of the Invar tape. Here, precision is the extent of repeatability (divergence) of repeated Invar tape reading at the same magnetic ring, and accuracy is the extent of closeness to the true depth of the ring.
The first test is to assess the precision of readings. This test was conducted in a laboratory at the National Cheng Kung University (NCKU), Tainan, Taiwan. As shown in Figure 4a, the probe was tied to a robotic arm that is part of the tri-axis testing platform in the laboratory. The robotic arm moved the probe until the resulting electric current reading (Section 2.1.2) indicates that the probe is at the magnetic ring center.
HUNG ET AL.  When this happened, an Invar tape reading was made. A reading is a measured vertical distance from a reference point at the tri-axis platform to the magnetic ring center. As shown in Figure 4b, the horizontal (x) and vertical coordinates (z) of the five dots have been determined before the test. The robotic arm moved the probe to locations near the five dots to determine the center of the magnetic ring. At each of the five dots, we collected five readings and then computed the standard deviation of the readings. Because the ring is at a level surface, the readings (25 in total) around the five dots should remain close. That is, the repeatability of the readings around the five dots indicate (1) how close the magnetic center is reached by the probe, and (2) how close the orientation of the ring is to a horizontal plane. The result shows that the standard error of the 25 readings is 0.4 mm.
We also conducted a field test to determine the repeatability (precision) of the Invar tape readings associated with the probe. The test was at the MLCW well at the Nan Gwang Elementary School (NGES), and the repeatability test was conducted to a 300 m depth. The test result is shown in Table 1. From Table 1, we infer that the NCKU laboratory test resulted in slightly better repeatability than the field test at NGES. The cause of the larger repeatability (poorer precision) from the field test is that the signals transmitted from the broadcasting device can be disturbed by temperature variations, water in the pipe and other factors that did not occur at the NCKU laboratory (in a dry condition up to 10 m).
The second test is to determine the accuracy of Invar readings and was conducted in the same NCKU laboratory. A device simulating the PVC pipe of a MLCW with 10 magnetic rings was installed every 1 m along the pipe ( Figure 5). The vertical distances of the rings along the pipe were determined to better than 1 mm before the test, and are regarded as "true" depths that can be used to determine whether depths from Invar tape readings of the MLCW are consistent with the true depths defined by the positions at the magnetic 10 rings.
At each of the 10 rings, five Invar tape readings were made to determine a mean reading. The difference between this mean value and the true depth was then computed. The standard deviation of the 10 differences is 0.8 mm (Table 1), which indicates the overall accuracy of an Invar tape reading. Note this accuracy (0.8 mm) was obtained to a depth of 10 m. An earlier work (Hung et al., 2012) indicates that this accuracy is applicable to depths of 300 m because the Invar tape has been calibrated against systematic errors.

Use of Compaction Measurements for Estimating Safe Groundwater Levels
A MLCW can measure compaction at multiple depth intervals with magnetic rings. If such measurements for a range of depths are collected over a sufficiently long time (>1 year), they may be useful in detecting elastic (recoverable) and inelastic (nonrecoverable) compactions due to changes in groundwater level (called GWL below for short) of natural or anthropogenic origins. GWL is the depth to groundwater in a well, relative to sea level. Elastic compaction is relatively small, mostly seasonal (often associated with HUNG ET AL.

Case
Test

Table 1 The Precision and Accuracy (the Result Column) Based on the Multilayer Compaction Well (MLCW) Measurements at the National Cheng Kung University (NCKU) Laboratory and the Precision From a Field Test
rainfall cycles) and is recoverable (Bell et al., 2008;Ezquerro et al., 2014). In contrast, inelastic compaction is permanent and is mostly due to excessive groundwater pumping that lowers hydraulic head to levels below a critical head threshold. Identifying this critical head can be approximated using simultaneous measurements of hydraulic head (or in general GWL) and compaction but is more difficult if residual compaction is an important component of total compaction. Residual compaction occurs in thick aquitards as heads in the thick aquitards equilibrate with heads in the surrounding aquifers (Terzaghi, 1925). Depending on the thickness and the vertical hydraulic diffusivity of a thick aquitard, fluid-pressure equilibration lags behind pressure (or hydraulic head) declines in the adjacent aquifers; associated compaction can require decades or centuries to approach completion. Thus, if the aquifer head declines below the previous lowest level for a relatively short period, the preconsolidation head in the aquitard is not necessarily reset to the new low value (Phillips et al., 2003).
Here, we summarize the theory for determining this limit and how compaction measurements from MLCW can provide the needed information, based on the works by Terzaghi (1925), Riley (1969), Burbey (2001), Sneed and Galloway (2000), Pope and Burbey (2004), Liu and Helm (2008), and Galloway and Burbey (2011), among other authors. First, the effective stress  e of a sediment layer can be expressed as where  T is the total stress and p is the fluid or pore water pressure. In a confined aquifer where total stress remains constant, such as F2 over the CRAF (see Section 5), the change in effective stress is due to the change in pore water pressure. That is, where Δ is change, ρ w is fluid (water) density, g is gravity and h is GWL (hydraulic head). Thus, a negative change in GWL (meaning a lowered GWL) will lead to a positive change (increase) in effective stress, which in turn creates a positive strain, τ, as where ΔB and B 0 are the change in thickness, or new distance between magnetic rings measured by MLCW, and the initial thickness, or original distance between magnetic rings, respectively. Note that ΔB is positive for compaction and is negative for expansion (in case of increasing pore pressures from rising GWL). When compaction occurs within a layer, the vertical distance measured by the two magnetic rings at the two ends of the layer is shorter than the original distance. Land subsidence measured by precision leveling or GPS at a MLCW site is the cumulative compaction measured by all magnetic rings, plus additional compaction occurring below the depth of the deepest ring (300 m in this study).
A hydraulic head (GWL) variation (Δh) and its resulting strain (Equation 3) can be used to determine the skeletal specific storage (S sk ) as which generally varies by one or two orders of magnitude depending on whether the effective stress exceeds a stress threshold. The greatest historical effective stress imposed on the aquifer system, sometimes the result of the lowest GWL, is the critical stress threshold, or the "preconsolidation stress," and the corresponding (lowest) GWL is the "critical head" (Leake & Prudic, 1991). When the effective stress is less than the preconsolidation stress, the compaction is elastic. Conversely, when the effective stress is greater than the preconsolidation stress, the compaction is inelastic. The preconsolidation stress can be used to distinguish between two different skeletal specific storages, corresponding to elastic and inelastic compaction as where  max is the preconsolidation stress, and S ske and S skv are the skeletal specific storages in the cases of elastic and inelastic compaction, respectively.
In this study, we define the safe GWL (relative to mean sea level) as the lowest GWL at which elastic compaction occurs. If GWLs decline below the safe GWL, inelastic compaction will occur, causing irreversible compaction and resultant land subsidence. Using paired measurements of strain and GWLs below the safe level, we can estimate S skv using the resulting slope (Equation 4) from a stress-strain diagram (see Section 5). Furthermore, we can use paired measurements of strain and GWLs above the safe level to estimate S ske using the resulting slope (Equation 4) from a stress-strain diagram. Examples of estimating skeletal specific storages using paired MLCW and GWL measurements are given in Section 5.

The Choushui River Alluvial Fan and Its Three Subfans
We tested the new MLCW in the southern part (Yunlin County) of the CRAF, which is the most important agricultural area on the west coast of Taiwan ( Figure 6). The CRAF is a flood plain created mainly by the Choushui River. It is bordered by Wu River in the north, Beigang Creek in the south, and the Douliu Hills in the east. The CRAF ends in the Taiwan Strait to the west and covers an area about 2,000 km 2 . Figure 7 shows the hydrogeological conceptual model of the CRAF. The entire formation comprising the CRAF can be divided into four aquifers (F1-F4) and four aquitards (T1-T4). The thickness of the sediments in the CRAF ranges from 750 to 3,000 m (Lin et al., 1992). According to the CGS Chen & Jiang, 1999), the CRAF can be divided into three subfans: proximal fan (fan top), middle fan (fan center), and distal fan (fan tail). The fan top spreads westward from Yunlin's foothills and the sediment is mostly composed of gravel and coarse sand. The fan center is a transition zone between the fan top and the fan tail, and is mostly composed of interlayers of sand and mud. The fan tail is mostly composed of interlaced layers of fine sand and mud due to the long distance from the sediment source (Chen & Jiang, 1999).
Because there are some highly compressible sediments in the fan center and fan tail, withdrawing groundwater without sufficient recharge will decrease the water level, reducing the pore pressure and increasing the effective stress beyond the preconsolidation stress, inevitably leading to land subsidence. Figure 8 shows the annual rates of land subsidence from 2017 to 2018 in Yunlin from precision leveling over a network of leveling benchmarks (Yang et al., 2019). Tuku and Yuanchang are two townships that experience the most severe land subsidence, with rates exceeding 5 cm/year. To assess the causes of land subsidence in Yunlin, 22 MLCWs ( Figure 8) were installed. These MLCWs have been equipped with magnetic rings, allowing measurements of compaction at different depth intervals using the new device presented in Section 2. Before the introduction of the automatic reading device described in this study (see Section 2.2), readings at the magnetic rings were made manually.

Groundwater Pumping Wells and Monitoring Wells
Because the current demand of surface water in Yunlin has exceeded the supply from rivers, here most of the water needed for agricultural, industrial, and domestic uses is extracted from groundwater. According to a survey conducted by the Ministry of Economic Affairs (Hsu & Lin, 2018), there are 125,699 pumping wells in Yunlin, of which 115,397 wells are for agricultural use, 356 wells are for industrial use, 4,398 wells are for domestic use (private and municipal drinking water), 5,102 wells are for aquaculture use, and the remaining 446 wells are for other uses. Figure 9 HUNG ET AL.
10.1029/2020WR028194 9 of 19  shows the distribution of these wells, which have different screen depths. The use of different screen depths permits withdrawal of groundwater from different aquifers. The depths of the wells for agricultural use are mostly shallow (the depths range from 0 to 100 m), and these wells withdraw groundwater from F1 and F2 (less energy for pumping is the deciding factor for depth). The depths for industrial and domestic uses are greater and the wells pump groundwater from F3 and F4 (Figure 7; less expense for processing deep groundwater is the deciding factor). As such, within a small area where several wells for different purposes are constructed, land subsidence (at land surface) is a result of compaction at different depths. A MLCW can be used to measure such depth-dependent compaction to identify the depth(s) of the aquifer system that contribute the most to the subsidence, and identify the depths, where mitigation measures would be effective for controlling land subsidence.
In the land subsidence-affected region of Yunlin, a multiple-monitoring system has been installed. The system consists of continuous GPS stations, benchmarks for precision leveling, MLCWs, and groundwater monitoring wells (Figure 8; benchmarks are not shown). The technique of interferometric synthetic aperture radar (InSAR) has also been used to monitor land subsidence here (Yang et al., 2019).

Analyses of Compaction in the Proximal Fan, Middle Fan, and Distal Fan of the CRAF for Safe Groundwater Extraction
Section 5 analyzes the compaction measurements and their relationships to GW-level changes and rainfall in the proximal, middle and distal fans of the CRAF. Conclusions about maximum stresses and safe GWLs (Section 3) are drawn from these compaction-GWL analyses. Table 2 summarizes the soil types, GWL changes, compaction, and other information for the three fans. The three subsections of Section 5 show the results for the three fans.
The GWL amplitudes and slopes in Table 2   compaction or expansion in response to GWL changes because this is an issue that cannot be addressed by our GWL and MLCW measurements in this study. A sample paper on such time delays is Hung et al. (2012).

Proximal Fan
Our analysis begins with data at site JNES located in the proximal fan. Figure 10a shows the cumulative compactions from the starting epoch (12/15/2014) to five ending epochs at JNES. The cumulative compaction at a depth in Figure 10a is the total compaction from 300 m to the depth (this definition is applicable to the compactions in Figure 12a and Figure 14a below). Figure 10b shows the time-depth diagram of compaction (TDDC) at JNES over the depth range of 0-300 m. A compaction in Figure 10b is the shortening between two successive magnetic rings (this definition is applicable to the compactions in Figures 12b  and 14b below). According to our field investigation, the evident compaction from depth 0 to 10 m at JNES was caused by constructions of walls and a security building, rather than groundwater pumping. These man-made structures create loads that compress the soil around them. Figure 11a shows the cumulative compaction at JNES for F1-F4, GWLs at two different depths at Huxi, and rainfall at LNJS (Figure 8). Compaction in Figure 11a is the change in the distance between two magnetic rings that bracket each aquifer.
GWLs at Huxi-1 and Huxi-2 fluctuate seasonally between 26 and 31 m above sea level but are fairly stable annually. The long-term rate of GWL change here is slightly positive (meaning increasing GWLs) ( Table 2). The amplitude of the annual oscillation of GWL at Huxi-2 is larger than that at Huxi-1. Anomalous lows occurred in the dry seasons (from November to the next May) of 2015, 2017, and 2018. In general, May is the month with the lowest GWL in the proximal fan, where GWL rose soon after rainfall. Figure 11a shows that when ignoring the upper 10 m, less than 1 cm of compaction occurred within the four aquifers at JNES. Here, the cumulative compaction in F2 oscillated with GWLs. Figure 11c (Table 2), the compaction here is expected to be elastic and this expectation is consistent with the result of our compaction measurements. In summary, the current groundwater pumping in the proximal fan does not result in land subsidence because the major soil types are sand and gravel where inelastic compaction is negligible. Consequently, any water level is a safe level. As such, Table 2 shows "na" in place of S skv and safe GWL at JNES. Figures 12a and 12b show respectively the cumulative compactions and the TDDC at site STES located in the middle fan. Figure 12b shows that the amount of compaction was largest in the depth interval 0-50 m, intermediate in the depth interval 50-200 m, and smallest in the depth interval 200-300 m. Figure 13a shows the cumulative compaction and the rainfall at STES, and GWLs at STES (depth 134 m), and HLES (depth 36 m, 5 km to STES; Figure 8).

Middle Fan
Both GWLs in F1 and F2 show significant seasonal variations, but with different declining mean annual rates. The rate in F1 is more than twice that in F2 (−0.18 vs. −0.07 m/year, Table 2). The lowest GWLs occurred in the dry seasons of 2015, 2017, and 2018. Like the proximal fan, the middle fan experienced the lowest GWL in May and GWL rose soon after rain.
The main compaction at STES occurs in the clay layer at the depths of 10-30 m ( Figure 12). Furthermore, there is a thick sand layer between depth 30 and 40 m. The permeability of this sand layer is high and thus this layer is the favorite layer to pump groundwater for irrigation. This heavy pumping of groundwater generated excessive stresses at the depths of 10-30 m, resulting in the major compaction seen in Figure 13a.
At STES, relatively large GWL changes occurred in F2. The GWLs in F2 over January-April in 2015, 2017, and 2018 were especially low. During these months, there were significant compaction in F2. By contrast, the compaction in F3 and F4 during these months were small.
A comparison between GWLs and rainfall (Figures 13a and 13b) suggests that GWLs increased significantly with rainfall at STES. During the dry seasons, GWLs dropped rapidly, causing large compaction. The largest GWL drop occurred at STES (Figure 13a). Figure 13a shows that the variations in compaction are affected by the variations in GWL.
HUNG ET AL.   Figure 13c shows that when the GWLs were above the −8 m level, the slopes of the lines fitting stresses and strains (green dashed), corresponding to the elastic compaction are similar. Using the theory in Section 3, we determined the skeletal specific storages S ske from the inverses of the slopes. Furthermore, when GWLs were below −8 m, the compaction become inelastic. In this case, we fit a line to the GWLs < −8 m and strains to determine the slope of the stress-strain relationship (red dashed line in Figure 13c); the inverse of this slope is the inelastic specific storage S skv . Table 2 shows the specific storages at STES for both the elastic and inelastic compactions, indicating that S skv is one order of magnitude larger than S ske . This shows that when the GWLs are below −8 m, the resulting compaction is irreversible and is much larger per unit of head change than the compaction when the GWLs are above −8 m. In general, the compactions are elastic when the GWLs are above −8 m (relative to mean sea level), below which the compaction becomes inelastic. Therefore, we infer that −8 m is the safe GWL for a safe groundwater yield at aquifer 2 around STES.

Distal Fan
Our last compaction analysis is in the distal fan using data at site YWJS. Figures 14a and 14b show respectively the cumulative compactions at YWJS and the resulting TDDC, indicating evident compaction occurred from depth 0-300 m, with the major compaction occurring in the depth interval 50-250 m. Figure 15a shows the cumulative compaction in F1-F4 and GW-level changes at two depths, and Figure 15b shows rainfall at station YWJS. The GWLs at YWJS (F2) ranged from −15 to −36 m, with a large decline rate of −0.36 m/year and an annual amplitude of 5.98 m. The range, rate, and amplitude of GWL changes in F3 at YWJS are similar to those in F2 (Table 2). Like those in the proximal and middle fans, GWLs in the distal fan were low in the dry seasons of 2015, 2017, and 2018 and the lowest occurred in May. The GWL also rose immediately after rain.
F2 (depths 55-160 m) at YWJS experienced the largest compaction, followed by F3, 4, and 1. Figure 15a indicates that during January-April in 2015, 2017, and 2018, there were evident drops in GWL, and the average compaction in 2018 was the largest. while the compaction in 2016 was elastic. From Figure 15c, we infer that −28 m is the approximate safe GWL around YWJS to prevent land subsidence from happening.

Toward Land Subsidence-Free Groundwater Use
The MLCW device in this study can measure stratum compaction at 25 depths with an accuracy and precision of 1 mm. Along with GWL data, measurements from MLCWs can be used to detect the layers of major compaction, how compaction is related to GWL changes and whether the compaction is inelastic or elastic. The analyses of the stress-strain relationships in Section 5 can offer critical information on the occurrence and evolution of compaction, and can be used to estimate the safe GWL. This information can be used by governmental agencies to prevent or mitigate land subsidence.
Compared to the middle and distal fans, our paired MLCW and groundwater-level measurements in the proximal fan indicate that the compaction in the proximal fan was elastic. This can be largely explained by geological conditions alone. The proximal fan contains coarse-grained materials (such as gravel and coarse sand); therefore, inelastic compaction of sands and gravels is negligible. Additionally, groundwater is easily replenished because of the high rainfall infiltration rate for the coarse-grained materials. In addition, the proximal fan groundwater is for commercial use, and the duration of commercial pumping is short. As such, the GWL in the proximal fan did not decline significantly during dry seasons. By contrast, the groundwater use from the middle and distal fans are dominated by rice plantations that need a large amount of water. The main rice growing periods are over January-June and July-December. The rice yield and rice price in January-June are higher than those in July-December. In addition, the typhoon season occurs during July-October, during which heavy rains are frequent, making the farmland in Yunlin prone to flooding and property damages. As such, the area of rice cultivation and the groundwater demand in January-June over the CRAF are twice those in July-December. Because the first period (January-June) overlaps with the dry season (January-May), groundwater is pumped to compensate for insufficient surface water for rice HUNG ET AL. irrigation. The pumping lowers GWLs to induce aquifer-system compaction and land subsidence, which in turn can be measured by MLCWs. Thus, MLCW measurements in January-May can help to understand the amount of pumped groundwater in a dry season to prevent land subsidence.
In Yunlin (in the southern CRAF), there are more than 120,000 pumping wells, resulting in land subsidence with the maximum annual rates exceeding 5 cm/year (Figure 8). Although this number of wells is well estimated, the actual pumping volume is not, so it is unclear how much water that these wells can pump without causing subsidence. Therefore, the government cannot formulate effective policies to manage the groundwater withdrawal of these wells. This difficulty in estimating safe yield of groundwater may be overcome by using MLCW observations, along with GWL observations. For example, the combination of these two kinds of observations indicate that the safe GWL at STES (middle fan of the CRAF) is about −8 m (relative to sea level), and that at YWJS (distal fan) is about −28 m. As shown in Section 5, the MLCW and GWL observations can also be used to estimate skeletal specific storages associated with elastic and inelastic compactions. The MLCW installation technique and the method for interpreting the MLCW observations over the CRAF can be applied to other parts of the world for monitoring multiple-layer compaction. In addition, as a preventive measure in a severe land subsidence area, surface water can be regularly used to recharge groundwater to reduce the chances of GWLs dropping below the safe level. This final goal is subsidence-free groundwater use.

Safe Groundwater Level and Sustainable Groundwater Development
The safe GWLs identified in the middle and distal fans (Sections 5.2 and 5.3) can be related to the sustainability of a developed aquifer system. One important parameter to understand when considering the sustainability of an aquifer system is the safe groundwater yield, which can be estimated by an annual water balance analysis using a three-dimensional numerical model based on the physics of groundwater flow, in situ hydrogeological parameters, GWL observations, and volumes of groundwater withdrawal. However, factors such as pumping and dynamic development of groundwater discharge can introduce large uncertainties to the resulting safe yield (Zhou, 2009). Over the CRAF, the difficulty in achieving a reliable estimate of safe groundwater yield is due to the following two problems. First, over the land subsidence-affected area, there are many private groundwater wells, where the depths of pumping wells and amounts of groundwater withdrawal are not precisely known. Second, the recharge of groundwater by the infiltration (input) of surface water has not been accurately estimated. These two problems make it impossible to carry out an accurate water balance analysis.
Instead of estimating a safe groundwater yield, MLCW measurements can be used to estimate the water extraction capacity of an aquifer system without introducing the risk of severe land subsidence. Specifically, MLCW measurements can be used to identify the safe GWL, above which the groundwater-pumping-induced compaction is elastic (recoverable). This safe level can be used to estimate the maximum safe volume of groundwater extraction, beyond which inelastic (nonrecoverable) compaction (land subsidence) will occur. Such MLCW measurements can be used to tighten the control of groundwater withdrawal at the main compaction layers (such as aquifer 1 near Sioutan) and can be dynamically fed to a three-dimensional numerical model to achieve a best equitable, safest use of a groundwater system. Furthermore, groundwater users can be advised to spread out the pumping time for the compaction layers, for example, by switching to nighttime pumping in some wells to avoid inelastic deformation due to the rapid drop of the water table beyond the safe level.

Application to the Safety Maintenance of Taiwan High Speed Rail
In central Taiwan near Yunlin, one particularly important application of MLCW is the safety maintenance of THSR (Figure 8). This application can be illustrated using the MLCW measurements at station STES. Since STES is near the Yunlin section of THSR, the MLCW measurements here can be used to assess the risk of THSR caused by land subsidence. A long section of THSR adopts a viaduct structure supported by friction piles to the depth of 70 m (Hwang et al., 2008). Therefore, when the soil above 70 m that holds the friction piles is compressed, the bearing capacity of the piles will be reduced because of the downward forces created by the soil compression (negative frictions). According to the design standard of the Japanese National Railway Construction and Construction (and followed by THSR), if the annual soil compaction rate at a THSR foundation pile exceeds 2 cm/year, a potential negative friction exerting on the pile must be considered in regular safety maintenances. If the compaction rate exceeds 4 cm/year, a 100% compaction-induced negative friction should be determined to assist the safety inspection of the pile. Therefore, if a long-line structure such as THSR passes through an area where its shallow soil layer is susceptible to compaction, a MLCW can be installed to determine the depth-dependent compaction rates to safeguard the structure. This is the advice (installing a MLCW around suspected piles) for the THSR section around STES because Figure 12a shows major compaction in the shallow soil layers around STES. Furthermore, two or more closely spaced MLCWs can provide data for monitoring differential settlements of THSR, which can also be detected by precision leveling, high-quality GPS measurements or high-resolution (spatially)/ high-precision InSAR analyses (see Section 1).

Conclusion and Prospect of a Fully Automated MLCW
One of the causes of the land subsidence in the CRAF is stratum compaction resulting from the over extraction of groundwater. Geodetic methods such as InSAR, leveling and GPS can be used to detect surface deformations only, rather than stratum compaction, which can be measured by a sensor such as MLCW presented in this study.
We show the details on how a MLCW is installed and how readings are recorded and interpreted. The MLCW measurements have been used to estimate skeletal specific storages and to infer the safe GWLs at aquifers in the middle and distal fans of the CRAF. Furthermore, it is possible to use MLCW measurements to construct a land subsidence model (Hung et al., 2012) to forecast land subsidence.
The MLCW measurements in this study were collected once a month and indicate compaction in F1 and F2 in the CRAF. In order to increase the MLCW sampling rate and improve its efficiency, we are developing a fully automated version of MLCW for continuous measurements (Hung et al., 2020), which can be used to analyze the real-time characteristics and amplitudes of compaction. The current sampling interval of this automated version is 10 min. Together with GWL measurements and a land subsidence model, continuous MLCW data can be used to infer the influence of depth-dependent pumping on stratum compaction and construct time-frequency spectra of land subsidence to provide a real-time reference for mitigating over-pumping and land subsidence in a groundwater-stressed aquifer.

Data Availability Statement
All the MLCW, groundwater level and rainfall data used in this study are freely available at http://space. cv.nctu.edu.tw/publications/#data.