Gravity Terrain Effect of the Seafloor Topography in Taiwan

1 Energy and Environment Research Laboratories, Industrial Technology Research Institute, Hsinchu, Taiwan, ROC 2 Institute of Geophysics, National Central University, Chung-Li, Taiwan, ROC * Corresponding author address: Dr. Lun-Tao Tong, Energy and Environment Research Laboratories, Industrial Technology Research Institute, Hsinchu, Taiwan, ROC; E-mail: tong@itri.org.tw doi: 10.3319/TAO.2007.18.4.699(T) Gravity terrain correction is used to compensate for the gravitational effects of the topography residual to the Bouguer plate. The seafloor topography off the eastern offshore of Taiwan is extremely rugged, and the depth of the sea bottom could be greater than 5000 m. In order to evaluate the terrain effect caused by the seafloor topography, a modern computer algorithm is used to calculate the terrain correction based on the digital elevation model (DEM). Based on the results of this study, the terrain effect caused by the seafloor topography makes a significant contribution to the total terrain effect in every region of Taiwan. In the southern coastal area of Taiwan, the terrain effect caused by the seafloor topography is approximately six times that caused by land topography. Therefore, the terrain effect caused by the seafloor topography cannot be neglected in the reduction of gravity data measured on land. The terrain correction should be calculated from the DEM of both land topography and bathymetry. The optimum correction distance for different regions in Taiwan is different. In eastern and southern Taiwan, the proper correction distance is 300 km, while in western and central Taiwan, it is 100 and 200 km, respectively. The accuracy of terrain correction over the rugged topography of Taiwan can be improved by compensating for the gravitational effect caused by the seafloor topography. The complete Bouguer gravity map of eastern Taiwan is reproduced based on the total terrain correction introduced in Terr. Atmos. Ocean. Sci., Vol. 18, No. 4, October 2007 700 this paper. The lineation of the gravity anomaly coincides with the known fault traces and provides the possible position of extension in those areas that are difficult to trace on the surface. In this paper, the southern extension of the Chihshang fault is well delineated in the vicinity of Luye from the gravity anomaly map and confirmed by a ground resistivity survey. (


INTRODUCTION
When ground gravity surveys are conducted in areas with varying topographies, a simple plate Bouguer correction must be compensated by adding two components: one that accounts for the topographic masses unaccounted for above the station and another imaginary component that accounts for any excess mass existing below the station (Nettleton 1976).In other words, for ground surveys, the terrain correction is always positive for a planar model and its value depends on the topographic relief in the vicinity of the station.Terrain correction is used to compensate for the gravitational effects of the topography residual to the Bouguer plate.Therefore, the terrain correction used must be consistent with the Bouguer model (LaFehr 1998;Vanicek et al. 2001).If the terrain correction due to the topography residual to the Bouguer plate is not applied, its gravitational signal contaminates the derived gravity anomaly, thereby leading to a possible misinterpretation of the subsurface structure.
Historically, gravity surveys have been classified into three types-ground surveys, marine surveys, and airborne surveys.In each type of survey, the terrain correction is calculated differently.Figure 1 shows the conceptual mode of terrain correction for ground and marine surveys.For ground surveys, terrain correction involves the summation of the gravitational attractions from parts A and B. The density of A is equal to that of B. For marine gravity data, the Bouguer correction consists of the replacement of the water column with rock material.Due to the variations in bathymetry, it is likely that this Bouguer correction might either neglect the actual seafloor relief (B part) or count it twice (A part).The terrain correction varies according to the adjacent seafloor relief, that is, sometimes it is added, while sometimes it is subtracted (Nettleton 1976).
A complete gravity survey of Taiwan was conducted by Yen et al. (1990).They also calculated the terrain corrections of 603 island-wide gravity stations based on a 1-km spacing land topographic grid of Taiwan (Yen et al. 1994).Tsai (2001) evaluated the use of a 40-m land digital elevation model (DEM) grid to calculate terrain corrections.Hwang et al. (2003) developed an algorithm to calculate terrain correction using Gaussian Quadrature based on the 80-m DEM grid.A gravity map of the Taiwan-Luzon region was compiled by Hsu et al. (1998).However, in the previous studies (Yen et al. 1994;Tsai 2001;Hwang et al. 2003;Lin 2004), only the land DEM was used for calculating the terrain correction.The assumption of terrain correction in ground gravity surveys as described in traditional textbooks might not be applicable to Taiwan because the sea bottom off the eastern coast is at a great depth.The gravitational effect caused by the water body should be compensated.This situation occurs when marine gravity surveys are conducted close to the coast.The traditional terrain correction for marine gravity surveys does not consider the effect caused by the mass above sea level.However, within some distance, the mass (mountains) above sea level will decrease the gravitational attraction measured in an offshore area.Therefore, land topography also needs to be incorporated when conducting terrain correction for marine gravity surveys within the correction distance.
The depth of the seafloor beyond 45 km from the eastern coastline of Taiwan is greater than 5000 m, and the density of seawater is less than that of rock material.Will the huge volume of seawater surrounding Taiwan produce a significant gravitational effect when compensating the gravitational effect of the topography residual to the Bouguer plate?The objective of this paper is to evaluate the terrain effect of the seafloor topography on the basis of land topographic and bathymetric DEMs.

COMPUTATION OF TERRAIN CORRECTION
Terrain correction is the most time-consuming calculation in the reduction of gravity data.Historically, terrain corrections were computed using Hammer (1939) charts at each station.However, terrain corrections can now be computed efficiently from the regular grid of a DEM (Kane 1962;Cogbill 1990).Nowadays, there have been considerable enhancements in the capabilities of laptop computers; with digital terrain data and computers, terrain corrections can be calculated in a matter of minutes.In this paper, terrain corrections are calculated using a combination of the methods described by Kane (1962) and Nagy (1966).The DEM data is resampled to a grid mesh centered on the station for which the correction is to be calculated.Kane (1962) suggested a calculation based on three zones, namely, the near zone, intermediate zone, and far zone.Various approaches for calculating the gravitational attraction of each zone are described below.
In the near zone, that is, 0 to 1 cell from the center, the terrain correction is calculated from the effects of four sloping triangular sections that describe the surface between the gravity station and the elevation at each diagonal corner.For each triangular section, the terrain correction is calculated by using the formula given below (Kane 1962): where g is the gravitational attraction; ρ T , the terrain density; θ , the horizontal angle of the triangular section; G, the gravitational constant; H, the difference between the station elevation and the average elevation of the diagonal corner; and R, the grid spacing.
The range of the intermediate zone is 2 to 8 cells from the station.The terrain effect is calculated for each cell by using the flat-topped square prism approach proposed by Nagy (1966).For each prism, the terrain correction is calculated using equation ( 2) as follows: where g is the vertical component of the attraction; ρ T , the terrain density; G, the gravitational constant; and r, the distance between a unit mass and the station.The region that extends beyond 8 cells is the far zone.The calculation of the terrain effect for this zone is based on the approximation of an annular ring segment to a square prism, as described by Kane (1962).The gravitational attraction is calculated from equation (3) as follows: where g is the gravitational attraction; ρ T , the terrain density; A, the length of the horizontal side of the prism; R 1 , the radius of the inner circle of the annular ring; R 2 , the radius of the outer circle of the annular ring; and H, the height of the annular ring or prism.The total terrain correction at each station is the summation of the local and regional terrain corrections.Both these corrections can be calculated from the DEM.A precise DEM surrounding the station is used to calculate the local terrain correction from zero to a certain distance, this distance is called the inner distance.A coarse DEM is then applied to calculate the terrain correction for the region that extends significantly beyond the inner distance.The distance to which the regional correction should be calculated is called the outer distance.In practical computations, the calculation of the regional correction is the most computationally expensive component of the calculation.In this paper, the DEM simulated using 500-and 1000-m grids (Fig. 2) provided by the National Center for Ocean Research, Ocean Data Bank, Taiwan is used as the local and regional DEMs, respectively to evaluate the terrain effect caused by the land and seafloor topographies.The 40-m DEM grid is used as the local DEM when the real gravity data is processed.
The computation of terrain correction for a station located on land can be divided into two parts.The terrain correction derived from the topography above sea level can be calculated by using traditional methods.The terrain correction derived from the seafloor topography can be calculated by assuming that the water body is replaced by rock material and then subtracting the gravitational attraction caused by the water body.In this paper, the density of seawater used for calculating the terrain correction from the seafloor topography is 1.03 g cm -3 .Finally, the total terrain correction is the summation of the terrain corrections derived from the land topography and the seafloor topography.

TERRAIN EFFECT BASED ON FULL TOPOGRAPHY
Using equation ( 1), (2), and (3), the calculation of terrain correction can be simplified as shown below.
Where G is the gravitational constant, ρ T is the terrain density, and T(s) is a function of the space coordinates that represent topographic relief.The term G T s ⋅ ( ) is called the terrain correction factor; its unit is mgal (g cm -3 ) -1 .The correction factor can be calculated based on the DEM only at an arbitrary point considered to be a virtual gravity station, following which, the correction grid can be created.The correction grid can be reused to calculate the terrain correction by extracting an interpolated value from the correction grid and then multiplying it with the specified terrain density.Thus, the correction grid only needs to be computed once.This is considerably time saving when conducting terrain correction in real gravity reduction.
Figure 2 shows the topography around Taiwan at various distances from the coast.It is evident that the water depth in the eastern offshore region of Taiwan is deeper than 5000 m, even extending to 1000 km away from the coastline.While the water depth in the western offshore region of Taiwan is shallow, the average depth of water in this area is approximately 200 m only.Four DEMs are used to evaluate the terrain effect of the seafloor topography-the Taiwan-land-only, land-only, seafloor-only, and full terrain models.For the Taiwan-land-only model, only the DEM of Taiwan land is used.For the land-only model, the DEM includes Taiwan land and the surrounding lands.For the seafloor-only model, only the bathymetric DEM is used.The DEM used for the full terrain model includes the land and seafloor topography.Since the 40-m DEM grid of Taiwan does not include bathymetric topography, the 500-m DEM grid is used as the local DEM and the 1000-m DEM grid is used as the regional DEM in order to compare the terrain effects of various models and to save computation time.The maximum correction distance is 1000 km, which is the maximum available distance based on the 1000-m DEM grid.
Figure 3a shows the correction grid of the land-only model.The correction factors in the western coastal plain are less than 1 mgal (g cm -3 ) -1 ; however, the correction factors in the Central Range are larger than 4 mgal (g cm -3 ) -1 .The highest correction factor is 24.22 mgal (g cm -3 ) -1 , which is obtained from the top of the highest mountain in Taiwan.Figure 3b shows the correction grid of the seafloor-only model.It is clearly evident that the high correction factor is mainly distributed over the eastern region of Taiwan.This phenomenon may be related to the effect of the seafloor topography.The maximum correction factor is 7.34 mgal (g cm -3 ) -1 , and it cannot be neglected in the terrain correction.Figure 3c shows the correction grid of the full terrain model.The main features of Figs.3a and c are nearly the same with the exception of the eastern and southern regions of Taiwan.The maximum correction factor is 26.94 mgal (g cm -3 ) -1 , which is larger than the correction factor of the land-only model.Fig. 2. 1000-m DEM used for calculating the regional terrain correction.

OPTIMUM CORRECTION DISTANCE
The correction distance for terrain correction should ideally extend to infinity.To reduce the computation time, it is necessary to determine a reasonable distance beyond which the terrain effect is negligible.This distance depends on the relief of the topography and the details of anomalies under investigation.For gravity data reduction in Taiwan, Yen et al. (1994) suggested that the optimum correction distance is 100 km, and Hwang et al. (2003) suggested that the best inner and outer radii of terrain correction computation were 20 and 200 km, respectively.
Figure 4 shows the correction factor versus the correction distance for each model at six virtual gravity stations around Taiwan.The positions of these stations are shown in Fig. 3. Station A is located in the eastern coastal area of Taiwan, while stations B and C are located in the southern coastal area.Station D and E are located on the northern and southwestern plain.Station F is located on top of the Central Range.It is difficult to differentiate between the terrain effects of the land-only and Taiwan-land-only models (Fig. 4).The terrain effect caused by mainland China can be neglected.Over a short distance, the terrain effect caused by the land topography dominates the total terrain correction.As shown in Fig. 4, the terrain effect caused by the seafloor topography becomes a major component of the total terrain correction beyond distances of 23, 180, and 400 km at stations A, D, and E, respectively.The terrain effect caused by the seafloor topography for stations A, D, and E is approximately 6, 4, and 1.5 times the terrain effect caused by the land topography, respectively.The terrain effects caused by the seafloor topography and land topography are at nearly the same level for station B. For stations C and F, the terrain effect caused by the seafloor topography is approximately 1/3 and 1/5 of the terrain effect caused by the land topography, respectively.It is evident that the terrain effect caused by the seafloor topography is an important part of the total terrain correction.
Figure 5 shows the ratio of the terrain correction factor from the seafloor-only model to that from the Taiwan-land-only model.It is evident that both the Huatung Longitudinal Valley and the Hengchun Peninsula yield a high ratio; therefore, the terrain effect caused by the seafloor topography is dominant in those areas and cannot be neglected.Figure 6 shows the increasing percentage of the correction factors versus the correction distance for each model at six virtual stations.For all the stations other than station F, the percentage of the correction factor reaches 90% within 100 km for the land-only and Taiwanland-only models.These findings are coincident with the result obtained by Yen et al. (1994).However, the optimum correction distance should be defined as the distance at which the percentage of the correction factor reaches 90% of the total terrain correction.Thus, as shown in Fig. 6, the optimum correction distance for stations A, B, C, and F is 325, 125, 225, and 200 km, respectively.
As shown in Fig. 6, for stations D and E, the percentage of the correction factor reaches 90% at 625 and 100 km for the full terrain model and land-only model, respectively.It is believed that the land topography mainly contributes to the terrain effect within 100 km, while beyond a distance of 100 km, the terrain effect caused by the seafloor topography becomes dominant.However, the total terrain correction for stations D and E is only 1.41 and 1.64 mgal, respectively; these values are calculated for a rock density of 2.57 g cm -3 and a water density of 1.03 g cm -3 .Thus, it is suggested that if one mgal of the total terrain correction is obtained, the practical correction distance is 150 and 95 km for stations D and E, respectively.

APPLICATION
Hu and Chen (1986) had published the gravity map of eastern Taiwan without applying terrain correction.At that time, only local terrain correction could be applied in some areas of Taiwan due to the lack of DEM and appropriate computer algorithms.In the eastern area of Taiwan especially, no gravity map with full terrain correction has been published.It is well known that the topography of the land and seafloor near this region is very rugged.To demon- strate the use of the full digital terrain model in terrain correction, the gravity data acquired in eastern Taiwan (Hu and Chen 1986) is used.Figure 7a shows the Bouguer gravity map without any terrain correction.Only a few major characteristics can be recognized in Fig. 7a.The 40-m DEM grid is used for calculating the local terrain effect within a distance of 6.5 km.The topographic effect beyond this distance is calculated by using the 1000-m DEM grid, and the outer distance is set at 1000 km.
Figure 7b shows the terrain corrections of the survey area calculated for a rock density of 2.57 g cm -3 as suggested by Yen et al. (1990) and a seawater density of 1.03 g cm -3 .The average correction is approximately 17.2 mgal.Figure 7c shows the complete Bouguer gravity map of eastern Taiwan.
Figures 8a and b show the topography and residual gravity anomaly of the southern part of the Huatung Longitudinal Valley.An interpretation of the gravity data is not the objective of this paper; however, some obvious characteristics that appear in the residual anomaly can be observed.The fault traces are superposed on the map for comparison.The gravity pattern of the Central Range is quite different from that of the Coastal Range; this indicates that the geological structures are different.The lineation of the gravity anomaly coincides with the known fault traces and provides the possible position of extensions in those areas where outcrop information is lacking.
As shown in Fig. 8b, the lineation of gravity anomaly and the Chihshang fault line published from the surface geological survey (Yu et al. 1994) are coincident.Due to the lack of clear evidence on the surface, it is difficult to determine the extension of this fault to the south of Chihshang.The lineation of gravity anomaly might provide an understanding of the extension of this fault.Figure 9 shows the tilt derivative gravity anomaly map in the vicinity of Luye.The southern extension of the Chihshang fault can be estimated by analyzing the distribution of the lineation, as shown in Fig. 9.Both seismic and resistivity profiles were obtained to confirm the position of the fault estimated from the gravity map. Figure 10 shows the resistivity image of profile CS-JY01E located in the northeastern direction of Luye.The low-resistivity zone shown in Fig. 10 may be related to the Lichi Formation and the high-resistivity zone may be related to the presence of alluvium.The position of the Chinshang fault is clearly delineated from the lateral discontinuity of the resistivity, as shown in Fig. 10.

CONCLUSION
In this paper, the terrain correction factor of Taiwan is calculated based on a full DEM.A 500-m regional correction grid of Taiwan has been compiled in this paper.This grid can be reused to calculate the regional terrain correction by extracting an interpolated value from the correction grid and then multiplying it with the specified terrain density.The depth of the seafloor in the eastern offshore region of Taiwan is greater than 5000 m.The density of seawater is less than that of rock material.The gravitational effect caused by the seafloor topography should be compensated for in the terrain correction.From the results of this paper, the terrain effect caused by the seafloor topography makes a significant contribution to the total terrain effect in every region of Taiwan.In the southern coastal area of Taiwan, the terrain effect caused by the seafloor topography is approximately six times that caused by land topography.Therefore, the terrain effect caused by the seafloor topography should not be neglected in the reduction of gravity data measured on land.The optimum correction distance for terrain correction for different regions in Taiwan is different.In eastern and southern Taiwan, the proper correction distance is suggested to be 300 km, while in western and central Taiwan, it is 100 and 200 km, respectively.To improve the precision and reduce computation time, the calculation of terrain correction can be divided into several parts; the total terrain correction is the summation of the terrain correction of each part.The slope of the terrain surrounding the gravity station within 170 m can be measured in situ and used to calculate the terrain correction near the station.A precise local DEM such as the 40-m DEM grid can be used to calculate the local correction between 170 m and 6.5 km; a coarse 500-or 1000-m DEM grid can be used to calculate the regional correction beyond 6.5 km.
The complete Bouguer gravity map of eastern Taiwan is reproduced based on the total terrain correction introduced in this paper.The lineation of the gravity anomaly coincides with the known fault traces and provides the possible position of extensions in those areas that are difficult to trace on the surface.In this paper, the southern extension of the Chihshang fault is well delineated in the vicinity of Luye from the gravity anomaly map and confirmed by a ground resistivity survey.
Acknowledgements I thank Drs. Christopher Bishop and Tai-Rong Guo for their helpful discussions.I also thank Mrs. Huey-Jen Wang for her assistance in the editing.

Fig. 1 .
Fig. 1.Illustration of Bouguer model for ground and marine survey.

Fig. 3 .
Fig. 3. Terrain correction grids of land-only model (a), seafloor-only model (b), and full terrain model (c).Six virtual gravity stations around Taiwan are labeled as A, B, C, D, E, and F.

Fig. 4 .
Fig. 4. Terrain correction factor versus the correction distance for the representative stations.The labels shown above represent the correction distance of the same correction factor from the land topography and seafloor.

Fig. 5 .
Fig. 5. Ratio of the terrain correction factor of the seafloor-only model to that of the Taiwan-land-only model.

Fig. 6 .
Fig. 6.Percentage of terrain correction factor versus the correction distance for the representative stations.The labels shown above represent the optimum correction distance for the land-only model and the full terrain model, respectively.

Fig. 7 .
Fig. 7. Shaded-relief gravity maps of eastern Taiwan.(a) Bouguer gravity without terrain correction; (b) terrain correction calculated for a rock density of 2.57 g cm -3 ; (c) complete Bouguer gravity with terrain correction.

Fig. 8 .
Fig. 8. Topography and residual gravity map of the southern part of the Huatung Longitudinal Valley.(a) Shaded-relief topography map; and (b) shadedrelief residual anomaly map.

Fig. 9 .
Fig. 9. Zoomed-in area of the black dashed box in Fig. 8 showing the shadedrelief map of tilt derivative gravity in the vicinity of Luye.

Fig. 10 .
Fig. 10.Resistivity image of profile CS-JY01E.This low-resistivity zone to the right may be related to the Lichi Formation and the high-resistivity zone to the left may be related to alluvium.The position of the fault could be clearly delineated from the lateral discontinuity of resistivity.