Seasonal Variations in Hydrological Influences on Gravity Measurements Using gPhones

Hydrological influences on a local gravity field may reach amplitudes on the order of 10 microgals. Since 2007, fifteen Microg LaCoste gPhones have been successively installed in gravity stations in China. The outputs from gPhones include ten data channels with second sampling such as raw gravity, corrected gravity, long level data and cross level data, ambient and sensor temperature, ambient and sensor pressure, and others. In this study, we select six stations in northwest China (GaoTai, LaSa, LanZhou, ShiQuanHe, WuShi, XiAn) and one station in the northeast (HaiLaEr). We have modeled the major tides (earth solid tide, ocean tide and pole tide), corrected for atmospheric loading effects using local measurements, fitted instrumental drift using segmental fitting based on the distinct characteristics of gravimeter drift, and ultimately obtained the monthly residual gravity with amplitudes of 10 ~ 20 microgals. We find that the results obtained by the gravimeter for those stations with stable conditions and no large disturbances are obviously correlated with hydrologic loading as modeled by the Global Land Data Assimilation System and Climate Prediction Center. We also notice that at some stations there are obvious phase lags with a period of three months or more between the residual gravity and the influence of hydrological loading. These large discrepancies may be associated with local hydrologic effects, local topography or some other complex tectonic movement and geodynamical mechanism, which were not considered in this paper.


