SGOTL: A Computer Program for Modeling High-Resolution, Height-Dependent Gravity Effect of Ocean Tide Loading

SGOTL, a computer package coded in FORTRAN, has been developed to model the gravity effect due to ocean tide loading (OTL), especially for a coastal station with large ocean tides. SGOTL uses a regional and a global tide model to account separately for near (inner) and far (outer) zone contributions, and optimizes an inner-zone region and grid interval for numerical convolution. Height dependent Green’s functions for Newtonian and elastic effects are developed. The coastline is defined by the full-resolution GMT shoreline, and optionally a digital elevation model (DEM). A case study using gravity observations at the Hsinchu superconducting gravity station and some offshore islands around the Taiwan Strait suggests that SGOTL outperforms some selected global OTL programs and achieves an accuracy of 0.1 μGal for 8 leading tidal constituents.


InTRODuCTIOn
Gravity effect due to ocean tide loading (OTL) is needed for correcting gravity measurements, assessing ocean tide models and validating the theory of tidal loading (Sato and Hanada 1984;Shum et al. 1997;Yamamoto et al. 2001;Baker and Bos 2003;Boy et al. 2004;Neumeyer et al. 2005;Hwang et al. 2009).There are many OTL programs in the earth science community which serve such needs, e.g., GOTIC (Sato and Hanada 1984), OLFG/OL-MPP (Scherneck 1991), SPOTL (Agnew 1996), NLOADF (Agnew 1997), CARGA (Bos and Baker 2005), GOTIC2 (Matsumoto et al. 2001) and g7.0 (g7 User's Manual 2007).With a typical spatial resolution of about 30' × 30' in latitude and longitude, these programs are most effective for global applications and may have a degraded accuracy and resolution when a gravity station is situated in a special location, e.g., a coastal area with complicated coastlines, high station elevation and large ocean tides.This was shown in a study of OTL corrections for superconducting gravity (SG) measurements at the Hsinchu station by Hwang et al. (2009), who found that some of the global OTL programs produce anomalously large values that cannot be used for gravity corrections.The Hsinchu SG station is about 8.6 km to the Taiwan Strait and situated at an elevation of 78 m.Earlier studies suggest that gravity effects and displacements due to OTL can be very large at offshore islands with large tidal amplitudes (Lambert et al. 1998;Yamamoto et al. 2001;Huang et al. 2008;Yeh et al. 2008;Hwang et al. 2009).
This study is motivated by the need for an accurate OTL program that can treat stations at high elevations near a coastline, anywhere in the world.The OTL computer program is named SGOTL (SG ocean tide loading) to emphasize the OTL accuracy required in many applications of superconducting gravimetry.For example, time-varying gravity changes detected by GRACE satellite gravimetry (Tapley et al. 2004) are typically at the μGal (10 -8 m s -2 ) level.The ground validation of such gravity changes by superconducting gravimetry will require an OTL model accurate to the same level.Other examples of high-precision applications of SG can be found in Crossley and Hinderer (2009).OTL corrections in such applications are critical to achieving optimal results.

The newtonian Effect
Here we summarize the theories and formulae for computing the gravity effect of OTL for SGOTL.There are several variants of such formulae, and the formulae presented in this paper are an extension of previous publications such as Farrell (1972), Melchoir (1983), Agnew (1997) and Bos and Baker (2005).In particular, the formulae of height-dependent loading of Green's function are presented for the first time.The gravity effect of OTL consists of two parts: (1) vertical attraction (acceleration) by the ocean tidal mass (the Newtonian effect), and (2) gravity change due to deformation caused by loading of the ocean tide mass (elastic effect).The Newtonian effect is the vertical gradient of the potential change caused by ocean tide mass.As shown in Fig. 1, at a gravity station p, the potential due to the total mass of ocean tides is where pq } : spherical angle between p and a running point at q(φ q , λ q ) sea located at with φ q , λ q being latitude and longitude h q : tidal height above the local mean sea surface w t : density of seawater G: Newtonian constant, which can be replaced with a combination of g and M (the mean acceleration and mass of the earth) R e : the mean radius of a sphere osculating the geoid near the station D: the entire domain of a spherical earth; see section 3 for the subdivision of D , and the differential surface area on the sphere is dS R d q e q 2 v = r p : geocentric distance of p, and is approximated by r p = R e + H p , with H p being the elevation of p (defined as the orthometric height above the local geoid) In the integral of Eq. ( 1), a "spherical approximation" has been used.This spherical approximation is necessary because the tidal height is counted from some "mean sea level" around the gravity station of interest.The underlying "sphere" will osculate the mean sea level or the local geoid.The radius of such a sphere, R e , is computed using the Gaussian formula given in Torge (1989, p.98), result-Fig. 1. Geometry showing a gravity station p around oceans for the evaluation of ocean tide (h q ) gravity effect.See text for the meanings of the symbols.
ing in the Gaussian mean radius.Because the distance of the point of the ocean tide ("q" in Fig. 1) to the geocenter is already approximated by the Gaussian mean radius, the station distance to the geocenter should be approximated in the same way.That is, the sum of the orthometric height and the Gaussian mean radius, r p , is regarded as the geocentric distance of the station.Such an approximation is also used in the height-dependent Green's function in section 2.2.The density of seawater can be regarded as a constant and we adopt the value of w t = 1030 kg m -3 .The Newtonian effect is then the negative vertical gradient of T expressed as where the dimensionless Green's function is Consideration of station height H p in the Newtonian effect is necessary because the station can be placed along a coastal, high-altitude location.To see the effect of a station height, we computed the Newtonian effects at the Hsinchu SG station (latitude = 24.7925°Nand longitude = 120.9855°E)and two proposed SG stations in Taiwan: One is at GPS station YMSM (latitude = 25.1657°N and, 121.5740°E) and the other at Lulin (latitude = 24.4710°Nand longitude = 120.8808°E) in the Central Range of Taiwan.The nearest distances to the sea for Hsinchu, YMSM and Lulin are 8.6, 10.1 and 80.0 km, respectively.Figure 2 shows the Newtonian effects (due to M 2 ocean tide using SGOTL) with respect to station heights.For Hsinchu, the Newtonian effect increases with the station height.The M 2 amplitude of the Newtonian effect at the Hsinchu SG station (station height: 78 m) is about 0.75 μGal from SGOTL (see section 4.1), which is close to the amplitude at sea level (0.70 μGal).If the Hsinchu station were to be at a nearby site with an elevation of 500 m, the error due to ignoring the station height in the modeled Newtonian effect (M 2 ) is about 1.60 -0.70 = 0.90 μGal, which is about the accuracy of the mean value from a 30-minute set of FG5 absolute gravity observations.The dependence on height at YMSM is similar to that in Hsinchu; the actual height of YMSM is 780 m, where the M 2 Newtonian effect is 2.4 μGal.If we assume the height of YMSM is zero, the M 2 Newtonian effect is 0.8 μGal, resulting in an error of 2.40 -0.80 = 1.60 μGal.On the other hand, the Lulin station is relatively distant (80.0 km) from the sea and the Newtonian effect (about 0.70 μGal) is insensitive to station height.In conclusion, because of the Newtonian effect, the height of a coastal gravity station must be taken into account in order for an OTL program to achieve a one-μGal accuracy.

The Elastic Effect
The loading of the ocean tide-induced mass deforms the earth's crust and then changes gravitational potential.The vertical acceleration due to such a change of potential is the loading gravity effect.Standard treatises such as Melchior (1983), Moritz and Muller (1987) and Farrell (1972) have presented detailed derivations of the loading gravity effect.The elastic effect is a convolution of a Green's function with the ocean mass re-distribution caused by tidal fluctuations of the sea surface.Like the Newtonian effect, in this paper we present a height-dependent Green's function for the elastic effect of ocean tide.According to Moritz and Muller (1987) and referring to Fig. 1, the Green's function for vertical displacement due to ocean tide loading is and the Green's function for the change of potential due to tide-induced deformation of the earth is where hn l and kn l are the loading Love numbers of degree n, and cos Pn p q } ^h is the Legendre function.The vertical displacement and the deformation-induced potential change can be experessed in the convolutions: If we approximate the gravity gradient, / g H 2 2 , by Moritz and Mueller 1987), the gravity change due to the vertical displacement and the deformation-induced potential change is where and + ^h .A similar formula for atmospheric loading effect depending on height can be found in Guo et al. (2004).Figure 3 shows the difference between the elastic effects (due to M 2 ocean tide) at sea level and at a given station height as a function of height for the Hsinchu, YMSM and Lulin stations.The distance of a station to the coastline will also change the elastic effect: the M 2 elastic effects at sea level are 3.1, 2.8 and 1.5 μGal for Hsinchu, YMSM and Lulin stations, respectively.The relative errors in the elastic effect at Hsinchu and YMSM due to discounting station height are about 0.27% and 0.9% of those at a station height of 500 m, which amount to few nanogals.Compared to the Newtonian effect, the elastic effect is therefore much less dependent on station height.For an island (remote from the sea) station such as Lulin, the relative error due to discounting station height is smaller than that at a coastal station such as Hsinchu.The results in Figs. 2 and 3 confirm that a remote gravity station is less affected by OTL than a coastal station, and this is consistent with other studies, e.g., Boy et al. (2003).Since the noise level of a SG instrument may be at the nanogal level at the frequency domain, for SG applications to such studies as free oscillations and tides, the station height cannot be neglected when computing corrections for the OTL gravity effects.
The combined Newtonian and elastic effects are called OTL for short, and can be expressed as where is the total Green's function.For the OTL modeling in this paper, discrete values of K T for a given station height were pre-computed at an 0.01°interval from } = 0° to 180°.The needed K T value at the numerical convolution of OTL was then linearly interpolated from the nearest two K T values.In this paper, we adopted the loading Love numbers of Farrell (1972) to compute the Green's function Gn } ^h in Eq. ( 11).User-defined loading Love numbers, for example, the PREM-based Love numbers (Mathews 2001), are also acceptable in our program, see Appendix A.

