應力場變化對於三維裂隙系統演化的影響：以台灣東北部為例

全文

(1)國立臺灣師範大學理學院地球科學系碩士論文 Department of Earth Sciences National Taiwan Normal University Master’s Thesis. 應力場變化對於三維裂隙系統演化的影響：以台灣東北部為例 Influences of stress variation on the evolution of 3D fracture system: An example in NE Taiwan. Pham Thi Huong 指導教授：葉恩肇博士 Advisor: En-Chao Yeh, Ph.D.. 中華民國 109 年 9 月 September 2020.

(2) ACKNOWLEDGEMENTS First of all, I would like to express my thank to my patient and supportive supervisor, Prof. En-Chao Yeh. Prof. En-Chao continuously provided encouragement and was always willing and enthusiastic to assist in any way he could throughout my research project. I am extremely grateful for our friendly chats and your support in my academic and my daily life in Taiwan. I am incredibly grateful for following people at Structural Geomechanics Laboratory of Dept. Earth Sciences of National Taiwan Normal University, without whom I would not have been able to complete this research. Special thanks for the helpful contribution from my friendly lab-mate Ping-Chuan Chen, Shin-Kuan Yang, Guan-Hao Chen, Sing-Jhih Liu, Wei-Lun Tsai, and Yung-An Liao and valuable suggestion from lab assistants Yung-Tien Chiu and Xiao-Jun Peng. My sincere thanks also go to all Professors and staff of the Department of Earth Sciences of National Taiwan Normal University. My completion of this project could not have been accomplished without their support. I would also like to extend my special thanks to Dr. Chih-Hsiang Yeh and my friend Kelly from National Chiao Tung University for their useful guidance in my methodology to carry out my research. Finally, I would like to express deep and sincere gratitude to my parents for their love, caring, and sacrifices for educating and preparing me for my future. I am very thankful to all my Vietnamese and Taiwanese classmates for their love, understanding, and support.. i.

(3) ABSTRACT When nucleation of fractures and reactivation of pre-existing fractures was supposedly triggered by earthquake events, the stress state characterized by the focal mechanism of earthquake events can manipulate the 3D fracture patterns. That 3D fracture pattern might present the distribution, density variation, and connectivity of the fracture system in spatial always received special attention from geologists for the considerations regarding the mining and exploration development, disaster forecasting, and even the evolution of mountain building processes. Northeast Taiwan provided a special opportunity to understand the evolution of stress state from compressive subduction at the south to back-arc extension at north in the northeastern part of this active orogenic belt. The stress state in cross-section generated by inverted the focal mechanism data, which was collected from 1991 to 2016. The certain stress regime of Northeast Taiwan primarily changes from reverse faulting at the south to strike-slip faulting at the north, which approximately conducts the transformation of the orientation of the nucleated fractures network system. Even the presence of different in-situ stress regimes can cause reactivation of pre-existing in later. This region basically provides to the important situation to evaluate the evolution of fracture patterns with the stress field, contemporaneously. Therefore, our study aims to seeking, identifying, investigating the geometric characteristics of fracture planes of NE Taiwan by analyzing Digital Elevation Models (DEM). From DEM, the recent development of GIS and remote sensing technology can support to access and collect data of the fracture system in the inaccessible region.. ii.

(4) Our results show that in three domains, we recognized 187 surface fractures by applied the cross-cutting relationship between foliation and fracture planes and also calculated the fracture attitudes. In general, the spatial distribution of mean striking of fracture plane changed from N-S in the Domain III via NE-SW in the Domain II to ENE-WSW orientation in the Domain I. However, to get better results, only 120 fractures planes with R2 ≥ 0.85 and length ≥ 1000 meters are compatible to predicted fractures generated from focal mechanism analysis. After evaluating the fracture instabilities, we realized that nucleated surface fractures in DEM could be correlated with the applied stress field in the evolution of the tectonic setting from the south to the north along the coastline area at NE Taiwan. The number of nucleated surface fractures of Domain I is highest (64%) in comparison with other domains. Additionally, we examine the slip tendency analysis to measure the reactivation of nucleated surface fracture. The results confirmed the important roles of the extensional force of back-arc Okinawa Trough opening to the reactivated fracture in NE Taiwan. The examination result of slip tendency presented 54% nucleated fracture planes of Domain I could be reactivated under the stress state. Therefore, mapping and classifying fracture planes provided valuable information to the 3D fracture map, which should be able to further implication in terms of the landslide, groundwater, and scientific research issues. Keywords: Fractures Instability Analysis, 3D Fracture Mapping, Stress Inversion Analysis.. iii.

(5) CONTENT ACKNOWLEDGEMENTS .............................................................................................. i ABSTRACT ..................................................................................................................... ii CONTENT ...................................................................................................................... iv LIST OF TABLES ......................................................................................................... vii LIST OF ABBREVIATIONS ....................................................................................... viii LIST OF FIGURES ......................................................................................................... ix INTRODUCTION ............................................................................... 1 Motivation .......................................................................................................... 1 1.1.1 Tectonics of NE Taiwan .......................................................................... 1 1.1.2 The model of fracture development at coastline area in NE Taiwan ...... 5 1.1.3 The reliability of digital elevation data in the geological examination ... 7 Purposes ............................................................................................................. 8 GEOLOGICAL BACKGROUND ...................................................... 9 Tectonic setting .................................................................................................. 9 2.1.1 Taiwan mountain belt .............................................................................. 9 2.1.2 Kinematics in NE Taiwan ..................................................................... 12 Geological Background ................................................................................... 16 2.2.1 Pre-Tertiary metamorphic rock ............................................................. 19 2.2.2 Tertiary metamorphic rock and sediment deposit ................................. 20 METHODOLOGY ............................................................................ 21 Geological definitions. ..................................................................................... 22 3.1.1 Foliation ................................................................................................. 22 3.1.2 Fracture .................................................................................................. 22 3.1.3 Stress tensor ........................................................................................... 23 iv.

(6) Methods ........................................................................................................... 23 3.2.1 DEM Database....................................................................................... 23 3.2.2 Red Relief Image Model (RRIM) .......................................................... 24 3.2.3 Lineament Map ...................................................................................... 29 3.2.4 Extraction of Fracture Attitude .............................................................. 30 3.2.5 Stress Inversion and Instability Analysis .............................................. 34 RESULTS.......................................................................................... 39 Identification of fracture and foliation ............................................................. 39 4.1.1 Attitude of foliation planes .................................................................... 39 4.1.2 Fracture plane mapping ......................................................................... 51 Analysis result of focal mechanism data ......................................................... 61 4.2.1 Distribution of focal mechanism. .......................................................... 61 4.2.2 Stress inversion ...................................................................................... 65 Fracture Instability ........................................................................................... 70 4.3.1 Threshold value of fault instabilities ..................................................... 70 4.3.2 Evolution of fracture system in NE Taiwan .......................................... 76 DISCUSSIONS ................................................................................. 79 Mapping limitations of DEM in metamorphic area ......................................... 80 Stress inversion results and stress evolution path in NE Taiwan .................... 83 Spatial distribution of surface fractures ........................................................... 86 Evaluation of reactivated fractures in NE Taiwan ........................................... 91 5.4.1 Slip tendency analysis ........................................................................... 92 5.4.2 The reactivated fractures in the tectonic evolution scenery .................. 96 CONCLUSIONS ............................................................................. 101 REFERENCES ............................................................................................................. 103. v.

(7) Appendix A Surface fractures orientation of domain I ................................................ 112 Appendix B Surface fractures orientation of domain II ............................................... 113 Appendix C Surface fractures orientation of domain III .............................................. 115 Appendix D Related parameters of surface fractures at domain I ................................ 117 Appendix E Related parameters of surface fractures at domain II ............................... 118 Appendix F Related parameters of surface fractures at domain III .............................. 120 Summary of questions/answers and suggestions .......................................................... 123. vi.

