Measuring Aquifer Specific Yields With Absolute Gravimetry: Result in the Choushui River Alluvial Fan and Mingchu Basin, Central Taiwan

Accurate and densely covered specific yields (Sy) are essential for estimating the storage capacity of a groundwater reservoir. A cross‐well pumping test can determine Sy, but its high cost often makes it unsuitable for sampling high‐resolution Sy. The gravity‐based method (GBM) based on gravity changes near existing groundwater wells may outperform cross‐well pumping tests in estimating Sy. We established 10 gravity sites close to groundwater wells in the aquifer‐rich Choushui River Alluvial Fan and Mingchu Basin in central Taiwan and measured gravity changes with two FG5 gravimeters from 2012 to 2017. Thirty‐nine Sy values with formal errors are determined by natural rises and falls in gravity and groundwater level. The representative Sy values (0.04 to 0.29) from GBM are in general consistent with those from cross‐well pumping tests (0.03 to 0.24). Repeated groundwater level changes over similar depth range at different times serve as revisit tests, showing that GBM can reproduce a reliable Sy value at a given site and depth. Soil moisture and compaction data show that the effects of gravity changes originating from unsaturated zones and deep aquifer layers are minor. Using the cylinder model for aquifers with limited lateral extents, we assess the validity of the Bouguer model by quantifying gravity differences and relative Sy differences originating from the model assumption. Improvements in environmental resilience and transportability achieved by recent atomic gravimeters may increase the potential of GBM to replace or supplement cross‐well pumping tests in densifying Sy point densities for an improved groundwater resource management.


Introduction
Specific yield (S y ) is an essential hydrogeological parameter for establishing groundwater models and assessing underground resources in an unconfined aquifer. S y is typically estimated by cross-well pumping test based on the hydraulic method (Boulton, 1963;Moench, 1995;Neuman, 1975). A S y value from a cross-well pumping test only represents the aquifer's storage capacity in the area between the pumping well and the observation well in the test. If the aquifer's hydrogeological property is heterogeneous, it can be challenging for the test to estimate a correct S y value. In addition, a S y value from a cross-well pumping test may be corrupted by large uncertainties. The uncertainties are caused by the fitting error of type-curve, inappropriate boundary conditions, storage effect, and well design (Chen et al., 1999;Kollet & Zlotnik, 2005;Mao et al., 2013;Raghavan, 2004;Sánchez-Vila et al., 1999;Yeh & Huang, 2009), among other reasons. In addition, the hydraulic method can be expensive and time-consuming because of the requirements of well construction and pumping. The high cost and inefficiency can often result in insufficiently sampled S y values over an aquifer, leading to a poorly estimated storage capacity of an aquifer. Thus, an alternative method for an accurate, fast, and low-cost estimation of area-averaged S y value is highly needed.
One alternative method is the gravity-based method (GBM), which uses gravity changes and groundwater mass changes to infer S y . Over an unconfined aquifer, groundwater changes between dry and wet seasons can induce significant gravity changes, depending on the magnitudes of S y and groundwater level oscillation. Following reasons show that GBM can be a good alternative to and can even outperform the hydraulic method: (1) Gravity changes reflect groundwater mass changes in the underlying aquifer over a lateral extent of several hundred meters, providing an effective S y that is more ergodic compared to the S y derived from the hydraulic method. (2) Gravity values are collected near an existing groundwater monitoring well, thus requiring no fees for constructing wells needed in the hydraulic method. (3) With the contemporary absolute gravimeter (AG) technology (but see the discussion in section 6.3), typically it takes 12 hr to measure a gravity value with a microgal accuracy; two such gravity values and two coincident groundwater level measurements (this means two gravity surveys) at two different times can be used to estimate an effective S y . Accordingly, gravity changes measured by AGs and relative gravimeters (RGs) over time spans of several months or years have been successfully used for estimating S y (Gehman et al., 2009;Howle et al., 2003;Pfeffer et al., 2011;Pool, 2008;Pool & Eychaner, 1995). In addition, a previous study by Chen et al. (2018) shows that GBM can determine S y using short-time gravity changes such as those happening during cross-well pumping tests and those induced by post-rain groundwater level recessions. Furthermore, superconducting gravimeters (SGs), which are the most precise RGs, can be used to measure gravity changes. Sample studies using SG over unsaturated zones are Kazama and Okubo (2009), Krause et al. (2009), Longuevergne et al. (2009), Kennedy et al. (2016), and Watlet et al. (2020). However, a SG is an observatory type of gravimeter and not designed for high mobility, making it expensive and time-consuming to measure gravity changes at multiple stations using a SG.
As surface water in Taiwan is in short supply due to uneven rainfalls, reservoir siltations, and growing water demands in a number of industries, there is an urgent need of the conjunctive use of groundwater and surface water in Taiwan. The Choushui River Alluvial Fan (CRAF) and Mingchu Basin (MB), both shown in Figures 1a-1c, contain two of the most important aquifer systems in central Taiwan for a potential conjunctive groundwater-surface water use and for supplying primary irrigation water. However, in the past several decades, central CRAF has experienced land subsidence due to over-pumping (Hwang et al., 2010). In response to the increasing groundwater demand, the Water Resource Agency and Central Geological Survey of Taiwan have been supporting several studies to investigate the groundwater resource over CRAF and MB. These studies, including pumping tests, hydraulic surveys, and development of efficient groundwater management schemes, have resulted in useful datasets for assessing the groundwater potential of CRAF and MB, and most importantly hydraulic S y values for this study.
Although S y values determined by GBM have been published in a number of studies, such S y values have not been fully compared with those from cross-well pumping tests in terms of aquifer storages implied by these two kinds of S y values. In addition, it has not been shown whether a gravity-based S y value can be reproduced using gravity measurements at different times over the same layer of an aquifer. That is, a revisit test of gravity-based S y has not been conducted. Also, there are little discussions about the impact of the gravity approximation of the Bouguer plate model on gravity-based S y estimations (see section 3.2). To address these issues, we have collected gravity data by AGs near 10 groundwater monitoring wells from 2012 to 2017 in the fan-top area of CRAF and MB. Near these 10 wells, S y values from cross-well pumping tests (hydraulic method) are available for comparison with the gravity-based S y values. At selected sites, we will also determine S y values over different phases of groundwater level changes (high-low and low-high) as revisit tests to show the repeatability and reliability of GBM in estimating S y . This paper aims to confirm that GBM can be widely used with high confidence for realistic estimates of aquifer storage coefficients needed in groundwater resource assessment.

