A Study on the Compatibility of 3-D Seismic Velocity Structures with Gravity Data of Taiwan

The Bouguer anomaly of Taiwan has been revised in this study based on more accurate terrain data provided by the Taiwanese Digital Terrain Model compiled by the Taiwan Forestry Bureau. Three seismic velocity models, those determined by Rau and Wu (1995), Kim et al. (2005), and Wu et al. (2007) respectively, were selected for our study. We converted their velocity models to density models using the relationship between P-wave velocity and rock density proposed by Ludwig et al. (1970) and Barton (1986), and then calculated their corresponding gravity anomalies. According to the correlation coefficient between the Bouguer anomalies calculated from the velocity models and the revised Bouguer anomalies, the Kim et al. model was more compatible with gravity data than the other two velocity models. The differences between the revised gravity anomaly and the calculated gravity anomalies trend toward positive values at elevations higher than 2000 m. This indicates that the velocities at the shallower depths beneath the mountainous area of the three models are overdetermined, i.e., higher than the real velocities. This ratiocination implies that the crustal thickness beneath the Central Range is less than 55 km which was obtained from the velocity models.


InTroDuCTIon
Taiwan is located at a complex juncture between the Eurasian and Philippine Sea plates as shown on the simplified geological map (Fig. 1).The mountains in Taiwan are very young, geologically speaking, formed as a result of the collision between an island arc system and the Asian continental margin (Wu 1978;Ho 1982;Tsai 1986).Its orogeny began about 5 mybp (Teng 1987) and is continuing vigorously.The geological structures of Taiwan trend mainly in a NNE-SSW direction.The Longitudinal Valley (LV) is considered the suture zone of the Eurasian and the Philippine Sea plates and separates the Coastal Range from the Central Range to the west.West of the Central Range is the Foothills and then the Coastal Plain.In northern Taiwan, the Central Range is composed of two ranges, the Backbone Range in the east and the Hsuehshan Range in the west.The southern Central Range is a single range.
In order to improve the understanding of the tectonics of Taiwan, an island-wide gravity survey was carried out from 1980 to 1987.The survey focused especially on the mostly-inaccessible mountain range where few gravity measurements were made before 1980 (Chang and Hu 1981).In all, 603 gravity stations were surveyed wherein 308 of these stations are located at elevations of 500 m or higher.This survey had notably improved the observed coverage and collected gravity data to construct a Bouguer anomaly map (Fig. 2a) after applying latitude, Free-air, Bouguer, and terrain corrections.It is well-known that terrain correction, especially in areas of steep, rugged and erratic slope topography, is the most important factor in calculating the Bouguer anomaly.Since a significant portion of our gravity stations are in mountainous areas, the terrain correction must be carefully estimated.The terrain corrections of all gravity stations were jointly calculated using two methods in Taiwan (Yen et al. 1995).The Hammer method (Hammer 1939) was used for distances up to 6.5 km from the station, while the line mass integral method (Nozaki 1981) was Terr. Atmos. Ocean. Sci., Vol. 21, No. 6, 897-904, December 2010 used for distances greater than 6.5 km.The average elevation within a single compartment for the Hammer method was estimated from the topographic contours within it.For the line mass integral method, the topographic database was obtained from a topographic map with a scale of 1 : 50000 by reading the elevations at the corners of the 1 km × 1 km grids.However, due to limited resolution in topographic data and an "intuitive" averaging within each Hammer's chart compartment, terrain corrections have inherent errors.The sources of error in these terrain calculations may arise from (1) discrepancies between the actual terrain and the flat surface obtained by the average elevation, especially in the innermost zone; (2) neglect of the effects from both upward attraction (hills) and lack of downward attraction (valleys) due to the average elevation within each compartment; and (3) poor estimation of the average elevation of all sectors (Ketelaar 1987;Telford et al. 1990;Herrera-Barrientos and Fernandez 1991).In principle, the above errors can be reduced to arbitrarily low values by constructing digital terrain data and decreasing the digitization interval.
In the past two decades, seismic velocity models in the Taiwan region have been constructed by many researchers.Their results show that the 3-D velocity structures of Taiwan are not quite consistent.In this study, we will evaluate which of those models is most compatible with the gravity data.