(8) LIST OF TABLES Table 2.1: Lithological units in the study area. .............................................................. 18 Table 3.1: The summary of rock mechanics data. .......................................................... 36 Table 4.1: The foliations attitudes calculated from the DEM data at Domain I............. 44 Table 4.2: Foliation attitude from the geological map ................................................... 44 Table 4.3: The distribution of foliation attitudes calculated from DEM data of Domain II. .................................................................................................................. 46 Table 4.4: The distribution of foliation attitudes calculated from DEM data of Domain III. ................................................................................................................ 50 Table 4.5: The distribution foliation attitudes of field data analyzed from the geological map. ............................................................................................................. 50 Table 4.6: The principal stress orientation and standard deviation. ............................... 70 Table 4.7: Threshold values of surface fracture instability examination. ...................... 76 Table 5.1: Threshold values of slip tendency in this study. ........................................... 96. vii.

(9) LIST OF ABBREVIATIONS CeR. Central Range. DEM. Digital Elevation Model. EP. Eurasian Plate. GPS. Global Positioning System. LVF. Longitudinal Valley Fault. LZV. Luzon Volcanic Arc. NE. Northeast. NF. Normal Fault. PSP. Philippine Sea Plate. RF. Reverse Fault. RRIM. Red Relief Image Model. SS. Strike-slip Fault. viii.

(10) LIST OF FIGURES Figure 1.1: (A) Tectonic map of Taiwan showing tectonic units. (B) Distribution of focal mechanisms (169 EQ events, M>4). ............................................................. 4 Figure 1.2: A proposed evolution model of fracture patterns in NE Taiwan. .................... 6 Figure 2.1: 3D model of the tectonic setting of Taiwan.. ................................................. 11 Figure 2.2: The kinematics of Taiwan island from stress tensor inversion.. .................... 14 Figure 2.3: The evolution of the stress state in NE Taiwan. ............................................. 15 Figure 2.4: (A) Geological units of Taiwan, (B) The Geological map of NE Taiwan ..... 17 Figure 3.1: Flowchart of the method................................................................................. 21 Figure 3.2: Elements of topographic openness map. ........................................................ 26 Figure 3.3: Red Relief Image Model. ............................................................................... 28 Figure 3.4: Illustration of plane attitude from LiDAR. ..................................................... 32 Figure 3.5: Demonstration of fracture planes from lineaments in the 3D view ............... 32 Figure 3.6: The calculation of fracture and foliation plane attitude in lower hemisphere projection.. ................................................................................................... 33 Figure 3.7: Fault instability in Mohr ‘s diagram. .............................................................. 38 Figure 4.1: Location map of geologic features in the study area. ..................................... 40 Figure 4.2: Foliation traces and foliation planes in the 3D view. ..................................... 41 Figure 4.3: Comparison of foliation attitude between DTM-derivation and field data at the location 1................................................................................................ 43 Figure 4.4: Comparison of foliation attitude between DTM-derivation and field data at the location 2................................................................................................ 47 Figure 4.5: Comparison of foliation attitude between DTM-derivation and field data at the location 3................................................................................................ 49. ix.

(11) Figure 4.6: Characteristics of the fracture plane in the 3D view of Arc Sense................. 52 Figure 4.7: Spatial distribution of surface fracture planes in DEM data.. ........................ 54 Figure 4.8: Classification of the fracture planes based on R2 value.. ............................... 56 Figure 4.9: Relationship of subsurface rupture length with the magnitude of earthquake events.. ......................................................................................................... 58 Figure 4.10: Distribution of subsurface rupture length of fracture planes in study area. . 59 Figure 4.11: Classification of fracture planes based on the horizontal distance. .............. 60 Figure 4.12: Distribution of focal mechanism of NE Taiwan (M>2). .............................. 63 Figure 4.13: Focal mechanism data distribution along cross-section AA’. ...................... 64 Figure 4.14: Focal mechanism distribution of B-B’ cross-section. .................................. 65 Figure 4.15: Map of stress domains in the study.. ............................................................ 66 Figure 4.16: Exhumation rate of NE Taiwan area ............................................................ 66 Figure 4.17: Stress inversion results with different depths. .............................................. 67 Figure 4.18: The reasonable stress state of each domain for fracture instability.............. 69 Figure 4.19: The threshold of fault instability at Domain III. .......................................... 71 Figure 4.20: The threshold of fault instability at Domain II.. ........................................... 73 Figure 4.21: The threshold of fault instability at Domain I. ............................................. 75 Figure 4.22: The evolution of surface fracture patterns in NE Taiwan. ........................... 78 Figure 5.1: The effective surface process on fracture plane mapping. ............................. 82 Figure 5.2: The current kinematic model of NE Taiwan .................................................. 86 Figure 5.3: Distribution of mean strike of surface fractures from southern to northern part of NE Taiwan coastline. .............................................................................. 88 Figure 5.4: Nucleated fracture at NE Taiwan under the prevailing stress regime at each cell from southern to northern part of NE Taiwan. ..................................... 90. x.

(12) Figure 5.5: A proposed evolution model of reactivated fracture patterns in NE Taiwan. 91 Figure 5.6: The simulation model of slip tendency in applied stress state ....................... 92 Figure 5.7: Threshold of slip tendency at Domain I ......................................................... 93 Figure 5.8: Threshold of slip tendency at Domain II ........................................................ 94 Figure 5.9: Threshold of slip tendency at domain III ....................................................... 95 Figure 5.10: Reactivation of surface fractures in the scenario of tectonic evolution in NE Taiwan. ........................................................................................................ 98 Figure 6.1: Characteristics of foliation and fracture in the 3D environment. A. Foliation; B. Fractures. ............................................................................................... 127 Figure 6.2: One surface fracture are parallel with the lithology boundary in the study area. ................................................................................................................... 128. xi.

(13) INTRODUCTION Motivation 1.1.1 Tectonics of NE Taiwan Stress state derived from the focal mechanism of earthquake events might govern the attitude distribution of surface fractures. The focal mechanism solution was on analyzed the first motion of seismic faults in the slip plane (Morris et al., 1996). Each fault type is developed in a distinctive stress state. For example, the reverse fault is commonly correlated with the compressional stress regime, whereas the normal fault is primarily indicated the extensional stress regime. The evolution of the stress state highly manipulated in the connection and distribution of fracture system, even might be useful to identify the fracture generation. Therefore, the stress field played important roles in constructing the formation of nucleated surface fracture and sparking the reactivation process of fractures. We realized that Taiwan orogen provides us an unique opportunity to study the evolution of structural characteristics at the transition area between of Philippine Sea Plate (PSP) and Eurasian Plate (EP). NW-SE collision at east Taiwan and N-S subduction at NE Taiwan between two plates implied the change of different environments from the collision, post-collision collapse at east Taiwan, and subduction to back-arc extension at NE Taiwan (Teng, 1990; Shyu et al., 2005). In eastern Taiwan, PSP collides with EP since about 5 Ma (Ho, 1986; Teng, 1990). Conversely, the topography dramatically dropped from Hualien to Ilan Plain might present the collapse of NE Taiwan orogen due to PSP subduct beneath EP in the NE Taiwan (Teng, 1996) and the back-arc opening of southern Okinawa Trough examined by Sibuet et al. (1987). 1.

(14) and Nishimura et al. (2004). Thus, northeastern Taiwan has been highly involved with different stages of the tectonic evolution of Taiwan orogen. Numerous previous studies (e.g., Angelier et al., 1990; Clift et al., 2008; Wu et al., 2009; Kang et al., 2015) already discussed the tectonic evolution from the northern segment of Longitudinal Valley Fault (LVF) to the northernmost of Taiwan island. Based on their analyses, the local crustal deformation of Taiwan might be experienced through: (1) The ongoing oblique collision of NE-SW trending Chinese continental margin of EP with the mostly N-S trending Luzon Volcanic Arc (LVA) of PSP since 5 Ma (Ho, 1986; Teng, 1990; Shyu et al., 2005) formed the NE-trending Taiwan orogen. (2) Offshore of NE Taiwan has been influenced by the Ryukyu subduction system where the EP overrides the PSP along the Ryukyu Trench with the rate approximately 40 mm/yr in NE orientation (Lallemand and Liu, 1998; Nakamura, 2004) and (3) The back-arc opening in N170○-180○E direction with the rates of 50 mm/yr developed by the active rifting of the Okinawa Trough behind the Ryukyu arc (Nishimura et al., 2004; Sibuet et al., 1987) (Fig 1.1A). The consequences of the arc-continental collision, subduction, and rifting have been suggested to lead the first-order crustal deformation pattern in NE Taiwan. Besides that, Lai et al. (2009) documented some earthquake events with the moderate-sizes (e.g., 2005 Mw 5.4, 2015 Mw 5.7, and 2016 Mw 5.6 events) at the Ilan plain. They are basically consistent with the strike-slip faulting motion (Lai et al., 2009). And the M = 7.3 earthquake event happened on the Milun Fault at the part of northern LVF segments, with two following M 7.0 earthquake sequences occurred in other segments of LVF (Chen et al., 2008). They have been questioned whether the manipulation of the focal mechanism from earthquake events motivated and modified the deformation patterns at. 2.

