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
25
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).
26
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
27 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.
28
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).
29 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.
30 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
31
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.
32
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.
33
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.
34 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
35
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 1 and 3 Mohr circle. Thus, the stress states compatible with a known were calculated by:
𝜎
1= 𝑘
1+ 𝑘
2𝜎
2= 𝑘Φ + 𝑘
2 (Eqn. 3.1)𝜎
3= 𝑘
2The effective stress was calculated by applying the assumption of Mohr’s circle tangent with the Mohr-Coulomb failure criteria. Therefore, the notable magnitude of
36
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-0002-V5PN_SNFD2017). 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).
37
∆𝜏 = 𝜏 − 𝜏𝐶 (Eqn. 3.2) Where 𝜏𝐶 = 𝐶 + 𝜇(𝜎𝑛− 𝑝) (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 DEM-derivation surface fracture.
38
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).
39
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
40
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.
41
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) Cutting-relationship 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.
42 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.
43
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.
44
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
45 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.
46
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
47
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
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