The Choushui River Alluvial Fan and Mingchu Basin for the Gravity-Based S y Experiment
The study area (Figure 1c) is over the upper CRAF and MB, located in the middle reach of Choushui River. The geological map (Figure 1d) indicates that the upper CRAF and MB are largely composed of gravel and sand. Furthermore, the upper unconfined aquifer can be regarded as a single aquifer system (without major aquitards) that covers an area of about 200 km 2 . Major land subsidence (opposite of vertical displacement of land) occurs in central CRAF, with a largest subsidence rate of about 7 cm/year. A recent analysis of land subsidence over CRAF is given by Yang et al. (2019). The average annual precipitations in the catchment  Tables 1 and 2) are from cross-well pumping tests at sites close to the newly established gravity sites (yellow stars; the numbers in the brackets are the representative S y values from GBM in Table 2 and Appendix A). (b) The locations of CRAF and MB in Taiwan. (c) A hydrogeological profile parallel to Choushui River (orange line, from Chukuo to Haifeng). (d) Distribution of soil layers along the profile.

10.1029/2020WR027261
Water Resources Research area of Choushui River and CRAF are about 4,000 and 1,500 mm, respectively. Summer precipitations contribute mostly to the recharges of the aquifers in CRAF and MB. On average, the groundwater level here peaks from May to September. The average depths of groundwater level are about 10 to 20 m, and the annual groundwater level changes range from 0.5 to 8 m in the unconfined aquifers of the upper CRAF and MB.
The Water Resource Agency and Central Geological Survey of Taiwan have conducted more than 15 cross-well pumping tests over CRAF and MB in the past 25 years (Central Geological Survey, Taiwan, 2014; Figure 1 shows 12 of these 15 values). Table 1 shows the coordinates of 7 pumping test sites, the screen ranges of the pumping wells, and the S y values estimated from the hydraulic method near the primary recharge areas (Central Geological Survey, Taiwan, 2014) of CRAF and MB (one can use the coordinates to locate the 7 sites in Figure 1a). The current pumping test results show that S y values in CRAF and MB range from 0.03 to 0.24 and hydraulic conductivities range from 0.02 to 152 m/day. The Water Resource Agency also constructed 224 groundwater monitoring wells in CRAF and MB over different aquifer systems and depths. In addition, 33 multilayer compaction monitoring wells in central CRAF were constructed to detect land subsidence at different depths (Hung et al., 2017).
Because of a large number of evenly distributed wells and the stable recharge by Choushui River, CRAF and MB are two ideal places for experimenting with the methods of impulse flow and hydraulic tomography (Kennedy et al., 2017;Tsai et al., 2017). Furthermore, the shallow groundwater levels, large seasonal groundwater level variations, and abundant cross-well pumping test sites make CRAF and MB two ideal places for experimenting with the gravity-based S y determination.

Models of Groundwater Volume Change and Gravity Change for the Gravity-Based Method
The method of gravity-based S y estimation considers gravity change as a function of groundwater volume change in a porous medium. First, in an unconfined aquifer, the portion of a void pore in the unit soil texture can be expressed as where n is the porosity, S y is the specific yield, and S r is specific retention (Bear, 1979) and is also the portion of undrainable water during the decline of groundwater table. As shown in Figure 2, the groundwater level at time T 1 is Z 1 , changed to Z 2 at a later time T 2 . The groundwater level at T 1 is higher than T 2 . Assume that S y is a constant. The resulting groundwater volume change (Δv) is where v T1 − T2 is the total volume of the soil (including soil skeleton and groundwater in the medium). In reality, S y can vary within the underlying aquifer and hence is not a constant. The schematic is shown in Figure 2 with a gravimeter and different groundwater levels below the gravimeter. The blue-shaded region in Figure 2 corresponds to the actual change in groundwater volume from T 1 to T 2 . This volume change induces a gravity change that can be detected by the gravimeter. In the local planar coordinate system (x is to east, y is to north, and z is down), the gravity change (Δg) induced by the groundwater volume change is the sum of all vertical attraction components ( Figure 2):

Water Resources Research
where dv ¼ dxdydz : differential volume change G : gravitational constant (6.67 × 10 −11 m 3 kg −1 s −2 ) ρ w : the water density (¼ 1,000 kg/m 3 ) r : the distance between the gravimeter and the center of the differential volume x, y, z : the three coordinate components of the differential volume x 0 , y 0 , z 0 : the coordinate components of the gravimeter S y (x, y, z) : specific yield It can be challenging to map the three-dimensional function of S y (x, y, z) in Equation 3 considering aquifer heterogeneity and the high cost of conducting pumping tests over different depths, despite some successes reported in the literature Yeh & Liu, 2000;Zhu & Yeh, 2005). Furthermore, the groundwater level surfaces at T 1 and T 2 are functions of (x, y), which may be inferred using the electrical resistivity imaging method Rayner et al., 2007). As a result, data for S y (x, y, z) and groundwater level changes in Equation 3 are hard to obtain. This difficulty in acquiring data for Equation 3 may be alleviated by the fact that groundwater levels can be nearly flat or follow roughly the terrain in an unconfined aquifer. This allows the use of a cylinder (Figure 2, brown plate) to approximate the groundwater volume change between times T 1 and T 2.
To use the cylinder approximation, one must measure the levels Z 1 and Z 2 at T 1 and T 2 ( Figure 2) relating to the groundwater level functions in Equation 3. With a constant S y within this cylinder, the gravity change corresponding to this cylindrical volume of groundwater is (Heiskanen & Moritz, 1967) Figure 2. Real and approximate groundwater volume changes from time T 1 to T 2 , with a monitoring well and a gravimeter to detect gravity changes. The blue-shaded area shows the real aquifer volume change. The brown plate shows the cylindrical approximation of the volume change with a fixed radius for the plate (R p ) and with a thickness equal to the groundwater level change from depth Z 1 to Z 2 . The mesh-red plate shows the Bouguer plate when R p → ∞. Inserted on the upper left is a typical time series of water table including measurements Z 1 and Z 2 at T 1 and T 2 (T 2 is later than T 1 ). The arrow indicates the direction of groundwater movement from T 1 to T 2 .