(15) NE Taiwan, and, whether the combination of tectonic setting and earthquake events can generate the existence of a fracture pattern at NE Taiwan. Previous seismic investigations practically indicated the transitional pattern of a stress state and seismicity in NE Taiwan (Lin et al., 2004; Wu et al., 2007; Wu et al., 2008; Chou &Chen, 2012 and Huang et al., 2012). In particular, Wu et al. (2008) illustrated that in the northeastern seismic zone, the intermediate-depth earthquakes are associated with the PSP-EP subduction system, whereas shallow events follow the opening of the Okinawa Trough and collision between EP and LVA. In Figure 1.1B, based on the results of relocated hypocenters and focal mechanism first motion, Huang et al. (2012) presented the evolution of stress state in NE Taiwan. In the northern part, normal faulting focal mechanisms (blue balls) are dominant, while 2 clusters of mostly left-lateral strike-slip focal mechanisms (light green balls) are discovered in the transition area (southern Ilan Plain and Nanao area). In the southern part, reverse faulting focal mechanisms are predominantly presented at the Hualien domain with 2 NNE-SSW trending clusters (red balls). Therefore, we supposed that the evolution of the in-situ stress pattern could govern and construct the surface fractures pattern in NE Taiwan by different stages.. 3.

(16) Figure 1.1: (A) Tectonic map of Taiwan showing tectonic units. (B) Distribution of focal mechanisms (169 EQ events, M>4). (A) Hengchung ridge in blue, Luzon volcanic arc in light orange, Ryukyu trench is in the light purple area. The red arrow indicated the predicted plate motion vector from NUVEL-1. (B) The kinematic model of the tectonic evolution in NE Taiwan was proposed by Huang et al., (2012) based on the distribution of focal mechanism. The study is located inside the black rectangle.. 4.

(17) 1.1.2 The model of fracture development at coastline area in NE Taiwan The stress state is an important controlling factor for the instability behavior of fractures (Sauber & Molnia, 2004; Sutinen et al., 2009 and Vavrycuk et al., 2015). The slip vector of fractures might correspond with the distribution of the stress state, such as paleo-stress analysis comes from fault-slip data (Bott, 1959; Delaney et al., 1986). The structural characterization, network system, and assessment of the faults and fractures are the keys to analyzes and demonstrate in-situ stress. According to the evolution of the tectonic setting from south to north along the coastline area at NE Taiwan described in Section 1.1.1, we divided in coastline area into three different Domains I, II & III, respectively (Fig 1.2). In this figure, we proposed the simulated fractures evolution model of the coastline area. In the early stage, the whole coastline area supposedly experienced the compressional environment; hence it produced the fault and fracture slipped in the applied reverse faulting regime (red area in Fig 1.2A). The result of the flipping of subduction polarity between PSP and EP at the secondary stage, which is led the post-orogenic collapse from Hualien to Ilan Plain, has been conducted the presence of strike-slip faulting in the Domain I and II (Fig 1.2B). In the late stage, the opening of southwestern Okinawa Trough since 2 Ma years ago become an important driving force in NE Taiwan. At that time, it dedicated the normal faulting regime in northern Taiwan and promoted the construction of fracture patterns in Domain I (Fig 1.2C). However, pre-existing fractures generated by the paleo-stress state still under considering the reactivation by the current stress state. Therefore, the estimation stress state based on the fracture investigation is crucial. Thus, this study aims to mapping the fracture pattern in different domains and examining it with a regional stress state.. 5.

(18) Figure 1.2: A proposed evolution model of fracture patterns in NE Taiwan. (A) First stage: The reverse faulting regime (red area) created a fracture pattern, which is corresponded to the collision time. (B) 2nd stage: Before the opening of Okinawa Trough, the flip of the subduction system PSP-EP in the northeastern most of Taiwan triggered the strike-slip faulting regime in Domain II and I (light green area). It developed fractures slipped and dilated in this stress state. (C) Final stage: The opening of the Southern Okinawa Trough was influenced the construction of fracture in normal faulting regime and punctured to NE Taiwan to form the Ilan Plain.. 6.

(19) 1.1.3 The reliability of digital elevation data in the geological examination Fracture and fault can be indicated by lineaments, which is expressed in the topography. Additionally, lineaments are usually associated with the fault, the linear zone of fracturing, folding, as well as linear chains of some geological features. Quantitative approaches to reveal and classify fractures and faults by examining lineaments from remotely sensed images have been repeatedly consideration (Wilson, 1941; Lattman, 1958). Vinogradova and Eremin (1971) are the first researchers who operated photogrammetric procedures to reveal and classify fractures, as well as to measure their strike and dip angle. In current time, accompanied with the progressive development of remote sensing research, geologists started to employ some standard techniques of digital elevation model to identify the fractures, such as a 3D topographic model constructed from aerial or satellite images (Campagna and Levandowski, 1991; Morris, 1991; Kukowski et al., 2001, 2008). Moreover, in recent years, the Digital Elevation Model (DEM) process can enhance “the contrast” of topography, which allows us to identify lineaments and fractures exposed on the surface. It has been providing the visualization of geological lineaments and their relationship to nearby geometric and geologic features. Especially, in the 3D model investigation, those geologic features exposed on the surface can be identified and classified by high-resolution DEM data. Further important applications of DEM data are the characterization of the fault and fracture system (Masoud and Koike, 2011a&b) and the construction of fracture plane from interpreted lineaments (Koike et al., 1998). The previous studies used DEM data to identify the dip direction of the bedding plane in sedimentary rock and examine equivalent structure geology (Cracknell et al., 2014).. 7.

(20) In this study, we examined and recognized surface fractures at NE Taiwan from 5meter DEM data. The lithology of NE Taiwan is dominated by metamorphic rock. In the counting, 2/3 rocks are schist and related schist. Fracture planes might be able to distinguish from foliation planes by cross-cutting relationships, which are present in the 3D DEM environment.. Purposes The main purpose of this study is to correlate the existing or pre-existing fracture pattern on the surface with applied stress state inferred from the underground focal mechanism. The contribution of the 3D nucleated fracture map construction in this study aim to achieve steps (1) to investigate the significant characteristics to distinguish the fracture and foliation in the 3D DEM environment by cross-cutting relationship because we examined fracture planes in metamorphic rock; (2) produces the 3D nucleated fracture map with reliable, determined strike and dip of fracture planes from DEM data and constructs the cross-section of stress state by inverted focal mechanism data; (3) evaluate fault instability of reasonable surface fractures in an applied stress state to construct the 3D nucleated fracture map at NE Taiwan.. 8.

(21) GEOLOGICAL BACKGROUND Tectonic setting 2.1.1 Taiwan mountain belt The Taiwan orogenic belt was developed by the ongoing oblique collision between EP and PSP since the late Cenozoic (Ho, 1986; Teng, 1990). It is the consequence of the collision between mostly N-S trending of the Luzon volcanic arc and NE-SW trending of EP continental margin, which mainly constructed the current high topography of this island (Chi, 1981, Suppe 1981). According to Yu et al. (1999), PSP overrode northwestward in EP along the Manila trench with the rate 80 mm/ year. And this collision had begun in northern Taiwan a few million years before just now happened in the southern Taiwan island (Huang et al., 2012). And the suture zone of this collision is proposed as the Longitudinal Valley in southeastern Taiwan (Shyu et al., 2005). Besides that, the flipping of subduction polarity between PSP and EP caused by slab break-off (Teng et al., 1996; Teng et al., 2000) promoted PSP subducted under the EP resulting in the Ryuku Trench and Okinawa Trough (Yu et al., 1997; Wu et al., 2010) in the late Cenozoic (Fig 2.1). This subduction system, accompanied by the volcanic arc-continent collision, has been proposed to cause the post-orogenic collapse in NE Taiwan (Teng et al., 2000) and furthermore to form the Ryukyu trench and develop the NW-SE opening of the back-arc basin Okinawa Trough in northeastern Taiwan (Fig 2.1).. 9.

