TAO, Vol.4, No.l, 141-154, March 1993 Numerical Simulations of Collision Process in Eastern Taiwan

In the process of tectonic convergence between the plates, the Longitudi­ nal Valley in eastern Taiwan is located on the suture zone, under the influence of the Philippine Sea plate carrying volcanic ridges in its north-western front. It has been colliding with the Eurasian plate at 7cm/y in the direction of about 315°. T he area is thus possessed of its unique tectonic products and motions, and has been investigated in this study. Left strike-slip motions along the fault (in suture zone) have been implanted into the finite element elastic models incorporating with the split nodes. Prediction of tectonic evo· lution in eastern Taiwan therefore can be made by comparing outputs with observed paleo-stress/strain patterns shown on the island. T he results indicate that left-slip motion exerting on the fault plane favors shear-releasing in the oblique convergent process, such that the suture zone might turn eastward at the southern end to avoid stress accumulating in its vicinity. In addition, the process will encourage the Ryukyu trench to bend clockwise relative to its converging direction at a distance of 50-100 km from the collision margin.


INTRODUCT�ON
Taiwan is in the convergent zone between the Philippine Sea plate and the Eurasian plate.The tectonic features shown in that region are believed to result from these two plates having collided obliquely with each other beneath the Longitudinal Valley at a convergent rate of about ?cm/year moving from southeast to northwest (Figure 1) (Karig. 1973;Seno.1977; Tsai, 1978Tsai, , 1986;;Barrier & Angelier. 1986).
Off southern Taiwan, the South China Sea oceanic plate has been subducted wider the Philippine Sea plate along the Manila trench and is widemeath the Hengchun ridge where the plate is bent; however, the trench disappears after pass ing Lat.21°N northward.To the north an associated eastward dipping Benioff zone exists but becomes hard to be recognized off Taitwig offshore around 23°N (Tsai, 1986).This indicates that convergence between the Philippine Sea plate and the Eurasian plate in the vicinity of Taiwan is dominated by sub duction process south of 21°N, and the process gradually transits to collision toward 23°N (Big, 1972;Bowin et aL, 1978; Barr ier & Angelier, 1986; Shyu & Chen, 1991; Wu et al.,  1991; Huang et al., 1992).
In eastern Taiwan a collision front of the Philippine Sea plate called the C.oastal Range, 10 km wide and 130 km long , begins at Taitung City (around 23°N) and extends northward to Hualien City.To the west of the Coastal Range, a fold-and thrust belt in front of the Eurasian plate has formed the Central Range (Suppe, 1981; Davis et al., 1983; Dahlen et  al., 1984) with a suture zone called the Longitudinal Valley in between (Ho, 1975(Ho, , 1986)).
Along the eastern margin of the Longitudinal Valley, a buried tectonic boundary features aseisroic left-lateral strike-slip fault (Bowin et al. 1978;Barrier & Angelier, 1986;Yu et al. 1991).This boundary was originated by a collision process between the Eurasian continental plate and the Philippine Sea plate with some isolated Luzon paleo-volcanic ridges located in its front (Chai, 1972;Karig, 1973;Rowlett & Kelleher, 1976;Lin & Tsai, 1981;Tsai, 1986;Barrier & Angelier, 1986).As a result, the Coastal Range has been fonned with an average uplift rate of 5mm/year under compression (Peng et al., 1977).Barrier & Angelier (1986) figured that there is about 80 percent of horizontal tectonic forces acting as compressive agent in uplifting mountains through overthrust motion on fault planes.
The purpose of lhis study is to investigate the effect of left-lateral strike-slip motion, ei ther in its stress dissipating consequence, or in controlling strain process along plate boundary in eastern Taiwan under oblique convergence.