Amplitude and Phase of OTL
A brief description of the computational principle used in SGOTL is presented below.A constituent of ocean tide and its gravity effect depend on both space and time.They share the same frequency but have different phases.In the SGOTL computer program, the amplitude and phase of a given OTL constituent at a given location are computed.The actual OTL values are then computed from the amplitude and phase.A tidal height h can be expanded into a time-dependent Fourier series as (Foreman and Henry 1979) where , A j h z m ^h and , g j h z m ^h are the location-dependent amplitude and phase, f t j ^h and u t j ^h are the nodal modulation terms, and j is the tidal constituent with V t j ^h being its astronomical argument.Given a value of V t j 0 ^h for j at a reference epoch t 0 , its value at any time t can be computed as , where w j is the frequency of the j tidal constituent.In the following development, we will use V t V t u t j j j = + ^^ĥ h h for simplicity.Equation ( 14) can be expressed as where , , , , , Substituting h in Eq. ( 15) into the integral of Eq. ( 12) and using the total Green's function in Eq. ( 13) leads to Thus, OTL is a function of time and geographic location, and can be expanded into a Fourier series in exactly the same way as ocean tide: where the amplitude A j G and g j G phase can be computed as

nuMERICAL COnvOLuTIOnS Of OTL AnD PROGRAM SuMMARy
The integration in Eq. ( 18) was carried out numerically to obtain the amplitudes and phases of OTL at any location specified by latitude, longitude and station height (orthometric height).The integration uses a regional tide model for the inner zone effect and a global tide model for the outer zone effect.The inner zone is defined as a rectangular box bordered by two meridians (west and east) and two parallels (north and south).The outer zone covers the entire globe from 0° to 360°E longitude and -90° to 90° latitude, excluding the inner zone.The grid interval of the outer zone is always equal to the nominal spatial resolution of the adopted tide model, i.e., 0.5° for CSR4.0 and NAO.99b.The land and ocean are based on the full-resolution coastline of Generic Mapping Tool (GMT) determined with the command "grdlandmask" (Wessel and Smith 1999), and optionally a digital elevation model (DEM).One source for the DEM is the Shuttle Radar Topography Mission (SRTM; http://www2.jpl.nasa.gov/srtm/).The outputs of grdlandmask and the reading program of the given DEM are files containing 0 for land and 1 for ocean.The implementation of SGOTL and its auxiliary programs with examples are given in Appendix A.
The numerical integration was carried out using the Gauss quadrature method (Press et al. 1989).This method has been used for a rigorous computation of terrain effect, and proved to be a highly efficient and accurate numerical integration method (Hwang et al. 2003).We experimented with several grid intervals and box sizes for the inner zone.In all, nine cases were tested at the Hsinchu SG station, as listed in Table 1.It was found that the OTL program from Case 2 agrees the best with the SG result at Hsinchu (section 4.2).To see the program accuracy and computational efficiency relative to those of Case 2, Fig. 4 compares the error ratios and computing time ratios of all cases.A time ratio is the ratio between the CPU time of a case and that of Case 2. An error ratio is defined in the same way, but for the errors of the OTL program for the M 2 constituent.It turns out that only Cases 3 and 6 are marginally faster than Case 2; Cases 3 and 6 result in degraded program accuracies (the error ratios are greater than 1).Since Case 2 delivers the best accuracy with only a minor increase in computational time, we used the box size and grid interval of Case 2 for SGOTL.
The computer program developed in this paper is termed SGOTL for further discussion.In summary, the SGOTL program features: (1) Height-dependent Green's functions for the Newtonian and elastic effects; (2) Gauss quadrature for numerical convolution; (3) Loading Love numbers from Farrell (1972); user-defined Love numbers are acceptable; (4) A regional and global tide model for the inner zone and the outer zone effects; (5) High-resolution coastline from GMT and/or a DEM; (6) Expected accuracies of better than 0.1-μGal in amplitude and 10° in phase for any tidal constituent.

Program Development
Here we demonstrate an example of OTL modeling by SGOTL around the Taiwan Strait.The regional tide model by Hu et al. (2010) and the NAO.99b tide model were used for the inner zone effect and the outer zone effect, respectively.NAO.99b is a global model that is shown to be the most accurate in the western Pacific region, partly due to the assimilation of tidal records in the model (Hwang et al. 2009).The amplitudes and phases of the tide models of Hu et al. (2010) models are given on a 5' × 5' grid.The coastline of Taiwan was defined by an 80 × 80 m DEM, provided by Ministry of the Interior, Taiwan.This DEM model was derived by aerial photogrammetry in 2007, with a decimeter accuracy for coastal areas of Taiwan.The new Taiwan DEM grid and its documentation are available by application at https://dtm.gps.moi.gov.tw/dtm/dtm/index.aspx.Outside of Taiwan, the full-resolution coastline of GMT was used.The batch job to produce the OTL gravity effects is given in Appendix A.
To see the pattern of OTL variation, we generated the amplitudes and phases of an M 2 OTL at sea level around the Taiwan Strait on a 5' × 5' grid by program SGOTLG (see Appendix A); their contours are given in Fig. 5. SGOTLG is a special version of SGOTL that can generate amplitudes and phases on a regular grid.Figure 5 shows a band of amplitude highs of M 2 OTL along the southeast China coasts at latitudes > 23.5°.On land, the ocean tide is zero by definition, but the OTL gravity effect can differ from zero.In general, OTL is the largest near the coast and decreases toward inland.As such, over Taiwan, the OTL gravity effect is the smallest along the Central Range.Figure 5 also shows an interesting phenomenon; the phase contours of M 2 OTL peak at a zone near Kinmen and form concentric circles around this offshore island.The phases show that both the waves of ocean tide and OTL gravity effect propagate from the Pacific Ocean towards the coasts of southeast China and Taiwan.However, the two waves respond to the bottom topography and coastlines differently, leading to different phase distributions.

Model Assessment by SG and AG Gravity Data
The model results of SGOTL (for the amplitude, phase and time series, see below) around the Taiwan Strait were assessed by superconducting gravimeter (SG) measurements at the Hsinchu SG station and absolute gravity measurements at offshore island stations around the Taiwan Strait.The distribution of these stations is given in Fig. 5.The Hsinchu SG station is now included in the SG network of the Global Geodynamics Project (GGP, http://www.eas.slu.edu/GGP/ggphome.html)(Hwang et al. 2009).In addition, the Hsinchu SG station is about 8.6 km from a coastal tide gauge (TG) station.The SG and TG data are sampled at 1-sec (SG) and 6-minute intervals, respectively.The raw 1-s SG data were filtered and decimated to one-minute interval for extracting the OTL constituents below.The amplitudes and phases of 8 leading OTL constituents from the one-minute SG data were computed using the method from Boy et al. (2003).These are called "observed" amplitudes and phases and will be used to validate the "modeled" am-

3°4°5°3
plitudes and phases from SGOTL and other programs.The absolute gravity (AG) data for the case study were collected at some island stations and the Hsinchu SG station.A Micro-g FG5 absolute gravimeter (serial number: 224) was used for the data collection.The distances between the gravity stations to the nearest coasts range from a few hundred of meters (Penghu) to a few kilometers (Hsinchu).Table 2 lists the record lengths and distances to sea.The AG measurements were collected at a set 30-minute interval with each set consisting of 50 drops.All environmental corrections, except OTL, were applied.
In addition to SGOTL, several popular OTL programs were also assessed in this paper, and they are GOTIC2 (Matsumoto et al. 2001), BS (Bos and Scherneck 2009; see also http://froste.oso.chalmers.se/loading//index.html),and the g7.0 program provided by the FG5 vendor.The pro-grams g7.0 and BS use FES2004 and NAO.99b ocean tides, respectively, to compute OTL without splitting the computations into an inner zone and an outer zone.The numerical approach for SGOTL is similar to that for GOTIC2, and GOTIC2 uses tide models NAO.99Jb and NAO.99b for the inner and outer zones, respectively.Table 3 compares observed amplitudes and phases of 8 leading OTL constituents with those from SGOTL and other programs at the Hsinchu SG station.For all 8 constituents, the amplitudes from SGOTL match the observed amplitudes to 0.1 μGal.However, the deviations in phase between the SGOTL result and the SG observations can be up to 5° (O 1 ).GOTIC2 performs equally well as SGOTL in this comparison.Other two programs (g7.0 and BS) produce slightly worse results.
In Table 3, we also performed an experiment to replace the DEM-defined coastline of Taiwan by the GMT-defined coastline in SGOTL.The GMT-defined coastline of Taiwan yields amplitudes and phases at the Hsinchu SG station that are less consistent with the observed SG data.The differences in amplitude and phase caused by different coastlines can be up to 10° (N 2 ) and 0.19 μGal (K 1 ).This highlights the value of a good DEM in modeling OTL.
Because the AG gravity records were relatively short (less than 6 days), an indirect accuracy assessment of the OTL programs, based on the FG5 set scatters, was made.A set scatter is roughly equal to the standard deviation of the single measurements of a FG5 in a set and is a descriptor of AG data quality.Table 4 compares the set scatters of FG5 measurements with/without the OTL corrections applied.For all stations, applying M 2 OTL corrections has decreased the set scatters.However, GOTIC2 produces very large Newtonian gravity effects at the four island AG stations, and applying OTL corrections by GOTIC2 increases the set scatters significantly.We suspect that the coastlines used in GOTIC2 around Taiwan result in such anomalous OTL values.In fact, GOTIC2 may lack accurate shorelines for small islands such as Penghu, Kinmen and Matzu, so that this model may create spurious Newtonian gravity effects here.However, the shoreline database of GOTIC2 over large islands such as Taiwan is sufficiently good to produce useful results.On average, SGOTL produces the least set scatters among the programs in Table 4.The results in Tables 3 and  4 conclude that SGOTL outperforms GOTIC2, g7.0 and BS at the Hsinchu SG station and the island stations, and is the optimal program for OTL corrections around Taiwan and possibly the southeast China coastal areas.

COnCLuSIOnS
This paper presents the theory and the numerical method of OTL modeling in the computer package SGOTL.The height-dependent Green's functions for the Newtonian and elastic effects were introduced for numerical convolutions.For an island or coastal station, failure to consider station height may result in an error in OTL at one μGal level (section 2).The case study around Taiwan using SG and AG data shows that SGOTL performs the best when compared to some selected OTL programs.In addition to the rigor-

Fig. 2 .
Fig. 2. The amplitude of the Newtonian effect (due to M 2 ocean tide) as a function of station height at Hsinchu, YMSM and Lulin SG stations.
11) is Green's function of the elastic gravity effect.The heightdependent effect is contained in the R R H e e p + ^h term.If H p ≈ 0, then R R H 1 e e p .

Fig. 3 .
Fig. 3. Differences in M 2 OTL between sea level and non-zero station height at Hsinchu, YMSM and Lulin SG stations.

Fig. 4 .
Fig. 4. The time ratio and error ratio of a case relative to Case 2.

Fig. 5 .
Fig. 5.The amplitudes and phases of M 2 OTL at sea level.Star (SG) and circles (AG) show gravity stations.

Table 1 .
Matrix of cases for OTL numerical convolution, with row and column being respectively box size and grid interval of inner zone contribution.

Table 2 .
Starting time and record length of absolute gravity measurements at one superconducting gravimeter (SG) at Hsinchu station and four absolute gravimeters (AGs) at offshore station.

Table 3 .
Amplitudes and phases of 8 leading OTL constituents from SG observations and models at Hsinchu SG station.NOTE:1 amplitude in μGal and phase in degree. 2 SGOTL result with Taiwan coastline from DEM.3SGOTL result with Taiwan coastline from GMT. a amplitude in μGal.b local phase lag in degree and positive.theories and numerical methods used in SGOTL, a high-resolution tide model and a detailed coastline defined by a DEM around Taiwan may add to the improved result of SGOTL over other programs.According to GGP, some proposed SG stations may be installed in coastal areas.The experience with SGOTL in this paper will be useful for the OTL modeling in this case; users of SGOTL can optimize the OTL corrections by incorporating a user-defined regional tide model, coastline data set, and even a new set of Green's functions in SGOTL around the region of interest (see the sample batch job in Appendix A). ous