(22) According to previous works, the crustal thickness decreases from > 33km to < 18km affiliated with the topographic drop from northern Taiwan to southern Ryukyu (Wu, 1978; Lee et al.,1980; Yen et al., 1995). Moreover, the formation of Okinawa Trough was associated with the Pliocene-Quaternary andesites of orogen strata. Extensional structures resulted from actively numerous submarine volcanoes located in the southwest Okinawa Trough (Lee and Wang, 1988) cut across the compressional structure by the andesites. Therefore, the Taiwan mountain belt involves the kinematic evolution of the different orogenic processes that evidential for the transition of the EP continental to PSP oceanic subduction from south to north of NE Taiwan (Malavieille et al., 2019).. 10.

(23) Figure 2.1: 3D model of the tectonic setting of Taiwan. In southwestern Taiwan, the oceanic crust of the South China Sea (2) subducted eastward under the PSP and created LVZ (4) along the transition area China continental margin. In eastern Taiwan, the northwestward movement of the western edge of PSP (3) collided with LVZ and overrode the EP (5) in the south, whereas it subducted beneath EP in the north along the Ryukyu arc (6) of NE Taiwan. (Malaivielle et al., 2019).. 11.

(24) 2.1.2 Kinematics in NE Taiwan NE Taiwan, where is experienced in different tectonic environments from the collision in the south to extension in the north (Fig 2.2), indicated that the stress regime has evolved from the southern to the northern part (Wu et al., 2008). The collision primarily constructed orogen and explained the existence of the thrust sheet, fold, and other compressional structures. Teng et al. (2000) supposed that the Luzon arc in the early Pliocene overrode the continental margin rapidly during the movement of PSP-EP collision system. As the eastward subduction of the continent margin, the collision orogen grows to the maximum in central Taiwan. In the north, the tectonic setting getting more complicated. Sibuet and Hsu (1997) supposed that the post-orogenic collapse of NE Taiwan related to the westward propagation of the PSP-EP collision system via the PSP-EP subduction system and back-arc spreading of the southern Okinawa Trough. And, the flipping of the Ryukyu subduction polarity system caused the extensional process, which is highly connected with the post-orogenic collapse. The break-down of NE Taiwan orogen is referenced by the reduction of topography in Ilan Plain and the eastern of oceanic crust of southern Okinawa Trough. The fast rifting caused the southward bending of the mountain belt in nearly E-W striking (Hsu et al., 2009). Therefore, numerous extensional and compressional focal mechanisms were recorded and addressed in northeastern Taiwan orogen (Wu et al., 2008; Huang et al., 2012) (Fig 2.2). The compressional stress regime (thrust and strike-slip faulting) is dominant in most of Taiwan (Yeh et al., 2014) (Fig 2.2). However, local tectonic variation in NE Taiwan induced the different distribution of the stress regime. In particular, Wu et al. (2008) and Wu et al. (2010) suggested the focal mechanism range from thrust and. 12.

(25) strike-slip faulting in the Hualien to normal faulting in the transition area and Ilan Plain. Especially, Wu et al. (2008) indicated that the trend of the tensile axis near Ilan Plain is close to NW-SE trending, whereas the NW-SE orientation of the compression axis is dominant in the southern part of the Heping area (Fig 2.2). Therefore, it is significant to evaluate the evolution of the stress pattern in NE Taiwan. Huang et al. (2012) determined regional stress patterns at NE Taiwan from the focal mechanism in NE Taiwan by relocated EQ events with magnitude >4.0 (depth <30km) (Fig 2.3). Particularly, six different stress groups were identified from 169 focal mechanisms. Within them, three groups are associated with our studies. Firstly, the thrust faulting regime with NW-SE principal stress (. 1). is dominant in the Hualien area. (Domain III), which is displayed by the red circle. In the Heping area (Domain II), Huang et al. (2012) found the two different clusters of strike-slip faulting regimes, the first strike-slip faulting in the southern Ilan Plain (Domain I) has consistent with the opening direction of Okinawa Trough, whereas. 1. 3. orientation. orientation of the. second strike-slip faulting regime in the Heping area might have corresponded with arccontinent collision. The light green square stress regime (Fig 2.3) of the Heping area represented the earthquake events with an NNW-SSE orientation of WSW orientation of. 1.. 3. and an ENE-. Meanwhile, in southern Ilan Plain, the light green triangular. stress regime (Fig 2.3) has a trans-tension stress pattern.. 13.

(26) Figure 2.2: The kinematics of Taiwan island from stress tensor inversion. Surface projection of the 1 and 3 axes (left) and the different principal axes in the equal-area projection of the lower hemisphere (Wu et al., 2008).. 14.

(27) Figure 2.3: The evolution of the stress state in NE Taiwan. The result of the focal mechanism analysis from Huang et al. (2012) emphasized the stress state evolution in Northeastern Taiwan. In this figure, the normal-faulting focal mechanism cluster illustrated by the blue dots is dominant in Ilan Plain. High density of reverse-faulting cluster was discovered in Hualien county with the red dots and the strike-slip faulting cluster highlighted by green dots located in the central part of study area.. 15.

(28) Geological Background Taiwan orogen can be classified by six narrow geologic units parallel to the long axis of the island (Fig 2.4A). From west to the east, they comprised of the Coastal Plain (CP), the Western Foothill (WF), the Hsueshan Range (HR), the Central Range (CeR); the Longitudinal Valley (LV), and the Coastal Range (CoR). Inside, CeR contains the Tananao Schist (TC) (Ho, 1988). The Tananao Schist is the oldest metamorphic basement of Taiwan, experienced several phases of orogenic deformation, magmatism, and metamorphism (Wang, 1966; Ernst et al., 1983; Stanley et al., 1981). Especially, the higher metamorphic grade that occurred in the northeastern part of Tananao schist might be caused by the strong influences of the tectonic evolution of Taiwan. Furthermore, Tananao Schist can be divided into two metamorphic belts: the western Tailuko belt and the eastern Yuli belt (Ho, 1986). This research addressed the northeastern Tailuko belt (red box, Fig 2.4B). NE Tailuko belt is mainly composed of Pre-Tertiary metamorphic rock such as pelitic schist, gneiss, migmatite, metamorphosed limestone, greenschist, and amphibolite. The detail will be described below.. 16.

(29) Figure 2.4: (A) Geological units of Taiwan (Ho, 1988), (B) The Geological map of NE Taiwan (Geological Central Survey, 1995).. 17.

(30) LATE PALEOZOIC TO MESOZOIC. LATE MESOZOIC. CENOZOIC. Table 2.1: Lithological units in the study area (Central Geological Survey, 1995). Stratigraphical division. Lithology. Nansuao formation. Slate; Metastone and slate in alternation. Terrance deposit. Gravel; sand and mud. Alluvial deposit. Gravel; sand and mud. Suao Formation. Slate and argillite with a thin bed of metastone. Yuantoushan gneiss. Gneiss or metamorphosed granodiorites. Metabasite. Metabasite. Kangyang gneiss. Green-greenish granite and gneiss. Dekeli gneiss. Green-greenish granite and gneiss. Fanpaochienshan gneiss. Gneiss. Amphibolite. Amphibolite. Paijang Schist. Greenschist; siliceous schist and marble interlayer. Happen Marble. Thick – bedded marble. Fongshushan amphibolite. Amphibolite and horblende schist. Tungao Schist. Marble; Amphibolite; Graphite schist and quartzmica schist. Wuta Schist. Marble and chlorite schist; Chlorite schist; Quarztmica schist. Chiuchiu Marble. Thick-bedded marble, dolomite and chlorite schist intercalation. Nanaoling Schist. Quarzt schist; paragneiss. Kuyan schist. Phyllite; micaceous schist and quartz to mica schist. 18.

