The Sub-Crustal Stress Field in the Taiwan Region

We investigate the sub-crustal stress in the Taiwan region. A tectonic configuration in this region is dominated by a collision between the Philippine oceanic plate and the Eurasian continental margin. The horizontal components of the sub-crustal stress are computed based on the modified Runcorn’s formulae in terms of the stress function with a subsequent numerical differentiation. This modification increases the (degree-dependent) convergence domain of the asymptotically-convergent series and consequently allows evaluating the stress components to a spectral resolution, which is compatible with currently available global crustal models. Moreover, the solution to the Vening Meinesz-Moritz’s (VMM) inverse isostasy problem is explicitly incorporated in the stress function definition. The sub-crustal stress is then computed for a variable Moho geometry, instead of assuming only a constant Moho depth. The regional results reveal that the Philippine plate subduction underneath the Eurasian continental margin generates the shear sub-crustal stress along the Ryukyu Trench. Some stress anomalies associated with this subduction are also detected along both sides of the Okinawa Trough. A tensional stress along this divergent tectonic plate boundary is attributed to a back-arc rifting. The sub-crustal stress, which is generated by a (reverse) subduction of the Eurasian plate under the Philippine plate, propagates along both sides of the Luzon (volcanic) Arc. This stress field has a prevailing compressional pattern.


INTRODUCTION
According to Runcorn's (1964) definition the subcrustal stress (induced by mantle convection) can be determined from the gravity field.Runcorn (1967) simplified the Navier-Stokes' equations to derive horizontal coordinate components of the sub-crustal stress based on assuming a two-layered Earth model.Liu (1977) applied Runcorn's formulae to calculate the convection pattern and the stress system under the African tectonic plate.He later conducted a similar study for the sub-crustal stress under the Eurasian plate (Liu 1978) and presented a theory for the sub-crustal stress concentration and its relation with a seismogenic model for the Tangshan Earthquake (Liu 1979).Liu (1980) also studied the relation between the intra-plate volcanism and the stress generated by mantle convection.McNutt (1980) implemented the regional gravity field for studying the stress field within the crust and upper mantle.Ricard et al. (1984) investigated a possible link between the lithospheric stress and the geoid undulations.Pick and Charvátová-Jakubcová (1988) modified Runcorn's formulae to reduce the far-zone gravity anomaly contribution and the geoid undulations for local applications.Pick (1993) presented the closed-form formulae for the integral kernel in terms of the gravity anomaly and the geoid height to model the sub-crustal stress.Eshagh (2014) derived expressions for determining the sub-crustal stress using gravity gradiometry measurements.
A spectral resolution of the sub-crustal stress field, computed according to Runcorn's formulae, is restricted up to a spherical harmonic degree of 25 due to the divergence of an asymptotically-convergent series above this degree.Runcorn (1967) presented sub-crustal stress results using spherical harmonics only up to a degree of 8. Liu (1977) later limited his result up to a degree of 25 in order to avoid the divergence problem.In a more recent study Pick (1993) proposed a solution based on assuming that the sub-crustal stress is computed at sea level, instead of beneath the Earth's crust.The convergence is then guaranteed, but the results might be unrealistic.Eshagh and Tenzer (2015) proposed a procedure based on utilizing the stress function with a subsequent numerical differentiation in order to improve a spectral resolution of the sub-crustal stress field.They demonstrated that the spectral expression of the stress function is convergent (at least) up to a degree of 180.This resolution is thus comparable with the currently available global crustal models (e.g., CRUST1.0).
Another disadvantage of applying Runcorn's formulae is that it assumes only a constant Moho depth.This assumption might be sufficient in regional studies with relatively small variations in Moho depths.In global and regional studies with a variable Moho geometry, however, this assumption can yield some errors.To overcome this problem, Eshagh and Tenzer (2015) incorporated the Vening Meinesz-Moritz's (VMM) inverse isostasy problem (Vening Meinesz 1931;Moritz 1990;Sjöberg 2009) to define the stress function.The stress field is then computed for a variable Moho geometry.
In this study we apply the procedure proposed by Eshagh and Tenzer (2015) to determine the sub-crustal stress in the Taiwan region.Since the stress function is defined in a frequency domain, detailed regional gravity (Hwang et al. 2014) and crustal structure datasets cannot be used.Instead, we used global models to compile the regional stress field.The sub-crustal stress and crustal thickness computation methodologies are reviewed in section 2. The sub-crustal stress regional model is presented in section 3. A tectonic interpretation of the regional stress field is discussed in section 4. The summary and concluding remarks are given in section 5.

THEORETICAL MODEL
The method for a gravimetric determination of the sub-crustal stress combines Runcorn's formulae (of solving the Navier-Stokes' problem) and the solution to the VMM isostatic model.The sub-crustal stress and crustal thickness are then determined simultaneously according to the expressions recapitulated in this section.Eshagh and Tenzer (2015) presented the expressions for determining the meridional and prime-vertical components S x and S y of the sub-crustal stress in the following form

