Delineating a Volcanic Aquifer Using Groundwater-induced Gravity Changes in the Tatun Volcano Group, Taiwan

The Tatun Volcanic Group (TVG) is an active volcano that could cause volcanic hazards in northern Taiwan. The latest phreatic eruption of the TVG occurred some 6000 years ago. Understanding the state of groundwater around the TVG can be a crucial step towards effectively assessing the risk of phreatic explosion by providing information about the sources of groundwater and the media it ows. We measured gravity changes at a superconducting gravity station and several groundwater-sensitive sites to examine the way the groundwater altered the gravity values around the TVG. Groundwater-induced gravity changes are simulated by two hydrological models (A and B). Both models show coherent seasonal variations in groundwater level and gravity value in the center of the TVG (Chintiengang). This coherence indicates inter-connected porous media for free groundwater ows below Chintiengang. However, inconsistencies between the modeled and observed gravity changes occurred in the eastern part of the TVG, suggesting here highly heterogeneous formations with fractures and barriers may exist below Chihsinshan and Dayoukeng. The gravity consistencies and inconsistencies between the observations and the models are used to delineate a volcanic aquifer, which can provide additional information for assessing the probability of a potential phreatic eruption over the TVG. outow (recharge and discharge) in a unit volume is equal to the ow change in a unit time. Rrecharges and discharges in Model A are constrained by the groundwater levels at the monitoring wells. At the junction of land and river, the groundwater levels are equal to the river levels with a steady state. The method of minimum curvature (Wessel et al., 2013) was used to interpolate the observed groundwater levels and river levels onto a 40 m The gridded groundwater levels were used to determine the gravity changes at the gravity sites. for the is bounded by the surface of the gridded levels to the homogeneous