Water Resources Research
where R p is the cylinder's radius, which in general is unknown because of a lack of knowledge in aquifer boundaries. This uncertainty in R p prompts the use of the so-called Bouguer plate model (R p → ∞, the mesh-red pate in Figure 2; Heiskanen & Moritz, 1967): In Equation 5, the cylinder's radius used in Equation 4 has vanished because of the assumption of R p → ∞. When applying the Bouguer plate model (called Bouguer model for short below) for modeling groundwater volume change, S y represents an average specific yield within the model area. This average S y from GBM is different from a S y determined by a cross-well pumping test, which produces well-dependent results that change with time and pumping test setups. Using the density of water and constants for G and π, we can simplify Equation 5 to obtain the following Bouguer model for S y estimation within the range of Z 1 to Z 2 : where Δg m is the gravity change from two successive gravity surveys associated with the groundwater level Z 1 and Z 2 between the two surveys. The Bouguer model is currently the benchmark for S y estimation and the only plausible method for GBM because the actual 3-D extent of a given aquifer can never be known. Therefore, we have used the Bouguer model (Equation 6) to estimate gravity-based S y values presented in section 5. In addition, Equation 6 suggests that GBM can produce repeatable S y values if the ratios between gravity changes and specific ranges |Z 2 − Z 1 | remain nearly constant, irrespective of groundwater level fluctuations from 1 year to another. Also, if we increase the frequency of gravity surveys at a given site, GBM can provide a vertical profile of S y values at the site . We will show these advantages of GBM in section 5.

Limits of the Bouguer Model in S y Estimation
The Bouguer model in Equation 5 is valid for S y estimation only if the lateral extent of a given aquifer is much larger than the groundwater depth (Pool, 2008). However, it is not possible to quantify an adequate ratio between the lateral extent and the depth in the Bouguer model (Equation 5) because this model contains no parameter for the lateral extent. On the other hand, the cylinder model in Equation 4 is a function of aquifer radius (R p , equivalent to lateral extent) and depths (Z 1 and Z 2 ) and thus can be used to quantify the ratio. Note that if the cylinder model fits the geometry of water volume change of an aquifer, then R p can be used to estimate the minimal areal (horizontal) coverage of the aquifer around a gravity site, which is πR 2 p .
Using the cylinder model, below we show formulae for computing the gravity difference and relative S y difference when the Bouguer model (Equations 5 and 6) is used for an aquifer with a finite lateral extent. These differences are based on an initial groundwater depth of Z 1 at T1, groundwater level change of Z 2 − Z 1 , and a radius of R p .

Gravity difference
The gravity difference is defined as the gravity change from Equation 5, minus the gravity change from Equation 4: 2. Relative S y difference The relative S y difference is obtained by first differencing S y values computed by Equations 4 and 5 and dividing the result by the S y value from Equation 5 (Δg C ¼Δg B ¼ observed gravity change): 10.1029/2020WR027261

Water Resources Research
The interpretations of D and P are as follows. If we use the Bouguer model to estimate a S y for an aquifer with a finite radius of R p , the resulting difference in gravity is D and the resulting relative difference in S y is P. In fact, D in Equation 7 is not the true gravity difference because the cylinder model still cannot truly model groundwater volume change. In addition, the relative S y difference (P) depends only on R p , Z 1 , and Z 2 and is independent of S y .
Figures 3a and 3b show how D and P vary with R p (1 to 1,000 m) and the initial depth Z 1 (10 to 100 m) under a fixed change of |Z 2 − Z 1 | ¼ 1 m and S y ¼ 0.2. Note that P is not affected by S y and one can obtain a new gravity difference by simply multiplying the difference in Figure 3a by the ratio between this new S y and 0.2. From Figures 3a and 3b, we can draw the following two general rules concerning the use of the Bouguer model: 1. For the same initial depth (Z 1 ) and groundwater level change |Z 2 − Z 1 |, the gravity model difference (Equation 7) and relative S y difference (Equation 8) increase with decreasing R p . 2. For the same R p and groundwater level change |Z 2 − Z 1 |, the gravity model difference and relative S y difference increase with the initial depth (Z 1 ).
The result from Figure 3b suggests that if the lateral extent (R p ) of a given aquifer is about 10 times larger than the initial depth (Z 1 ), the resulting relative S y difference is below 10% for |Z 2 − Z 1 | ¼ 1 m. As a numerical example for a shallow aquifer, consider the case in which Z 1 ¼ 10 m. In this case, R p must exceed 100 m to have a relative S y difference smaller than 10%. Another example is for a deep aquifer with Z 1 ¼ 100 m. In this case, R p must exceed 1,000 m to have a relative S y difference smaller than 10%. The mesh-green shaded parts in Figures 3a and 3b show the typical gravity differences and relative S y differences in the upper part of CRAF and MB, where the initial depths (Z 1 ) of the studied aquifers (Table 1) are smaller than 20 m. Because the lateral extents of these aquifers are far greater than 200 m, the use of the Bouguer model here will cause relative errors (P) that are smaller than 10%. In section 6.2, we will discuss the issue of lateral aquifer coverage using observed groundwater level changes at multiple wells around a gravity site. In the following ). The dashed line shows the 10% relative difference. These model differences are for six initial depths of groundwater level (Z 1 ) ranging from 10 to 100 m. We assume S y ¼ 0.2 in (a) and groundwater level change |Z 2 − Z 1 | ¼ 1 m in (a) and (b). The mesh-green shaded areas show the typical model differences in the upper part of CRAF and MB for Z 1 < 20 m.

10.1029/2020WR027261
Water Resources Research development, we will use the term "lateral extent-depth ratio" when discussing the validity of the Bouguer model. This term is defined as the ratio of R p to the mean of Z 1 and Z 2 .