_
Huchon et al. (1986) and Lee (1986) calculated stress trajectories in Taiwan region under collision process with the a finite element models.In their models the slip component along the Longitudinal Valley Fault were excluded in the calculations; hence, the estimates made by their models would poorly represent the strain behavior occurri ng in the vicinity of the convergent region where left-lateral strike-slip motion has been involved.
In this two-dimensional (2-D) finite element model, different types of left-lateral strike slip motion along the convergent boundary were implanted by incorporating some split nodes with assumed left strike-slip values (the split-node method) to simulate discontinuous motions across slip fault plane.This split-node technique was applied to calculate the stress variation and shown to be valid in simulating crustal deformation in the focal region following the earthquake triggered by the slip process on fault plane (Yeh & Niu, 1991).Prediction of tectonic evolution involving deformation of the plates in the vicinity of the Longitudinal Valley hence can be made by comparing the results given by the model with the measurements of geodetic deformation on the island.

Assumptions
(1) This 2-D model is based on plane-stress calculation: In this case the vertical component (which points down into the Earth) of stress is zero or is comparatively small to the hor izontal (which lies on the Earth surface) stresses exerted by tectonic driving or resisted forces resulting from the plates' motion.Under this condition, it favors inducing thrust or normal faulting in response to strain process.Nevertheless, strike-slip component of the strain along the plate boundary will be included in the model (by split nodes method to be discussed later), then plane-strain assumption is not suitable.
(2) Force balance in finite element analysis follows the Equilibrium Equations, whose body force is neglected: The equations are shown as erij,j = 0, where i,j indicates x,ycomponent, respectively; er is the stress.
(3) The plates have been assumed to be isotropic with elasticity as the rheology: The constitutive relationship in plane-stress condition is thus given as Uxy where E is the Young's modulus given as 70GPa for the basalt or sandstone (Turcotte & Schubert, 1982); v is the poissons ratio as 0.3 in the model.f xx , er xx represents the strain and stress in +x direction.(4) Free-slip condition is applied onto the bottom and surface of the plates; in addition, it is frictionless as plate subducts into the mantle along trench type boundary: Although this assumption is not satisfied at seismic time scale of less than tens of year, it is acceptable at time scale greater than 100 years employed in this model.
(5) Geometry of the model is designed according to the tectonic model shown in figure 1: Northern boundary of the Philippine Sea plate in the model follows the swface trace of the Ryukyu trench, and it intersects with the trace of the Longitudinal Valley Fault, representing collision boundary between the plates near Hualien city at an angle of 45°.
In southern extension of the Longitudinal Valley offshore of Taitung city, it is assumed that the convergence is the act of free-slip following the trace of North Luzon Trough.Hence, deformation occurred in Taiwan region is simply controlled by compression and slip processes acting across the Longitudinal Valley Fault.With this assumption, in the region of Taiwan the act of collision can be simulated in terms of pushing a 200 km wide stick (with the Coastal Range in the front) against the Eurasian plate.