Introduction
The Tatun Volcano Group (TVG) is located in northern Taiwan and originated from the active plate convergence in the western Paci c ring of res. Because the TVG neighbors Taipei City and New Taipei City, its eruption will threaten the lives of about 7 million people in the two cities, and potentially damage the two nuclear power plants at the northern anks of the TVG, creating major economic and social crises in Taiwan (Konstantinou, 2015). The latest phreatic eruptions of the TVG occurred about 6000 years ago (Belousov et al., 2010). Recent studies have shown that the TVG is a potentially active volcano (Lin et al., 2005b;Murase et al., 2014). Data from a new seismic network  over the TVG showed repeated occurrences of local seismic swarms associated with uid migrations. In addition, Lin (2016Lin ( , 2020 detected magma chambers below the TVG using the seismic S-wave shadows and P-waves delay. Hydrothermal reservoirs in the eastern limb of the TVG (Fig. 1a) have been detected using volcanoearthquake signals (Lin et al., 2005b), precise leveling survey (Murase et al., 2014), audiomagnetotellurics (AMT; Komori et al., 2014). The variations of fumarolic gas compositions indicated highly hydrothermal activity below the TVG (Lee et al., 2008).
Gravity observations have been used to assess the states of volcanoes around the world. For example, Rymer and Brown (1989) used time-lapse gravity changes as a precursor, on the ground that ascending magma before an eruption can lead to mass changes. Using gravimetric and geodetic measurements, Battaglia et al. (2006) showed that the uid migration in the Campi Flegrei caldera hydrothermal system was the cause of ground deformation and geological unrest. Kazama et al. (2015) showed that gravity changes can originate from both magma and non-magma sources. Depending on the location, the largest contributor of gravity change around a volcano can be hydrological variation, which can overwhelm magma-induced gravity changes around a volcano by many orders of magnitude. Using the absolute gravity measurements in 2004-2007 around the TVG, Mouyen et al. (2016) suggested the existence of a tube-like structure that allows a large-scaled hydrothermal uid circulation below the TVG.
A gravity-based study of the TVG began with the installation of an AG site at a continuous GPS station, called YMSG, in the TVG in 2004 (Kao et al., 2017). In 2012, a superconducting gravimeter (SG, serial number T49) was installed at YMSG. Prior to 2012, 4 new absolute gravity sites and 27 relative gravity sites along 5 hiking trails of the TVG were constructed (Fig. 1a). In 2012, four gravity surveys were carried out and the surveys covered the approximate volcanic regime of the TVG, including Chihsinshan, Chintiengang, Dayoukeng and Dingshan in an area of about 10 km 2 (Fig. 1a). In the gravity surveys, one absolute gravimeter (AG) and two relative gravimeters (RG) were also used. In contrast to the continuous, one-second measurements by the SG/T49 meter, the gravity measurements by the AG and RG meters are time-lapsed with a 3-month measuring interval.
In most gravity studies of volcanoes, gravity effects originating from hydrological changes are estimated by models that do not consider 3-D ows of groundwater, which are largely affected by local hydrogeological settings. In principle, a 3-D hydrological model of an aquifer below a volcano is able to show how groundwater transports, as in an alluvial plain (Singhal and Gupta, 2010). However, over a volcano the hydraulic conductivities of volcanic materials can vary widely, at a range of nine orders of magnitude. For example, the hydraulic conductivities of lava ows and cinder beds are high and those of ash beds, intrusive dikes and sills are much lower (Fetter, 2001). Thus, validating a 3-D model for groundwater ows over a volcanic region can be challenging.
Gravity changes can be a type of data for validating a 3-D groundwater ow model because a gravimeter can sense mass changes due to groundwater re-distributions from the model. An example in Taiwan is the investigation of the active state of the Hsinchu Fault using gravity changes from a superconducting gravimeter (Lien et al., 2014). In the TVG, the sources of gravity changes are seasonal hydrological changes, transitional uid migrations (Mouyen et al., 2016) and secular plate motions (Kao et al., 2017).
In typical solid earth applications of gravimetry, gravity changes of hydrological origins (called the hydrological gravity effect) are regarded as noises and are removed from the raw gravity measurements before solid-earth applications. At a gravity site, the hydrological gravity effect is the sum of the effects over the unsaturated zone and the saturated zone below the site. Typically, the effect from the unsaturated zone is determined and removed using soil moisture measurements and the methods for this effect can vary widely, from the simple Bouguer plate model (Torge, 1989) to a sophisticated approach such as those presented by Crossley et al. (1998) and Kazama et al. (2015). Likewise, the gravity effect from the saturated zone can be modeled by a simple model such as Bouguer model (Torge, 1989) or a 3-D model such as the method for the Hsinchu superconducting gravity station (Lien et al., 2014).
In this study, the hydrological gravity effects from the unsaturated zone in the TVG will be still regarded as noises and removed. However, the effects from the saturated zone will be modeled by two approaches that are more comprehensive in 3-D extent than the Bouguer plate model. Such modeled gravity effects are then compared with the observed gravity changes from the four gravity surveys in 2012 to infer the potential boundaries of an aquifer around the TVG. In the following development, we will describe in detail how our gravity measurements were collected and processed, and how two groundwater ow models are constructed to model gravity changes for inferring the aquifer boundaries. Because groundwater ow is important for understanding volcanic unrests (Jasim et al., 2018), this result from this study will bene t the understanding of the aquifer distribution in the TVG for volcanic hazard modeling in northern Taiwan.

Gravity data collections and corrections
In this study, the time-lapse gravity measurements were collected by one FG5 absolute gravimeter (AG) (serial number 231) and two CG5 relative gravimeters (RG) (serial numbers 050800136 and 050800137) in April, June, September, and November, 2012. These gravity measurements in the four seasons were expected to sense the TVG hydrological changes in the wet and dry seasons. The gravity sites were installed along ve hiking trails, covering the region over 25.13°N − 25.19°N and 121.525°E-121.65°E. The purpose of the AG measurements is twofold: (1) obtaining gravity changes associated with groundwater level changes, and (2) constraining the RG measurements in the relative gravity adjustments (see below). The RG measurements along a trail were collected in a double-run (back and forth) survey to eliminate drifts in the relative gravimeters and to enhance the point gravity accuracies. The 4 AG sites (YAG1-3 and YMSG) are co-located with 4 RG sites (Fig. 1b).
Following previous experiences made in studies using high-precision relative gravimetry, at each RG site and for each survey we determined a mean reading every 90 seconds using a sampling frequency of 6 Hz for the CG5 gravimeters (Debeglia and Dupont, 2002; Scintrex Limited, 2010; Mouyen et al., 2013). In total, 15 mean gravity and temperature readings were acquired at one site occupation. On average, it took about 30 minutes for such measurements at a site. From the 15 readings, we selected 8 most stable readings using the following criteria: (1) neglecting the rst two readings, (2) choosing continuous readings in which the gravity changes of the last 3 readings are less than 5 µGal, and the trend of the continuous readings is less than 5 µGal/reading, (3) selecting the readings when the temperature variations of the sensor chamber of CG5 are within ± 7 mK, and the temperature changes between two successive readings are less than 0.1 mK, and (4) selecting the readings when the ranges of tilt variations in X and Y directions are within ± 20". Note that, the tilt-induced gravity errors can be neglected when the X and Y tilt variations are less than ± 3" (Merlet et al., 2008). Criterion 3 is to use only readings when the chamber temperature is stabilized to cause little variation in reading (Fores et al., 2017).
The temporal gravity effects due to solid earth tide in all gravity measurements were removed using the model values of ETGAB (Wenzel, 2002;Dehant et al., 1999). The Newtonian effect and the loading effect of ocean tides (called OTL correction) were removed using the FES2004 tide model (Letellier et al., 2004) and the SGOTL software (Hwang and Huang, 2012), which considers the elevation of a gravity site. We also removed the effects of barometric pressure and polar motions using a SG-based gravity-atmosphere admittance coe cient (see below) and the earth rotation parameters from the International Earth Rotation Service and Reference Systems. Finally, the gravity effects at the unsaturated zone were removed using the soil moisture measurements at YMSG and the Bouguer plate method given in Torge (1989).

Adjustments of relative gravity measurements
For each of the four gravity surveys, the RG and AG measurements were least-squares adjusted using the weighted constraint method of Hwang et al. (2002) by holding xed the gravity values at the ve AG gravity sites. The outliers in the RG measurements were detected and then deleted in the adjustments. The histograms of the residuals of the RG gravity measurements survey follow the normal distribution, suggesting that there were no systematic errors in the adjusted gravity values.
To quantify our RG data quality, two separate adjustments of the RG measurements were made: one adjustment used RG measurements based on all 15 readings at one site occupation, and the other based on the selected 8 readings (Sect. 2.1). The adjustment results are shown in Table 1, which shows that the standard deviations (STDs) of the gravity values in the case of using the selected 8 readings are smaller than those using all 15 readings. The accuracy improvements due to the use of the 8 selected reading for the four surveys are 11.1%, 53.8%, 10.0%, and 20.0%, respectively. The reductions in standard errors suggest that our data selection strategy based on the ve criteria (Sect. 2.1) is effective. Table 1 Standard errors of point gravity values (in µGal) from relative gravity adjustments using all (15)  In addition, we found that applying the OTL corrections using SGOTL has improved the point gravity accuracies in the adjustments by 1.4% to 6.2/%, compared to the cases without the OTL corrections. The improved accuracies are most pronounced at the sites with elevations > 1000 m, suggesting that the altitude of a gravity site should be considered when computing the OTL effect. This conclusion is consistent with that made in Hwang and Huang (2012).
Another contributor to the overall improved RG point gravity accuracy is the SG-based admittance coe cient for atmospheric pressure corrections. We experimented with two gravity-atmosphere admittance coe cients: one is the standard value of -0.30 µGal/hPa (Torge, 1989) and the other is -0.35 µGal/hPa from the analysis of the SG49 gravity and barometric records (see Fig. 1 for the SG site, YMSG). Table 2 shows the averaged standard errors of the adjusted point gravity values. In the second survey, the difference between the use of the two coe cients ( -0.30 and − 0.35) is 1 µGal.
This is because the gravity-pressure admittance coe cient under large atmospheric pressure changes can be signi cantly different from the standard admittance coe cient (Hwang et al., 2009).
For example, during the second survey, the air pressures around Taiwan were affected by Typhoon Doksuri on June 28 and 29, 2012. In this case, the use of an admittance coe cient of -0.35 µGal/hPa improves the point gravity accuracy in the second survey. Hence, SG measurements can be used to determine an optimal gravity-pressure admittance during a eld work to calibrate the atmospheric pressure corrections. Although the SG-based gravity-pressure admittance coe cients improved the gravity accuracy only in the second survey, we believe such coe cients can improve atmospheric pressure corrections in future gravity data collected in the TVG. 4th 8 8 Figure 2 shows the gravity changes at all the gravity sites in the TVG in April, July, September and December, 2012. A gravity change is the difference between an observed gravity value and the mean of the four observations from the four surveys. The gravity changes in Fig. 2 will be used to infer potential aquifer boundaries around the TVG using Models A and B (see below).
Gravity changes in Fig. 2 were correlated with rainfalls. As an example, Fig. 3 shows the gravity changes at sites Y01, Y02 and Y03 (see Fig. 1a), rainfalls and in ltrations at two weather stations WZS and SD.
The rainfalls (red and blue bars) in Fig. 3 were converted to in ltrations (red and blue lines) using an in ltration coe cient method (Patil et al., 2019; Tarka et al., 2017) as follows: (1) multiplying rainfall with an in ltration coe cient to obtain initial in ltration, and (2) determining the nal in ltration by moving average with a time lag. The observed gravity changes at Y01 to Y03 are coherent with the variations in in ltration. This coherence is due to the fact that in ltrated rainfall recharges groundwater to create gravity changes as sensed by the gravimeters. The gravity change in September was negative, indicating groundwater discharges were more than rainfall recharges in the dry season. In general, at Y01, Y02 and Y03 the gravity changes correspond well to the groundwater level uctuations.

Groundwater level, rain and soil moisture data
The area covered by our two hydrological models is 10 km10 km and the center is at YMSG, as shown in Fig. 1a. One continuous groundwater level (TB23, Fig. 1a) automatically recorded water levels at 15minute intervals. TB23 is located 500 m southwest of YMSG. In addition, from April to December 2012 the water levels at 22 wells were manually recorded once a month. At YMSG, a rain gauge and soil moisture sensors were installed. The sensors recorded soil moistures at the depths of 10, 30, 50, 70, and 100 cm every 15 minutes. The Central Weather Bureau (CWB) and Tatun Rain Gauge Network installed additional 13 rain gauge stations that provided data for this study.

Model A: a continuous ow model
We assume that the formation under the TVG area contains porous media for groundwater ows, which introduce gravity changes that can be used to infer the borders of the media. Our inference employs two models, namely Models A and B. Model A constructs a layer of groundwater-level changes by interpolations from the observed groundwater levels at the 23 monitoring wells and the river levels (the colored stars and lines in Fig. 4a). Figure 4b shows the gridded water levels in April 2012. A water level is the elevation relative to the mean sea level at Keelung, which is the origin (zero point) of elevation in Taiwan. Then, the groundwater-level changes were converted to gravity changes using the method of volumetric integration described in Lien et al. (2014). Note that the gravity contribution of the unsaturated zone has been removed from the gravity observations using the method (Sect. 2.1) and hence is not considered in the volumetric integration.
In Model A, we assume that the sum of in ow and out ow (recharge and discharge) in a unit volume is equal to the ow change in a unit time. Rrecharges and discharges in Model A are constrained by the groundwater levels at the monitoring wells. At the junction of land and river, the groundwater levels are equal to the river levels with a steady state. The method of minimum curvature (Wessel et al., 2013) was used to interpolate the observed groundwater levels and river levels onto a 40 m × 40 m grid. The gridded groundwater levels were used to determine the gravity changes at the gravity sites. The vertical domain for the gravity computation is bounded by the surface of the gridded groundwater levels (top) and the surface corresponding to the lowest water level of the TVG (bottom). We assume that each vertical column of a grid is isotropic and homogeneous with a constant speci c yield. The water storage capacity of a water-bearing body is proportional to , and indicate the ease of groundwater movement. Proper values of these hydrogeological parameters will result in modeled groundwater levels that match the observed levels, which serve as the boundary conditions (values) for the partial differential equation in Eq. 1. Table 3 shows the parameter settings for Model B in the TVG. Figure 5a shows the concept for Model B, which considers one aquifer that is bounded by the top and bottom boundaries at the ground surface and at the sea level, respectively, and also bounded by the two faults (Shanchiao Fault and Kanchiao Fault) to the aquifer's western and eastern sides. The boundary conditions for hydraulic values are the same as those for Model A, i.e., groundwater levels at monitoring wells and rivers. In addition, the rainfalls at the 14 stations (Fig. 1a) are fed to Model B as recharges (Fig. 5a).  As shown in Table 3, the number of hydraulic conductivity zones is set to 7 according to the hydrogeology map from the Central Geological Survey of Taiwan (Fig. 1b). The hydraulic conductivities, speci c yields and porosities begin with some initial values and are nalized when the modeled groundwater levels from Model B best t the groundwater observations (boundary conditions) in the least-squares sense. Finally, Model B uses the same integration method as that for Model A to compute the gravity changes due to groundwater level changes. The two faults in Fig. 5a  In a mountainous area like the TVG, the variations in groundwater level and the hydraulic head gradients can be relatively large compared to those in an alluvial fan in western Taiwan. However, in the TVG there are only few borehole wells that provide the needed hydrogeological parameters to address such large variations. Hence, in Model B we assume one layer of materials covering the unsaturated and saturated zones from the ground level to the sea level to avoid discontinuities in water level at the 40 m×40 m grids. The bottom of the layer is at the sea level and is a no ow setting that functions as a low permeability boundary. To estimate recharges in the TVG, the study area was divided into 14 sub-areas according to 14 geographic distributions of the rainfall stations ( Fig. 1a and Fig. 5a) using the Thiessen polygons method. Each sub-area was assigned with several runoff coe cients. The runoff is the product of rainfall and runoff coe cient. The runoff coe cients in one sub-area were given according to mountain slopes and aspects (windward and leeward sides of the TVG). The recharge of groundwater in each unit area is the difference between rainfall and runoff.

Model B: a dynamic hydrological model
The boundary conditions from the groundwater wells and river levels constrain Model B in a different way than that for Model A. In Model B, the river levels were held xed, while the observed well levels were used to validate the modeled water levels. The geological formations of the TVG are composed of lava ows, tuff breccia, sandstone and shale (Fig. 1b). The lava ows from the latest eruption are composed of andesite on the top layer, followed by tuff breccias, Mushan formation, and Wuchishan formation. Mushan formation is composed of alternations of sandstone, shale and intercalated coal seams. The compositions of Wuchishan formation are moderate to coarse gravels and sandstone interbedding with thin carbonaceous shales. The initial values for hydraulic conductivity, speci c storage, speci c yield, and porosity were from a geological map released by the Central Geological Survey of Taiwan (http://www.moeacgs.gov.tw/; a modi ed version is shown in Fig. 1b).
The hydraulic gradients can be high in a mountainous area like the TVG because of steep slopes. We assume that the vertical ow of groundwater is higher than the horizontal one. Hence the initial vertical hydraulic conductivity of volcanic rock is 10 times more than the horizontal one. Then, we adjusted these parameters by comparing the modeled groundwater levels with the observations at the wells. The adjustments were repeated until the differences between the modeled and the observed groundwater levels are less than 0.5 m for all wells. Note that this tolerable difference of 0.5 m is relatively large because the TVG is over a rugged terrain with higher elevations, and the grid size is a relatively large 100 m in Model B. If the grid size is smaller than 100 m and/or the terrain is at, the tolerable difference should be smaller than 0.5 m.

Result Of The Modeled Groundwater Level Changes And Gravity Changes
The results of Models A and B show seasonal groundwater level variations, which are largely coherent with the variations in rainfall. As examples, Fig. 4b and Fig. 5b show the modeled groundwater levels from Models A and B in April 2012. In general, Model A produces very localized groundwater variations, while Model B leads to regional groundwater variations. According to the CWB precipitation data in 2012 (Appendix A), the monsoons in May and June resulted in heavy rainfalls over the TVG, where it was dry from October through November. Figures 6 and 7 show the groundwater level changes relative to the mean levels from April to December 2012. Figure 6 (Model A) shows a dry pattern (negative water level changes) from September to November, and a wet pattern (positive water level changes) from April to June. The seasonal pattern is consistent with the rainfall pattern (see Appendix A). Due to a lack of monthly river level data, a constant river level for a given month was used in Models A and B. Hence, the monitoring wells data dominate the modeled groundwater level changes in Model A. The in uence area of one monitoring well is local and depends on the method of interpolation. Figure 6 shows that the modeled groundwater levels are constant in the eastern TVG because there are no continuous water level data here, with only constant river boundaries.
Model B uses the governing groundwater ow equation (Eq. 1) to produce regional groundwater levels (hydraulic heads). faults is the watershed of the rivers in the TVG, where groundwater is mainly recharged by rains and stream ows (Fig. 1a). Figure 5b and Fig. 7 show sharp differences in the groundwater levels at the two sides of each of the faults.
Both Model A and Model B result in seasonal gravity changes corresponding to the modeled groundwater variations. Figure 8 and Fig. 9 show the gravity changes induced by the modeled groundwater variations during the times of the four gravity surveys. Figure 8 shows that Model A-derived gravity changes in the eastern TVG are small (few µGal). Near Dayoukeng and Dingshan (see Fig. 1a for their locations), the groundwater variations are nearly zero, but not gravity changes. This suggests that the gravity changes at these two points were created by the gravitational effects of water level changes at locations not nearby Dayoukeng and Dingshan. Figure 9 shows distinct patterns of gravity changes at Chishinshan (Y20), Chintiengang (YMSG), Dayoukeng (Y08) and Dingshan (Y03), which will be discussed in Sect. 5.

Gravity changes associated with seasonal groundwater changes and hydrothermal uid migrations
The gravity changes from the observations and from the two models are mostly seasonal with the extreme values occurring in June and October (see Appendix B for the time series of gravity changes). We carried out an empirical orthogonal function (EOF) analysis to identify the leading variabilities of the gravity changes in space and time. Figure 10 shows the time variations (left column) of the rst three modes from the observed gravity changes and the spatial distributions of the rst-mode variabilities (right column; see also Appendix C for the spatial distributions of the second-and third-mode variabilities). The rst three modes of EOF explain 100% of the variabilities of the gravity changes. For both the observed and modeled gravity changes, the rst modes explain up to 78.6% of the total variabilities. Because the observed gravity changes have been corrected for temporal gravity changes (see Sect. 2.1) and are largely affected by the local hydrology of the TVG, the rst-mode variabilities should be of hydrological origins. Because groundwater level data were used to constrain Models A and B, the rst-mode variabilities from the two models are largely due to groundwater level changes. The second-mode variabilities of the observed gravity changes explain about 20 % of the total variabilities, followed by the third-mode variabilities.
The EOF analysis of the gravity changes enhances the spatial patterns of the gravity variabilities. From the gravity observations, the rst-mode variability near Chishinshan is the largest (near 1) and the variability near Dayoukeng is the smallest (near − 1). The variabilities at these two locations from Models A and B are different from the variabilities from the gravity observations. This suggests that Models A and B cannot adequately model the gravity changes here. These inconsistencies in the rst-mode variabilities between the observed and modeled gravity changes over Chishinshan and Dayoukeng suggest that the hydrogeological properties here are far more complicated than the current settings used in Models A and B, which cannot produce gravity changes that match the observations. Figure 11 shows the differences between the observed and modeled gravity changes at the gravity sites. Four gravity sites are located at Dayoukeng, namely, Y06, Y07, Y08 and Y09, where relatively large positive gravity changes were observed in September (see Appendix B). However, there was little rainfall in September (Appendix A). Because the modeled gravity changes are based purely on rainfalls, the positive gravity differences at Dayoukeng in September (Fig. 11) were probably caused by the mass surplus not predicted by the model. This mass surplus may be related to the hydrothermal uid beneath the TVG, which cannot be detected in our hydrological models. Evidences for such hydrothermal uid were described in Murase  The differences in Fig. 11 suggests that, in September the uid migrated from a deep hydrothermal reservoir under Chishinshan to a shallow reservoir under Dayoukeng (Fig. 1a and Fig. 12). On the other hand, in April the uid below Dayoukeng migrated to the deep hydrothermal reservoir, causing negative differences (mass de ciency) between the observed and the modeled gravity changes (April, Fig. 11).

The aquifer boundaries from gravity consistency and inconsistency
Because the predicted gravity changes from Models A and B are largely due to hydrological effects, the discrepancies between the modeled and the observed gravity changes indicate that (1) the hydrogeological parameters for the uid ow media in two models are imperfect at the sites of large discrepancies, (2) additional sources in uencing the temporal gravity corrections for the raw gravity observations (Sect. 2.1) should be considered. Figure 11 shows that at the gravity sites around Chintiengang, i.e., YMSG, Y05, Y10, Y12, Y13, Y14, Y15, and Y24, the gravity differences are small (few µGal), suggesting that both Models A and B perform adequately in predicting gravity changes with the current settings in the models. On the other hand, at the sites near Chishinshan (Y17, Y18, Y19, Y20, Y21, and Y23) and Dayoukeng (Y06, Y07, Y08, and Y09), the discrepancies between the observed and the modeled gravity changes are much larger, indicating that the two hydrological models are inadequate here.
Model A is constrained by the groundwater level observations at the 23 wells and at the rivers without using the groundwater owing dynamics (Eq. 1). The modeled groundwater tables are predicted from interpolations to the observed groundwater levels. The groundwater layer is assumed to be homogeneous and contains porous media without low permeable/boundary layer in vertical or horizontal directions. Thus, the sites where the gravity changes from Model A and the observations are consistent should be located in a region with porous or permeable formations without groundwater ow barriers.
Model B uses the continuous groundwater ow equations (Eq. 1) with inputs from rainfall recharges and optimized hydraulic conductivities, speci c yields, and porosities (Table 3) to predict groundwater level changes. With these settings, Model B assumes that the underlying media below the study area (Fig. 1b) are porous and resemble an alluvium formation. As shown in Fig. 11, this assumption and parameter settings are valid only over the gravity sites where the observed and modeled gravity changes are consistent. This consistency occurs largely at the gravity sites over Chingtiengang, where the formation contains inter-connected porous media for groundwater to ow as in an alluvial aquifer. Th extent of the inter-connection is determined by the hydraulic conductivities in Table 3, which range from 0.0005 to 0.02 m/day. Moreover, Model B assumes that there is a low permeable layer under Chishinshan and Chintiengang, following the concept of Komori et al. (2014). This assumption also contributes to the consistency between the observed and modeled gravity changes over Chingtiengang.
Using the gravity consistencies and inconsistencies found in this study and the ndings in previous studies, we present a possible con guration of a hydrogeological scheme along pro le A-B (see Fig. 1b) in Fig. 12. The aquifuge in Fig. 12 is based on Komori et al. (2014) and it prevents hydrothermal uids below the sea level from owing to the formation below Chintiengang. This aquifuge also does not allow the groundwater below Chingtiengang to leak into the deep hydrothermal reservoir identi ed by Murase et al. (2014).
As mentioned above, Model A and Model B are expected to simulate groundwater level variations over strata containing continuous media without boundaries (low permeable materials). The observed changes in water levels are the boundary values for the two hydrological models. However, if the underlying strata contain discrete cracks, the observed groundwater levels will be in favor of a scenario in which groundwater level variations correspond to a very local fractured block, instead of a regional one.
That is, a large difference between the observed and modeled gravity change at a site implies that the materials below the site may contain volcanic conduits or inter-connected cracks, rather than a continuous porous formation. Figure 11 show that, at most of the gravity sites over Chishinshan, the differences (red circles) between the observed and modeled gravity changes (observation minus model) are positive in the rst (April), second (July) and fourth (November) gravity surveys. This implies that groundwater levels over Chishinshan are underestimated. Because the hydraulic conductivities and storage coe cients in Model B are assumed to be constants along the vertical direction at each grid, the underestimated groundwater levels suggest that the hydraulic conductivities and storage coe cients in Table 3 are not suitable for Chishinshan, but are still optimal over Chingtiengang. Another interpretation is that the governing equation of MODFLOW (Eq. 1) is simply not realistic for the formations over Chishinshan, which may contain discrete volcanic conduits and fracture materials. As such, MODFLOW cannot produce reasonable groundwater ows here.
The differences between the observed and modeled gravity changes could be also caused by different recording times of the gravity and groundwater level measurements. For example, the negative gravity differences in September 2012 around Chishinshan (Fig. 11) might be caused by a sudden rise of groundwater level due to the heavy rainfall brought by typhoon Jelawat in the TVG region. This sudden water level rise was recorded in our well data and led to large gravity changes from Model A and B. However, the gravity survey in September 2012 was carried out before typhoon Jelawat, so that the Jelawat-induced water mass surplus was not re ected in the gravity observations, leading to the large negative gravity differences (observed minus modeled gravity changes) around Chishinshan in Fig. 12.
As a nal remark, the Model B result suggests that there is a volcanic aquifer below the superconducting gravity station (YMSG near Chintiengang) bounded by an aquifuge at the sea level and by unknown lateral boundaries (Fig. 12). We hypothesize that this aquifer may partly supply source water to replenish the deep hydrothermal reservoir 2 km below Chishinshan (Murase et al, 2014). If the aquifuge is ruptured by earthquakes, groundwater in this aquifer will ow to the deep hydrothermal reservoir or even to the magma below Chishinshan. This rupture may greatly increase the probability of a phreatic eruption over the TVG. It is recommended that the interaction between groundwater and magma over the TVG should be investigated, based on recent works such as Jasim et al. (2018) and Hsieh and Ingebritsen (2019).

Conclusion
This study uses the time-varying gravity measurements in 2012 to delineate an aquifer around the TVG by hydrological modeling and EOF analysis. The gravity changes from the RG measurements in the four surveys were carefully processed to achieve a mean standard error of point gravity at about 8 µGal, which is small enough to sense the seasonal hydrological signal associated with the aquifer.
Our 3-D hydrological models (A and B) show porous media below the SG station (Chintiengang) around the TVG, where the hydraulic conductivities are optimized by the observed gravity changes through the modeled groundwater levels by Model B (Table 3). The consistency between the observed and modeled gravity changes over Chintiengang implies that the underlying strata contain inter-connected porous media above a low permeability layer (aquifuge). On the other hand, the sites with inconsistent observed and modeled gravity changes over Chishinshan imply that the subsurface layers at these sites may contain cracks and fractures that are not considered in Model B.
The aquifer identi ed in this study may pose a potential risk of phreatic eruption in the TVG. The risk can occur if the aquifuge below the aquifer is damaged by weathering and by external forces such as earthquakes. Groundwater leaking through a damaged part of the aquifuge can increase the water volume at the deeper parts of the TVG and thereby increase the water vapor pressure to raise the probability of an eruption. How long the aquifuge can sustain damages or whether leaking groundwater can cause a phreatic eruption can be potential subjects for future studies.
Appendix A: The TVG rainfall data in 2012   Groundwater levels for and from Model B. (a) Topography (gray background) and two major faults bounding the aquifer, water levels at groundwater wells (black stars) and rivers (colored lines), and  Differences between the observed and modeled gravity changes. The left column is for Model A and the right column is for Model B. Positive and negative differences are shown in red and blue circles, respectively. Yellow circles correspond to the sites where the gravity differences are small. Model B produces a better t than Model A over gravity sites near YMSG (the location around SG49). This better t implies the porous media around YMSG from Model B.

Figure 12
Possible hydrogeological structure along pro le A-B (Fig. 1b). The pink solid circles are the two hydrothermal reservoirs in Fig. 1a. The porous media form an aquifer below the SG49 site (YMSG). The aquifer's aquifuge and the hydrogeological settings for Model B result in gravity changes that are consistent with the observed gravity changes.

Supplementary Files
This is a list of supplementary les associated with this preprint. Click to download.