Correlation between the Bouguer gravity anomaly and the TAIGER tomography of the Taiwan region

The Bouguer gravity anomaly derived from observed gravity data and calculated from a 3-D P-wave velocity model were used to investigate the compatibility between the two to understand the crust structure of the Taiwan region. The seismic velocity model determined by Kuo-Chen et al. (2012) was used for our study. We converted the velocity model to a density model using the relationship between Pwave velocity and rock density proposed by Brocher (2005), and then calculated the corresponding gravity anomaly. The differences between observed gravity anomaly and calculated gravity anomaly in vicinity of the dense TAIGER seismic stations are in general small. To discuss the anomaly discrepancy between shallow and deep structure, we used the upward and downward continuation method to separate the gravitational signal into shallow and deep effects for the comparison of gravitational effect of the 3-D velocity model. A conspicuous gravity low which is lower than the observed Bouguer gravity anomaly occurred beneath the main edge of the Central Range, as shown in the calculated deep-structure gravitational map. This indicates that the Moho depth estimated by the seismic tomography is deeper than that estimated by the gravity data. The negative anomaly location differences resulted from deep-structure effects suggest that the locations of crustal thickening estimated by gravity and seismic tomography are different. Article history: Received 4 August 2017 Revised 28 February 2018 Accepted 1 March 2018