The reVISeD bouGuer AnomAly AnAlySIS
The digital terrain data of Taiwan compiled by the Taiwan Forestry Bureau became available in 1998.This data set was retrieved from topographic maps with a 1 : 5000 scale, grid spacing of 40 m and average elevation accuracy of 1 m.It is worth recalculating the terrain corrections using this data set.In this study, we did the recalculations to generate a revised Bouguer anomaly map as shown in Fig. 2b.The isogals in Fig. 2b do not reflect the topographic signatures clearly and display more smoothly than in Fig. 2a, particularly in the rugged mountain range.This means that the error in the terrain correction over the mountainous area has been minimized using the digital terrain data.
Figure 3 is a plot of the differences between the new and previous terrain corrections with respect to the station elevations.The figure shows that at stations located at elevations of 1000 m or higher, the differences trend toward the negative and larger magnitudes.A negative value means that the topographic effect was overestimated in the previous correction process.The maximum correction of 93 mgal (Yen et al. 1994) is at Yushan, the highest peak (3952 m) in Taiwan, while the correction is 80 mgal after recalculation.The previous terrain corrections obviously overestimated the topographic effects at the stations in the mountainous area.More specifically, the isogals in Fig. 2b with a generally NNE-SSW trend are similar to those in the previous map (Fig. 2a).They are also consistent with the structural trends of the island.Bouguer anomalies over mountainous areas, however, are commonly smaller than in the previous one.

bouGuer AnomAlIeS ConVerTeD from
The 3-D VeloCITy STruCTureS Seismic tomographic analysis is the popular method used to study subsurface structures in terms of seismic velocities.However, due to the limited numbers and coverage of seismic stations in the study areas, velocity structures may not be quite resolved.In this case, additional information may be used to compensate for incomplete data.Gravity data is one example of such information.
In the Taiwan region, seismic tomographic studies have been carried out by many researchers (e.g., Roecker et al. 1987;Rau and Wu 1995;Ma et al. 1996;Chen and Shin 1998;Kim et al. 2005;Wu et al. 2007).Roecker et al. (1987) conducted their study by using P wave arrival time data observed by the Taiwan Telemetered Seismographic Network (TTSN).The TTSN consisted of 25 stations equipped with vertical component, short-period seismometers (Liu and Tsai 1978).This network had less coverage of stations in the Taiwan region.In 1991, the TTSN was incorporated into the Central Weather Bureau Seismic Network (CWBSN).
The CWBSN consists of 71 telemetered stations equipped with three-component S13 seismometers to collect abundant high-quality and high-resolution seismological data.Since then, many seismic tomographic results have been obtained using that available earthquake data.Using a similar database observed from the CWBSN, however, the crustal struc-ture images obtained by Rau and Wu (1995) and Ma et al. (1996) are generally not consistent with each other.Kim et al. (2005) studied the 3-D P and S wave velocity structures by jointly using data sets from the CWBSN and two temporary seismic arrays in Hualien and Pingtung (Chen 1995(Chen , 1998)).More recently, Wu et al. ( 2007) conducted a seismic   The 3-D velocity structures of Taiwan obtained by the above-mentioned authors are not quite consistent.In order to find out which model would be more suitable to describe the surface structures in the Taiwan region, we evaluated them using gravity data.First, we chose three models for comparison, including Rau's (Rau and Wu 1995), Kim's (Kim et al. 2005) and Wu's model (Wu et al. 2007).Kim's and Wu's models were obtained using more data than the others.Rau's model and Ma's model were basically derived using similar databases.We chose Rau's model because Rau and Wu (1995) imaged the lithospheric structures beneath Taiwan and found a very significant variation of crustal thickness across the entire island.Ma et al. (1996) constrained to crustal depth and therefore did not resolve the deeper structures beneath Taiwan.
In this study, we use a simple strategy that remaining the geometry of each velocity model the same as well constrain and only converting the velocity to the density.In other words, we converted the velocity models to their corresponding density models for the three models chosen, and adopted a linear interpolation and extrapolation between the different size density blocks to 5 km × 5 km × 5 km density blocks.Each converted density model of this study is specified by a box that extends 250 km in the E-W direction, 400 km in the N-S direction and 100 km in depth for the forward the gravity anomaly calculation.In this study, an empirically determined density-velocity relationship (Ludwig et al. 1970;Barton 1986) was applied to the conversion.After the transfer process, we obtained the density models with respect to their corresponding velocity models.
Then, we calculated gravity anomalies for the three models respectively.The gravitational effect of each block was taken as where G is the gravitational constant and t D is the density contrast of each block.We summed up the gravitational effect from each block g D to obtain the Bouguer gravity anomaly of each gravity station.Thus, the Bouguer anomaly distributions were obtained for the three models.The correlation coefficients of each station between the revised Bouguer anomaly and the calculated Bouguer anomaly for the different models were obtained, respectively.

