Factors Affecting Transpression Folding as Inferred from Numerical Analysis

Three-dimensional folding behavior of an elasto-viscous layer induced by transpression, and biaxial compression on the boundary is studied. The factors studied include the strain rate, elastic competency contrast, and viscous competency contrast. It was found that the influence of elastic competency contrast fades as the strain rate slows until finally folding behavior is completely controlled by viscous behavior. If biaxial compression is applied, induced folding along each principal stress direction is dependent. The resulting wavelength in any particular direction will not be affected by initial shortening applied in another direction. The final fold form is the supposition of fold forms in two individual directions.
For transpression, the resulting fold axis may deviate from being perpendicular to the direction of boundary displacement. As the strain rate decreases, wavelength increases. The fold axis can be curved at slow strain rates. The directions of principal stress and the fold axis are nearly perpendicular to one another in most cases. However, for medium strain rates (έ=10^(-13) ~ 10^(-14) sec^(-1)), the angle of deviation from the mutual perpendicular direction may reach 20 degrees. This deviation possibly accounts for the orientation of slaty cleavage not being perpendicular to the fold axis. Accordingly, interpretations or assertions are made based on the results of numerical analyses for the folding structure of southeastern Taiwan.


INTRODUCTION
1.1 Two-Dimensional Buckle Folding Biot (1959Biot ( , 1961) ) opened a window allowing us to look into the underlying mechanism affecting the waveform of buckle folding for elastic or viscous competent layers embedded in a matrix.For a single elastic layer embedded in spring-typed matrix, the governing equation has the following form: where Y, E, I, K, and F are the vertical displacement along the x-axis, the Young's modulus and the moment of inertia of the layer, the stiffness of the spring, and the horizontal force required to buckle the competent layer, respectively.The solution of Eq. ( 1) indicates that a dominant wavelength of the fold exists and can be expressed in terms of the elastic competency contrast as (Biot 1961;Currie et al. 1962): where ( ) l cr e is the dominant wavelength, E is the Young's modulus of the layer and E o stands for the Young's modulus of the matrix.On the other hand, if the layer is a viscous material, the dominant wavelength can be expressed in terms of the viscous competency as: where ( ) l cr v is the dominant wavelength, h is the thickness of the layer, η is the viscosity of the layer and η o the viscosity of the matrix.As revealed by Eqs. ( 2) and (3), elastic folding and viscous folding are affected by the elastic competency contrast R e and the viscous competency contrast R v , respectively.These two competency contrasts are defined as: Based on Biot's solution regarding an elastic layer in a spring-typed matrix, in which a sinusoidal waveform with a single wavelength is found (Fig. 1b), a general solution for the same governing equation can be latterly obtained (Jeng et al. 2001); in addition, two other types of sinusoidal waveforms, one with a single frequency and attenuating wave amplitude (Fig. 1c) and the other with two frequencies (Fig. 1a), are also possible.The type of waveform depends on the level of the lateral force during buckle folding.A critical force F cr is defined as: When F F cr > , the sinusoidal waveform has two wavelengths and the relationship between the applied force F and the resultant wavelength (l) can be obtained as (Karman and Biot 1940): When F F cr = , the waveform is a sinusoidal wave with a single wavelength l cr found as: Finally, when F F cr < , the waveform is a single-frequency sinusoidal wave with the attenuating amplitude, and the relation of F with the yielded wavelength l can be expressed as: The illustration of the typical F-l relation for the above-mentioned equations is shown in Fig. 2a.
If the layer has a finite length, L, the wave number, n, can be calculated based on the wavelength as n = L / l and the F-l curve can be transformed into a F-n curve, as shown in Fig. 2b.
Recent studies show that the strain rate of global lateral shortening is another dominant factor affecting waveforms (Zhang et al. 1996(Zhang et al. , 2000;;Schmalholz andPodladchikov 1999, 2000;Jeng et al. 2002).Based on the results of numerical analyses, it was found that (Jeng et al. 2002): (1) A relatively faster strain rate leads to a greater boundary force, and vice versa.
(2) For two-stage shortening at different strain rates, the final waveform depends on either the earlier strain rate inducing buckling or the later strain rate applied in the post-buckle stage of deformation.If the later strain rate is relatively fast, the final waveform will be similar to the one yielded by the fast strain rate alone; a substantial amount of elastic energy can be accumulated during the subsequent fast deformation.On the other hand, if the later strain rate is relatively slow, the earlier waveform is retained and further amplified during the slower post-buckle deformation.
(3) The dominant wavelength ratio ( λ d ) proposed by Schmalholz and Podladchikov (1999)  serves as a good measure to determine the elastic or viscous state of buckle folding and the relative speed of the strain rate.λ d is defined as: Where G is the elastic shear modulus of the competent layer and ε is the lateral strain rate.It has been found that the elasto-viscous model tends to have elastic or viscous folding behavior for λ d < 1 and λ d > 1, respectively (Schmalholz and Podladchikov 2000).
For a further detailed review on the development of research on the two-dimensional folding, please refer to Jeng et al. (2002).

Three-Dimensional Buckle Folding
Observations of oblique convergence of plates reveal that a set of en-échelon folds is a typical result of transpression (Harland 1971;Sanderson and Marchini 1984;Sylvester 1988).The three-dimensional fold produced by transpression is in fact a folding induced by two-Fig.2. The F-l relationship for a layer embedded in spring (modified after Karman and Biot 1940).(a) The variation of the lateral force while buckling and the resulting wave length; (b) The variation of the lateral force while buckling and the resulting wave number.For F F cr > , the developed wave- form is comprised of two frequencies (or wavelengths).For F F cr = , the waveform has a single frequency and constant amplitude.For F F cr < , the waveform also has a single frequency yet with attenuated amplitude.R e is the elastic competence contrast, defined as E /E o where E and E o are the Young's modulus of the competent layer and the matrix, respectively.dimensional boundary movements, as shown in Fig. 3a, in which the competent layer is subjected to both compression and transcurrent boundary displacements.In particular, the set of en-échelon folds in southeastern Taiwan is induced by the oblique convergence of the Pacific Sea Plate toward the Eurasian Continental Plate, as shown in Fig. 4 (Lu and Malavieille 1994;Lu et al. 1995Lu et al. , 2001Lu et al. , 2002)).The underlying mechanism affecting the three-dimensional folding awaits further study.

Scope of This Paper
Owing to the complexity of 3-D folding, analytical solutions for 3-D folding are difficult to obtain.With this in mind, numerical simulation for 3-D folding is adopted in this research.Before applying the 3-D numerical model, established in this research, to the 3-D folding, it was first necessary to test the adequateness of the model.Therefore, the 3-D model was initially used to simulate 2-D folding until it matched 2-D solutions.On this basis, the 3-D model was then used to simulate 3-D folding.For two-dimensional folding, the strain rates discussed in the aforementioned literature range from 10 -6 to 10 -14 sec -1 .This research further extends the slow range of the strain rate from 10 -14 to 10 -16 sec -1 to highlight the effects of viscous competency contrast R v .
For three-dimensional folding, factors including the strain rate, elastic competency contrast (R e ), and viscous competency contrast (R v ) as well as their influences on the fold form induced by the three-dimensional transpression are of interest and are to be explored in this paper.Understanding of three-dimensional folding behavior can be helpful when interpreting the evolution of the tectonic structure in southeastern Taiwan.

Model configuration
For numerical simulation of two-dimensional folding, the geometry of the analyzed model and boundary constraints are identical to those adopted by Zhang et al. (1996Zhang et al. ( , 2000)), Mancktelow (1999), andJeng et al. (2002), allowing for a direct comparison.The competent layer has a length of 99 h, where h is the thickness of the layer (Fig. 5b).The matrix has a thickness of 33 h on both sides of the layer and a "roller" constraint (free slip in horizontal direction and no vertical displacement) is imposed along the peripheral boundary of the model.The adopted material model, the Maxwell elasto-viscous model, is also identical to the material model adopted by the above-mentioned researches.The properties of the Maxwell model are listed in Table 1.
For three-dimensional folding, a three-dimensional model is constructed by more than 12000 three-dimensional elements, as depicted by Fig. 6.The competent layer is extended into a square layer with geometry similar to the two-dimensional model (Fig. 6).To enable a stable and realistic numerical outcome, proper boundary constraints have to be applied; therefore, infinite elements are applied to the surfaces B and B ' , as indicated in Fig. 3a, and form a threedimensional model, as shown in Fig. 6.Since the three-dimensional elements have a greater number of nodes and degree of freedom than two-dimensional elements do, the mesh sizes of the elements in the three-dimensional model are unavoidably coarser than those of the twodimensional model.

Boundary Conditions
The end-rotation perturbation method (Jeng et al. 2002) was adopted to induce a buckling fold following a particular distance of initial global shortening of the layer (Fig. 5a).By employing this end-rotation, a moment is induced and applied to the shortened layer, which accounts for the subsequent buckle folding.By doing this, the initial geometric perturbation is no longer needed and hence it will not affect the resulting waveform.
For the three-dimensional folding, two types of boundary conditions were studied.To simulate transpression, both normal displacement and transcurrent displacement are applied on the surfaces A and A ', forming an inclination angle α , as indicated in Fig. 3a.
Another type of three-dimensional boundary condition may exist, the so-called "biaxial compression" condition, as illustrated in Fig. 3b.Such a type of boundary condition is of interest  Table 1.Input material properties of the basic case.
especially when observation is to be made on the direction of principal directions of displacement.
Since the boundary displacements are imposed along the 1-direction and 2-direction, as shown in Fig. 3, the influence of the boundary compression in each direction on the resulting fold forms can be analyzed.
In addition to the above-mentioned two types of boundary conditions, the three-dimensional model can also be used to simulate two-dimensional folding by imposing a roller-condition on all surfaces between B and B ' so that the displacements of these two surfaces along the 2-direction can be inhibited and enable a plane-strain mode of deformation on the plane formed by the 1-axis and the 2-axis.

Verification of 2D and 3D Simulation
In the previous work done by Jeng et al. (2002), two-dimensional simulation of folding was compared to the works done by Mancktlow (1999) and Zhang et al. (1996Zhang et al. ( , 2000) ) and identical results were obtained.
To test the adequateness of the three-dimensional model in this research, this three-dimensional model is imposed on a plane-strain condition so as to simulate two-dimensional folding.Accordingly, the resulting wave numbers upon various two-dimensional initial shortenings (or boundary forces) were then compared to the two-dimensional δ -n curve, as shown in Fig. 7.In general, the resulting wave numbers are close to those predicted from two-dimensional buckle folding except for the fact that discrepancies occur in the higher frequencies of the Type A waveforms (Fig. 1a).These discrepancies possibly result from the coarser mesh size of the three-dimension model.Overall, this three-dimensional model is recognized as acceptable for further simulations of three-dimensional folding.

Influence of R e and R v -Two-Dimensional Folding
The influences of R e on the elastic F-l curve and the F-n curve are illustrated in Figs.2a and b, respectively.A greater R e leads to a smaller F cr and a greater l cr .
To identify the influence of strain rate ε and R v on the resulting waveform, simulations were made with ε varying from 10 -10 to 10 -16 sec -1 and R v varying from 100 to 10000 under a constant R e and η o of 100 and 10 13 or 10 14 MPa-sec, respectively.The resulting wavelengths l and the corresponding λ d for each case are summarized in Table 2.It was found that: (1) When the strain rate is relatively fast, e.g., ε ranging from 10 -10 to 10 -12 sec -1 with the associated λ d greater or much greater than 1, the resulting waveforms conform to the elastic F-l curve for all R v magnitudes tested (ranging from 100 to 10000), as shown in Fig. 8a.This indicates that R v controls the lateral force exerted on the boundary at the moment of folding and R e controls the waveform.This behavior was also observed by Jeng et al. (2002).
(2) When the strain rate progressively slows, with ε decreasing from 10 -13 to 10 -16 sec -1 , the yielded waveforms gradually diverge from the relationship of the elastic F-l curve, as shown in Figs.8b and c.This indicates that elastic folding fades away upon relatively slower strain rates.(3) For extremely slow strain rates, ε = 10 -15 ~ 10 -16 sec -1 with corresponding λ d being less than 0.2, the resulting waveform is no longer related to the elastic F-l curve (Fig. 8c) and shows that a greater R v results in a greater wavelength.This phenomenon agrees with the prediction of Eq. ( 3) and indicates that viscous behavior totally controls the resulting waveform for extremely slow strain rates.
The factors ε and R v also affect the strain energy stored in the competent layer and the surrounding matrix, as shown in Fig. 9: (1) As shown in Figs.9a and b, the strain energy stored in the matrix increases upon more global shortenings during the post-buckling stage.As global shortening increases from 3.5% to 30%, the amount of stored strain energy increases from about 1.5 to 2 orders of magnitude.
Meanwhile, a greater strain rate results in a greater strain energy stored in the matrix.For most cases shown in Figs.9a and b with ε = 10 -16 ~ 10 -11 sec -1 , the increase in strain energy is less than 0.5 orders of magnitude.
(2) On the other hand, when looking into the strain energy stored in the competent layer (Fig. 9c), a slower strain rate also results in smaller strain energy.Due to an increase in global shortening, the increase in strain energy is less than one order, compared to 1.5 ~ 2 orders in most cases.Due to the increase in ε from 10 -16 to 10 -11 sec -1 , the increase in strain energy in the competent layer is more than 2 orders, compared to less than 0.5 orders.

Biaxial Compression of the Elastic Layer
When the elastic layer is subjected to biaxial displacement with initial shortenings in the 1-direction (applied on surfaces A and A '; referred as δ 1 ) being a constant of 2% and in the 2-direction (applied on surfaces B and B ' ; referred as δ 2 ) being 2%, 1%, 0.5%, -0.5%, -1%, -2%, and -4%, the resulting F-n relationships are shown in Figs.10a to b.
Given a constant δ 1 , the varying δ 2 affects the boundary force exerted on surfaces A and A '; however, the F-n relationships along the 1-direction meet well with the F-n curve obtained from the two-dimensional elastic folding, as shown in Fig. 10a.Namely, the folding Table 2. Summary of the resulting wavelengths and the corresponding λ d numbers.behavior along 1-direction is not affected by δ 2 applied on surfaces B and B ' except for the fact that δ 2 the boundary force on surfaces A and A ' .On the other hand, the F-n relationships along the 2-direction shown in Fig. 10b also meet well with the F-n curve and are not affected by δ 1 .Remarkably, for the extension δ 2 , the layer will not buckle along the direction where extension force is applied (the 2-direction) since buckling folding will only occur under compressive conditions.Therefore, when the layer is subject to biaxial folding, it can be concluded that folding behavior is independent in each direction and not affected by initial boundary displacement in the other direction.
Consequently, the resulting fold-forms are a combination of the fold-forms in the 1-and 2-direction, as shown in Figs 11a to f.When there is no folding along the 2-direction, resulting from extension δ 2 (being negative), there are only foldings along 1-direction, as shown in Figs.11d, e, and f.On the other hand, when there are also foldings along the 2-direction, resulting from compressive δ 2 , the supposition of fold forms in the 1-direction and the 2-direc- tion leads to the three-dimensional fold forms, as shown in Figs.11a to c.

Transpression of the Elastic Layer
When the competent layer is subjected to transpression, for a constant shortening along the 1-direction of 1%, the variation of the fold forms with α angles is shown in Fig. 12.When α increases (a more transcurrent boundary displacement), the fold axis rotates and the wave- length changes.Although the fold axis is not perpendicular to the direction of the boundary's initial displacement, it tends to be perpendicular to the maximum principal stress (in compression), as illustrated in Fig. 12.
Since the transpression boundary condition differs from the biaxial compression boundary condition, it is of interest to know whether the folding behaviors induced by these two types of boundary conditions are similar.Therefore, the magnitudes of the maximum and the minimum principal stresses were read from the transpression simulation (Fig. 12); these magnitudes  were applied to the biaxial compression, which resulted in fold forms, as shown in Fig. 13.A comparison of the elastic fold forms induced by similar principal stresses but two different types of boundary conditions (Fig. 13) reveals similarities in the waveforms.However, discrepancies exist for the resulting fold forms.It can be asserted that the folding behavior inferred from the previous sections is still applicable to this transpression boundary condition.Nevertheless, to best simulate the fold form induced by oblique convergence, the transpressiontype boundary condition should be adopted.

Transpression of the Elasto-Viscous Layer
Subject to transpression, the variation of fold forms with varying α and ε is illustrated in Fig. 14. Figure 14 shows that significant changes of the fold form, including the direction of the fold axis and the wavelength, are induced by variation in strain rates.A smaller strain rate Comparison of folds induced by a varying α and strain rates (R v =100).The angle α is defined in Fig. 3a.tends to increase the wavelength, which is similar to the variation of two-dimension folding.Remarkably, in addition to the rotation of the fold, the fold axis is curved when viscosity is involved in the folding behavior.

Geological and Tectonic Setting of Taiwan
The Taiwan Mountain Belt is one of the most active and youngest mountain belts on Earth's surface.It has resulted from an oblique convergence between the Eurasian plate and the Philippine Sea plate (Biq 1971(Biq , 1972;;Chai 1972;Suppe 1981;Fig. 4).The convergence has propagated accretion of the Luzon arc southward and subsequently resulted in continental growth of the Asian continent.The current convergence velocity has been estimated at 8.2 cm yr -1 in the N54°W direction (Yu et al. 1997).The Luzon arc trends N10°W, whereas the Central Range trends N16°E.As it is located in the mountain belt, the southeastern Central Range of Taiwan serves as an adequate site for studying the early history of the structural evolution.
Southeastern Taiwan is underlain by a thick series of Miocene accretionary wedge deposits, including dark gray argillites and flysch deposits with occasional interbeds of gray compact sandstone and disseminated marly nodules (Hu et al. 1981).They were metamorphosed under the prehnite-pumpellyite to lower greenschist facies during the Taiwan mountain building processes (Chen et al. 1983) and exhibit intense layer parallel shearing with or without slaty cleavage (Pelletier and Hu 1984).Considering the fact that the actual Luzon arc is located to the east of the study area, these rocks could have already been exposed and exhumed before the arc-continental collision during such processes.
The topographic map of southeastern Taiwan (inset A of Fig. 4) may give one the impression that the NE en-échelon pattern of rivers and ridges is well developed (inset B of Fig. 4, Lu et al. 2002).This pattern seems to be the simplest structure in Taiwan.This particular topography is mainly controlled by regional structures of folding, thrusting left-lateral strike-slip faulting, and erosion (inset C of Fig. 4) resulting from transpressional tectonics (inset D of Fig. 4).
The accretion of landmass at the southern end of the Coastal Range in combination with oblique convergent tectonics and indentation of the forearc basement has provided a specific left-lateral right step transpressional environment in the southeastern Central Range.

Interpretation of Transpression Folding
A direct comparison of the simulated transpression folding (Fig. 14) to the folding configuration of southeastern Taiwan (Fig. 4) is not acceptable owing to at least the following reasons: (a) The actual boundary of the oblique convergence, especially at the time when folding was initiated, is not the current costal line shown in Fig. 4; (b) The actual boundary condition additionally involved subduction of the Philippine Sea plate, which is not considered at this moment; (c) The actual thickness of the folding layer is not known since the current topographic features are the combinational results of the continuous folding and erosion; and (d) The actual parameters of the layer and the matrix, especially at the time of folding initiation, are not necessarily adopted in this research.
Therefore, a direct comparison may be misleading.However, the understanding of the three-dimensional folding inferred from this study is still useful for interpreting the following issues: (a) It has been found that the orientation of slaty cleavage deviated from the direction perpendicular to the fold axis.To answer this puzzle, a hypothesis of multi-stage oblique convergence was proposed (Lu et al. 2001(Lu et al. , 2002)).However, with the results shown in Figs. 14 and 15, one stage oblique convergence and folding are possible, as long as the strain rates are medium (ranging from 10 -13 to 10 -14 sec -1 ).In addition, the mean crustal strain rate induced by the convergence of plates was found to be within the order of 10 -14 per second (Price 1975;Yu et al. 1997).Therefore, the inferred strain rate seems to be realistic.
Remarkably, the folding induced by this strain rate is within the transition of pure elastic folding and pure vicious folding, within which the fold axis can be curved as mentioned when discussing Figs. 14 and 15.(b) If the orientation of the oblique convergence boundary is parallel or sub-parallel to the current coastline, in conjunction with the fact that the fold axis is almost west to east in its orientation, the convergence direction is most likely pointed northerly since in all cases the fold axis is parallel or sub-parallel to the direction of convergence.The direction of convergence seems to be nearly northern-pointed at the time of the triggering folding and is in a N54°W direction currently (Lu et al. 2002;Chang et al. 2003;Chen 2006).If such an inference is correct, it is reasonable to assert that a counter-clockwise rotation of the convergence direction has occurred.This assertion is partly supported by the fact that a series of left strike-slips and en-échelon faults occurred after the folding (cutting the folds) with the orientation shown in Fig. 4 (inset D).

CONCLUSION
Three-dimensional folding behavior was studied.For fast strain rates, folding is dominated by elastic buckling.For slower strain rates, the influence of elastic competency contrast fades and the importance of viscous competency increases.Since the mean tectonic strain has medium strain rates, it can be expected that most of the tectonic, transpression foldings are a mixed result of elastic and viscous foldings.The fold axis can deviate from the perpendicular direction to the shortening direction, and most importantly to the principal stress direction.The latter phenomenon possibly accounts for the deviation of the cleavage orientation from being perpendicular to the fold axis.Overall, the strain rate tends to be the key factor affecting the resulting fold forms since it determines the rotation of the fold axis, i.e., whether the fold axis will be curved as well as the contribution of the elastic and the viscous competency contrast to the fold form.
Yet, if simulation of the actual three-dimensional folding induced by oblique convergence is desired, there are still more issues to be resolved, e.g., the subduction-type boundary condition, the effect of basement friction, and the effect of gravity, etc. Remarkably, the oblique convergence is actually applied on one side of the layer, rather than two sides.This research explored fundamental behavior, which hopefully may serve as a foundation for further study.

)Fig. 1 .
Fig. 1.Schematic illustration of the three types of waveforms obtained from numerical analyses.(a) Type A waveform comprises two frequencies, yielded when F F cr > ; (b) Type B waveform is a single-frequency-wave with a constant amplitude over the full length of the competent layer (when F F cr = ); and (c) Type C waveform is a single-frequency-wave characterized with an amplitude attenuation away from the perturbed end (when F F cr < ).