INTRODUCTION
As shown on the geological map presented in Fig. 1, Taiwan is located at a complex juncture between the Eurasian and Philippine Sea plates.The Ryukyu trench in the east is where the Philippine Sea Plate (PSP) subducts northward beneath the Eurasian Plate (EP); to the south the South China Sea plate (SCSP), a marginal oceanic subplate of EP, subducts eastward beneath the PSP at the Manila trench.The geological structures of Taiwan trend mainly in a NNE-SSW direction.The surface geology of Taiwan is dominated by Tertiary rocks, except on the east side of the Central Range, where pre-Tertiary metamorphic complexes are exposed.The Longitudinal Valley (LV) is considered the suture zone of the Eurasian and Philip-pine Sea plates, and it separates the Coastal Range from the Central Range to the west.West of the Central Range is the Foothills and then the Coastal Plain.The orogeny is considered started since ca. 4 -6 Ma (Suppe 1984;Teng 1990;Liu et al. 2001) and is currently very active with a convergence rate about 80 mm yr -1 in a direction of N54°W (Yu et al. 1997).The tectonic activities in Taiwan region is vigorous, resulting in numerous active faults, rapid crustal deformation (Yu et al. 1997), and abundant seismicity (Wu et al. 2008) due to this high rate and the junction location between two subductions.
To better understand the complicated structure beneath Taiwan, gravity studies, in addition to seismic and geological studies, provide useful geophysical information for studying shallow depths.Gravity observation is one priority method that is used to explore the subsurface structure.As this is an easy, fast, and low-cost method, dense gravity observation stations are distributed across the whole of Taiwan, including mountainous areas instead.In order to improve the understanding of the tectonics of Taiwan, an island-wide gravity survey was carried out from 1980 to 1987 (Yen and Hsieh 2010).A newly revised Bouguer anomaly map was published by Yen and Hsieh (2010) to examine the compatibility of 3-D seismic velocity structures with the gravity Terr.Atmos.Ocean. Sci., Vol. 29, No. 5, 473-483, October 2018 data of Taiwan.The differences between the revised gravity anomaly and the calculated gravity anomalies tend 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 overestimated.
Due to the non-unique solution of gravity inversion, we need the appropriate seismic velocity model as the initial model.The initial density model is reconstructed from the three dimensional tomographic velocity model using an empirical relation between velocity and density (Cheng 2004;Zhang et al. 2004;Hsieh and Yen 2016).Furthermore, the joint inversion of different geophysical datasets is an effective means to reduce the non-uniqueness and improve the reliability of geophysical inversion (Masson et al. 2012;Li et al. 2014).If a substructure model can be consistent with the two physical observations, it increases the credibility of the model.Before the joint inversion, we should examine the correlation between the original velocity model and the observed gravity data.Because the intensive seismic stations in the TAIGER project were sets in the mountains, increasing the resolution of the shallow layer velocity model.Therefore, we mainly study on the compatibility of TAIGER seismic velocity structures, Taiwan.TTT.KWR.2012(https:// ds.iris.edu/ds/products/emc-taiwantttkwr2012/),with gravity data.Follow the previous work of Yen and Hsieh (2010), we calculated the gravitational effect based on the Taiwan.TTT.KWR.2012velocity model (Kuo-Chen et al. 2012;Wu et al. 2014) derived from the TAIGER project and compared it with the observed gravity value in this study.Then, we divided the gravitational effect of velocity model into shallow and deep parts, the results of which are discussed.By comparing shallow and deep observations and calculating gravity, we can analyse the differences between different geophysical observations.According to the different distribution of gravitational isogal map, we can discuss different geophysical observations to analyse the differences in crustal thickening.

The Observed Bouguer Anomaly Analysis
Gravity measurements and analysis are particularly useful as a reconnaissance tool for a large tectonic region.A high-resolution regional gravity map over the Taiwan area (21.5 -25.5°N and 119.8 -122.2°E)can help us understand its tectonic signature.In this study, the gravity data were compiled from two sources: (1) 603 island-wide gravity stations that had been established by the Institute of Earth Sciences, Academia Sinica (Yen et al. 1995) since 1980 and (2) 3195 gravity stations that were chosen from the Ministry of the Interior (MOI).The processing flow is based on the data processing of Yen and Hsieh (2010).This gravity survey was well covered and collected gravity data to construct a Bouguer anomaly map after applying latitude, Free-air, Bouguer, and terrain corrections.The data were evenly resampled to 1076 points and integrated to obtain a Bouguer anomaly map (Fig. 2a).
The zero value of contour is on the Central Range main edge in the Bouguer anomaly map.On the western of Taiwan, the gravity anomalies show negative values.The lowest value is located in the midwestern Taiwan, reaching -60 mgal.Other negative gravity anomalies in the southwest Taiwan, such as Kaohsiung-Pingtung gravity low area, is reaching to -40 mgal which is caused by covering with thick sediment.The positive anomalies of the higher gradient zone are in the eastern part of the Central Range, Longitudinal Valley and the Coastal Range.The NW-SE trend of the isogal contour feature in the northern extremity of Taiwan, nearly perpendicular to the main struc tural trend, is similar to the structural trend.
Here we apply a mathematical algorithm for isolating gravity anomaly sources at depth.The purpose of the method is to identify the effect of the observed gravity that is harmonic above a given depth, h.We can treat this function as a gravitational effect below depth h.The main idea of the algorithm is to eliminate sources from the Earth's surface to a prescribed depth by means of upward and downward continuation.The algorithm allows separation of the effects of shallow and deep objects and extracts the gravitational effect of sources located in horizontal layers between given depths h1 and h2 (Vasin et al. 1996;Martyshko and Prutkin 2003).
The algorithm is based on upward and downward continuation (Vasin et al. 1996;Martyshko and Prutkin 2003).Firstly, we continued the data upward to height h to eliminate the effect of the sources in the near-surface layer.This operation caused main errors in the vicinity of the boundary of the area.The gravity data coverage should be extended further than the range of the target area.We performed upward continuation of the residual field by means of the following formula: , , , , Secondly, we continued the obtained function downward to depth h, that is, to the distance 2h in a downward direction.The problem of downward continuation is an ill-posed linear inverse problem; therefore, regularization must be performed.To this end, we applied a formula similar to Eq. (1): , , , , where h1 = 2h.We treated the formula as an integral equation: the right-hand side was given, and we had to determine the unknown function U(x, y, -h).
Although the influence of shallow sources is diminished after upward continuation, the remaining sources are still present in the field; therefore, we performed downward continuation through the sources.Then, the function was continued upward back to the surface (h = 0), and we separated the effects of deep (below depth h) and shallow sources.Using different h values of upward and downward continuation, we can extract the effect between h1 and h2.After subtracting this field from the total effect of deep sources, we obtained the contribution of the shallow structure.

The Computed Bouguer Anomaly Analysis
Seismic tomographic analysis is the popular method used to study subsurface structures in terms of seismic velocities.Since 1911, the modernization programme of the seismic network instrumentation is beginning.With growing stations and records of Taiwan region, it has been a number of tomographic studies conducted in past decades (Roecker et al. 1987;Rau and Wu 1995;Ma et al. 1996;Kim et al. 2005;Wang et al. 2006Wang et al. , 2009;;Wu et al. 2007Wu et al. , 2009a)).The geological processes underlying the collision and its relation to subduction are the primary targets of the Taiwan Integrated Geodynamic Research (TAIGER) project.Newly acquired passive and active source data on land, supplemented by ocean bottom as well as permanent seismic network data, were used to derive a new 3-D tomographic velocity model (Kuo-Chen et al. 2012).In particular, in the TAIGER project, the intensive seismic array was set up along four transects in the mountain area.
Since density anomaly is one dominant factor responsible for the corresponding seismic velocity anomaly and the fact that the gravity inversion is highly non-unique, seismic models are commonly used as an important calibration reference for the gravity interpretation.In this study, we used a simple strategy that retained the well-constrained geometry of 3-D velocity model and only converted the velocity to the density.In other words, we converted the velocity model to the corresponding density model, and linear interpolation and extrapolation were used to convert the different size density blocks to 3-km × 3-km × 3-km density blocks.The converted density model was specified by a box extending 252 km in the E-W direction, 402 km in the N-S direction, and 102 km in depth for the forward gravity anomaly calculation.In this study, an empirically determined densityvelocity relationship (Brocher 2005) was applied to the conversion.The following polynomial regression, fit to the selected t and V p values presented by Ludwig et al. (1970), is valid for V p values between 1.5 and 8.5 km sec -1 : . . .
After the transfer process, we obtained the density models with respect to their corresponding velocity models.Then, we calculated the gravity anomalies for the velocity models.The gravitational effect (g ij ) at the ith gravity station from a certain density block, j, can be represented as where j t is the density contrast of a certain block, j, G is the gravitational constant, Z j is the depth of certain block, j, and R ij is the distance between the gravity station, i, and a certain block, j.The forward Bouguer anomaly (g i ) created by the initial density model at the ith gravity station is where n is the total number of model blocks, G is the gravitational constant, and t D is the density contrast of each block.We summed the gravitational effect from each block, Δg to obtain the Bouguer gravity anomaly for each gravity station.Thus, the Bouguer anomaly distributions were obtained for the velocity models.The residual of each station, based on the difference between the observed Bouguer anomaly and the calculated Bouguer anomaly, was obtained.The calculated Bouguer anomaly is shown in Fig. 2b.