DISCuSSIon
As Fig. 2b shows, the revised Bouguer anomaly distributions trend generally NNE-SSW, in accordance with the structural trends of the island.The general trend of the Bouguer anomalies increases from west to east across the island, with gradients being much higher in the east than in the west.Negative anomalies cover a major part of the island; positive anomalies dominate in the eastern part of the Central Range, the Coastal Range and the northern extremity of the island.Conspicuous gravity lows are found in west central (with a minimum value of -63 mgal) and southwestern Taiwan (about -45 mgal), over the Tertiary and Quaternary sedimentary basins.The high mountains on the island lie in a zone of the high Bouguer gradient, more or less parallel to the structure strike.The highest anomaly value reaches 110 mgal along the east coast.The contact between the Philippine Sea plate and the Eurasian plate causes the Bouguer anomaly to rise sharply.Rau and Wu (1995) found a significant variation of crustal thickness across the entire island and suggested that the estimated crustal thickness beneath central Taiwan is about 55 -65 km, with a P-wave velocity of 7.5 km s -1 .The Bouguer anomaly distribution calculated from Rau's model (Fig. 4a) appears generally more complex and erratic than the revised Bouguer anomaly map (Fig. 2b) and is almost inconsistent with structural trends.Negative anomalies cover the western and northern parts of the island.A conspicuous low is located in southwestern Taiwan, with a steep gradient and a minimum value of -120 mgal.Although a gravity low is also found in southwestern Taiwan in Fig. 2b, its magnitude of anomaly is only -45 mgal.A striking gravity low, encircling the west central Taiwan, is found in Fig. 2b.But this gravity low is not shown in Fig. 4a.Two minor gravity lows appear in the western coast instead.A discrepancy between the revised and calculated gravity anomaly was observed in eastern Taiwan.The positive anomalies of the steep gradient zone and the highest anomaly value of 180 mgal are evident only in the central part of eastern Taiwan.From this gravity high area to the south, the Bouguer anomaly pattern without steep gradient implies that the lateral velocity change is not very obvious.This somewhat incongruent result does not reflect the Longitudinal Valley, which is believed to be a suture zone between the Eurasian and Philippine Sea plates.The histogram of the Bouguer anomaly differences between that calculated from Rau's model and the revised data tabulated by the number of gravity stations is shown in Fig. 4b.The differences range between -110 to 125 mgal.The correlation coefficient between the revised and calculated gravity anomaly station by station is 0.57.

Kim's model
Kim's model (Kim et al. 2005) suggested that in the western Coastal Plain and Western Foothills, the depth of the Moho is around 35 km.The Moho deepens gradually eastward and reaches a maximum depth of 55 km beneath the eastern Central Range.Figure 5a shows the corresponding Bouguer anomaly distribution of this model.Basically, the distribution resembles the revised Bouguer anomaly distribution (Fig. 2b).A good general agreement with the isogal distribution is reached, especially in eastern and southwestern Taiwan.More specifically, negative gravity anomalies also cover a major part of Taiwan.Two noticeable gravity lows are also in west central and southwestern Taiwan.The magnitude of the southwestern gravity low is about -100 mgal, which is -45 mgal lower than that in the revised Bouguer anomaly map.The positive anomalies of the higher gradient zone in the eastern part of the Central Range, Longitudinal Valley and the Coastal Range also follow the pattern of the revised Bouguer anomaly map.The NW-SE trend of the isogal contour feature in the northern extremity of Taiwan, nearly perpendicular to the main structural trend, is also similar.Figure 5b shows the histogram of the Bouguer anomaly differences between the calculated from Kim's model and the revised data tabulated by the number of gravity stations.Their differences are more centralized between 0 to 40 mgal.The correlation coefficient station by station is 0.68.

Wu's model
Wu's model (Wu et al. 2007) proposed that the Moho interface may reach down to 60 km beneath the Central Range.The Bouguer anomalies calculated using this model exhibit a disordered and localized anomaly distribution (Fig. 6a).Generally speaking, the Bouguer anomaly distribution does not correspond with the structural features.A clearly visible gravity low with an anomaly value smaller than -100 mgal is in southwestern Taiwan.Northward, the positive gravity anomaly covers in west central and northern Taiwan (north of 24°), with the exception of the northeastern tip of Taiwan.In contrast, the revised Bouguer anomaly map (Fig. 2b) shows that negative anomalies dominate in west central and northern Taiwan, except for at the northern extremity of the island.In the eastern part of Taiwan, the steep gradient zone, isogals trend NNE in consonance with the structural trend of the Longitudinal Valley and the highest anomaly value of 110 mgal is not found in Fig. 6a.This is a very remarkable discrepancy to note on the revised Bouguer anomaly map.The histogram of the Bouguer anomaly differences between those from Wu's model and the revised data tabulated by the number of gravity stations is shown in Fig. 6b.Although the lowest difference reaches to -130 mgal, the number of the positive variation is greater than that of the negative variation.It is astonishing that the correlation coefficient is 0.36.Basically, Wu et al. (2007) added the TSMIP data set to improve source-station path coverage tremendously and provided much better constraints and resolution in velocity structure determination.Oddly enough, the Bouguer anomaly from this model is less compatible than those from Rau's and Kim's models.