Fig. 4 .
Fig. 4. Geological setting of Taiwan and the studied area.(a) Main stratigraphic units of Taiwan (modified after Ho 1986; Biq 1989; Lu and Malavieille 1994) CC: Chiuchih Fault; CK: Chukou Fault; CYL: Chenyulanchi Fault; FR: Frontal Range Thrust; FT: Foothill Thrust; LO: Loshan Fault; LS: Lishan Fault; MTL: Matalanchi Fault; PKH: Peikang basement high; TL: Tili Fault; YS: Yushan; and HS: Hsuehshan.(b) A Topographic map shows the studied area and the well-aligned en-échelon structure pattern in the southeastern part of the Central Range; B Regional structures; C Perspective view of B; and D simplified interpretation of transpressional setting.

Fig. 5 .
Fig. 5. Illustration of the two-dimensional model used for numerical analysis.(a)Schematic illustration of the proposed end-rotation method adopted in numerical simulation.(b) Schematic setup of analyzed model.The resultant force on the competent layer (F) and on the matrix (F 0 ) can be computed for each node.Either a pure elastic material model or an elasto-viscous material model, for both layer and matrix, can be selected for analysis.

Fig. 6 .
Fig. 6.Illustration of the three-dimensional model used for numerical analysis.The peripherals of the competent layer and the matrix are surrounded by elements with the boundary extending to infinity, the so-called "infinite element".The global setup of the 3-D model is illustrated in Part (a).Part (b) shows the close-up of the competent layer, the matrix and the infinite elements.

Fig. 7 .
Fig. 7. Comparison of the δ -n curve with the resulting δ -n relationships, ob- tained from the three-dimensional model simulating the two-dimensional plane-strain folding.δ is the amount of initial shortening caused by lateral compression at the moment of buckling.For elastic material, a greater δ corresponds to a greater later compression force.The line indicates the developed waveform being comprised of two frequencies (or wave-lengths); and the line indicates the waveform also having a single frequency yet with attenuated amplitude.

Fig. 13 .
Fig. 13.Comparison of the fold form induced by transpression and the biaxial compression boundary conditions.The principal stresses of the fold are set to be similar in magnitude.

Fig
Fig. 14.Comparison of folds induced by a varying

Fig. 15 .
Fig. 15.Comparison directions of principal stresses of folds induced by a varying α and strain rates (R v =100).