(31) 2.2.1 Pre-Tertiary metamorphic rock The wide-spread lithological units in the basement can basically be classified into several complex series of geological features as listed below: Kuyam schist belongs to Nanao Peak Schist (124-125.5K) located in the western Central Cross-Island Highway. Nanaoling Schist incorporated with Tungao Schist became Donao Schist at 115122.6 K along Central Cross-Island highway. In the Suao-Nanao area, they are lenticular bodies that interbedded with other kinds of metamorphic rocks. Chiuchiu Marble is accompanied by the Happen marble located in the SuaoHualien highway, north of the Heping area. They are lenticular metamorphosed limestone, which intercalated in the schistose rock (Ho, 1988). It is characterized by massive and thick-bedded and ranged from fine to coarse-grain. Paijang schist was featured by lenticular greenschist, siliceous schist, and marble alternation. The greenschist are derived from mafic volcanic flows or pyroclastic rocks. Especially, the chlorite schist of Wuta schist formation is easily accompanied by the greenschist, and it has fine to medium-grain. Amphibolite is elongated lenticular bodies that well-foliated rock distributed in the northern Nanao area. Its foliation plane is generally parallel to that of surrounding schist. Six gneiss bodies were discovered by Yen (1954) contained two large gneiss bodies in the north and others contacted with metamorphosed limestone in the south. Granitic texture, which explained for the intrusion of the granite mass into the pervasive greenschist facies metamorphism. (Fuh, 1962) is characterized as the orthogneiss.. 19.

(32) 2.2.2 Tertiary metamorphic rock and sediment deposit Nansuao and Suao formation located in the north of the study area. The rocks contained black to dark grey argillite, slate, and phyllite, and occasional interbeds with dark gray compact sandstone are characterized by the low metamorphic grade. Gravel, sand, and mud were deposited as alluvial and terrace that were distributed in the stream area.. 20.

(33) METHODOLOGY The methodology of our research was briefly summarized in Fig.3.1. There are three major parts: (1) Using the 5-meter Digital Elevation Model to identify surface fracture planes and calculate their attitudes; (2) Operating the stress inversion program to obtain the stress state at NE Taiwan from focal mechanisms data of 1991-2016; and (3) the correlation of surface fracture pattern with a given stress state can help us creating 3D Fracture Nucleation Map. Several geological definitions will be clarified in the first part of this chapter.. Figure 3.1: Flowchart of the method. 21.

(34) Geological definitions 3.1.1 Foliation In general, foliation is a geological term defined for any planar fabric or curve structure (Fossen, 2016). There are two types of the foliation: primary foliation was formed during the formation of magmatic rock or deposition of sedimentary rocks, and secondary foliation is tectonic foliation responded to the deformation process. Tectonic foliation is the penetratively structural fabric defined by zones of different grain size and flattened object or recrystallized tabular grain with the uniform orientation (Fossen, 2016). Therefore, the planar, parallel, triangular faces observed from the topography of metamorphic terrain is characterized as the existence of specified sheet texture in the 3D environment of Arc Sense. Under the surface processes, the foliated rock may easily split when the force overcomes its cohesion.. 3.1.2 Fracture Fracture in geology is a common term composed of structure features such as faults, joints, veins, dikes, and sills (Schultz and Fossen, 2008; Peacock et al., 2016). However, in this study, we proposed that the tectonic fracture is composed of only faults and joins, and they can reveal in digital elevation data. The study of Pollard and Aydin (1988) aimed to distinguish the joint and fault by the opening model and shear model. ▪. Extension fracture: where the fracture walls move apart, opening model, and may contain joints or faults.. 22.

(35) ▪. Shearing-mode fractures: a general term for discontinuity features along the wall rocks move sub-parallel to the plane of the fracture. These include fault and fault zones.. 3.1.3 Stress tensor Stress tensors are a specialized mathematic equation that can be used to define the Earth’s forces. Individual components in a stress tensor are tractions acting perpendicular or parallel to three orthogonal planes. The normal axes to the three orthogonal planes define a Cartesian coordinate system. The stress tensor has nine components. Each of them has an orientation and a magnitude. Three of these components are normal stresses perpendicular to plan, whereas the other six are shear stresses, and their traction is applied along the plane in a particular direction.. Methods 3.2.1 DEM Database In recent years, spaceborne satellite images, high attitude aircraft images, and their products, such as DEM, are highly efficient dataset able to characterize geologic features on the landform. Furthermore, it also is a reasonable source of information for the detection of crustal deformation patterns caused by faults and fractures (Wladis, 1999; Morelli and Piana, 2006). Because geology fieldwork is a costly, time-consuming, and sometimes dangerous undertaking, thus any beneficial techniques supported tool for effortless and convenient fieldwork should be considered. Therefore, in our research, we. 23.

(36) investigated the DEM of 5m raster cell data to analyze the surface fracture pattern along the coastline area of NE Taiwan under the 3D software environment. The surface patterns of the fracture plane can be recognized by geological lineaments. Therefore, we try to mapping lineaments through 2D view of the Red Relief Image Model (RRIM). Finally, 12536 lineaments were identified and picked through RRIM.. 3.2.2 Red Relief Image Model (RRIM) RRIM is a new conception of the geomorphology analysis that took a new approach of examining surface to classify geological structures. Compared with other topographic visualization methods such as contour map and shaded relief, RRIM eliminates the dependency of incident light direction and shows the 3D image without shadow (Chiba et al., 2008; Chiba and Hasi, 2016). Two landform element layers resulted in the RRIM (Chiba et al., 2008) are the topographic openness map and slope map. We used the Relief Visualization Toolbox (RVT 1.3) to create the openness map and slope map from DEM data. The RRIM was further created in Arc Map by overlap the topographic openness map and slope map. The adjustment of contrast, brightness, and transparency parameter of each layer in Arc Map software is required to clearly display RRIM.. 3.2.2.1 Topographic Openness Map The topographic openness map, including positive and negative openness maps, is an essential component to produces RRIM. Openness expressed the degree of dominance or enclosure of location on an irregular surface. The positive openness map 24.

(37) refers to a calculation of concave-downward zenith angles. Yokoyama et al. (2002) concluded that positive openness represented the topographic ridges. On the others, negative openness referred to an evaluation with concave-upward nadir angles and characterized the valley (Fig 3.2A). Therefore, the concept of openness strongly emphasizes the dominance of concavities and convexities on the landscape. The zenith and nadir angles were computed the subtraction 90○ to slope angles and calculated in eight different azimuth directions (D= 0; 45; 90; 135; 180; 225; 270; 315). The positive and negative values of the openness map were obtained by averaging angle for each direction. The openness map is designed by the gray tone to enhance its optimal detail and contrast. Positive openness takes a higher value of gray tones and generally encountered an expanse of terrain by one or several elevated relief features such as ridges (Fig 3.2B), while negative openness expressed the degree closure of the lower location as valley, river, or crater by elevated surroundings (Fig 3.2C) (Yokoyama et al., 2002). To successfully display the positive and negative map into ArcMap software, we need to adjust the contrast, brightness, and transparency parameters of each map. And the adjustments are continuously repeated until we observe clearly the topography. In summary, openness parameters highly support to highlight the different terrain features such as ridge, crests, gullies, or valleys (Fig 3.2D).. 25.

(38) Figure 3.2: Elements of topographic openness map. (A) Using the nadir point to calculate the negative openness and zenith point to identify the positive value on openness map analysis (modified from Yokoyama et al., 2002). (B) Negative openness map. (C) Positive openness map. (D) Diagram of the concave and convex topography based on the concept of topographic openness map. 26.

(39) 3.2.2.2 Slope Map The slope map was computed as the ratio of relief to irregular terrain surface and can be roughly identified topographic features of the area. If the openness map displays the gray tone in RRIM, the slope map layer is usually generated by a red color pattern because the red color has the richest tones for human eyes (Chiba et al., 2008). In similar, the contrast, brightness, and transparency parameters of the slope map also need to adjust. The slope map should be overlap with the topographic openness map and repeatedly adjust the parameter until the clearly observe topography (Fig 3.3). The yellow lineaments were mapped along the gullies and primarily indicated for the surface fracture traces in the 2D view of ArcMap software.. 27.