Model Fomulation
The mesh of the model is shown in figure 2. The left half of the rectangular region represents a portion of the Eurasian plate where the Taiwan fold-thrust belt is overlaid.The rest of the portion of the mesh represents the Philippine Sea plate.
There are 140 serendipity elements in the mesh of the model, with 8 nodes (4 at the corners, other 4 at midpoints of the sides) in each element (Zienkiewicz, 1977).The model applies displacement method with two unknowns u and v, representing the displacement value in x and y components respectively, at each node.Therefore, there are 16 degrees of freedom or polynomials in each element and they merge to solve respective (u,v) pair of displacements at each node.Finally stress values at each node can be obtained by the defining strain-stress relationship.In figure 2, the mesh contains 475 nodes.
In this model the "split node method" is applied (Melosh & Raefsky, 1981).This method gives extra-force terms exerted by discontinuous displacements across the build-in slip line, where split nodes are inserted.For instance, for implanting left-lateral strike-slip motion along collision boundary of the plates (as shown in figure 2 A disadvantage of using the "split node" method is that those discontinuous slip vectors can only be as input parameters in solving strain/stress pattern in response to faulting process; however, with the technique the trend and amount of slip motion occurring along the fault plane could be predicted by finding minimum strain energy exerted in deformed regions through feeding different combinations of inputs to split nodes.

Experiment I: Continuous domain without slip component
This experiment demonstrates strain/stress behavior around the Taiwan region while the Eurasian plate and the Philippine Sea plate have been interlocking together along the points of contact in eastern Taiwan during the convergence.
Figure 3a shows horizontal displacement at each node as the Philippine Sea plate moves northwest 14 cm relative to the Eurasian plate in a two-year period.In southern Ta iwan (south of the contact comer near Taitung City at 23°N) the Eurasian plate is dragged by the convergence under locked condition along the plate boundary in the Longitudinal Valley, by showing that the vectors make as large as 50° clockwise rotations pointing toward the convergent zone.This phenomenon will encourage collision process propagating southward ,  'i \ H'\'q":"t'l'"ii�\• q \ \ , as well as stress accwnulation near the southern convergent comer.Figure 4a shows the shear stress occurring in that region is about twice the magnitude of the northern portion of the collision zone.
However, in the rest of the portion of the Taiwan region (to the north of 23 °N) displace ment vectors rotate clockwise as much as 20° from the relative converging direction between the Eurasian plate and the Philippine Sea plate, which shows agreement with the paleo-motion occurring in northern Taiwan obtained from geomagnetic data (Lee et al., 1991a(Lee et al., , 1991b)).Magnitude of horizontal displacement near the convergent zone is about half the amount of the convergent motion between the Eurasian plate and the Philippine Sea plate (fi gure 3a).
In this experiment convergent motion begins to retard within the Philippine Sea plate at a distance of 100-200 km from its collision front.As a result, relative motion between the Coastal Range and the Central Range across the Longitudinal Valley becomes small relative to the plate convergent rate (or 7crn/y).This is not the situation observed in eastern Taiwan, that Yu et al. (1991) estimated .The present-day Philippine Sea plate's motion relative to the Eurasian plate across the Longitudinal Valley is about half of the convergent rate.This is because the slip process was excluded in this experiment, and it is expected that the strike-slip process exerted along the Longitudinal Valley Fault shall prevent the retardation of horizontal velocity in the Philippine Sea plate under the convergence.
The direction of compressional stress in and around Taiwan is nearly parallel to the .relative motion between the Philippine Sea plate and the Eurasian plate (figure 5a).This phenomenon was supported by earthquake focal mechanisms and paleosttess study in the Tai wan region (Angelier et al., 1986;Yeh et al., 1991 ).However, there are two characteristics at variance with other studies: (1) In southern Taiwan at the south end of the convergent zone it acts as a radiant center; the principal axis lies in an east-west direction (at azimuth of 105° /285°) and changes its bearings gradually clockwise until compressional axis is in accordance with the convergent direction.It indicates that the signifi cant locking effect exerts strongly at south convergent comer across the contact boundary.Eventually, along the Longitudinal Valley Fault going northward, the compressional stress gradually changes its direction from about 0° (or parallel to) to about 15° clockwise with the convergent direction of motion.This result shows different sttess pattern given by Angelier et al. (1986) and Huchon et al (1986).(2) Around the north end of the convergent zone offshore of the Hualien area, the results do not show the complex stress pattern determined by Yeh et al. (1991).
According to figure 5a, when compression continues such that shear stress exceeds the yield stress, fracturing develops along the zone of convergence in the Longitudinal Valley because the fracture always follow the line intersecting with the axis of maximwn com pression at an angle of 30°-45° (Turcotte & Schubert, 1982).This experiment supports the break tending to start at the south end of the convergent boundary because shear stress is at maximum there.

Experiment II: Implanting a Left-slip Convergent Boundary
Shear sttess caused by compressive motion can be dissipated efficiently through slip motion upon fault plane.The Philippine Sea plate has been obliquely pushing Taiwan island; hence, it is reasonable to suggest that there exists a slip component on the contact plane along the Longitudinal Valley fault during the collision process.As a result, the strain/stress pattern shown in Taiwan region shall not totally be the same as demonstrated in the Experiment I, which did not include stress dissipation effect contributed by the strike-slip process acting upon the collision boundary.In this experiment a left-lateral strike-slip motion is implanted across the contact bound ary, with split-node method.Relative magnitude of strike-slip motion across the Longitudinal Valley Fault is equal to the amount of the component of the convergent vector that is parallel to the strike of the Fault line.
Figure 3b shows displacement vectors under the influence of left-lateral strike-slip mo tion upon the Longitudinal Valley Fault with other conditions the same as Experiment I.At each accommodated split node, there are two vectors representing relative motion on either side of the strike-slip fault (see right box shown in fi gure 3b).Two horizontal displacements across the fault therefore represent the joining effect derived from the slip and the strain processes wider the collision.
Displacement field shows the following characteristics: (1) Relative motion between the Coastal Range and the Central Range across the narrow Longitudinal Valley becomes comparable to the half of the plate convergent rate (the rate is 7cm/y) under the effect of strike-slip process, which matches the situation estimated in eastern Taiwan (Yu et al., 1991).
(2) East of the convergent zone the Philippine Sea plate turns into the direction parallel to the strike of the contact boundary without losing its magnitude of motion compared with the plate convergent rate under the effect of strike-slip process.
(3) In the central and northern Taiwan region (to the north of 23°N) displacement vectors still pose about 20° -30° clockwise rotation relative to the convergent direction between the Eurasian plate and the Philippine Sea plate.It indicates that the strike-slip process does not affect the direction of motion in the Taiwan region obtained in Experiment I; the influence is primarily on the magnitude of the motion.' Hence, judging from the displacement field, it is not surprising that shear stress devel oping around the collision region will be dissipated by the slip process.In figure 4b it shows as much as 50% of shear that is dissipated except around the corner of the southern end of the collision zone.
The direction of the compressional principal stress is almost identical with that shown in Experiment I (figure Sb).except in the two regions: (1) In the northern portion of the Philippine Sea plate at a distance of about 50-100 km from the collision front near the Hualien area where there exists significant axes rotation; however, this phenomenon can be supported by earthquake focal mechanisms (Yeh et al., 1991); (2) fu southern Taiwan the direction of principal stress does not exactly match the observation (Angelier et al., 1986), indicating that (i) following the southern extension of the Longitudinal Valley the plate under the ocean may be in the process of collision; although this phenomena does not affect the prediction made by this model to the region around the Longitudinal Valley and northern Taiwan; as well as (ii) the strike-slip direction along the southern portion of Longitudinal Valley Fault is at variance with the north-south trend that given in the model.

Experiment III: A Curve Left-slip Convergent Boundary
Stress dissipation at the southern contact comer was not efficient under the slip process conducted in the Experiment II.This Experiment will try to curve the convergent boundary at its southern end to southeast direction, and hope that high shear stress anomaly which developed in the vicinity of southern contact corner can be diminished.
Three left-slip vectors belonging to the split nodes located on the elements 107 & 115 (figure 2) turn 5°, 10°, and 20° toward the east (without changing their slip magnitudes applied in the experiment II), respectively (figure 3c).In other words, south trace of the Longitudinal Valley Fault begins to turn southeast at a distance of about 30-40 kilometers north of its comer.Under the circumstances, shear stress in the convergent zone is efficiently dissipated (figure 4c).The stress anomaly with a smaller magnitude transits to the location where the strike-slip fault begins to curve.The northward migration of the east-west trended compres sional principal axis (figure Sc) compares with the stress pattern obtained in Experiment II (fi gure Sb), indicating resistance of slip motion originated in that location.It is reasonable to suggest that this stress anomaly can be totally removed by slightly adjusting the bend of the fault.
According to this experiment, it is predicted that the collision boundary tends to turn southeast at its southern comer, but the bowidary will go straight into the sea at its northern end.Tracking of the surface trace of the Longitudinal Valley on a map indicates this sug gestion.In addition, this experiment shows agreement with geodetic measurements in the region of the convergent zone of eastern Taiwan (Yu & Lee, 1986) (figure 6), implying the validity of this mode.The similarities between them are as follows: (1) In the central portion of the Longitudinal Valley Fault, two horizontal principal strains show they are both in compresion with E2 (the second principal strain) pointing normal to the convergent direction.This phenomenon encourages giving larger strain in the direction normal to the convergent plate surface (or in the case of giving higher uplifting in mountain-building).
(2) The magnitude of compressional strain in the convergent direction tends to retard from the central to southern portion of the fault.In the case of being without existence of a curved strike-slip phenomena, the maximum compressional strain shall occur near the southern end of the Longitudinal Valley Fault (as shown in figure 6, _Experiment I case).
ll!ill geodetic measurement experiment 1 • computed principal strains Yu& Lee,1986 Fig. 6.Computed horizontal principal strains in eastern Ta iwan.Arrow length is proportional to the strain magnitude (with the meeting arrow pair for com pressional strain, the departing arrow pair for extensional strain).They are compared with principal strains obtained by geodetic measurements (3) They all show similar compressional directions for the principal strains, except at the southern end of the fault.This variance may be caused by the influence of being a boundary node in the finite element calculations.This experiment suggests (1) collision process has been experienced in eastern Taiwan and may extend southward toward under the sea, (2) left-lateral strike-slip component is an important motion across the Longitudinal Valley fault in dissipating the strain/stress.(3) the trend of the Longitudinal Valley Fault may bend southeastward at its south end.

CONCLUSION
1.In a time scale longer than tens of years, the buried boundary in the Longitudinal valley poses a left-lateral strike-slip motion triggered by oblique convergence between the Philippine Sea plate and the Eurasian plate.This motion has a tendency to dissipate the shear efficiently caused by collision process acting in Taiwan region.2. Longitudinal Valley shall go straight northward into the sea, and it sustains less strain than at its southern corner.Nevertheless, at the southern comer of the convergent zone a strike-slip fault shall tum southeast to avoid accumulating too much of the shear stress in that area.3.There is a tendency that stress is easily accumulated by both sides of the southern portion of the collision boundary.It will sustain highest compressive force acting perpendicular to the strike of the contact fault representing resistance of the slip motion acting in that region.4.This model suggests that the north margin of the Philippine Sea plate tends to move in a direction parallel to the collision line at a distance of 50-100 km away from the collision front.This is also supported by the fact that the Ryukyu trench may tum north off the east coast of Taiwan.In addition, this model does not support that the existence of right-lateral strike-slip fault found at the location where the trench curves because of there being too little shear stress developing there as suggested by the experiments of the model.

Fig. 1 .
Fig. 1.Schematic diagram showing collision between Philippine Sea plate and Eurasian plate and plate tectonic setting of Taiwan (modified from Barrier & Angelier, 1986).
), the nodes shared by the elements 111-115 (belonging to western boundary of the Philippine Sea plate) and the elements 103-107 (belonging to eastern boundary of the Eurasian plate) are called split nodes when extra-slip vectors representing strike-slip component across the fault are inserted in the model.Thus in each split element there are 6 extra-force terms given by those extra-slip vectors on the fault line for matrix solving to get the solution (Melosh & Raefsky, 1981).
Fig. 2. Geometry and applying te ctonic framework of the finite element model.Model area has been divided into 140 elements, representing two homo geneous tectonic domains: Philippine Sea plate on the left and Eurasian Plate on the right.Two types of plate boundary are shown: collision boundary (Longitudinal Valley Fault denoted by double arrows indicating its left-slip behavior) and subduction boundary (Ryukyu trench and north Luzon trough with triangles on the overlying side).The area occupying the overlying region in the vicinity of the boundary is vacant in order to simulate free-slip condition along subduction boundary in geological time scale.

Fig. 3 .
Fig. 3. Computed relative horizontal displacements at finite element nodes of (a) experiment l, (b) experiment 2, and (c) experiment 3 as the Philippine Sea plate moves at a distance of 14.14 cm parallel to its upper or lower boundary, and western boundary of Euranian plate in the model is fixed .Arr ows point in the directions of motion under the influence of the convergence.The diagram shown in the respect ive box on the right in dicates the displacements across the Longitudinal Valley Fault in each experim�nt, where the double line outlines the assumed rupture zone of the fault and the single line indicates no slip exened along the fault.

Fig. 4 .
Fig. 4. Contours of computed deviatoric stresses of (a) experiment 1, (b) experw iment 2, and (c) experiment 3. Boundary conditions are the same with that described in figure 3. Unit of contours is in Pascal.