Absolute gravity change in Taiwan: Present result of geodynamic process investigation

Gravity values at 24 sites over 2004 2016 measured with absolute gravimeters are used to study geodynamic processes in Taiwan. We model rain-induced gravity effects and other temporal effects of non-geodynamic origins to obtain residual gravity, which cannot be fully explained by GPS-derived vertical displacements. We explain the gravity changes associated with deposited debris, earthquake, volcanism and Moho deepening. Gravity changes of 53.37 and 23.38 μGal near Sinwulyu and Laonong Rivers are caused by typhoon Morakot, leading to estimated volumes of 6.0 × 105 and 3.6 × 105 m3 in deposited debris. The observed co-seismic gravity change near the epicenter of the M 6.9 Pingtung earthquake (26 December 2006) is 3.12 ± 0.99 μGal, consistent with a dislocation-based gravity change at the μGal level, thereby supplying a gravity constraint on the modeled fault parameters. The AG record at the Tatun Volcano Group is the longest, but large temporal gravity effects here has led to a current gravity signal-to-noise ratio of less than one, which cannot convince a sinking magma chamber, but supply an error bound for gravity detections of long-term or transient magma movements. The gravity values at Ludao and Lanyu decline steadily at the rates of -2.20 and -0.50 μGal yr-1, consistent with the expected magma states of the two extinct volcanoes. The gravity rates at an uplifting site in central Taiwan and three subsiding sites in eastern Taiwan are negative, and are potentially caused by Moho deepening at a rate of -3.34 cm yr-1 and a combined Moho deepening and plate subduction at the rates of -0.18, -2.03, and -1.34 cm yr-1. Article history: Received 10 February 2017 Revised 12 June 2017 Accepted 13 June 2017