Summary
Figure 7 plots the Bouguer anomaly differences between the values derived from the above three models and their respective revised Bouguer anomalies with respect to the station elevations.The figure shows that the differences seem to be positive when the elevations are higher than 2000 m.In other words, Bouguer anomalies derived from the models are higher than those from observations at elevations higher than 2000 m.This may imply that velocities at the shallower depths beneath the mountainous area of those three models were overestimated compared to the real velocities.
It is plausible that the proposed crustal thickness of 55 km and more underneath the Central Range in previous studies was obtained commonly from seismic data.In comparison with other older orogenic zones, Zhao et al. (1994) modeled the deep structures across Japan subduction zone using local, regional, and teleseismic data.They obtained the crustal thickness of 36 km for the Japan Kurile arc (Hokkaido) and 40 km for the Japan trench (Honshu), respec-tively.Davey et al. (1998) conducted a geophysical study in the Southern Alps of New Zealand.A Moho depth reaching 45 km was determined.In contrast, the state of disequilibrium is most probably related to the fact that the mountains are so young that the crust has not yet had time to respond (Yen et al. 1998).The perplexed problem, the younger orogeny accomplished with the thicker crustal thickness (reaching to 55 km or more) from the different tomographic results, will be unraveled from above ratiocination.Therefore, the crustal thickness beneath the Central Range should be less than 55 km.

ConCluSIonS
In this study, we recalculated the terrain corrections using digital terrain data of Taiwan retrieved from topographic maps by the Taiwan Forestry Bureau with a 1 : 5000 scale, grid spacing of 40 m and average elevation accuracy of 1 m to obtain a revised Bouguer anomaly map (Fig. 2b).The isogals of the map do not reflect the topographic signatures concisely and display more smoothly than before (Fig. 2a), particularly in the rugged mountain range.This means that the error in terrain correction over the mountainous area was minimized by using the digital terrain data.
Seismic velocity models in Taiwan region have been carried out by many researchers in the past two decades.Their results have shown that the 3-D velocity structures of Taiwan are not quite consistent.In order to find out which model would be more suitable to describe the surface structures in the Taiwan region, we evaluated them using gravity data.Three of the seismic velocity models were chosen for our study.We converted the three velocity models into density models.Then, their corresponding gravity Bouguer anomalies were calculated.Comparing the Bouguer anomaly distribution and the correlation coefficient between the revised and calculated anomalies from the three velocity models, the velocity model of Kim et al. (2005) was more compatible than other models.The Bouguer anomaly differences between the calculated from the three models and revised gravity data indicated that velocities at the shallower depths beneath the mountainous area were overestimated compared to the real velocities.It may conclude that the crustal thickness beneath the Central Range should be less than 55 km which was obtained commonly from seismic studies.

Fig. 3 .
Fig. 3. Plot of the differences between the renewal and previous corrected values with respect to the station elevations.
using the P and S wave arrival times from the CWBSN and the recorded seismograms from the Taiwan Strong Motion Instrumentation Program (TSMIP).

Fig. 4 .
Fig. 4. (a) The Bouguer anomaly distribution calculated from Rau's model.(b) The histogram of the Bouguer anomaly differences between calculated Rau's model and the revised data tabulated by the number of gravity stations.

Fig. 5 .
Fig. 5. (a) The Bouguer anomaly distribution calculated from Kim's model.(b) The histogram of the Bouguer anomaly differences between calculated Kim's model and the revised data tabulated by the number of gravity stations.

Fig. 6 .
Fig. 6.(a) The Bouguer anomaly distribution calculated from Wu's model.(b) The histogram of the Bouguer anomaly differences between calculated Wu's model and the revised data tabulated by the number of gravity stations.

Fig. 7 .
Fig. 7. Plots of the Bouguer anomaly differences derived from (a) Rau's, (b) Kim's and (c) Wu's velocity models and the revised Bouguer anomaly, respectively, with regard to the station elevations.