Sub-Crustal Stress
where S n is the stress function of degree n, Y n is the spherical harmonic function of degree n, and N is the maximum degree of a spherical harmonic expansion.A horizontal position is defined by spherical co-latitude i and longitude m .The stress function S n in Eq. ( 1) reads where 9.81 g , m s -2 is the Earth mean surface gravity, and  Sjöberg (2009) presented an approximate solution (up to the second-order term) to the VMM problem for finding the Moho depth D in the following form

Moho Depths
( , ) ( , ) sin where } is the spherical distance between the computation and integration points, v is the unit sphere, and sin d d d and the (nominal) Moho-depth spherical functions (D 0 ) n are defined as follows where G = 6.674 × 10 -11 m 3 kg -1 s -2 is the Newton's gravitational constant, g i 0 is the normal compensation attraction, , n 0 d is the Kronecker's delta, c t is the density distribution function, and H is the information on the solid topography (i.e., topographic heights on land and bathymetric depths offshore).The parameters : , , ..., , are given by The normal compensation attraction g i 0 in Eq. ( 5) is computed from (Sjöberg 2009) The Laplace's harmonics of the product H c t are defined using the following integral convolution The density distribution function c t in Eq. ( 8) reads where c t is the crustal density, and w t is the seawater density.
The computation of the third constituent on the righthand side of Eq. ( 3) is practically realized by limiting the surface integration area over the near zone while disregarding the distant zone contribution.The solution of the integral equation in Eq. ( 3) then gives the Moho depth without applying an iterative procedure.Moreover, the singularity for 0 " } is solved analytically (Sjöberg 2009).

RESULTS
The numerical models given in section 2 were applied to determine the crustal thickness and the respective sub-crustal stress field on a 5 × 5 arc-min spherical grid in the Taiwan study area (between the parallels of 17 and 30 arc-deg northern latitudes and the meridians of 115 and 127 arc-deg eastern longitudes).The refined Bouguer gravity anomalies at the surface points were calculated with a spectral resolution complete to a spherical harmonic degree of 180 (which corresponds to a half-wavelength of 1 arc-deg, or about 100 km on the Equator).The same spectral resolution was used to determine the Moho depths and the subcrustal stress components.The tectonic setup of this study area is reviewed and regional models of the refined Bouguer gravity anomalies, Moho depths and sub-crustal stress are presented in this section.

Regional Tectonic Setup
The most pronounced tectonic feature in the Taiwan region is the Philippine oceanic plate collision with the Eurasian continental margin (Fig. 1).The Eurasian continental margin in this region consists of the Sunda, Yangtze, and Okinawa sub-plates.Taiwan is directly located at the convergent tectonic boundary between the Yangtze sub-plate and the Philippine plate.The present convergence rate there is about 8 cm per year (Seno 1977;Seno et al. 1993;Yu et al. 1997).This tectonic plate boundary is rather complex, because it comprises an active orogenic belt formed by the collision between the Luzon (volcanic) Arc and the Eurasian continental margin (Teng 2011) and two subduction zones with reverse polarities.
To the east of Taiwan, the Philippine plate is subducted underneath the Eurasian plate along the Ryukyu Trench.The volcanism and back-arc rifting associated with this oceanic subduction occur at the Ryukyu Arc and the Okinawa Trough, respectively.The Okinawa Trough is an active, initial back-arc rifting basin (evolving from arc type) formed by the continental crustal extension.
To the south of Taiwan the Eurasian plate is subducted by the overriding Philippine plate along the Manila Trench.This oblique collision formed the Luzon Arc (e.g., Defant et al. 1989) which is a part of the Philippine Mobile Belt separated geologically and tectonically from the rest of the Philippine plate.This tectonic margin is characterized by Philippine plate subduction beneath the Philippine Mobile Belt along the Philippine Trench and the East Luzon Trench.

Regional Gravity Field
The coefficients of the global gravitational model GO-CO03S (Mayer-Gürr et al. 2012) were used to generate the free-air gravity anomalies.The normal gravity component Fig. 1.Regional tectonic setup.Tectonic boundaries (blue -subduction zone, green -transform zone, red -divergent zone, cyan -convergent zone) were modified according to Bird (2003).
was computed according to the GRS80 parameters (Moritz 2000).The refined Bouguer gravity anomalies were computed according to the following expression  8).The density distribution function in Eq. ( 9) was evaluated for the continental upper crustal density of 2670 kg m -3 and the seawater density of 1027 kg m -3 .The topographic and bathymetric gravity corrections are shown in Fig. 2 and the statistics are summarized in Table 1.The (step-wise) corrected gravity anomalies are presented in Fig. 3 with the statistical results given in Table 2.The refined Bouguer gravity anomalies are positive everywhere with maxima over marine areas and minima over land (see Fig. 3c).