(40) Figure 3.3: Red Relief Image Model. Yellow lineaments presented in the RRIM are emphasized the convex topography. Lineaments correspond to the straight or slightly curved edge of convex topography and might be affected by the geological features, such as join, fault or cleavage, etc. (Koike et al., 1998).. 28.

(41) 3.2.3 Lineament Map Geological lineaments can be associated with the rupture of Earth’s surface from a different origin such as structure, lithology, surface processes, etc. They commonly demonstrate a particular type of fracturing (Hobbs, 1904; Wise et al.,1985; Twidale et al., 2007). The topographic expression of the rupture area with the offset of the surface, differential erosion of juxtaposed units, or erosion of damaged rock could help us to classify the geological structures such as faults, fault zones, or joints, etc. Therefore, mapping lineaments is a necessary step for interpreting and examining the fracture plane. Typically, lineaments displayed by straight or slightly curved components can be observed and recognized through high-resolution remote sensing data. The identification of lineaments distribution in the DEM is strongly dependent upon the effect of illumination such as azimuth, imagery declination angle, and view direction (Koike et al., 1998). However, RRIM enhances lineament quantity by the independent illumination for any viewing direction (Chiba et al., 2008). The lineament map of NE Taiwan was generated by hand-picking the lineaments through RRIM. In this study, we recognized and picked lineaments along gullies, straight river banks, or broken ridge. It basically displays in RRIM by yellow dash lines (Fig 3.3). In total, we mapped 12536 lineaments in the coastline area of NE Taiwan. The fracture planes can be identified based on the distribution of lineaments in the 3D view.. 29.

(42) 3.2.4 Extraction of Fracture Attitude RRIM overcomes the shortcoming of the illumination problem to identify the lineaments, which might be considered as a trace of the fracture plane. Fracture planes are three-dimensional and are commonly represented in the spatial by attitude (Fig 3.4). In the field investigation, the attitude of fractures is examined on the surface by using the clinometer. The strike and dip direction of the fracture surface are commonly measured. In recent years, the development of the reformative resolution of digital elevation data was able to obtain the strike and dip of the sedimentary bedding (Yeh et al., 2014). It follows the same traditional principle method that the dip angle was calculated by measuring the angle between bedding and horizontal plane (Fig 3.5). They defined that the sedimentary bedding planes in the LiDAR are measurement triangles. Each measurement triangle contains more than three points to compute a regression plane. The strike and dip of the bedding planes are acquired from the attitude of the regression planes (Fig 3.6). In this study, we tried to use the DEM with 5 meters of resolution to conduct the 3D fracture map. We also assume that fracture planes in the spatial are shaped as the triangle and should be fitted with the regression planes. At least three measurement points defined the triangle, but more than three points are needed to obtain the smooth and correctly triangle. To conduct the 3D fracture map, we floated the lineaments map into the 3D environment of Arc Sense. After that, hand-picking the triangle planes (S), which are defined by the group of green dots (Fig 3.5) in 3D environments. The attitude. 30.

(43) of fracture planes was calculated by measuring both the intersection and the inclined angle between a triangle (S) and the horizontal plane (H) (Fig 3.5). In the TM2 coordinate system, the x and y gave the longitude and latitude of the point, while z is the elevation above a given datum. The regression plane can be represented by z= ax+ by+ c. To derive the strike and dip of the measurement triangle, we take the cross product of the normal measurement vector and the normal horizontal plane vector (Fig 3.6). Calculating the gradient of two-variable function can reveal the dip direction of any point on the surface: 𝑣⃗𝑑𝑖𝑝 = 𝛻𝑧 = (𝜕𝑧⁄𝜕𝑥 )𝑖⃗ + (𝜕𝑧⁄𝜕𝑦)𝑗⃗ , where 𝑣⃗𝑑𝑖𝑝 is the fracture plane dip direction and 𝑖⃗, ⃗⃗𝑗 are the unit vector of x, y-direction. The dip angle was calculated by the following equation: 𝑆𝑑𝑖𝑝 = ‖𝑢 ⃗⃗𝑑𝑖𝑝 ‖ = ‖𝑣⃗𝑑𝑖𝑝 ⁄‖𝑣⃗𝑑𝑖𝑝 ‖‖. The strike vector can be calculated at any point over the regression surface using the following equation: 𝑣 ⃗𝑠𝑡𝑟𝑖𝑘𝑒 = 𝑛⃗𝑑𝑖𝑝 × 𝑛⃗𝑧𝑠𝑢𝑟𝑓𝑎𝑐𝑒 , where 𝑣⃗𝑠𝑡𝑟𝑖𝑘𝑒 is the strike direction, 𝑛⃗⃗𝑑𝑖𝑝 is the dip unit vector, and 𝑛⃗⃗𝑧𝑠𝑢𝑟𝑓𝑎𝑐𝑒 is the unit vector of the regression surface in z. The angle between the north and strike direction is considered as the strike value.. 31.

(44) Figure 3.4: Illustration of plane attitude from LiDAR. (A) The bedding planes. (B) The spatial vector calculation for bedding planes A and B. Two triangles are the planes regressed from the black dots (Yeh et al., 2014).. Figure 3.5: Demonstration of fracture planes from lineaments in the 3D view of Arc Sense software. The fracture plane was defined as the yellow triangle. The strike and dip angle calculated from the result of the intersection between the fracture plane (S) and the light green horizontal plane (H). Many green measurement points were conducted on the fracture plane based on several lines in the 3D environment.. 32.

(45) Figure 3.6: The calculation of fracture and foliation plane attitude in lower hemisphere projection. The dark green imaginary fracture plane is plotted in the lower hemisphere. The red, dark blue, yellow, and purple arrow indicated for the strike, dip direction, normal measurement triangle, horizontal vector. They are supported to calculate the fracture plane attitude.. 33.

(46) 3.2.5 Stress Inversion and Instability Analysis 3.2.5.1 Focal mechanism data In this study, we used 295 focal mechanisms data from earthquake events in Taiwan from 1990 to 2016 (Wu et al., 2008) for producing the cross-section of the stress state of NE Taiwan because the focal mechanism of earthquake events displays the stress regime patterns resulted from regional tectonic setting. Additionally, Scholz (2002) indicated that the aspects of the stress field might be associated with the fracture possesses in the Earth’s crust. Also, the stress fields might construct the new fracture from the intact materials or change the kinematic behavior from pre-existing fractures. For this study, we calculated the fracture instability to evaluate the development of nucleated fractures from stress evolution. The approach of reactivation potential at any existing fractures will discuss in Chapter 5 by examining slip tendency. In detail, to model the stress state of three different domains in NE Taiwan, we analyzed the movement of the seismic fault of the focal mechanism. This method solves the maximum shear stress direction on the fault, which was assumed by Wallace (1951) and Bott (1959), to determine the stress state of NE Taiwan. Furthermore, we try to calculate the three direction cosine of principal stress and relative size ( ) in the Mohr’s circle by applying the reduced stress tensor method.. 3.2.5.2 Reduced stress tensor Reduced stress tensor was first introduced by Bott. (1959), which applied to reconstruct three principal stresses axes and stress ratio ( ) from the conjugate faults system in the isotropic media. The conjugate faults are mainly caused by the consequence of deformation-induced by the stress field. The fault movement along the 34.

(47) fault plane, when earthquake happens, is assumed to be induced by stress tensor. The direction and sense of slip vector that occurred in the fault plane are assumed to be the same with those maximum resolved shear stresses. Hence, whether we knew the orientation and senses of fault movement on the fault plane, reduce stress tensor can be instituted. To solve reduced stress tensor, four unknown independent parameters should be obtained. Therefore, we primarily calculate the eigenvalue and eigenvector in order to determine the orientation of principal stress axes and stress ratio.. 3.2.5.3 Effective stress calculation In order to analyze seismic fault instability, we need to obtain the effective stress. The reduced stress tensor included only the direction of the principal stresses and the stress ratio ( ). And the stress ratio has governed the shape of Mohr’s circle but remains Mohr’s circle size and position unknown because it depends on the unknown quantities of differential stress and means stress magnitude. To determine the exact position of Mohr’s circle, Lisle and Srivastava (2004) assumed that the k1 and k2 are unknow parameters related to the absolute size and position of the Mohr’s circle, and the frictional sliding envelope is tangential to the states compatible with a known. 1 and. 3. Mohr circle. Thus, the stress. were calculated by:. 𝜎1 = 𝑘1 + 𝑘2 𝜎2 = 𝑘Φ + 𝑘2. (Eqn. 3.1). 𝜎3 = 𝑘2 The effective stress was calculated by applying the assumption of Mohr’s circle tangent with the Mohr-Coulomb failure criteria. Therefore, the notable magnitude of. 35.