Standard Errors of Gravity Differences and Gravity-Based S y
Each of the gravity-based S y values is associated with a standard error (the 68% confidence region), which is estimated as follows. The uncertainty of a gravity-based S y is caused by the groundwater model error and the errors in the gravity and groundwater level measurements. The total error of a gravity measurement is a combination of the errors in the gravimeter and environmental corrections. The groundwater model error is not considered in our result (section 5) but is discussed in section 6.2. Potential environmental error sources in absolute gravity measurements, such as those originating from seismic and volcanic activities, are shown in Van Camp et al. (2017), who also compared the sizes of such errors with the sizes of natural temporal gravity variations. Furthermore, Kao et al. (2017) listed all the sources of errors in the absolute gravity measurements collected in Taiwan over 2006-2016. In general, the total error of a measured gravity change (σ Δg m ) can be express as is the sum of the squared standard errors (error variances) from all sources at two successive epochs T 1 and T 2 (Figure 2), respectively. The typical standard error of σ Δg m using an FG5 gravimeter is 2.5 μgal and that of groundwater level change (σ Δh ) using a commercial piezometer is 1 cm. Applying error propagation to Equation 6, we can express the standard error of S y as where Δh is groundwater level change (Z 2 − Z 1 ). In Equation 10, the second error term due to σ Δh can be neglected because the precision of a commercial piezometer is normally better than 0.1% (relative precision), that is, σ Δh /Δh < 0.1% and the resulting error (the second term) is much smaller than that due to the gravity measurement error (the first term). Note that Equation 10 does not consider the error caused by the imperfect relationship between gravity change and groundwater water change. Both the Bouguer and cylinder models cannot perfectly describe this relationship, but the cylinder model can estimate the impacts of lateral extent and aquifer depth on an estimated specific yield.

Absolute Gravity Data Collection and Processing
Establishing a suitable site for absolute gravity measurement is much easier and cheaper than drilling a well or multiple wells needed for a cross-well pumping test. According to Equation 6, measurements of gravity and groundwater level change are both needed for estimating S y by GBM. To satisfy such measurement requirements, we inspected many existing wells maintained by the Water Resource Agency over CRAF and MB. We then decided whether a candidate gravity site near a given existing monitoring well should be established. For a good quality of gravity measurements, a permanent room with a power supply and free from ground vibration and large changes in humidity and temperature is preferred. In addition, we prefer candidate sites that are close to cross-well pumping test sites and close to existing monitoring wells where significant groundwater level changes (>1 m) have occurred between wet and dry seasons. Using these criteria and the information in Table 1, we established 10 gravity sites in our study area ( Figure 1 and Table 2). These gravity sites are mostly located at permanent, quiet rooms at elementary or junior high schools, and the gravity pillars are over solid foundations.
At these 10 gravity sites, we measured absolute gravity values in the wet and dry seasons that varied with the natural groundwater level changes from 2012 to 2017. One gravity survey session lasted 12 to 16 hr, depending on the performance of the gravimeter and other factors. A gravity survey session contains several gravity sets (values). A set of gravity value was obtained from 100 drops (measurements) over 0.5 hr, and all such set gravity values within a survey session were averaged to produce a mean gravity value on a date closest to the mean time of these set measurements. These gravity measurements were collected by two FG5 AGs 10.1029/2020WR027261