Regional Moho Model
The VMM Moho depths were determined from the DTM2006.0 solid topography coefficients and the GO-CO03S gravity anomalies according to the expressions in Eqs. ( 3) -( 6).The normal compensation attraction in Eq. ( 5) was computed according to the expression in Eq. ( 7) using the CRUST1.0Moho depths (Laske et al. 2012).The Moho density contrast of 480 kg m -3 was adopted from the Preliminary Reference Earth Model (Dziewonski and Anderson 1981).The regional Moho model is shown in Fig. 4. The VMM Moho depths vary from 9.1 -27.3 km.The Moho geometry resembles a regional tectonic configuration.The most pronounced feature is the contrast between the oceanic and continental crustal structures along the continental margins.According to our results the Philippine plate crustal thickness is mostly less than 14 km.The Sunda plate crustal thickness is typically between 14 -17 km.Crustal deepening is seen along the continental margins.The Moho depth maxima in Central Taiwan and China locally exceed 25 km.
In principle the gravimetric Moho models determined based on adopted isostatic theory more or less disagree with the regional Moho models obtained from detailed seismic surveys for several reasons.One of the limiting factors is a relatively low resolution to which the isostatic models can reproduce the actual mass balance.Bagherbandi (2011) demonstrated that the isostatic equilibrium can be attained using a variable Moho geometry only at a long-wavelength spectrum approximately up to a spherical harmonic degree of 60 (which corresponds to a half-wavelength of 3 arc-deg, or about 330 km on the Equator).Similarly, Kaban et al. (1999Kaban et al. ( , 2004) ) concluded that the validity of isostatic models is limited mostly to spatial wavelengths not shorter than 500 km.Different estimates were given in earlier studies.Zhong (1997) reported on long wavelengths larger than 800 km, Sjöberg (1998)

Gravity corrections Min (mGal) Max (mGal) Mean (mGal) STD (mGal)
g  Table 2. Statistics of the (step-wise) corrected gravity anomalies: the GOCO03S gravity anomalies δg, the topography-corrected gravity anomalies δg T , and the topography-corrected and bathymetry-stripped gravity anomalies δg TB (i.e., the refined Bouguer gravity anomalies).underestimate the crustal thickness along the continent-tocontinent and ocean-to-continent collision zones.This is also evident from our result.As seen in Fig. 4 the maximum Moho depths in Central Taiwan reach about 25 km.The results of recent seismic studies, on the other hand, indicate that the Moho depths there locally exceed even 45 km (cf.e.g., Wang et al. 2010;Hsu et al. 2011).Despite these discrepancies between regional gravimetric and seismic models, our results (presented in the sub-section 3.4) show that the sub-crustal stress field can be detected realistically based on using the isostatic Moho model.This is explained by the fact that most of the stress energy is concentrated at the low-frequency part of a harmonic spectrum.

Regional Sub-Crustal Stress
The coefficients of the CRUST1.0Moho-depth and the DTM2006.0 solid topography were used to compute the stress function according to Eq. ( 2).The numerical differentiation was then applied to compute the horizontal subcrustal stress components.The regional model of the stress field is shown in Fig. 5 and statistics of its components are given in Table 3.The stress distribution closely resembles a regional tectonic configuration.The inter-plate stress field distinctively marks most of the tectonic plate boundaries with the maximum intensity detected at the oceanic subduction zone.In contrast the intra-plate stress is typically much smaller or completely absent.A more detailed interpretation of the stress field in the context of regional tectonic setup is given in the next section.

DISCUSSION
The maximum sub-crustal stress intensity (at the study area) was detected along the Ryukyu trench (see Fig. 5b), where the Philippine oceanic plate is subducted underneath the Eurasian continental margin (consisting there of the Okinawa and Yangtze sub-plates).This stress anomaly has two distinctive maxima, which indicate the subducted slab location.Along the Yangtze and Philippine sub-plate collision the subduction occurs to the east of Taiwan.Along the collision between the Okinawa and Philippine sub-plates the subduction along the Ryukyu trench extents northwards under the Okinawa sub-plate.The stress intensity significantly attenuates along a transform zone (with prevailing horizontal tectonic motions) in the proximity of a triple junction (of the Philippine plate and the Okinawa and Yangtze sub-plates) and in central Taiwan (see Fig. 5a).Moreover, the stress vector orientation indicates the presence of a lateral crustal extrusion in parts of northern and southern Taiwan, while a compressional stress prevails in the central orogenic belt of Taiwan.
The Philippine plate subduction underneath the Okinawa sub-plate generates additional stress field which is almost symmetrically distributed on both sides of the Okinawa Trough.The opposite orientation and increasing stress vector intensity can be explained by a tensional stress due to a back-arc rifting along this divergent tectonic plate boundary.In contrast, the Ryukyu Arc presents no pronounced stress anomalies.
To the South of Taiwan the stress anomalies are distributed on both sides of the Luzon Arc, with anomalies almost  completely absent in the central part.The stress vectors on the side of the Manila Trench have a prevailing eastward orientation, while the vectors on the side of the Philippine Trench are oriented mostly in the opposite direction.This suggests that the Sunda sub-plate subduction underneath the Philippine plate mainly generates a compressional (rather than shear) stress field.Some minor stress anomalies can also be recognized along the tectonic margin of the Yangtze and Sunda subplates.Interestingly, their spatial pattern closely resembles locations of transform, convergent and divergent sections (compare Figs. 1 and 5a).Convergent zones coincide with the maximum stress anomalies.The stress intensity decreases at transform zones and reaches minima at divergent zones (cf.Fig. 1).