(48) principal stresses compatible can be calculated with the rock mechanics data and stress ratio and orientations. The average frictional angle and rock cohesion for the failure criteria basically were calculated from the rock strength database with area weighting are 49.5○, 7.2MPa, respectively (Table.3.1). Besides that, the effective stress is considered. Pore pressure (p) static increases with depth. The rock density is 2.7 g/cm3. Consequently, we can estimate the instability of fracture varies as a function of its orientation concerning the magnitude of the principal stress. Table 3.1. The summary of rock mechanics data (SNFD-ITRI-TR2016-0002V5PN_SNFD2017). Area. Friction angle. Cohesion. Tensile strength. (m2). (o). (MPa). (MPa). 4306201. 57.1. 13.8. 8.41. Schist. 118651907. 48. 4.11. 8.38. Gneiss. 13532320. 56.8. 32.63. 10.44. 136490428. 49.16. 7.24. 8.58. Lithology Marble. 3.2.5.4 Fracture Instability The differently oriented faults have different susceptibility to be active and be individually unstable in the given stress fields. To evaluate the fault instability, Vavrycuk et al. (2015) primarily determined the unstable fault plane from the focal planes by quantified with the Mohr-Coulomb failure criterion. According to this criterion, shear traction (𝜏) on an activated fault must exceed a critical value (𝜏𝐶 ) (Eq. 3.2) (Beeler et al., 2000; Scholz, 2002), which is calculated from cohesion C, fault friction 𝜇, compressive normal stress (𝜎𝑛 ), and pore pressure (p) (Eq. 3.3).. 36.

(49) ∆𝜏 = 𝜏 − 𝜏𝐶 Where. (Eqn. 3.2). 𝜏𝐶 = 𝐶 + 𝜇(𝜎𝑛 − 𝑝). (Eqn. 3.3). They demonstrated that the distribution of fault planes inside the unstable area (red area in Fig 3.7) corresponded with the principal fault plane of the focal nodal plane. The concept of instability (I) was introduced, and fault instability of all fault orientation can be defined in the range from 0 to 1 (Fig 3.7) by the following formula (Vavrycuket al., 2015).. 𝐼=. 𝜏 − 𝜇(𝜎 − 𝜎1 ) 𝜏𝑐 − 𝜇(𝜎𝑐 − 𝜎1 ). (Eqn. 3.4). Where 𝜏𝐶 , 𝜎𝐶 is the shear, and effective normal stresses along the principal fault plane from focal mechanism data and 𝜏, 𝜎 is shear, normal stress along the DEMderivation surface fracture.. 37.

(50) Figure 3.7: Fault instability in Mohr ‘s diagram. The red area marked unstable areas of principal fault planes from focal mechanisms, and blue dots characterized the principal fault instability I=1. The red dot marks an arbitrarily oriented fault with instability I. τ, σ are the shear and effective normal stresses, respectively. σ1, σ2, and σ3 are the effective principal stresses. (Modified from Vavrycuk et al., 2015).. 38.

(51) RESULTS Identification of fracture and foliation On a large scale, the consequence of the deformation structure and weathering process might support the significant modification of landform by several activities such as uplift, crustal thinning, surficial erosion, or deposition process. These activities can occur in a variety of deformation styles, such as ductile, ductile-brittle, or brittle modes. The cutting-relationship among individual fracture planes or between fracture and foliation is considered as a priority to characterize the fracture and foliation planes in the metamorphic area. The cross-cutting relationship also provides the relative age information of fracture sequences. Therefore, four locations were picked up at different domains (Fig 4.1) to demonstrate the relationship between fracture and foliation plane and its characteristics.. 4.1.1 Attitude of foliation planes In the 3D view of Arc Sense, several foliation planes were exhibited in Figure 4.2A. And, purple lineaments (Fig 4.2B) were identified and mapped through RRIM in the same area (Fig 3.3). The foliation planes were highlighted by the small dark blue triangle planes in the DEM data (Fig 4.2C). The 3D parallel, planar, triangle faces are characterized as foliation planes in 3D view. Figure 4.2D indicated that the features of planar morphologic geometry of foliation (purple lines) in the ridge where was truncated by fracturing (red and green arrows). Fracture planes with NW-SE strike and west-dipping (red arrows) and ENE-WSW strike fracture planes (green arrows) cut and broke the foliation, and foliation became the discontinuously geological features (Fig 4.2D). Although in detailed observation, the fracture surfaces are gently unsmooth due. 39.

(52) to the influence of surface process, the fracture plane still is able to be mapping. In summary, the foliation and fracture planes can be observed by using the DEM with GIS relevant applications. Furthermore, their attitudes have the potential to calculate if they are clearly exposed on the surface by one plane. The calculation of foliation attitude was following a similar approach to fracture planes. The results were examined in three representative locations for three domains of this study area. The foliation characteristics in the 3D environment are shaped by parallel triangle faces, and, in comparison, the results of the foliation attitude calculation from DEM were similar to attitudes with field data from the geological map (Fig 4.2E). Therefore, the research methodology of this study is the potential to apply in the calculation of fracture attitudes.. Figure 4.1: Location map of geologic features in the study area.. 40.

(53) Figure 4.2: Foliation traces and foliation planes in the 3D view. (A) 5m DEM illustrated the regional foliation aspects in schist rock of the NW Heping area. (B) Purple lines were mapped as the foliation traces. (C) Regional foliation planes. (D) Cuttingrelationship between foliation system and fracture system. (E) The diagram showed that each stacked triangle reflected the foliation plane and the lower hemisphere projection presents a similar attitude between foliations from DEM-derivation and foliation from fieldwork. 41.

(54) 4.1.1.1 Domain I The systematic parallel solid purple lines in NE-SW trending are interpreted as the trace of the foliation plane and have been traced along the gullies of the Dongao greenschist area in the 3D environment of DEM data (Fig 4.3A). Foliations are gently west-dipping with a dip angle of 48○-58○ (Table 4.1). It can continually extend several hundred meters if the surface process does not modify the landform. The attitude of foliation from nearly field data (Table 4.2), which are marked by green dots (Fig 4.3A), was selected to compare with our foliation attitude calculation results. The lower hemisphere projection in Figure 4.2B illustrated that the green great circles presented the foliation attitude derived from field data, whereas the purple great circles described the foliation plane from DEM (Fig 4.2B). It shows a strong similarity of attitude between field mapping and calculation. Also, poles to the plane displayed a similar distribution of foliation in this area between field data and DEM-derivation. In this domain, the red arrows indicated the traces of NE-SW striking fracture planes that cut the NW-SE purple foliations. According to the cutting-relationship between the fracture traces and foliation, we can distinguish fractures from foliation.. 42.

(55) Figure 4.3: Comparison of foliation attitude between DTM-derivation and field data at the location 1. (A) The NW- SE purple lines indicated the foliation plane and the green dots referred for the existing foliation at the geological map. (B) Comparison results of foliation attitude between our interpretation from DEM (dark green dots) and field foliation data from Geological map (purple dots) in the lower hemisphere projection.. 43.

(56) Table 4.1: The foliations attitudes calculated from the DEM data at Domain I.. No.. Strike. Dip. Long. Lat. R2. 13. N. 76. W. 48. W. 331747. 2714714. 1. 14. N. 83. W. 49. W. 331641. 2715032. 0.99. 15. N. 78. W. 47. W. 331747. 2715133. 1. 16. N. 84. W. 49. W. 331904. 27213. 1. 17. N. 70. W. 52. W. 332842. 2713772. 0.98. 18. N. 77. W. 59. W. 332915. 2714096. 0.99. 19. N. 74. W. 58. W. 332944. 2714024. 0.99. Table 4.2: Foliation attitude from the geological map. No.. Dip Direction. Dip. 783. 009. 63. 785. 186. 39. 786. 189. 63. 820. 184. 79. 789. 193. 73. 787. 194. 64. 44.