RESULTS
The negative anomaly distribution in Fig. 2b calculated from 3-D velocity model derived from the TAIGER project covers the central parts of the island with a minimum value of -60 mgal.The pattern of the computed Bouguer anomaly distribution is generally consistent with the structural trend, except that there are some differences in the gravity low in the south-western Taiwan and the gravity high in the eastern Taiwan.
The residual distribution map (Fig. 3a) indicates that the residual anomaly is relatively small near seismic stations at the northern, central, southern, and N-S direction lines of TAIGER project.The histogram of residual distribution (Fig. 3b) ranges from -70 to 80 mgal but about 70% of the difference value are between -20 and 20 mgal at the 414 seismic stations of the TAIGER project.The computed Bouguer anomaly values are closer to the observed gravity values in the dense seismic station area.This result indicates that imaging of the crustal structure can be better resolved with dense seismic stations with active and passive sources.
The negative value of residual distribution map in western Taiwan is caused by the higher calculated gravity value compared with the lower observed gravity value, it indicated that the velocity model is faster than the real structure in western Taiwan.Conversely, the positive value of residual distribution in eastern Taiwan is caused by the lower calculated gravity value compared with the higher observed gravity value.It shows the velocity model is slower than the real structure in eastern Taiwan.
Forward computing the converted density grids between 0 -15 km depth and 15 -102 km depth of 3-D velocity model, the computed Bouguer gravity values can be decomposed based on the results obtained from the shallow (above 15-km depth) and deep (below 15-km depth) converted density models.A flow chart of the computed processing is presented in Fig. 4. First, we calculated the theoretical gravitation of the velocity model (step 1) and decomposed the shallow (step 3) and deep (step 4) calculated gravitation.Then, we decomposed the observed gravity (step 2) into shallow (step 5) and deep gravitational effects (step 6) after upward and downward continuation processing.