Water Resources Research
(serial nos. 224 and 231) in Taiwan. The measurement procedure and method for gravity uncertainty estimation are the same as those described in Chen et al. (2018) and Kao et al. (2017).
The groundwater level measurements at 1-hr intervals near a given gravity site were supplied by the Water Resource Agency after the gravity measurements at the site were collected (the latency can be up to 1 year). As such, we conducted gravity surveying without knowing the actual groundwater levels. As it happens, these latencies in gravity and groundwater data allow us to determine several S y values that represent different segments in an aquifer (see section 5) at a given gravity site. The water levels were measured by a commercial piezometer called GW200 that has a 0.1% relative precision, that is, 1-cm error over a 10-m range of water level variation. Because groundwater levels can vary with time, we averaged the hourly water levels to produce a daily mean that has a consistent definition with the mean gravity value measured by the FG5 gravimeters.
To obtain gravity changes for S y estimations, we removed the following nongroundwater gravity effects: (1) solid and ocean tides, polar motion, and atmospheric pressure change, and (2) land subsidence. In this paper, we used the ETGTAB solid earth tide model and the NAO.99JB ocean tide model. The polar motions were from International Earth Rotation and Reference system Service (IERS), and the gravity-atmospheric admittance was −0.35 μgal/hPa. These models and parameters were also used in our previous works (Chen et al., 2018;Hwang et al., 2009). Furthermore, we corrected the gravity effects due to significant vertical displacements at two gravity sites with continuous GPS stations (SJES and ELPS). We assume that 1-cm land subsidence induces a gravity change of 1.9 μgal (De Linage et al., 2007;Hwang et al., 2010;Van Camp et al., 2011). In addition, landslides and river debris deposits can also create large gravity changes in Taiwan (Mouyen et al., 2009(Mouyen et al., , 2014, but the 10 gravity sites are at locations free from such surface mass effects. Another source of gravity change not of groundwater origin (in a saturated zone) is soil moisture (in an unsaturated zone). To minimize this soil moisture effect in case of a rain event, we did not carry out a planned gravity survey immediately after the rain. Due to the budgetary constraint, in this study, we collected moisture data only at the gravity sites LHES, SJES, and JSES to only limited depths. In Appendix B, we show the gravity changes between two successive gravity surveys due to the soil-water-content (SWC) changes at these gravity sites and the expected uncertainties in estimating the soil moisture effects. The SWC-induced gravity changes are relatively small compared to the absolute gravity effects from SWCs. This can be explained by the difference between two successive gravity values expressed as (see also Equation 6 and Figure 2) where g 1 and g 2 are the raw gravity values and SM 1 and SM 2 are the gravity effects of SWCs at epochs T 1 and T 2 . The difference (SM 2 − SM 1 ) will be significantly smaller than the absolute value SM 1 or SM 2 , because SM 1 is similar to SM 2 ; see more details in Appendix B.

S y Values and Their Uncertainties From Gravity and Groundwater Measurements in Wet and Dry Seasons
A S y value is estimated using Equation 6 from the differences in gravity and groundwater level between two successive gravity surveys. Table A1 shows all the gravity changes, groundwater level changes, and the resulting S y values at the 10 gravity sites. There are 16 cases of rising water levels and 23 cases of falling levels from all successive gravity surveys, creating 39 gravity differences for S y estimations. Figure 4a shows all the S y values and the corresponding rises and falls of groundwater levels. The gravity-based S y values in Table A1

10.1029/2020WR027261
Water Resources Research range of depths between the two successive gravity surveys. For a given site, a S y value from a pair of gravity and groundwater level changes represent the mean value over the depth experiencing groundwater variations.
To maximize the signal-to-noise ratio and minimize the standard error of S y , we use the maximum groundwater level changes to estimate a representative S y for a given site. Figure 4b shows the representative S y values at the 10 sites and their standard errors (see section 5.2), along with the standard errors of the gravity changes. The representative S y value at a site can be compared with a given S y value from the hydraulic method at a nearby site if available. Table 2 shows the information of the representative S y values (by GBM) at the 10 sites and the S y values from the hydraulic method (Table 1; see also Figure 1a). As an example, there are six S y values from GBM at SLJH (Table A1). We select S y ¼ 0.14 ± 0.01 as the representative S y at SLJH because this S y is estimated from the largest groundwater level change (2.40 m). In Figure 4b, the representative S y values are derived from gravity changes between about 4 and 26 μgal and groundwater level changes between about 1 and 5 m. Figure 5 shows σ sy as a function of gravity uncertainties (σ Δg m Þ and groundwater level changes (Δh). The standard errors of the 39 gravity-based S y values (Table A1) are also shown in Figure 5 (circles). About 20 of the σ sy values are smaller than 0.03, as a result of using gravity differences with small gravity uncertainties and using large groundwater level changes. Figure 4b shows the σ sy values of the representative S y values at the 10 gravity sites. Except the relatively large σ sy ¼ 0.14 at CSES, the σ sy values of the representative S y values range from only 0.01 to 0.04. As shown in Chen et al. (2018), an investigation at a candidate gravity site should be made to ensure that the groundwater level changes at its nearby wells are sufficiently large to derive a gravity-based S y with a small standard error. Figure 5 suggests that the more the groundwater level difference, the higher the accuracy of the estimated S y value. For example, if the expected groundwater level change is 2.5 m between two successive gravity surveys and the acceptable σ sy is below 0.03, σ Δg m should be below 3 μgal.

The Repeatability of Gravity-Based S y
Because the groundwater levels varied over wide ranges and spanned over several years, we can examine the repeatability and reliability of S y values determined by GBM, as parts of our result. Equation 6 implies that if the aquifer is homogenous at a given site, then S y is a constant. This implies that the ratio between gravity  (Table A1). An arrow indicates the upward or downward movement of groundwater level between two successive gravity surveys that result in a S y . A site can have different S y values resulting from different depth variations of groundwater (for an arrow, the depth changes from one end to the other). (b) The best (representative) S y values with the smallest standard errors (the thickest arrows in (a)) from the maximum changes in gravity and groundwater level (Table A1). The vertical bars show the 68% confidence regions (one sigma) of the gravity changes.
change and groundwater level change should be a constant equal to S y / 0.0239, where S y is the specific yield (but subject to the condition that the Bouguer model is adequate, which is the case for the 10 gravity sites in this paper (see section 3.2). Figure 6 shows the scatter plots of gravity changes vs. groundwater level changes at SLJH, YLES, LHES, and ELPS, where the numbers of measured S y values from repeat gravity and groundwater level changes exceed 3. The squared correlation coefficient (R 2 ) at SLJH, YLES, LHES, and ELPS are 0.97, 0.66, 0.93, and 0.98, respectively. At YLES, the relatively low R 2 value of 0.66 is due to the negative S y value (−0.09) resulting from a negative gravity change of 4.6 μgal and a positive groundwater level change of 1.27 m in 2017 (Table A1). In fact, this anomalous S y value is the result of a low signalto-noise ratio in the gravity change measurement. The high R 2 values at SLJH, LHES, and ELPS show that GBM can produce coherent aquifer storage coefficients using gravity and groundwater measurements collected at different times.
Here we provide another assessment of GBM's repeatability using S y values determined from the rises and/or falls of groundwater levels over the same aquifer segments. The process of groundwater level rise and fall is similar to injecting into and pumping out of an aquifer. We expect that S y values from repeated injection or pumping tests in the same aquifer segment should be consistent. That is, S y should be reproduced (or repeatable) when several tests "revisit" the same part of the aquifer. Such repeatability criterion can be used to examine the adequacy of GBM in estimating S y . In other words, GBM should produce the same S y value irrespective of falling or rising groundwater levels or different times of gravity surveys, as long as the level changes occur over the same segment of the aquifer. Because our gravity experiment in CRAF and MB results in several S y values from GBM at each site, we can test this repeatability of GBM. However, we obtained the groundwater level data from the Water Resource Agency only after a latency of several months to a year. Thus, in practice, it is difficult to define exactly the same segment of an aquifer in a revisit test; this difficulty is true for both the gravity method and the hydraulic method. Using the groundwater levels coincident with gravity surveys, we determined the S y values corresponding to similar aquifer segments at each site in our revisit tests. This is explained by the criterion for selecting S y pairs (S y1 and S y2 ) at the lower right of Figure 7. This criterion is to ensure that a S y pair corresponds to a similar segment of an aquifer. As a result, we identify 10 qualified S y pairs at 3 sites.
In Figure 7, S y1 and S y2 are two independent gravity-based S y values at the same site. A diamond symbol shows a pair in which S y1 and S y2 are determined using groundwater levels moving in the same direction (both upward or downward). On the contrary, a cross symbol shows a pair in which S y1 and S y2 are determined using groundwater levels moving in the opposite directions. For example, the blue diamond at SJES in Figure 7 shows the pair S y1 ¼ 0.30 and S y2 ¼ 0.29 for which the groundwater levels both moved upward (26.38 to 28.21 m vs. 25.50 to 27.57 m). The green cross at LHES in Figure 7 shows the pair S y1 ¼ 0.08 and S y2 ¼ 0.09 for which the groundwater levels moved in the opposite directions (39.07 to 37.04 m vs. 37.04 to 39.98 m). In Table 3, we list the S y values for each selected pair and show the relative differences between S y1 and S y2 in the selected pairs. The relative differences range from 0.00 to 0.03, and the average difference is about 0.01. Table 3 indicates that the ratios between the signals (estimated S y values) and the standard errors of S y values are mostly larger than 3. In addition, we show the representative S y values from the gravimetric and hydraulic methods for comparisons at SLJH, SJES, and LHES. Note that the representative S y value (0.14) from GBM at LHES is not included in any of the selected S y pairs in Table 3, because the criterion for selecting a representative S y is different from that for the revisit test. Like Figure 6, Figure 7 and Table 3 suggest that GBM can almost reproduce the same S y value in the selected aquifer segment, irrespective of the times of the gravity measurements and the magnitudes of gravity changes. In Figure 5. The standard error of S y (σ sy , 68% confidence level) as a function of gravity difference standard error (σ Δg m ) and groundwater level change (Δh), with selected contours of S y standard errors. The standard errors (circles) of the 39 gravity-based S y values (Table A1) are at the positions of the circles.

10.1029/2020WR027261
Water Resources Research summary, the analyses using the results in Figures 6 and 7 indicate that GBM is a reliable method for measuring S y but subject to the condition that gravity errors are adequately modeled.

The Reliability of Gravity-Based S y
Using the information given in Table A1 and Figure 7, here we present three numerical examples at LHES, SLJH, and SJES to show the reliability of GBM for measuring S y . The first example is given at LHES, where we determined six S y values, 0.12, 0.08, 0.09, 0.07, 0.13, and 0.06, from seven pairs of gravity and groundwater level changes. The results show that the rising and falling groundwater levels at LHES, along with gravity changes, result in almost the same S y values at two segments of the aquifer below LHES. The S y values at shallower layers (34 to 37 m) are 0.12 and 0.13, but the S y values at deeper layers (37 to 40 m) are 0.08 and 0.09. Therefore, irrespective of the moving directions of the groundwater level, GBM produced repeatable (and hence reliable) S y values. This implies that GBM has the potential to determine depth-dependent S y values for a given site (Chen et al., 2019).  Table A1).

Water Resources Research
The second example is from the result of SLJH. We made seven gravity surveys to determine six S y values from 2015 to 2017 at SLJH. The groundwater level here varied from 171.02 to 173.60 m during our gravity surveys. The six S y values at SLJH are 0.10, 0.18, 0.09, 0.12, 0.10, and 0.14. In particular, the estimated S y values (0.09, 0.12, and 0.10) over the same depth range of 172.78 to 170.97 m are almost the same. These three repeated S y values are based on two falling groundwater levels and one rising level, which produce repeatable S y regardless of the moving directions of groundwater. Like LHES, at SLJH, the larger S y values of 0.14 to 0.18 occur in a shallower segment than the segment with the smaller S y values of 0.09 to 0.12.

The Effect of Land Subsidence on S y Estimation
The first subject of discussion is the effect of land subsidence on estimating S y using GBM. Land deformation can induce gravity changes (Torge, 1989), which should be removed before estimating S y using GBM. Among the 10 sites in this study, only SJES and ELPS experience significant land subsidence and are equipped with continuous GPS stations. The remaining gravity sites experience only small land subsidence rates, and they have no co-located GPS stations. SJES and ELPS are located in central CRAF, where subsidence is caused by massive groundwater pumping for irrigation and the rates of subsidence are about 2 cm/year. Because two successive gravity surveys spanned several months, we corrected for the land subsidence gravity effects using the vertical displacements from continuous GPS measurements.
The continuous GPS data at SJES show that the averaged vertical displacement rate is about −5 × 10 −3 cm/day in 2012 and −6 × 10 −3 cm/day in 2017. The cumulative subsidence here is about 14 cm from 2012 to 2017 (5 years). Furthermore, at ELPS, the averaged vertical displacement rate is about −6 × 10 −3 cm/day from 2015 to 2017. After removing the land subsidence-induced gravity changes (i.e., 1 cm of land subsidence is equal to 1.9 μgal of gravity change), we obtain S y ¼ 0.30 (before correction: Figure 7. Scatter plots of selected S y pairs (S y1 and S y2 ) determined from two different sets of gravity and water level changes. Inserted in the lower right is the criterion for selecting the S y pairs. A is the depth difference above the overlapped layer, and B is the depth difference below the overlapped layer. The summation of |A| and |B| must be smaller than 1.5 m. Gravity site SLJH has the most pairs, LHES has two, and SJES has only one. Water Resources Research 0.31) and 0.29 (before correction: 0.30) at SJES in 2012 and 2017, respectively. In addition, we obtained identically S y ¼ 0.10 at ELPS from two pairs of gravity and groundwater level changes in 2015 and 2017. Note that the two identical S y values (0.10) were obtained only after removing the land subsidence effects on the gravity measurements in 2015 (before correction: S y ¼ 0.12) and 2017 (before correction: S y ¼ 0.08). Thus, at ELPS, the corrections for land subsidence are about 20% of the corrected S y value (−0.02 or 0.02 relative to 0.10). At the two sites (SJES and ELPS) with land subsidence corrections, the removal of the land subsidence effects has improved the accuracies and consistencies of the estimated S y values.
It is known that S y will be changed if the void space in the corresponding aquifer is compressed. However, the estimated S y values at SJES and ELPS in different years are almost the same, suggesting that the S y values are not affected by compactions. Here we explain why this happens. We obtained compaction data at SJES and ELPS measured by multicompaction monitoring wells (Hung et al., 2017) near the two sites. Figure 8 shows the cumulative compactions at SJES and ELPS from 2007 to 2017. The major compaction at SJES occurs at depths from 125 to 300 m and that at ELPS from 37 to 300 m. The shallow-layer compactions (depths <15 m) contribute a subsidence rate of only 0.13 cm/year at SJES and only 0.07 cm/year at ELPS from 2007 to 2017. Because such shallow-layer compactions are small and the unconfined aquifers are over the shallow layers that experience the groundwater variations, land subsidence at deeper depths will not alter the water storage capacity here; that is, S y remains the same. Thus, the compaction data at ELPS and SJES suggest that S y values in an unconfined aquifer estimated by GBM at different times can remain unchanged despite the compactions occurring at the deeper layers of the aquifers.

The Validity of the Bouguer Model for S y Estimation
The theory in section 3.2 and the numerical examples in Figure 3b suggest that the lateral extent-depth ratio of an aquifer should be larger than 10 for the Bouguer model to be valid for S y estimation with a 10% error.
Here we use the representative S y values estimated by the Bouguer model (section 5.1 and Table 2) to show

10.1029/2020WR027261
Water Resources Research how relative S y differences vary with the approximate areal coverages implied by the cylinder model. Figure 9 shows the changes in S y for R p ¼ 5 to 500 m using Equation 4 with the gravity and water table data for the representative S y values at the 10 gravity sites. As expected, in Figure 9, S y values decrease with increasing R p and eventually converge to a stable value at each site. At the nine sites (other than LHES), the average depth of groundwater level is about 8 to 12 m. Thus, when R p exceeds 120 m at these sites, S y values experience little changes (the almost horizontal lines after R p ¼ 120 m in Figure 9). In contrast, the convergence distance of R p (the value beyond which S y becomes stable) is 200 m at LHES because here the average groundwater depth is 20 m. In these 10 examples, the convergences (P < 10%) in S y can be reached only after the lateral extent-depth ratios are greater than 10. In summary, if the real lateral extent of an aquifer is indeed small (say, 100 m), then the Bouguer model will lead to a relatively large error in the estimated S y . On the other hand, if the real lateral extent of an aquifer is indeed large (say, >500 m), the Bouguer model causes only a small error in S y .
Equation 3 shows the rigorous formula for determining S y value by GBM, but it is nearly impossible or can be costly to acquire a 3-D model of groundwater level changes needed in Equation 3. As stated above, since the lateral extent-depth ratio of the studied aquifers greater than 10, there is no significant difference (for P < 10%) in gravity-based S y between the uses of the Bouguer model and a complex groundwater model. As an example of showing groundwater-model-induced error in S y estimation, below we determine the difference between a S y value from 3-D groundwater changes (not employed at any sites) and the S y value from the Bouguer model at LHES.
The representative S y value at LHES from GBM is 0.12 (Table 2 and Figure 4b), which is derived from the measurements on 14 May 2015 and 9 November 2015. To establish a 3-D model of groundwater level changes around LHES, we collected measured groundwater levels at eight wells near LHES on these 2 days. Note that only groundwater level data at a well near LHES were used in gravity-based S y estimation in Table 2. Figure 10a shows contours of groundwater level changes between 14 May 2015 and 9 November 2015, interpolated from the daily groundwater level observations at the eight wells using Kriging. The model dimension in Figure 10a covers a 10-× 12-km area on a 100-× 100-m grid. The groundwater level changes range from 2.50 to 5.48 m in this area. Using Equation 3 and the level changes in Figure 10a, we obtained a , with the groundwater level change at the nearest well) results in a gravity change of 204.05 μgal. Thus, the gravity difference (0.6 μgal) is small compared to the gravity change. This example shows that there is no significant difference between the estimated S y values from the Bouguer model and that from the 3-D model of groundwater level changes. Figure 10b shows contours of the residual groundwater level changes, which are the differences between the original changes and a constant of 4.84 m. Here, 4.84 m is the groundwater level change at the nearest monitoring well, which is 0.1 km to LHES (the one located in the southwest of the gravimeter; see also Table 2). As shown in Figure 9, at LHES, the use of R p > 200 m results in the same S y value as that from the Bouguer model. This is because the gravity effect of groundwater level change decreases with squared distance; the residual level changes beyond 200 m (exclude the blue-shaded region in Figure 10b) create only small gravity effects on the gravity measurements at LHES (compared to the FG5 measurement accuracy of 2.5 μgal).
The results in Figure 10b indicate that the contribution of groundwater variations beyond 200 m is negligible and does not affect the S y value determined by GBM using the Bouguer model at LHES. This concept can be applied to any study site where the layer of groundwater level changes is a flat plate (Bouguer assumption) around a gravity site. In this study, the unconfined aquifers conform to this Bouguer assumption when determining S y , because the terrain over CRAF and MB is flat. Similarly, the gravity effect of soil moisture change close to a gravity site is more significant than those beyond R p . Note that this conclusion is under the condition that the studied aquifer (LHES in this example) has a radius of more than 200 m and the seasonal groundwater level variations are used for analysis. This conclusion is not applicable to the short-term groundwater level variations in a pumping test (Chen et al., 2018).

A Prospect of a Large-Scaled Estimation of S y Values Using Atomic Gravimetry
Our final discussion is on the impact of recent progress in atomic gravimetry on a large-scale gravity-based S y estimation. The precision of a terrestrial gravimeter limits the application of GBM to the estimation of aquifer storage coefficient. An AGlike FG5 uses repeat measurements to reduce measurement uncertainties. For example, in this study, we used FG5 gravity measurements collected in 12 to 16 hr to obtain an average absolute gravity value in a measuring session with a precision of several microgal levels. In addition to the errors from gravimeter measurements, groundwater level changes may also alter gravity values during the measurement times. The daily groundwater level fluctuations range from 0.05 to 0.2 m in the main recharge area of CRAF, which induce gravity changes of 0.4-1.7 μgal for S y ¼ 0.2. Thus, it is advised that the model error from groundwater level change should be considered if the gravity survey lasts more than 12 hr. As such, it is important that the time for acquiring a gravity value should be as short as possible to ensure that the gravity accuracy is at the microgal level.

Water Resources Research
This short-time gravity accuracy requirement may be met by a portable absolute atomic gravimeter. For example, Menoret et al. (2018) demonstrated a 1-μgal gravity accuracy from gravity measurements in a 1.5-hr session. An atomic gravimeter can even use the kinematic mode to acquire gravity values (Bidel et al., 2018;Wu et al., 2019). Such gravity measurement times are short enough to avoid the errors caused by groundwater movements in the saturated and unsaturated zones. The increasing portability in AGs can also provide more opportunity to acquire absolute gravity values in the field that are usually obtained by RGs. In short, the recent technological improvement in atomic gravimetry can greatly facilitate GBM to estimate the storage capacity of an aquifer system.
Here, we envision a scenario in which atomic absolute gravity measurements around selected, active irrigation wells over CRAF are used to estimate S y . If an irrigation well is actively in use, it should have provided a large amount of groundwater for irrigation and experienced large water table changes. Thus, any active well can be a candidate well. Figure 11a shows the distribution of 310,000 irrigation wells in CRAF. In the upper CRAF (the fan top), there are some 100,000 irrigation wells extracting groundwater in unconfined aquifers (Figure 1d) for irrigating a wide range of crops. The coordinates of every well have been measured by GPS, and these wells are mostly situated on the sides of paved roads with convenient accesses (readers may use the street maps of Google Earth to see these wells). The water tables of these wells can undergo large fluctuations between dry and wet seasons. If a microgal-level accuracy mobile AG is available (or a regular gravimeter needing only a short measurement time to achieve 1-μgal-level accuracy), one can use this gravimeter to collect the gravity value near a well in about 1 hr or so (with a gravity accuracy at 1 μgal) and use a level gauge to measure the water table (Figure 11). With the convenient road infrastructure in CRAF, one can finish 8 to 10 wells in a working day. In a month, gravity values and water tables at some 200 to 300 wells can be collected. This campaign of gravity and groundwater table measurement can be repeated in the wet/dry season (depending on the season of the first campaign). Using the method in section 3, we can estimate a large number of S y values with a relatively small budget compared to pumping tests, depending on the technological development of absolute gravimetry.

Conclusion
This study shows the results of measuring aquifer specific yields in two of the most aquifer-rich flood plains in central Taiwan using GBM. GBM requires only measurements of gravity changes and groundwater level changes around existing wells, without the need of constructing pumping wells and conducting pumping tests. This paper uses natural groundwater variations and the induced gravity changes to estimate multiple S y values for each of the 10 sites. Such S y values are repeatable and reliable, and represent the average aquifer storage capacities around the wells. Although land subsidence occurs around two of the gravity sites, it is properly removed and hence does not affect our result.
We use the cylinder model of groundwater level change to examine the validity of the Bouguer model for S y estimation for aquifers with finite lateral extents. The assessment result shows that one should be cautioned with the interpretation of a S y value from the Bouguer model, which can cause a relatively large error in an estimated S y value when the real aquifer radius is small. In the cases of CRAF and MB, the gravity model differences and relative S y differences are negligible because the lateral extents of the studied unconfined aquifers are more than several hundred meters counted from the gravity sites and the groundwater depths are mostly less than 20 m. In general, the gravity-based S y values from this study are consistent with those from the cross-well pumping tests.
The conventional, pumping-based hydraulic method took more than 25 years to obtain the 15 hydraulic S y values in the main recharge areas of CRAF and MB (all funded by Water Resource Agency and Central Geological Survey, Taiwan), owing mainly to the budget. In contrast, the GBM provides 39 S y values and their error estimates at the gravity 10 sites in only about 3 years (2015 to 2017). In fact, the use of 1 year of gravity and groundwater level measurements may result in more or less the same (close) 10 S y values. However, if one wishes to obtain a representative S y of an aquifer, one would need to acquire a priori knowledge about the times of the highest and lowest groundwater levels in a year to conduct two gravity surveys at the two times. In terms of spatial scale, such a representative S y value may better represent the storage capability than a S y value derived from a conventional cross-well pumping test. Because the current point density of S y values in CRAF and MB (and perhaps in most flood plains of the world) is still very low, it is expected that more resources and times are needed to obtain a sufficiently dense field of S y by pumping tests. GBM may serve the need of densifying S y values with a reduced budget and time if the technology of atomic gravimetry is fully developed. Table A1 shows all the gravity values, the gravity changes, and the groundwater level changes from two successive surveys, and the resulting S y values at the 10 gravity sites in this study (see). In total, 39 S y values are estimated.

Appendix B: Soil Moisture Gravity Effect and Discussion
This study focuses on measuring S y in an unconfined aquifer. Thus, whenever we have soil moisture measurements during gravity surveying, we removed gravity changes induced by soil moisture changes from the raw gravity data. The movements of water in the saturated and unsaturated zones are governed by different mechanisms. The groundwater movement in the saturated zone is a function of hydraulic gradient and aquifer parameters. Since the head dissipates quickly in the saturated aquifer, groundwater level changes measured at a well are the result of water balancing in the regional aquifer. In contrast, infiltrations in an unsaturated zone are controlled by soil tension and can result in different soil water contents (SWCs) at different depths. Moreover, a measurement of SWC by a sensor at a point only reflects the water content in a small and enclosed area near the sensor, especially in an urban area or a region with complex soil formations. Most soil moisture sensors are installed by hand-drilling at shallow depths, typically from 10 to 100 cm below the surface. Furthermore, because the upper CRAF and MB contain mostly gravels, it is difficult to place a sensor at a deep depth using hand-drilling instruments. As such, it is challenging to accurately measure water mass in the entire unsaturated zone needed in this study.

Water Resources Research
Because SWC can vary with depth, we measured SWCs at different depths. Because of the budget, we measured SWCs only at the gravity sites LHES, SJES, and JSES up to 1 m. As a result, there are only two to three SWC measurements at a given site ( Figure B1). For each of these sites, we averaged the SWCs over the different depths to obtain a representative daily SWC that is consistent with the daily mean gravity value and The SWC-induced gravity changes (Δg s ) are computed using where Δθ is the average change in SWC (%) and Δz is the depth range of the unsaturated zone, which in this study is the vertical distance from land surface to the lowest probe. The depths of the lowest probe at LHES, SJES, and JSES are 0.6, 0.6, and 1 m, respectively. The SWC-induced gravity changes from 13 June 2017 to 25 September 2017 about 0.01, 0.85, and 1.35 μgal at LHES, SJES, and JSES, respectively. The magnitudes of these SWC changes and the SWC-induced gravity changes are consistent with those based on continuous SWC measurements for correcting hydrological effects in SG measurements (Creutzfeldt et al., 2010;Kazama et al., 2012;Krause et al., 2009). Because of the study area (CRAF and MB) in a sub-tropical agricultural area, the land is productive throughout the entire year. As such, the soil moisture of the land is usually kept at a similar level for agricultural productions. Thus SM 1 and SM 2 at two successive epochs can be close, and the difference (SM 2 -SM 1 ) is much smaller than SM 1 or SM 2 . As a numerical example, if the gravity change induced by SWC is 1 μgal, the resulting error in S y is σ sy ¼ 0:02 over an area with 1-m groundwater level change based on Equation 10 and Figure 5. Further detailed discussions can be found in Pfeffer et al. (2011). Note that the uncertainty in SWC-induced gravity change can be very large (exceeding the expected uncertainty of an FG-5 measurement, about 2.5 μgal) when estimating such an effect using only sparse point measurements. In addition, SWC values from the deepest probe (0.6 or 1 m, in this paper) to the surface of the saturated zone can be hard to estimate and such SWC values can also cause a large uncertainty in SWC-induced gravity change. The SWC measurements collected at the three sites are point measurements outside of the buildings housing the AGs. These buildings may create the "umbrella effect" that result in different scenarios of soil moisture distributions around and below the buildings. In addition, there may be no SWC change right below a building as a result of the umbrella effect (Reich et al., 2019). Because of the shielding of the gravimeter buildings, the SWC-induced gravity changes can be much smaller than that computed from Equation B1 (Creutzfeldt et al., 2008;Kazama et al., 2012). In short, accurately modeling SWC-induced gravity changes remains a challenging issue when using gravity measurements for S y estimation.