InTRoduCTIon
Absolute gravimetry (AG) is a branch of gravimetry that uses absolute gravimeters to determine gravity values for scientific and engineering applications.Repeated AG measurements at a given site can be used to determine the site's gravity changes to study site-specific problems associated with mass change (Torge 1989).In Taiwan, AGdetermined gravity changes were not available until the establishment of the Laboratory of Geodesy and Geodynamics (LOGG, http://www.logg.org.tw/) in 2004 by the Ministry of the Interior (MOI), Taiwan.Since 2004, gravity changes at a number of geodetic control points were observed by the two absolute gravimeters of LOGG, mainly for collecting gravity data to assist the assessment of the stabilities of the control points.The first geodynamic study of AG-determined gravity changes was initiated in 2006, under a joint Taiwan-France project, called "Absolute Gravity for Taiwanese Orogeny" (AGTO; Masson et al. 2008).AGTO is a project dedicated to measuring gravity changes in Taiwan.Under AGTO, ten gravity sites were installed in a cross section of southern Taiwan.The gravity changes collected under AGTO over 2006 -2010 have been used to examine two crustal models in southern Taiwan associated with mountain Terr. Atmos. Ocean. Sci., Vol. 28, No. 6, 855-875, December 2017 buildings (Mouyen et al. 2009), and to investigate gravity changes due to debris flows caused by typhoon Morakot (Mouyen et al. 2013).
In addition to the ten AG sites installed under AGTO, 14 AG sites were established by MOI under the project "National Gravity Datum Service" (NGDS) beginning in 2004.Figure 1 shows the 24 AG sites that are distributed over the island of Taiwan and its offshore islands, along with the horizontal and vertical velocity fields determined by GPS (http://gps.earth.sinica.edu.tw/main.jsp).Typically, a gravity site is installed next to a continuous GPS site, because horizontal and vertical deformations derived from GPS measurements can aid the interpretation of gravity changes at the site.Hydrological effects will also contribute to gravity change and can be modeled using data from a rain gauge.Around the 24 gravity sites in Fig. 1, there are 29 rain gauges, deployed by the Central Weather Bureau (CWB), that can provide rainfall records to account for hydrological gravity effects.Since 2004, most of the AG sites under NGDS are re-visited annually, but some are re-visited as little as two times over 2004 -2016.Despite the varying frequencies and numbers of revisits, these records provide time-lapse gravity values that allow for investigations of a number of geodynamic processes in Taiwan.
The 24 AG sites (Fig. 1) are distributed over the 6 geological settings of Taiwan, extinct and dormant volcanoes (Lo et al. 1994;Lin et al. 2005;Konstantinou et al. 2007), areas affected by land subsidence (Hwang et al. 2010) and areas of active collisions between the Eurasian Plate and Philippine Sea Plate (Yu et al. 1997;Malavieille 2010;Kuo-Chen et al. 2012).At some AG sites, gravity values were determined before and after major typhoons and earthquakes, allowing to see how these events change gravity values and how such changes can feed back the models describing these events.In central Taiwan, precision leveling has detected vertical uplifts of up to 2 cm yr -1 , which cannot be explained fully by contemporary models of mountain building (> Million years) (Ching et al. 2011).The discrepancy between the geodetic result and the models of mountain building is partially due to the fact that the geodetic result spanned only a few years, while mountain building processes act at the million year time scale.Another source contributing to the discrepancy is Moho deepening that has occurred in tectonically active regions such as Tibet (Sun et al. 2009a).Despite the suggested thickening of the lower crust by Yamato et al. (2009), Moho deepening around Taiwan has not been identified.
The above geodynamic issues motivate a comprehensive investigation of the gravity changes measured at the 24 sites Fig. 1.Distributions of 24 absolute gravity (AG) sites (circles), along with their nearest GPS sites (squares), over six geological settings of Taiwan.Also shown are the GPS-derived horizontal rates (arrows, with error ellipses) at 317 sites.The vertical displacement rates from GPS are interpolated into an areal rate (color-shaded) to show the pattern of uplift (positive rate) and subsidence (negative rate) across Taiwan.The mean horizontal displacement rate of the Philippine Sea Plate relative to the Eurasian Plate is 8.2 cm yr -1 (Kuo- Chen et al. 2012).(Color online only) (Fig. 1) over 2004 -2016 in this paper.We will show how the absolute gravity data are collected, processed (section 2) and interpreted in connection with selected geodynamic processes (section 3).Some of our gravity investigations in this paper are presented for the first time in the literature, based on preliminary models.Through this study, we aim to point out the important role of absolute gravimetry in understanding some of the geodynamic processes in Taiwan.

RePeATed ABSoluTe gRAvITy oBSeRvATIonS: fRom RAw gRAvITy To
ReSIduAl gRAvITy

The geological Settings Around the Absolute gravity Sites
As stated in section 1, the 24 AG sites are distributed over the 6 geological settings of Taiwan.Table 1 shows the names of the AG sites, their nearest continuous GPS sites, rain gauge stations and the geological setting where a AG site is situated.Three of these AG sites are located at the offshore islands of Penghu (AG9), Lanyu (LYUG), and Ludao (AG1).Lanyu and Ludao are two extinct volcanic islands off southeastern Taiwan.AG9 and LYUG are located at the Penghu and Lanyu weather stations, and their distances to the ocean are 0.26 and 1.37 km, respectively.YMSG is a site located at the Tatun Volcano Group (TVG) and is situated between the Shanchiao Fault and the Taipei Basin.YMSG (established 2004) is co-located with a Superconducting Gravimeter (SG) site (Serial number #T49, established 2012).The latest volcanic activity of TVG occurred some 5000 years ago, so TVG is considered a dormant volcano (Konstantinou et al. 2007).The most likely causes of gravity change around TVG will be magma migration (Locke et al. 2003) and groundwater level change (see section 2.2).
The 10 AGTO sites are along an east-west transect across southern Taiwan and are selected mainly for measuring gravity changes associated with mountain buildings.In August 2009, the river courses around the eastern (AG3) and western (AG6) Central Range were clogged by debris triggered by typhoon Morakot (Mouyen et al. 2013); later the debris was removed for continued gravity observations.It turns out that AG3 and AG6 have the largest gravity changes among 24 sites (see section 3.1) as a result of the large mass increase due to debris deposited in nearby rivers.AG5 was the highest AG site, but it was destroyed by a landslide triggered by typhoon Morakot in August 2009, so here we collected only 3 AG measurements (2006 -2008).AG2a and AG2b are situated on the eastern and western sides of Coastal Range, and these two AG sites are in the Tuluanshan Formation.AG7 is located on the eastern side of Zuojhen Fault in western foothills.Zuojhen Fault is a left-lateral fault lying in the northwest-southeast direction across the township of Zuojhen.AG8 is located in the southwestern coastal plain, and the site is in the basement of a building belonging to the Department of Geomatics, National Cheng Kung University, Taiwan.
The purposes of the 14 NGDS sites are manifold, and the site selection is based on the guidance that gravity sites should be distributed evenly over the entire Taiwan and its offshore islands.In addition, the following NGDS sites also serve as gravity control points for relative gravity measurements (Hwang et al. 2014): WFSG, DSIG, HCHG, TCHG, KDNG, and FLNG.WFSG is located in the Wufenshan weather radar station in New Taipei City.DSIG is situated at the National Defense University in Taoyuan County, where the geological setting belongs to the lateritic terrace deposits.HCHG, co-located with Hsinchu SG (#T48; Hwang et al. 2009), is located in the Toukoshan Formation between alluvium and lateritic terrace deposits.HCHG is also near the Hsinchu Fault, which has been re-classified as an active fault of the second category (Chang et al. 1998;Lin et al. 2000;Hwang et al. 2009).This fault category implies that the Hsinchu Fault could be active in the last 100000 years.PKGG and TAES are located in Yunlin County and Chiayi County, where severe land subsidence exists (Hung et  2010).Earlier gravity measurements at TAES and PKGG were used for inferring land subsidence in the two counties and the result indicated that gravity changes here could originate from a number of sources.Therefore, the AG data at these two sites will not be analyzed in this paper because our objective is to study the natural geodynamic process rather than man-made land subsidence.Finally, SMLG, TLGG, FLNG, and YLIG are located in an area experiencing active mountain building and Moho deepening in central and eastern Taiwan (Lin and Ando 2004;Ching et al. 2011).

gravity uncertainties and Temporal gravity effects
From 2004 to 2016, we collected gravity measurements using two FG5 absolute gravimeters from Taiwan (FG5 #224 and #231) and one from France (FG5 #228).An FG5 gravimeter determines the vertical acceleration by measuring the free falling distances and times of a corner cube along its dropping chamber.Because of the high accuracies in distance and time measurements, an FG5 gravimeter can achieve one μGal gravity accuracy (Niebauer et al. 1995;Okubo et al. 1997) at a quiet site.However, gravity measurement accuracies can be degraded if the gravimeter is situated over a sedimentary plain, near oceans or an area of high human activity.For all the AGTO sites, we measured gravity values every November, with the exceptions that the gravity measurements at AG3, AG4, and AG5 were collected in March 2008.Collecting gravity data in the same month of the year could potentially remove common-mode temporal gravity effects, particularly the hydrology gravity effect, which is largely seasonal effect in Taiwan (Van Dam and Wahr 1998;Boy and Hinder 2006;Hwang et al. 2009).
At 14 NGDS sites, the AG measurements were collected by the FG5 gravimeters from Taiwan, but AG measurements at 10 AGTO sites were collected by both Taiwanese and French gravimeters.Specifically, at 10 AGTO sites, the gravity measurements in 2006, 2008, and 2010 were collected by FG5 #228, and the measurements in 2007 and 2009 by FG5 #224.For each measurement campaign using FG5 #228, we determined an offset between FG5 #224 and #228, which was then used to correct for the inter-gravimeter bias in the measurements.For each of the measuring campaigns and for each of the AG sites, we also determined a mean gravity gradient by a Graviton EG gravimeter, which was used to reduce the raw gravity measurements to a value at a common reference point and the reference point was one meter above the survey mark of the gravity site.
For all the AG sites, 12 and 24 hr are two typical lengths of measurement sessions.In a session, set gravity values based on 100 drops set -1 were determined, with a set lasting 30 min.In this study, we use the procedure recommended by Taylor and Kuyatt (1994)  A uncertainty is the standard uncertainty (error) of the mean gravity, computed by where n = 43 is the number of sets, and x i is the gravity at the i th set.Using the 43 sets of observations, we obtained m v = 0.5 μGal as shown in Table 2.This value is actually the standard deviation of the 43 gravity measurements without considering the Type B uncertainties.
The Type B uncertainties are associated with the uncertainties in the instrument parts of the FG5, and the uncertainties in the models for correcting temporal gravity effects.Table 2 lists all the contributors to the instrumental uncertainty and the site uncertainty.The final (total) Type B uncertainty is simply computed by where p is the number of contributors given in Table 2.For the case study in Table 2, we have p = 10 and 7 for the instrumental uncertainty and the site uncertainty.Using the standard uncertainties and degree of freedoms (DFs, DF = number of observations minus number of unknowns) by Welch-Satterthwaite model (Satterthwaite 1946;Welch 1947) to estimate the effective degree of freedom in Type B, we obtain the uncertainties of 1.9 and 1.5 μGal and the degree of freedoms of 85 and 18 for the instrument and the site, respectively.Combining m v , the instrumental uncertainty and the site uncertainty, we obtain the expanded uncertainty for the case in Table 2 + + = 2.5 μGal.An expanded uncertainty can vary with the number of sets [n in Eq. ( 1)], the condition of the FG5 at the measurement time, and the sophistication of the correction models for reducing the site uncertainty.For example, m v will decrease with increasing number of sets.On the other hand, if an FG5 gravimeter is well maintained and its instrument parts perform equally well at any time, the Type B uncertainty will remain a constant.The current correction models for the standard uncertainty have already achieved sub-μGal accuracies (Table 2), but there could be room for improvement.
In this paper, we modeled the gravity effects of nongeodynamic origins in several stages.In the first stage, a raw gravity value was corrected for the gravity effects of solid Earth tide (SET), ocean tide loading (OTL), atmosphere, and polar motion, resulting in the first-stage gravity value defined as: where g raw : raw gravity value relative to the mean of all gravity values g set : gravity effect due to SET effect g otl : gravity effect due to OTL effect g a : gravity effect due to atmospheric pressure change g p : gravity effect due to polar motion The modeling of SET, OTL, atmospheric pressure change and polar motion follows the procedure given in Hwang et al. (2009).
If gravity changes of geodynamic origins at a gravity site are sought, the vertical motion-induced gravity change at the site must be removed.The ratio between vertical motions and gravity changes vary, depending on the processes responsible for the ground motions.In this paper, we used a gradient combining the free-air gradient and the Bouguer gradient, i.e., Δg/Δh = -3.086(free-air gradient) + 0.1119 (Bouguer gradient) ~ -1.967 μGal cm -1 (De Linage et al. 2007).Considering the accuracies of vertical motions of GPS, we did not apply point-wise corrections of vertical motion, rather we computed the vertical motion-induced gravity changes and removed them from the first-stage gravity values.That is, the second-stage gravity value, g se , is determined by where g f is the first-stage gravity value defined in Eq. ( 3), Δg/Δh is the gravity-to-height ratio and V h o is the rate of vertical movement and Δt is the time difference relative to a reference epoch.

Rain-Induced gravity effect
After removing the vertical motion-induced mass change, the patterns of the second-stage gravity changes (absolute values) are correlated with the patterns of hydrological effects.In section 2.2, the g se in Eq. ( 4) still contains the hydrological gravity effect, which should be removed as: where g res is the residual gravity value after removing the hydrological effect, g h , and g normal is the mean gravity value.
In this paper, we modeled the hydrological effect, g h , at time i using rainfall data near a given AG site as (Crossley et al. 1998;Mouyen et al. 2016): where G is the gravitational constant (6.67 × 10 -11 m 3 kg -1 s -2 ), w t is the water density (1000 kg m -3 ), r j is the amount of rainfall (from a rainfall gauge) accumulated over time j, τ 1 and τ 2 are the charge and discharge times after rainfalls around the AG site.Equation ( 6) models the gravity effect caused by water flow in an infinite water layer (Bouguer slab model); see Crossley et al. (1998) for a detailed model description.A g h value in Eq. ( 6) is called the rain-induced gravity effect (RIGE).Here we show an example at Hsinchu AG/SG station, HCHG (Fig. 1).We experimented with τ 1 = 0.1 to 10 days and Typically, the RIGE at a gravity site is the largest in June, so that we measure AGTO data in November to avoid large uncertainties in modeling RIGE in June.In Taiwan, the hydrological gravity effects in the same month over different years are likely to have a similar magnitude.Therefore, this hydrological gravity effect can be reduced if repeat gravity measurements are collected in the same month.This is the case for gravity measurements at the gravity sites in the AGTO project: we collected data in November over 2006to 2010(except in March 2008)).In contrast, the gravity measurements at AG sites in the NGDS project will contain time-varying hydrological gravity effects that cannot be canceled, because here measurements were not collected in the same month.Therefore, hydrological corrections (in this paper, RIGE) for the gravity measurements at NGDS sites must be applied before interpreting the gravity changes in terms of geodynamic processes.
The RIGE values at 24 AG sites depend on a number of factors: topography is the leading factor, followed by geographic location.The hydrological correction can be very localized.For example, at AG site YMSG over TVG, we found that in the monsoonal season, a rain gauge at the windward side tends to receive more rains than a rain gauge at the leeward side.It is very important to decide the rain gauge station for use when correcting AG measurements.In the case YMSG, we used measuremrents from SG #T49 to decide the optimal rain gauge.In general, RIGE values in mountainous areas (central Taiwan) are larger than those over plains.Also, in the same season, RIGE values in eastern Taiwan are larger than those in western Taiwan.

gravity Changes for Investigating geodynamic Processes
Figure 3 shows the time-lapse residual gravity values [Eq.( 5)] at sites with more than two visits by FG5.Such residual gravity values are free from hydrological and vertical displacement effects within the uncertainties of the correction models.At AG sites DSIG, JSIG, LYUG, TCHG, TLGG, and YLIG, AG data were measured only twice, so it is not likely to obtain reliable gravity rates here without modeling RIGE values, which dominate the gravity changes and can reach 10 μGal.At some AG sites, because of the poor gravimeter performances and the effects of extreme rainfalls and typhoons, only part of the AG measurements were used to estimate the gravity rates; Table 3 shows the rates of gravity changes and their uncertainties based on the first-, second-stage, and residual gravity values determined by the least-squares method.In theory, the residual gravity rate g res o at a site is caused only by mass changes associated with one or more geodynamic processes and will be used for the investigations in section 3.In addition to the gravity rates, the gravity differences before and after a typhoon event and an earthquake event will be used to estimate parameters associated with these events.

Typhoon morakot-Induced large gravity Changes Associated with landslide and River deposits
During a typhoon, large gravity oscillations can occur because of the high-frequency changes of atmospheric pressure and rainfall conditions.A typhoon can also result in a large gravity-atmosphere admittance; from the gravity and pressure measurements at Hsinchu SG station, Hwang et al. (2009) showed that the average gravity-atmosphere admittance during typhoons over 2006 to 2010 is -0.48 μGal hPa -1 , compared to a typical value of -0.35 μGal hPa -1 .In particular, when typhoon Morakot was around Taiwan, the gravity-atmosphere admittance was -0.68 μGal hPa -1 and typhoon Morakot caused up to 3000 mm of rains in southern Taiwan from 6 to 10 August 2009.
Typhoon Morakot was a high-intensity typhoon -the strongest typhoon to hit Taiwan in the last 50 years, and triggered widespread landslides and debris flows in southern Taiwan (Chu et al. 2011;Lin et al. 2011), especially around AG3 and AG6.Short-term atmosphere-induced gravity changes, such as those due to typhoons, can only be predicted well with high-frequency atmospheric data (Petrov and Boy 2004;Chou et al. 2011).The three optimal AG sites for studying typhoon-induced gravity changes are AG3, AG5, and AG6.Because AG5 was destroyed by typhoon Morakot, here we show such gravity changes only at AG3 andAG6 over 2006 to 2008 (Fig. 3 2013) already reported large gravity changes (up to 285 μGal) in the Western Foothills, the Central Range and the Longitudinal Valley of Taiwan using measurements from absolute and relative gravimeters; such gravity changes were interpreted by a model of morphological change.In this paper, we re-examine gravity changes at AG3 and AG6 and present an extra interpretation of the changes before and after typhoon Morakot.
In the case study of Kim and Lowe (2004), they suggested that deposits of debris flows could occur over terrains with slopes of 3 -6°, and transportations could occur in the upper alluvial fan with slopes of 6 -26°.However, precise slopes allowing for debris flow deposit and transportation depend on many factors.Here a slope is defined as arctan (elevation change/distance change).Figure 4a shows the elevations around AG3 and AG6, indicating that these two sites are located over river sediments that are surrounded by mountains.The slope between AG2 and AG3 is about 5.6° (altitude difference divided by horizontal distance) and the slope is about 4.3° between AG6 and AG7.In-situ images around AG3 (Fig. 4b) and AG6 (Fig. 5b) show that the maximum thicknesses of deposits here can be up to 10 m (Fig. 4b) and 8 m (Fig. 5a).Satellite images from FORMOSAT-2 (http://geodac.rsh.ncku.edu.tw/files/11-1151-10403.php?Lang=zh-tw) indicate significant changes of the riverbed coverages before and after typhoon Morakot (Fig. 4c for AG3 and Fig. 5c for AG6).For comparison, the village around AG6 was partially covered by landslides triggered by typhoon Morakot, while the nearby Hsiaolin Village was completely covered by landslides that in total have a mudslide volume of 2.7 × 10 7 m 3 in an area of 2.5 km 2 (Wu et al. 2011).Furthermore, we adopted an average density of t = 2.0 g cm -3 for the debris flow deposit (from recent sediment), instead of the normal density of 2.67 g cm -3 (Yen et al. 1998).Assuming a mean distance of 8.5 km from Hsiaolin Village to AG6, we estimate that the gravity change due to the mudslide of Hsiaolin Village is 0.40 μGal using the 3D gravity method of Mouyen et al. (2013) and Kao et al. (2014).This gravity change is smaller than the uncertainty of FG5, and therefore we assume that the major gravity effect is within an effective radius of 0.5 km around AG3 (Fig. 4d) and AG6 (Fig. 5d).As stated earlier (Fig. 4), AG3 and AG6 are situated over flat areas along rivers, use of the simple Bouguer plate can generate sufficiently accurate gravity effects due to the river sediments.The point here is to highlight the value of gravity changes for erosion research in Taiwan.Specifically, gravity changes due to debris flows at AG3 and AG6 are approximated as where Δg b is the gravity change associated with a "partial" Bouguer plate defined by the ratio k a between the area of the deposited river course and the area within the effective radius (i.e., 0.5 km), t is the density of sediment, G is the gravitational constant, and H is the sediment thickness.The sources of the deposited materials around AG3 and AG6 are from the Central Range and Western Foothills.Around AG3 (Fig. 4d) and AG6 (Fig. 5d), the areas of deposited debris flows are 0.11 and 0.17 km 2 , respectively, and the landslide-covered areas are 0.02 and 0.11 km 2 , respectively.Using the area of the effective radius [π × (0.5 km) 2 = 0.785 km 2 ] and the area covered by debris flow, we estimated that k a = 0.14 (AG3) and 0.21 (AG6), which were used to estimate the thicknesses of sediments in the riverbed.Using Eq. ( 5) and the gravity changes of 53.37 μGal (AG3) and 23.38 μGal (AG6), we estimated that the thicknesses of sediments at AG3 and AG6 are 4.66 and 1.31 m, respectively.The two thicknesses result in 6.0 × 10 5 and 3.6 × 10 5 m 3 of volume changes in the riverbeds around AG3 and AG6.Such volume changes are close to but smaller than the volume changes of 3.0 × 10 6 and 3.  First, instead of measuring the gravity changes before and after the day of typhoon Morakot, we measured the gravity changes 3 months after typhoon Morakot.That is, at the time we did the gravity measurements (November 2009), the areas of erosion and landslide around AG3 and AG6 may have been modified.Second, the volume changes around AG3 and AG6 are partially contributed by landslides around the Laonong River and the Sinwulyu River.Third, the uncertainties of the RIGE model (section 2.3) can contribute to the discrepancies.Finally, after typhoon Morakot, dredging of deposits around AG6 was carried out because AG6 is near a popular tourist attraction.The dredging has decreased the gravity value at AG6 by 4.21 μGal from 2009 to 2010.However, no dredging of deposits was made near AG3.As such, continual deposits around AG3 have increased the gravity value at AG3 by 19.56 μGal from 2009 to 2010 (Fig. 3), while the gravity value at AG6 stopped to increase from 2010 onward due to the man-made dredging.
Another potential site for studying mass change due to typhoons is AG4.Here, the difference between the gravity values observed on 5 November 2008 and on 23 November 2009, is 11.59 μGal (the latter is smaller).Again, the difference was caused by typhoon Morakot.However, unlike sites AG3 and AG6, AG4 is situated at a summit of about 1500 m height.Satellite images before typhoon Morakot showed a clear slope on the eastern flank of this summit affected by landslides.Typhoon Morakot did not change the landscape dramatically, but decreased the gravity value at AG4.At AG4, it will require a different method of analysis than the one we used at AG3 and AG6 to understand how the mass in the original slope around AG4 has lost to typhoon Morakot (Mouyen et al. 2013).In addition, the large mass loss around AG4 implies that one has to be very cautious in interpreting gravity changes due to mountain building when using gravity data at a site like AG4 that is prone to erosion.

earthquake-Induced gravity Changes: example of the 2006 m w 7.0 Pingtung earthquake doublet
Co-seismic surface displacement and mass change occurring around the hypocenter of a major earthquake can be detected by sensors such as GPS and gravimeters.In Taiwan, the co-seismic surface displacements associated with major earthquakes have been studied (Lin and Ando 2004;Johnson et al. 2005;Wu et al. 2006).However, only a few studies on gravity detection of earthquake-induced mass change by gravimetry were reported (Yen et al. 1998;Lo and Hsu 2005) and no AG measurement was used in these studies.Here we present an example of gravity-detected mass change associated with the 26 December 2006 Pingtung earthquake doublet using the AG measurements at KDNG site.In this earthquake event, two M w > 6.8 mainshocks and thirty-seven M w > 3.0 aftershocks occurred around Ping-tung between 26 and 31 December 2006 (CWB seismological reports).The first mainshock occurred at 12:26:21 UTC (120.56°E,21.69°N).The moment magnitude is 6.9 and its focal depth is 44 km.The epicenter is located 37 km southwest of KDNG.After 8 min, the second mainshock (M w = 6.8) occurred at 120.42°E, 21.97°N with a 52-km depth, and its epicenter is 36 km to the north-northwest of the first mainshock and 38 km to KDNG.
We used the model of Sun et al. (2009b) to determine the deformations caused by the Pingtung earthquake based on the finite parameters of Yen et al. (2008).Table 3 shows the parameters associated with the focal mechanism of the Pingtung earthquake.We modeled the co-seismic gravity change, as seen from space, using (Sun et al. 2009b) where ( , , ) is gravity change, a, i, and { are the radius, latitude and longitude of the site of interest, g 0 denotes the gravity value on the Earth's surface, the variable k nm ij is the dislocation Love number of potential change with n and m being the degree and order, Y n m is the spherical harmonic function, U denotes the dislocation vector of the slip v i and normal components n j , and indices i, j represent different types of earthquakes when using Love numbers (Sun et al. 2009b): a strike-slip (superscripts ij = 12), a dip-slip (superscripts ij = 13 or 23), a vertical tensile fracturing (superscripts ij = 11 or 22), and a horizontal tensile fracturing (superscripts ij = 33).From the distribution of aftershocks, which plays an important role in revealing the mechanisms of the main earthquake (Wu et al. 2009), the first mainshock is a west-dipping normal fault (Fig. 6a) and the second mainshock is a strike-slip fault.The maximum co-seismic slips are 2.3 and 1.2 m, respectively (Yen et al. 2008), occurring near the hypocenter.A subsidence was detected by GPS at KDNM (Chen et al. 2008), and the co-seismic slip propagating directions of these two events were northward and down-dip direction on their fault planes.
Here we show how our gravity measurements (Fig. 6c) can benefit earthquake studies using the example of Pingtung earthquake.We used the AG measurements on 26 February 2006 and on 7 February 2007 and the interseismic gravity change rate (-1.35 μGal yr -1 ) to extrapolate the expected gravity values at the epochs right before and after this earthquake.This results in a co-seismic gravity change of 3.12 ± 0.99 μGal at KDNG due to this earthquake (Fig. 6c).The standard error of the gravity change (3.12 μGal) is derived from the measurement uncertainties and the time span of 11 months (26 February 2006 and7 February 2007) by error propagation.Using the fault parameters from the database of finitefault rupture models (Table 4) and the Love numbers from   Sun et al. (2006), we computed the gravity change and vertical displacement due to the Pingtung earthquake at KDNM (Fig. 6) (Note that KDNG is the AG site, and KDNM is the GPS site; they are in the same building but have different pillars).Our model runs result in a gravity change of 2.40 μGal, which differs from the observed value (3.12 μGal) by 0.72 μGal (KDNG).Note that the modeled vertical displacement of -1.20 cm is close to the GPS-observed value of -1.64 cm during the Pintung earthquake (Fig. 6d).The consistency (at the μGal level) between the modeled and observed gravity changes at KDNM during the Pintung earthquake shows that absolute gravity changes can provide a constraint for coseismic source modeling, in addition to GPS observations.As such, we recommend that AG measurements be collected in the earthquake-prone regions of Taiwan to provide strengthened constraints when estimating fault parameters in the event of a future major earthquake.

gravity Changes at dormant, Active, and extinct volcanoes in Taiwan
The AG sites YMSG, AG1, and LYUG are located over three volcanoes in Taiwan: TVG, Ludao, and Lanyu.Here we investigate the states of these volcanoes using gravity changes collected in this paper.Figures 7a and b show the free-air and Bouguer gravity anomalies (Hwang et al. 2014) around the three volcanoes and gravity changes from our AG measurements.First, we study TVG using gravity changes at site YMSG.TVG is a volcano group that has been classified as being active to dormant (see section 1).Two operational nuclear power plants are on the flanks of TVG.A magma chamber may still exist below TVG, based on four pieces of evidence: (1) compositional gases (Lee et al. 2008;Witt et al. 2008), (2) volcano-induced seismicity (Kim et al. 2005;Lin et al. 2005), (3) observational ground deformations (Chang et al. 2014), and (4) S-wave shadowing (Lin 2016).Presently, TVG mostly shows hydrothermal activities, whose intensities seemed to rise between 2004 and 2007 (Lee et al. 2008;Rontogianni et al. 2012;Mouyen et al. 2016).However, the existence of a TVG magma chamber is an active research subject.
Using the AG measurements at YMSG from 31 July 2007 to 14 December 2014, we obtained a residual gravity rate of -0.15 ± 0.51 μGal yr -1 .If the source of gravity change is entirely caused by the magma redistribution, the negative gravity rate will indicate a sinking magma chamber below TVG.However, the fact that the uncertainty of the gravity rate is higher than the rate indicates that the temporal gravity effects around TVG, especially the hydrological effect, may obscure the gravity signal originating from the magma redistribution.
Because the gravity uncertainty is larger than the gravity signal, here we carry out a preliminary estimate of the magnitude of magma sinking rate in order for the AG to detect it.Assume that a cylinder-shaped magma chamber is located at a depth H under the surface, with the difference in density between the crust and magma, t D , and radius R (Fig. 8).In Fig. 8, the vertical movement of the magma chamber is Δh.Furthermore, we assume that the gravitational attraction of the magma cylinder to YMSG is due to a point mass at the center of the cylinder, so that the acceleration at the magma cylinder -YMSG direction, g mc , can be expressed as where r is the distance from the center to YMSG.Using Eq. ( 10) and the parameters in Fig. 8, we can determine the gravity change (vertical component) sensed at YMSG as where x is the horizontal distance from the magma chamber to YMSG and H is the vertical distance between YMSG and the top of the chamber.Because Δh is much smaller than H, Δh can be ignored in the term (H + Δh).Under this assumption, the rate of gravity change is where h o is the rate of Δh.We assume that the depth of magma chamber is H = 8 km (Konstantinou et al. 2007), the radius of magma chamber is R = 2 km and the horizontal distance from magma chamber to YMSG is x = 1.3 km (the distance from the center of positive anomalies to YMSG in Fig. 7b).Assume that magma signal can be detected when the gravity signal-to-noise ratio is one (corresponding to a 68% confidence level).Under this assumption, and with g m o = -0.51μGal yr -1 , the magma sinking rate will be For a density contrast of 0.57 g cm -3 and g m o = -0.51μGal yr -1 , the rate is h o = -0.14mm yr -1 .That is, a magma-sinking rate must be larger than 0.14 mm yr -1 for our absolute gravimeters to detect the induced gravity signal.We expect to continue our model improvements of temporal gravity effects and extend AG measurements at YMSG to reduce the Next, we study Ludao and Lanyu, which are two islands in the Luzon Volcanic Arc (LVA) The LVA is the result of the continuous subduction of the Eurasia Plate into the Philippine Sea Plate (Yang et al. 1996).Due to the collision of the Eurasia Plate continental margin in the LVA, these two islands are moving westwards and the volcanisms around the two islands ended 1.5 and 0.5 -0.04 Ma, respectively (Yang et al. 1988;Lo et al. 1994).As shown in Figs.7d and g, the free-air gravity anomalies are positive around the two islands, with the largest values on the summits of the islands.Such patterns of gravity anomalies are typical for volcano islands created by the lithosphere's subduction into and collision with an island arc (Robert 2008).Furthermore, the vertical displacement rate relative to the origin of the ITRF is 0.29 cm yr -1 at LUDA (Ludao) and is 0.11 cm yr -1 at LANY (Lanyu), suggesting that the surface of these two extinct volcanic islets are rising.The GPS-derived vertical rates conform to the expected rising rates over an island arc such as the LVA.
One unique geological feature of Ludao and Lanyu is the seawater hot springs in the southeastern parts of the islands, implying ongoing geothermal activities similar to that around TVG.However, we found no study about the volcano states of Ludao and Lanyu linking to such geothermal activities.Using the AG measurements at AG1 and LYUG, we obtained a residual gravity rate of -2.20 μGal yr -1 for Ludao (Fig. 7f) and -0.50 μGal yr -1 for Lanyu (Fig. 7i).The negative gravity rates imply that the heavier materials that formed the volcanoes of Ludao and Lanyu are retreating to the mantle.More investigations are needed to understand the mechanism behind the negative gravity rates, and this is a subject of future study.

gravity Changes due To moho deepening and Plate Subduction
In this paper, the last two geodynamic processes to be investigated from gravity changes are Moho deepening and plate subduction.Moho deepening is the isostatic response to the vertical crustal extension of a plate, and has a long wavelength effect because it typically occurs at large depths.However, due to local mass disturbances, it is not realistic that we can use all 24 AG sites to study Moho deepening in Taiwan.In this paper, we believe the potential sites for studying Taiwan's Moho deepening are SMLG and some AGTO sites (AG2a, AG2b, AG3, AG4, AG5, and AG6), and sites for studying the combined effect of Moho deepening and plate subduction, situated in eastern Taiwan, are TLGG, FLNG, and YLIG.Here we exclude the AGTO sites where the AG measurements may have been contaminated by typhoon Morakot.
First, we focus on Moho deepening in central Taiwan using the gravity measurements at SMLG, where the first AG value was collected in 2005, the second in 2009 and the last in 2016.These three AG values show a steady decline at a rate of -1.54 ± 1.18 μGal yr -1 (see the second-stage gravity rate in Table 3) after correcting for the gravity effect of vertical displacements that have a rate of 0.55 cm yr -1 based on a nearby GPS site SUN1.The large uncertainty (1.18 μGal yr -1 ) of the estimated gravity rate (-1.54 μGal yr -1 ) is probably due to the hydrological gravity effect, because seasonal hydrological gravity effects may influence gravity rate estimations (Van Camp et al. 2010, 2016).After removing the hydrological effect using Eq. ( 6), we obtained a gravity rate of -0.98 ± 0.07 μGal yr -1 at SMLG (Fig. 9b), which indicates a mass loss.The erosion rates in Taiwan range from 1 -2 mm yr -1 around SMLG (Dadson et al. 2003), which contribute to some -0.1 to -0.2 μGal yr -1 of gravity change rates (depending on the densities of eroded materials) and can be neglected.According to Li et al. (2014), the Moho depth around SMLG is the deepest in Taiwan and is approximately 50 km (Fig. 9a).The Moho depths of Li et al. (2014) are derived from the contours of Vp = 7.50 km s -1 .A large Moho depth can be caused by a rapid Moho deepening.Therefore, Moho deepening is a candidate contributor to the mass loss detected by AG at SMLG.Note that Yen et al. (1998) has used gravity to study the isostasy of Taiwan, which is associated with Moho's static state in Taiwan, not Here, we postulate that a deepening of Moho in central Taiwan is ongoing and it is, therefore, possible to estimate the deepening rate using gravity changes as in Sun et al. (2009a) over Tibet.This deepening can be transient, but it is a plausible scenario to explain the negative gravity rates at SMLG, TLGG, FLNG, and YLIG.Note that we have only two AG observations at TLGG and YLIG.Following the method of Sun et al. (2009a), as a preliminary estimate we used a Bouguer plate model to compute the deepening rate by: b where t D is the difference in density between the mantle and crust and g res o is the residual gravity rate in Eq. ( 5).In this paper, we adopt t D = (3.3-2.6) g cm -3 .Using g res o = -0.98 μGal yr -1 , we obtained a deepening rate of -3.34 cm yr -1 at SMLG (b m o , a negative rate indicates deepening).Compared to the deepening rate of -2.30 cm yr -1 and vertical displace-ment rate of 0.12 cm yr -1 over Tibet (Sun et al. 2009a), the deepening rate of -3.34 cm yr -1 from AG and the vertical velocity rate of 0. The sites TLGG, FLNG, and YLIG in eastern Taiwan are situated on the same plate as SMLG.As such, the Moho deepening found at SMLG should be extended to the plate below TLGG, FLNG, and YLIG, where surface uplifts would have occurred.However, the vertical displacement rates from GPS at sites TLGG, FLNG, and YLIG are negative, indicating that these sites undergo surface and subsurface deformations that combine to yield negative rates of vertical motions.One source of vertical motion is plate subduction that can cause these sites to subside at rates larger than the uplifting rates caused by Moho deepening.In Fig. 10, we explain the role of plate subduction that causes the negative rates (or subsidence) at FLNG site; see also Johnson et al. (2005) and Ching et al. (2011).Using Eq. ( 14), we obtained the following combined rate (Moho deepening and plate subduction) at TLGG, FLNG, and YLIG: -0.18, -2.03, and -1.34 cm yr -1 .The largest rate occurred at FLNG, which is close to the spot of the largest Moho depth.The second largest rate occurred at YLIG, followed by TLGG.
Because of the convergence between the Eurasia Plate and the Philippine Sea Plate, the growth of mountain in Taiwan should be associated with the thickening of buoyant roots into the mantle (e.g., Ho 1986).However, Moho deepening in the Central Range has been proposed mainly based on the detections of varying Moho depths (e.g., Rau and Wu 1995;Yeh et al. 1998;Yen et al. 1998;Kuo-Chen et al. 2012;McIntosh et al. 2013;Zheng et al. 2013;Wu et al. 2014), the temporal variation of surface uplift rate/ exhumation rate (e.g., Gourley et al. 2007), or other geophysical observations associated with the 1999 M w 7.6 Chi-Chi, Taiwan earthquake (e.g., Hsu and Lo 2004;Shin et al. 2011).In addition, several 2D fault models have been used to simulate the surface deformation across the Taiwan mountain belt (e.g., Johnson et al. 2005;Ching et al. 2011).
In particular, the model constructed by Ching et al. (2011) indicates that not only the slips on major active faults but also other opening-mode sources are needed to account for the high uplifting rates derived from precision leveling and GPS across Taiwan's mountain belts.The opening-mode sources here may be the Moho deepening and plate sub-duction.The combined rates at TLGG, FLNG and YLIG in eastern Taiwan indicate that these sites are very close to the plate suture zone and the estimated surface motions are all undergoing subsidence, noting that the effect from Ryukyu subduction cannot be eliminated.In addition, because of the coupling on the plate interface of the Ryukyu subduction zone (Hsu et al. 2012), the materials close to the plate boundary (Wu et al. 2014) may be dragged into the deep of the Earth, resulting in the ground subsidence and the sinking of the Moho during the inter-seismic period.Despite the varying gravity rates at the four sites studied in this section, the rates are consistently negative and are indicative of losing mass below the surface where the depths of Moho are larger than elsewhere.Our gravity observations may fill the data gap in understanding Moho deepening in Taiwan, for which there is no direct evidence in previous studies.

ConCluSIon And PRoSPeCT
In this paper, we have used time-lapse AG measurements collected over 2004 to 2016 to study a number of geodynamic processes in the 6 geological settings of Taiwan.In total, AG measurements at 24 sites were used, but the numbers of measurements over 2004 to 2016 vary from one site to another.Compared to previous studies that used AG measurements in Taiwan, we applied vertical and rain-induced corrections to our AG measurements, thus avoiding seasonal gravity variations that can affect our interpretations of the proposed geodynamic processes.The geodynamic processes are associated with landslides and debris flows triggered by typhoon Morakot, Pintung earthquake, volcanoes at TVG, Ludao and Lanyu, and Moho deepening and plate subduction in central and eastern Taiwan.Except landslides and debris flows, for the first time our gravity change measurements provide reasonable estimates of the critical parameters for these processes.
We estimated the volumes of landslides and debris flows in southern Taiwan (AG3) and in eastern Taiwan (AG6) using gravity changes before and after typhoon Morakot.At AG3, the gravity-based estimated volume is lower than the official value because the effective area of gravity-based estimates is limited to a radius of 0.5 km.The variations of sediment thickness inferred from the gravity changes at AG3 and AG6 well explain the transient mass changes due to debris flows after typhoon Morakot.To improve the estimation of sediment thickness at the two sites, a detailed digital terrain model should be used and such a model may be obtained from Lidar and other sensors.
The fault slip distribution and the dislocation model at KDNG returned co-seismic gravity changes that are consistent with the AG measurements at the μGal level.This high consistency between the modeled and observed gravity changes suggests that AG measurements can be used to validate source parameters associated with earthquakes in Taiwan; more AG sites can be installed in areas prone to large earthquakes.We also study volcanoes in Taiwan using gravity changes.The long-term AG measurements at YMSG show little gravity variation around TVG, indicating a stable state of the proposed magma chamber.The estimated sinking rate of a potential magma below TVG is less than 0.1 cm yr -1 , which can be disputable due to its larger uncertainty caused by large hydrological gravity effects.The AG measurements at Ludao and Lanyu indicate that the two islands are extinct volcanoes, with gravity declining at the rates of 2.20 and -0.50 μGal yr -1 .Again, the estimated rates should be improved with more AG measurements.
The AG site SMLG, located at a spot in central Taiwan with the largest Moho depth and a large uplifting rate, shows monotonically declining gravity values over 2005 -2016 (11 years) that are highly indicative of Moho deepening in Taiwan.We also used three AG sites in eastern Taiwan to demonstrate the gravity effects of Taiwan's Moho deepening.Unlike SMLG, the three sites in eastern Taiwan are subsiding, thereby requiring an additional mechanism, i.e., plate subduction caused by the collision of the Eurasia Plate and Philippine Sea Plate, to explain the declining gravity values.
The result presented in this paper shows that gravimetry is an alternative tool to popular tools such as GPS and seismology for studying geodynamic processes in Taiwan, particularly those causing mass changes.Most of the methods in this paper for explaining the gravity rates are preliminary, but are sufficient to show the overall scenarios of the proposed processes.At the sites with few AG measurements and large hydrological gravity effects, the estimated gravity rates contain large uncertainties, which can be reduced as more gravity measurements are collected and improved correction models are developed.Therefore, we recommend continual AG data collections at the existing sites and at potentially expanded sites.
τ 2 = 10 to 60 days, and the resulting RIGE values were compared with the residual gravity values from SG #T48 between 1 April 2006 and 31 December 2010.Figure 2 shows the comparison of the RIGE and SGderived values (the gravity values in the second-stage).The three (τ 1 , τ 2 ) pairs in Fig. 2 result in correlation coefficients of -0.53, -0.86, and -0.77.Note that, because SG #T48 is in a tunnel 40 m below the surface, the correlation coefficients between the RIGE values and the second-stage SG gravity values are negative.The largest correlation coefficient (absolute value) of 0.86 led to the best pair of τ 1 = 2.2 and τ 2 = 41 days at HCHG.At other AG sites, we adjusted τ 1 and τ 2 values so that the resulting RIGE values best matched the measured gravity changes.It turns out the determination of the best (τ 1 , τ 2 ) pair at a given AG site can be very empirical and this is an issue needing further studies.
), in connection to typhoon Morakot.The residual gravity values collected in 2008 and 2009 differ by 53.37 μGal (AG3, value in 2009 minus value in 2008) and 23.38 μGal (AG6), which were due to the large landslides and debris flows along Sinwulyu and Laonong Rivers associated with typhoon Morakot.Mouyen et al. ( Fig. 4. (a) The terrain around sites AG3 and AG6, at altitudes of 383 and 411 m.(b) A 10 m-thick debris deposited near AG3 along Sinwulyu River, caused by typhoon Morakot.(c) Terrain change around AG3 shown in two images of satellite FORMOSAT-2.(d) The riverbed coverages around AG3 before (dark gray) and after (light gray) typhoon Morakot, and the dashed red circle shows the effective radius (0.5 km).(Color online only) Fig. 6.(a) and (b) Coseismic gravity changes and vertical displacements induced by the 2006 Pingtung earthquake doublet from the dislocation model of Sun et al. (2009b).(c) Gravity change from 26 February 2006 (before) to 7 February 2007 (after).The interseismic gravity rate at KDNG ( g res o in Table 2) is used to estimate the coseismic gravity change of 3.12 μGal on 26 December 2006.(d) The GPS-derived latitude, longitude and vertical displacement are 0.62, -2.05, and -2.10 cm at KDNM on 26 December 2006.(Color online only)

Fig. 8 .
Fig. 8.The geometry explaining gravity change due to magma chamber change.The gravitational attraction due to the change of height Δh of a cylindrical magma chamber with a radius of R. (Color online only) 55 cm yr -1 from GPS in central Taiwan are larger than those in Tibet.The different Moho deepening rates in Tibet and Taiwan are attributed to many factors, such as crustal age, size, and mountain building process.In summary, our time-varying gravity observations at SMLG shows the first evidence of Moho deepening in Taiwan, conforming to a deepening hypothesis presented by, e.g., Yamato et al. (2009), Wang et al. (2010), and Kuo-Chen et al. (2012).

Fig. 10 .
Fig. 10.Moho deepening (brown line) under SMLG and FLNG.FLNG is also affected by the subduction of the Phillipine Sea Plate: the uplift at FLNG (purple line) is offset by the plate suduction (blue dash line) to cause a negative vertical displacement rate (green line).The red circles are the resulting surface locations due to uplift only (SMLG), and due to a combined uplift and plate subduction (FLNG).The red triangles indicate the positions near the mantle estimated from the residual gravity values associated with Moho deeping [Eq.(14)].(Color online only)
. Absolute gravity, GPS and rain gauge sites over related geologic settings.

Table 2 .
Categories and types of error absolute gravity measurements.In the "Source" column, we illustrated the uncertainty values for the gravity data obtained at the LOGG gravity site in November 2008.

Table 3 .
4 × 10 6 m 3 given in the 2009 final report of Ministry of Economic Affairs (MEA; ISSN 0257.5442).An ongoing study uses 3D laser scanning and river stage detecting device to precisely estimate terrain changes and mass changes near AG6.Such measurements will result in more rigorous estimates of sediment volume changes.The sources of such discrepancies in the riverbed volume changes around AG3 and AG6 are analyzed below.Gravity rates at AGTO sites with at least three repeated measurements.
Note: *1: First-stage gravity rate: based on residual gravity values that are free from solid earth tide, ocean loading, air-pressure effect and polar motion effect.See Eq. (3).*2: Second-stage gravity rate: the rate from first-stage gravity change minus vertical motion-induced gravity change, See Eq. (4).*3: Residual gravity rate: final gravity rate without the effect of hydrological effects.See Eq. (5).*4: The AG was measured only once or twice at a site.*5: AG data in 2006 is not used.*6: AG data in 2005, 2006, and 2007 are not used.*7: AG data in 2009 and 2010 are not used.*8: AG5 was destroyed by typhoon Morakot in 2009 August.There is no AG data in 2009 and 2010.*9: AG data in 2006 is not used.

Table 4 .
Parameters for the finite-source rupture model in the two mainshocks of the Pintung earthquake, 26 December 2006.