Shallow and Deep Calculated Gravitational Effects of 3-D Velocity Model
To discuss the observed and calculated gravitational effects of the shallow and deep velocity models with a 15-km depth boundary, the similarities and differences be-tween seismic tomography and observed gravity need to be understood.
The calculated gravitational effect of shallow velocity model (Fig. 5a) exhibited a low value (-40 -0 mgal) in the Western Foothills and Western Plain and a high value (0 -40 mgal) in the Central Range and Coastal Range.This result is reasonable because the density of the shallow structure in the mountains is higher than that in the Western Plain.The variation of the calculated shallow gravitational effect is distributed near the sparse seismic stations.This result shows that the resolution area is limited in the shallow structure.
According to the calculated gravitational effect of the deep velocity model (Fig. 5b), the negative gravity anomalies is indicated deep crustal thickening.Conspicuous gravity lows occur beneath the main edge of the Central Range (with a minimum value of -50 mgal), that is, 10 mgal lower than the observed Bouguer gravity anomaly.

Shallow and Deep Observed Gravitational Effect
We use the upward and downward continuation method to separate the gravitational signal into shallow and deep effects to compare with the gravitational of velocity structures.In order to decrease the error at the boundary of the study area when using upward and downward continuation processing, we added an offshore gravity dataset to avoid this effect affecting our target research area and replaced the sea-water layer with continental rock (Hsieh et al. 2010).The revised Bouguer anomaly is shown in Fig. 6.Through the upward and downward continuation method, the revised Bouguer anomaly was separated into shallow and deep gravitational effects with a 15-km depth boundary.In the shallow observed gravitational effect map (Fig. 7a), low gravity areas, such as the Ilan, Hsinchu, Taichung-Puli, and Kaohsiung-Pingtung areas, are caused by sedimentary rock and coincide with the geological structure.
Figure 7b shows the deep observed gravitational effect.There are two low gravity areas in central and southwestern Taiwan.In addition to the presence of low density material, the long-wavelength negative gravity anomaly may be caused by a thicker crust.The location of the longwavelength negative gravity anomaly in observed shallow anomaly map is existed a distance of approximately 50 km west of the main ridge of the Central Range.Assuming a relatively constant crustal density, large positive Bouguer gravity anomalies correspond to thin crust and negative Bouguer gravity anomalies correspond to thick crust.From the viewpoint of gravity observations, the thickest crust occurs beneath the midwestern Taiwan.

DISCUSSION
A conspicuous gravity low occur beneath the main edge of the Central Range in the result of deep calculated gravitational effects (Fig. 5b), a long-wavelength of negative anomalies may be caused by a thickening of the crust.Furthermore, the calculated and observed negative anomalies location differences can be summarized by the differences in the location of the crustal thickening.Whether it is the result of observation or calculation of gravity, we can find that the distribution of shallow and deep anomaly has a different trend.This shows that the shallow and deep structure is quite different.
Following the maps of group velocity of fundamental mode Rayleigh waves, lower velocity regions are distributed in lower altitudes area with recent sediments, such as Western Coastal Plain, Pingtung area, Ilan plain, and Longitudinal Valley (Huang et al. 2012).The lower velocity regions are consistent with the shallow observed gravitational effect area.From the three-dimensional velocity and density structure models, low-density and low-velocity material are distributed in the shallow layer.This can be consistent with the extremely gravitational lows areas of shallow observed gravitational effect in the Pingtung area (Kim et al. 2005;Wu et al. 2007;Kuo-Chen et al. 2012;Huang et al. 2014;Hsieh and Yen 2016).According to the geomagnetic investigations in Pingtung plain, thick sediments or deeper basement are existed in Pingtung Valley (Yu and Tsai 1981).It is the reason for the extremely gravitational lows of Pingtung area.
Based on the geological interpretation, the high gravity area on the shallow observed gravitational effect map, such as the Peikang and Tai-Po areas, is caused by hard basement rock.The gravity high along the northern Eastern Central Range (ECR) results from the exposure of a high-density layer of the Tertiary metamorphic complex and gradually decreases from east to west.Along the southern ECR, the high velocity zone is distributed in shallow crustal rocks that cause the high gravity value on the shallow observed gravitational effect map (McIntosh et al. 2005).
Figure 7b shows the deep observed gravitational effect.In discussion the Moho depth variation in Taiwan from teleseismic receiver functions (Wang et al. 2010), it shows a comparison between Bouguer gravity anomaly map in Taiwan and Moho depth variation.The area of long wavelength negative effect of deep observed gravitational (Fig. 7b) is consistent with the position of the deeper Moho variation from teleseismic receiver functions.And we compare the Moho depth of the T5 line in ATSEE project (Kuo et al. 2016), which is more shallow than the Moho depth obtained by the seismic tomography (Fig. 8), and is more consistent with gravity data.It is plausible that the proposed crustal thickness of 55 km and more underneath the Central Range in previous studies (Rau and Wu 1995;Kim et al. 2005;Wu et al. 2007;Kuo-Chen et al. 2012;Li et al. 2014) was obtained commonly from seismic data.In this study, the results of the deep crustal structure, determined from two physical properties, are different.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), respectively.Kuo et al. (2016) use ray tracing method to model the crustal structures from the Wuyi-Yunkai orogen to the Taiwan orogeny.In the results, the obtained crustal thickness is range from 35 -40 km.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).Owing to the negative Bouguer gravity anomalies correspond to thick crust, we derived the calculated Bouguer gravity value from velocity model lower than the observed gravity value, we speculate that the crustal thickness beneath the Central Range should be less than 55 km.
We compared the minimum and maximum ranges of the calculated and observed gravitation values, the range of calculated gravitation is larger than the observed gravitation.This result indicates that the variation of the density model converted from the velocity model is relatively more intense than that of the observed gravity.This may imply that the velocity at the shallower depths beneath the mountainous area was overestimated compared to the real velocity, so a thicker crust is required to fit the travel time and resulting in a lower calculated gravitational effect.The negative anomaly location differences from deep effects suggest that the locations of crustal thickening estimated by gravity and seismic tomography are also different.