IntrodUctIon
Most geodetic observations show large annual or seasonal variations which are likely due to changes in continental water storage (Boy and Hinderer 2006).Water supply and storage in northwestern China are very important for local economic development; however, due to the complex processes of infiltration and redistribution to ground water (Harnisch and Harnisch 2006), there is little information on the variations of continental water storage and especially local ground water storage.At present several studies have shown that hydrological influences on a local gravity field may reach amplitudes on the order of 10 microgals and have obvious seasonal variations (Kroner et al. 2004).Most of these studies are based on Superconducting Gravimeters (SG) used in the GGP project.SGs have high sensitivity for long-term stability, but the number of these instruments worldwide is limited, so they can only detect hydrological effects at an isolated site.Another type of gravity instrument is the Microg LaCoste gPhone which have been used in many sites in China and have the capability of recording local or regional gravity variations.The gPhone meter has high precision for tide observations but performs poorly in non-tidal bands due to the high instrument drift.Many methods have been used to simulate the drift (Li et al. 1983;Chen et al. 2007).For this study, the least squares method was used to remove the effect of instrument drift.In the future, we will use absolute observations to calibrate the drift to obtain more data where the results are more scientifically valid.
At present, 39 gravity stations have been built in China and are distributed in 20 provinces (Fig. 1).Most Terr. Atmos. Ocean. Sci., Vol. 22, No. 2, 157-168, April 2011 of them are located in the southeast region where frequent earthquakes occur.The average distance between stations is now around 450 km, which will be much smaller when the Crustal Movement Observation Network of China II (CMONOC II) accomplished.During a project of Tenth Five-Year Plan Period (TFYPP) we have introduced fifteen Microg LaCoste gPhones (gMonitor User's Manual 2008; gPhone/PET Hardware Manual) which convey raw gravity data to the Gravity Network Center of China (GNCC) in Wuhan through Internet in real time.For this study, we select six stations in the northwest and one in the northeast (Fig. 1).The operating conditions of each station are shown in Table 1.
In this study, we show the processing of gravity data as well as our models of long-period gravity changes, including the effects of solid earth tide, ocean tidal loading, pole tide, instrumental drift, atmospheric loading and continental water storage loading.We also compare the gravity residuals with global continental water storage loading.

dAtA ProceSSInG And PrelIMInAry re-SUltS
The raw data of gPhone downloaded from GNCC database are first preprocessed.In order to quantify the seasonal variations induced by continental water storage changes  with gravity measurements, the factors must be modeled at seasonal timescales including (Wahr et al. 1998): solid earth tides, ocean tidal loading, pole tides and atmospheric loading effects.The amplitude of non-tidal ocean loading is less than 10 nms -2 (Zhou 2008), therefore non-tidal ocean loading effect is ignored in this study.The instrumental drift is removed by segmental polynomial fitting.We then calculate the continental water storage changes using global hydrology models, and compare estimates of storage with the gravity residuals.

Preprocessing
The second raw data are packed monthly and each data pack has a one-day lap with data from the adjacent day to provide continuous data after some filtering.Raw data are first changed to minute sampling data using a low-pass filter with a cut-off frequency of 720 cpd (cycle per day), and then corrected for undesired artifacts such as spikes, earthquakes, steps and gaps using the Tsoft software package (Vauterin 1998).Finally, those data are re-sampled to hourly values using a low-pass filter with a cut-off frequency of 12 cpd.

Solid earth tides
For each station considered, a multi-channel tidal analysis [program ANALYZE (Wenzel 1996)] is used to find the tidal parameters for 20 major tidal wave groups (not including the SA and SSA wave group due to the length of raw data) whose precision is represented by a standard deviation in the frequency domain and regression coefficient of the local air pressure in microgal/millibar.The local solid earth tides are modeled using tidal parameters and the results of tidal analysis are shown in Table 2.

Pole tide
The pole tide, induced by polar motion and lengthof-day variations, is modeled using daily Earth orientation parameters provided by the International Earth Rotation Service (ftp://hpiers.obspm.fr/iers/eop/eopc04/).The formula is specified in the IAGBN Absolute Observations Data Processing Standards (1992).The gravimetric factor is supposed to be equal to 1.16 (no frequency or latitude dependence).

Instrumental drift
The gravity sensor of the gPhone system consists of a low drift, metal, Zero Length Spring system whose characteristic is determined by the mechanics of elastic material (gMonitor User's Manual 2008).Instrumental drift is modeled perfectly by simulating mechanics of spring system taking into account the factors such as the elastic lag, creep deformation, elastic after-effect, age of spring, the temperature of inner sensor, and features of the stress loading on spring and so on.Due to its complexities we choose an easier but more effective method to fit instrumental drift.
Each gravity meter has unique characteristics induced by instrumental drift, wherein its varying trend will be largely changed after some disturbance (such as a power interrupt, re-start, lightning, etc.).For long time periods, if no disturbance occurs, the drift will be regular so that is able to perfectly fit using a simple polynomial.Those data influenced by serious disturbances need to be treated using different methods: if the time series become regular after a short period of interruption, for example in one week, those data are deleted; if the gravity data series take long time to come back to regular, special functions (such as polynomial or exponential function) should be used to fit the data.
The process of fitting instrumental drift is as follows.First, according to the characteristics of varying trends of gravity residuals time-series (gravity residuals after correc-tions of solid earth tide, ocean tidal loading and pole tide), the entire time-series are set into segments, then segmental polynomial fitting is applied to each segment, and finally combine together to compensate for long term meter drift.The polynomial fitting curve of HaiLaEr station is shown in Fig. 3 in detail.The segmental time interval, fit degree, polynomial fitting expression and residual standard deviation of all stations are shown in Table 3.

Atmospheric loading effects
According to section 2.2, the influence of varying air density in the atmosphere is eliminated by a linear regres- sion model using the local air pressure variations at the gravimeter site and a regression coefficient derived by the standard tidal analysis (shown in Table 2) (Crossley et al. 1995).More details of the air pressure model (e.g., regional air pressure distribution, deviations from the standard atmosphere, 3D global atmosphere model) are not included here.

Hydrology Models and continental Water Storage loading
The CPC model accounts for soil moisture in one layer with an effective holding capacity of 76 cm of water, which is constant over the space domain.The input consists of temperature and precipitation data.The output of the model is on global 0.5° × 0.5° soil moisture grids.The soil moisture data set is compared with in situ measurements, and it is found that the seasonal and annual variabilities are simulated reasonably well (Fan et al. 2003;Fan and van den Dool 2004;Schmidt et al. 2008).
The GLDAS model aims to produce high-resolution land surface states uncoupled from any atmospheric model (Rodell et al. 2004).In contrast to the CPC model, the GL-DAS model uses observation-derived surface parameters that can vary in time and higher-resolution data for surface parameters, such as vegetation class.The output is given as soil moisture in four layers, snow water equivalent and canopy water storage (Schmidt et al. 2008).The model output used in this study is based on the Noah land surface model (Chen et al. 1996) and is averaged monthly from the original 3-hour data.
We model the hydrologic loading as a thin-layer process acting on a spherical Earth surface, and calculate loading gravity effects using Green's Functions (Longman 1962;Farrell 1972).Note that the sign of the local contribution of modeling of the direct Newtonian attraction varies under different conditions: if the station is located below the ground, an excess of water mass will decrease the gravity and the opposite if the station is located at the surface (Boy and Hinderer 2006).In this study all stations are in caves so the situation is complex wherein a heavy rain occurs, for the first few hours or days the mass of water remains above the gravimeter so the gravity decreases, but as the water flows down to the underground, the opposite occurs (Matsumoto et al. 2000;Abe et al. 2006).Moisture from snow shows similar behavior before and after thawing.Due to the overburden thickness of caves is greater than 20 m, we consider all seven gravimeters as below the soil horizon and reverse the sign of the local contribution to simulate the effect.Hydrologic loading modeling in a more detailed manner would require precise topography for computing the direct attraction by nearby water masses, as well as local models such as precipitation data, underground data, and soil moisture at stations, which are not available for this study.
As an example, Figs. 2, 3 and 4 show the successive steps of data processing described above using the data series recorded with the gravimeter PET0030 from September 2007 to August 2009 at HaiLaEr station in Inner Mongolia province.The comparison between gravity residuals after the above corrections and modeled hydrology loading from the GLDAS and CPC models at other six stations are plotted in Fig. 5.
From Figs. 4 and 5, we find that comparing gravity residuals at seven stations mentioned above with hydrological loading effects from GLDAS and CPC models differs.For some stations there is a reasonable correlation between the loading and the residuals in both amplitude and phase (HaiLaEr, ShiQuanHe, GaoTai, LanZhou).Some stations have good agreement in amplitude but a minor phase lag (XiAn), or a good match in phase but the amplitude differs (WuShen).At LaSa station the estimated loading isn't coherent with the gravity observations.The gravity residual (in days) does not agree with the hydrology loading from GLDAS and CPC (both in months) and the former is not smooth as the latter (Fig. 5).

ShiQuanHe
This station is located in far western part of Tibet.The annual precipitation is about 69 mm.The nearest river is ShiQuan river, 6 km away from the station, the landscape around the station is farmland and there are few inhabitants in the area (Shang et al. 2006).Gravity data is collected automatically using the internet and thus data from ShiQuan-He have little man-made noise.As shown in Fig. 6, we find good agreement between the amplitude of gravity residuals and hydrological loading.Therefore the main component of seasonal variation results from hydrological influence.The poor agreement during July to December 2009 is related to  a large instrument interruption which caused a change in the trend of drift.The phase between gravity residuals and hydrological loading has a time lag of about 3 months.This may be associated with complex tectonic movements and geodynamical mechanism in the Tibetan Plateau, such as mass loss beneath the Tibetan Plateau (Sun et al. 2009) or the earthquake-prone areas along the Himalayas.

Gaotai
This station is located on the northern edge of the Tibetan plateau.There are many sources of disturbances.The  town of GaoTai is 5 km from the station; the town road is just 100 m away, along with several quarries scattered about approximately 1 km away.All the sources mentioned above lead to interruption of the normal measurements.Low-pass filter is applied to eliminate the noise in most of the sources which cause high frequency noise on gravity data.In Fig. 7, good agreement in the trend change of the gravity residuals and hydrological loading while a two or more month long phase lag is seen.

Hailaer
PET0030 is recording at HaiLaEr station with stable observation conditions and fewer man-made disturbances.This has an average annual precipitation of 352.2 mm and evaporation of 1250.0 ~ 1533.7 mm.This station is 0.5 km from HaiLaEr River.The seasonal variation of this river is distinct; the maximum flux is in summer and reached 7.02 billion m 3 in 1984; during its icing period in winter the flux decreases to 10 percent of the maximum.During January and February, the minimum flow of 1.22 billion m 3 was observed in 2002.The obvious seasonal changes of hydrology loading model of GLDAS provide validation for this application (Fig. 8).Comparing residuals with the hydrological contribution from GLDAS, a reasonable correlation in amplitude and phase, especially after May 2008 is found.This agreement is due to the stable condition of PET0030 whose drift is regular even after a severe earthquake in May 2008.The hydrologic loading plays a leading role in the seasonal changes of gravity measurements at this station.

XiAn
This station is located in the southern Ordos platform.The average annual precipitation is 550 ~ 880 mm and the average evaporation is 900 ~ 1150 mm.Due to the serious shortage of surface water, ground water has been excessively exploited which has caused 160 km 2 ground subsidence.In the left bottom frame of Fig. 5 an obvious phase lag between hydrology loading model from GLDAS and CPC can be seen.In Fig. 9 we plot all of the three kinds of data.It shows that the phase of gravity residuals agrees well with CPC model while it has a phase lag with GLDAS model, in that the model of GLDAS aims to land surface state while ground water is the leading factor on continental water storage variations at XiAn station.The amplitude of residuals is larger than the gravitational effects of hydrological loading in October, November and December 2008, which may be associated with the man-made preprocessing in the course of correcting several artifacts, or due to some local effects that are not considered in this study.

lanZhou
This station is one of the integrated geophysical observatories in China.It is located 1.5 km from the Yellow River, the second longest river in China.The nearest town is 5 km far away.LanZhou has average annual precipitation of 318.6 mm, groundwater level depth of about 5 ~ 30 m, and aquifer depth of more than 5 m.In Fig. 10, the variation of the residuals and hydrology loading is coincident to some extent, but a strong mismatch both in amplitude and in phase exists.The data for May, with a large decrease may be related to the WenChuan earthquake on 12 May 2008 and a series of aftershocks.The epicenter of the earthquake is only 557 km away from the station and both WenChuan and LanZhou are located in the South-North Seismic Zone.

WuShen
This station is in the southern TianShan seismic zone, Fig. 7.The gravity residuals and modeled hydrology loading from GLDAS model at GaoTai.where several earthquakes with magnitude 6.0 or greater have occurred since 1970.Broad field is around the station so there is no obvious man-made noise.Figure 11 shows that the variation of the residuals and hydrology loading are in good agreement but the amplitude of the gravity residual is about 4 to 6 uGal larger than the amplitude of hydrologic loading from GLDAS.A phase lag of about 3 months also exists.

laSa
Annual sunlight-hours in LaSa city exceeds 3000 hours, thus it is also called "the sunlight town."LaSa has an annual water surface evaporation of 2171.1 mm and annual precipitation of 472.5 mm (Fan et al. 2005).From the hydrologic loading of GLDAS (Fig. 12b) it can be seen that the amplitude of variations is only about 3 uGal.The gravimeter measurements began at this station in January 2008, and data are available for 15 months.In processing the data, we find the gravity residuals after tidal analysis using program ANALYZE have a different noise level before and after 12 May 2008 (Fig. 12a).The former contains obvious tidal components while the latter do not.We delete the data for the period of January to May 2008.Even so, the quality of data after May 2008 is not sufficient due to serious interruptions of the gravimeter, such as lightning strike on 23 June 2008, and several gaps induced by human disturbances.All of these factors might be the reason that there is a discrepancy between gravity residuals and hydrology loading from GLDAS (Fig. 12b).

conclUSIon
Research has shown that gravity observations present obvious annual or seasonal variations which are likely due to the changes of continental water storage (Wahr et al.  1998; Kroner et al. 2004;Boy and Hinderer 2006;Harnisch and Harnisch 2006).In this study, we have shown that the gravity series measured by gPhones contains the contribution from the hydrologic changes at seasonal timescales.We modeled three tides (earth solid tide, ocean tide and pole tide), and corrected for atmospheric loading effects using local measurements.We also fitted the instrumental drift using a segment fit approach based on distinct characteristic of gravimeter drift, and obtained monthly gravity residuals with an amplitude of 10 ~ 20 microgals.The gravimeter results from stations with stable observations and no large disturbances have an obvious correlation with the hydrologic loading modeled by CPC and GLDAS.We found that at some stations there were obvious phase lags with a period of three months or more between the gravity residuals and hydrology loading.Gravity measurements on land were strongly influenced by the local conditions such as local topography, geographical environment, climate, relative position of the gravimeter with respect to the ground and other (Kroner et al. 2004).The effects induced by local factors are difficult to model precisely; more research is needed.A more precise atmospheric pressure model both in local and global (2D and 3D) will be considered in future work.Instrument drift will be calibrated with absolute measurements in a more systematic approach.In the future, we will utilize additional 37 gravimeters to the Gravity Station Network of China, and the total number will reach approximately 60 in 2011.

Fig. 1 .
Fig. 1.The distribution map of gravity stations in China.

Fig. 4 .
Fig. 4. Gravity residuals after successive corrections at HaiLaEr station.From top to bottom, respectively: (1) station pressure measurements, (2) hourly residual after correction of solid earth tides, ocean tidal loading, pole tide, instrumental drift and atmospheric loading effect, and (3) daily residuals after the same correction as in (2), monthly hydrology loading from GLDAS and CPC.

Fig. 3 .
Fig. 3. Polynomial fitting of HaiLaEr station.From top to bottom, respectively: (1) gravity residuals (after correction of solid earth tides, ocean tidal loading and pole tide) and segmental polynomial fitting, and (2) residuals after the polynomial fit.The data series of around two years are used (September 2007 ~ August 2009), which was interrupted by an earthquake in May 2008.The varying trend of drift is changed obviously before and after this interruption.Therefore fitting the drift for the segmental time interval of 07.9.6 ~ 08.6.1 and 08.5.1 ~ 09.2.25 (the first two segments in the upper frame) respectively is reasonable.And each segment has about 30 days lap with adjacent one to remove border effect.On 26 February 2009 a power interrupt happened which lead a data gap and change of drift trend, so the single polynomial fitting is supplied in this time interval to simulate this segment.After the segmental fitting, RSD of gravity residual reaches to 4.44, 2.69, 3.86 uGal in the same order.

Fig. 5 .
Fig. 5.The gravity residuals of gPhone data and modeled hydrology loading from GLDAS and CPC model.

Fig. 6 .
Fig. 6.The gravity residuals and modeled hydrology loading from GLDAS model at ShiQuanHe.RMS of difference between gravity residuals and hydrology loading is also presented in the figure (below as the same).

Fig. 9 .
Fig. 9.The gravity residuals and modeled hydrology loading from GLDAS model at XiAn.

Table 1 .
The observation condition of stations.

Table 2 .
The results of tidal analysis of 7 stations.

Table 3 .
The instrumental drift fitting of gPhone.
* Each segment has about 30 days lap with adjacent one to remove border effect.** RSD: Residual Standard Deviation.