SUMMARY AND CONCLUDING REMARKS
We applied a gravimetric method to compute the subcrustal stress field in the Taiwan region with a spectral resolution completed to a spherical harmonic degree of 180.This resolution is about seven times more detailed than that obtained by applying Runcorn's formulae.This method allows modeling the sub-crustal stress more realistically by taking into consideration a variable Moho geometry.
The results revealed some general characteristics of the stress field that are valid over the regional Taiwan study area (but might not be applicable in other regions).The subcrustal stress is generated mainly by inter-plate tectonics, while intra-plate stress intensity is typically much smaller.The maximum (in magnitude) stress anomalies occur at the oceanic plate subduction zone underneath the continental plate.In contrast, the reverse subduction (of the continental plate under the oceanic plate) does not directly generate any substantial stress field.Along tectonic plate boundaries locations of maximum stress anomalies closely agree with convergent zones, while the stress intensity attenuates along transform zones and diminishes along divergent zones.The crustal extrusion can be depicted from the stress vector orientation, but not from their intensity.In other words this process is not seen in the stress anomalies map.
We have been able to identify the locations of three distinctive types of sub-crustal stress fields in the Taiwan region: (1) the Philippine plate subduction underneath the Eurasian plate generates shear stress along the Ryukyu Trench.
(2) This subduction also generates additional tensional stress along the Okinawa Trough due to a back-arc rifting.(3) The oblique collision of the Philippine plate with the Eurasian plate generates compressional stress on both sides of the Luzon Arc.
It can also be concluded that despite a limited spectral resolution (up to degree of 180) of the computed quantities, we have been able to realistically reproduce, for instance, sequences of the transform, convergent and divergent sections along the tectonic plate boundaries.This finding supports the feasibility of the applied numerical scheme for computing the stress field in regional studies.It is worth mentioning that the stress field resolution can be further improved using more detailed seismic Moho models and information on the crustal density structure.For this purpose, however, the expressions defined in a frequency domain must be derived in the closed analytical forms.The computation can then be realized in a spatial domain, while avoiding the divergence problem of the asymptotically-convergent series at higher degrees of the spherical harmonics.
value of the Moho density contrast.The parameter s = 1 -D/R is a ratio function of the Moho depth D and the Earth mean radius R = 6371 × 10 3 m.The first constituent in the parenthesis on the right-hand side of Eq. (2) describes the gravitational contributions of the topography and ocean density contrast in terms of the solidtopography spherical functions H n c t ^h .The second constituent specifies a spatial localization of the stress field computed beneath the Earth's crust by means of the (nominal) Moho-depth spherical functions (D 0 ) n .

Fig. 2 .
Fig. 2. Regional maps of the gravity corrections computed with a spectral resolution complete to degree 180 of spherical harmonics: (a) the topographic correction, and (b) the bathymetric stripping correction.

(Fig. 3 .
Fig. 3. Regional maps of the corrected gravity anomalies computed with a spectral resolution complete to degree 180 of spherical harmonics: (a) GOCO03S gravity anomalies, (b) the topography-corrected gravity anomalies, and (c) the topography-corrected and bathymetrystripped gravity anomalies.

Fig. 4 .
Fig. 4. Regional map of the VMM Moho depths (in km) computed with a spectral resolution complete to degree 180 of spherical harmonics.Red line indicates tectonic plate boundaries.

Fig. 5 .
Fig. 5. Regional map of the sub-crustal stress computed with a spectral resolution complete to degree 180 of spherical harmonics: (a) the sub-crustal stress vectors, and (b) the spatial distribution of the sub-crustal stress intensity (in MPa).Red line indicates tectonic plate boundaries.

Table 1 .
Statistics of the topographic g T and bathymetric g B gravity corrections.

Table 3 .
Statistics of the horizontal coordinate components of the sub-crustal stress field.