(57) 4.1.1.2 Domain II In the Heping area, we clearly observed the predominance of W-E striking lineaments. The representative lineament illustrated by purple, parallel lines (Fig 4.4A) supposedly indicates foliation traces located in the west flank of the ridge were cut by three distinctively fracture planes numbered by 1, 2, and 3 (Fig 4.4A). This is these fracture planes also cut across the ridges, and normally explained why we could clearly observe numerous discontinuous lineaments. Therefore, the planes cut the foliation traces, and together cut across the ridges on a large scale can consider as the fracture planes (Fig 4.4A). Field investigations are often difficult to perceive in highly forested metamorphic terrain. Although, only two field data were documented at two locations near the river area (Fig 4.1) and simply marked by green dots in the 3D scene of Arc Sense (Fig 4.4A). However, these valuable field data gave us an opportunity in foliation attitude comparisons. The attitudes calculation of foliations from DEM derivation are closely vertical and ranged from 62○ to 72○ (Table 4.3), and in the comparison, they are compatible with the foliation attitudes from field data (Fig 4.4B). N-S striking and west-dipping fractures numbered by 1, 2, and 3 (Fig 4.4A) were generally used to distinguish the foliation of metamorphic rock.. 45.

(58) Table 4.3: The distribution of foliation attitudes calculated from DEM data of Domain II.. No.. Strike. Dip. Long. Lat. R2. 1. N. 76. E. 72. W. 323948. 2704302. 0.84. 2. N. 72. E. 66. W. 323890. 2704421. 0.96. 3. N. 80. E. 60. W. 323806. 2704525. 0.95. 4. N. 72. E. 75. W. 323928. 2704176. 0.99. 5. N. 81. E. 71. W. 323784. 2704758. 0.8. 6. N. 73. E. 62. W. 325737. 2705467. 0.97. 46.

(59) Figure 4.4: Comparison of foliation attitude between DTM-derivation and field data at the location 2. (A) The NE-SW systematic purple parallel lines were identified as the traces of foliation in the Heping area, and foliation was cut by fracture surfaces numbered by 1, 2, and 3. (B) Stereonet showed the foliation attitude determined from DEM data and field data in the geological map.. 4.1.1.3 Domain III We identified the NE-SW striking and NW dipping foliation system at location 3 where it is greenschist (Fig 2.4) near the lithological transition area between marble and greenschist. The purple lines, which are displayed in the 3D environment of Arc Sense, are illustrated for the foliation traces (Fig 4.5A). The similarity of attitudes is observed in the lower hemisphere (Fig 4.5B).. 47.

(60) The dip angle of foliation planes with ranged from 41○-52○ was cut by one fracture plane. This fracture plane is composed of three planar surfaces numbered by 1, 2, and 3 and has attitude in W-E striking and south-dipping (Fig 4.5A). The fracture planes with strikes in NW-SE in this area were denoted by yellow arrows (Fig 4.5A). Those fracture planes are considered as the different fracture systems from previous W-E striking fracture planes. It could denote disparate deformation phases with the W-E fracture system due to the evolution of the tectonic setting.. 48.

(61) Figure 4.5: Comparison of foliation attitude between DTM-derivation and field data at the location 3. (A) The regional foliation with NW – SE striking (purple linear) was cut by a NE-SW trending fracture plane, which was distinguished by three planar planes numbered by 1, 2 and 3. (B) The great circle in stereonet presented a similar strike and dip angle of regional foliation calculated from DEM and the attitude of foliation in geological map.. 49.

(62) Table 4.4: The distribution of foliation attitudes calculated from DEM data of Domain III.. No.. Strike. Dip. Long. Lat. R2. 1. N. 48. E. 41. W. 308956. 2685017. 0.99. 2. N. 49. E. 41. W. 309033. 2684946. 0.98. 3. N. 33. E. 52. W. 309283. 2684879. 0.92. 4. N. 38. E. 56. W. 309182. 2684922. 0.98. Table 4.5: The distribution foliation attitudes of field data analyzed from the geological map.. No.. Dip Direction. Dip. Long. Lat. 386. 290. 28. 308413. 2684239. 387. 310. 36. 307750. 2684787. 382. 316. 30. 310948. 2686848. 383. 312. 28. 310235. 2686456. 50.

(63) 4.1.2 Fracture plane mapping 4.1.2.1 Interpretation of fracture plane The surface fracture set, in combination with the ductile foliation, define the regional morphotectonic lineament pattern. Therefore, Florinsky (2016) dedicated that the mapping morphological lineaments can be used to recognize geologic structures, especially the lineaments located along the straight river bank, across the ridge, etc., are able to present the plane of the fracture. The planes, when breaking the slope or scarp during the action of the tectonic or non-tectonic process supposedly are planar (Fig 4.6A) and exposed on the surface. In a 2D map, lineaments located along the plane supposedly is “traces” (Fig 4.6B). Tracking the “traces” in the 3D environment of Arc Sense can recognize the fracture plane and even evaluate the contribution and development of the fracture system in the 3D model (Figure 4.6C). Along lineaments, because the supposed fractures are represented the sensitive and risky area with the erosional surface process, the topography is easy modification and forms the gullies topography. Therefore, tectonic lineaments along the gullies structure implied the fracture planes. The paths of gullies are significant features to recognize the neotectonics fracture (Chen et al., 2016).. 51.

(64) Figure 4.6: Characteristics of the fracture plane in the 3D view of Arc Sense. (A) The footprint of the fracture plane on the surface significantly correlates with the distribution of lineaments. (B) Four lineaments (purple lines) along with the fracture footprint. (C) Four fracture planes are mapped based on the fracture trace in 3D view.. 52.

(65) 4.1.2.2 DEM-derived fracture plane In DEM mapping, 187 fracture planes were identified, mapped, and calculated their attitudes in three different domains. In the Domain I, 50 fracture planes was recognized with the orientation dominated in ENE-WSW striking, and the mean fracture plane was evaluated and displayed by the ENE-WSW striking great black circle in the lower hemisphere projection. 65 & 72 fractures in Domains II and III were recognized, respectively. The attitudes of fracture planes at Domain II were gently clockwise rotated to NE-SW compared with the attitude of N-S striking of fracture in Domain III. The evidence in the lower hemisphere illustrated 30○ differences in the strike of mean fracture attitude of Domain II and III (Fig 4.7).. 53.

(66) Figure 4.7: Spatial distribution of surface fracture planes in DEM data. 187 fracture planes were identified, mapped, and calculated their attitude in three different domains of NE Taiwan via used the 5m DEM data. 50, 65, and 72 fractures were mapped in Domain I, II, and III, respectively. The lower hemisphere projection represented the distribution of the fracture plane by the contour map. The mean fracture attitude was calculated and displayed by the mean pole (red dot) and the great mean circle (black curve lines) in different domains.. 54.

(67) 4.1.2.3 Restriction of Selected Fracture Plane The strike and dip of DEM-derived fracture planes are evaluated from each measurement triangle. Triangle selection in the study area is based on the following restrictions for better results of representative fracture plane results. Reliability of regression calculation During calculating the fracture attitude, the regression reliability (R2) values of the measured triangle can also be determined. Yeh et al. (2014) indicated that the selected triangles with high R2 have the goodness of fit for the regression plane. Therefore, the fracture planes were classified into ten groups based on different R2 values (Fig 4.8). The dark green to light green lines represented the fracture plane within lower R 2 (0.030.4), the yellow to orange line represented fracture planes of middle R 2 (0.41-0.8), fracture planes of high R2 (0.81-1) values were indicated by red lines. Totally, in three domains, 120 measurement triangles with R2 > 0.85 based on the above restrictions are used for the instability analysis.. 55.

(68) Figure 4.8: Classification of the fracture planes based on R2 value. 120 selected fracture planes (red lines) are reasonable for the input database of tendency analysis.. 56.