CONCLUSION
In this study, gravity observation plays the role of an examine tool by means of density-velocity relationship.In the area of dense seismic observation station, there is a small gravity value of difference.It indicates that imaging of the subsurface structure can be better resolved with dense seismic observation stations with active and passive sources.Analysis of the difference between the observed and calculated gravity values can be explained that the negative gravity value in western Taiwan is caused by the faster velocity model.Conversely, the positive gravity value in eastern Taiwan is due to the fact that the velocity model is slower than the real structure.
When discussing the observed and calculated gravitational effect of the shallow and deep model with a 15-km depth boundary, it is important to understand the similarities and differences between seismic tomography and observed gravity.The contours of shallow and deep effects of observed and calculated gravity are different, it indicates that the tectonic structure is also quite different.Shallow observed gravitational effect can be interpreted by geology and other geophysical observations.And deep observed gravitational effect shows the location of the crust thickening.In the mountainous region, the shallow calculated gravitational effects is larger than the shallow observed gravitational effects, it indicates that the shallow velocity model is probably overestimated.The study of Yen and Hsieh (2010) have proposed the same viewpoint, this study examines the 3-D velocity model that derived from dense seismic observation stations in mountain area to prove that this view is more credible.Overestimation of shallow model velocities also cause thickening of the crust that resolved by seismic tomography, so the deep calculated gravitational effects is lower than the deep observed gravitational effects in our results.From the viewpoint of gravity observation, the thickest crust of the island of Taiwan occurs beneath the midwestern Taiwan.
Although there are some differences in the effect of separating shallow and deep gravitational models, the pattern of calculated gravitational of 3-D velocity model is closer to the Bouguer gravity of Yen and Hsieh (2010).In conclusion, it is essential to combine dense seismic array data and gravity observations to uncover the complex tectonic structure in the Taiwan area.

Fig. 3 .
Fig. 3.The differences distribution of gravity anomaly.(a) The differences between the observed gravity and the Bouguer anomaly calculated using 3-D velocity model.(b) The histogram of residual between the observed gravity and the Bouguer anomaly calculated using the 3-D velocity model.The red dot indicates the location of the seismic stations using in TAIGER Project.The red bars show the differences at all gravity stations; blue bars show the differences at the seismic stations.Compared with the whole differences distribution (red bars), the differences in the locations of the seismic stations (blue bars) will be concentrated between -20 and 20 mgal.

Fig. 4 .Fig. 5 .Fig. 7 .
Fig. 4. Flow chart of data processing.With 15 km depth as the boundary, it is divided into the observed gravity part and the 3-D velocity model respectively for discussion.

Fig. 8 .
Fig. 8.Comparison of Moho relief in different velocity models with 7.8 km s -1 contour.Bouguer anomaly value (purple line) is referred to the left Y axis.Moho depth derived from ATSEE project (black line), Taiwan.TTT.KWR.2012model (red line), Kim's model (blue line), and Wu's model (green line) are refered to the right Y axis.The Mohr depth changes of T5 line obtained from the ATSEE project have similar trends to the gravity observations.