Aims & Scope

Journal of the Geological Society of Korea - Vol. 57 , No. 6
[Paper List] [Go to Volume List]

[ Article ]
Journal of the Geological Society of Korea - Vol. 57, No. 6, pp. 809-826
Abbreviation: J. Geol. Soc. Korea
ISSN: 0435-4036 (Print) 2288-7377 (Online)
Print publication date 31 Dec 2021
Received 25 Sep 2021 Revised 04 Nov 2021 Accepted 09 Nov 2021
DOI: https://doi.org/10.14770/jgsk.2021.57.6.809

이산요소법을 이용한 습곡-충상단층대 모형의 변형률 분석 방법으로써 유한차분법 및 유한요소법의 적용성 평가
안수정 ; 소병달‡
강원대학교 지구물리학과

Applicability of finite difference and finite element methods as a post-processing tool for strain analysis of fold-and-thrust belt model developed with the discrete element method
Soojung An ; Byung-Dal So‡
Department of Geophysics, Kangwon National University, Chuncheon 24341, Republic of Korea
Correspondence to : ‡+82-33-250-8582, E-mail: bdso@kangwon.ac.kr

Funding Information ▼
초록

본 연구는 불연속적인 입상 물질 모사에 널리 활용되는 이산요소법을 이용하여 취성(brittle)-소성(plastic) 변형이 광역적으로 발생하는 습곡-충상단층대 모형을 개발하였다. 본 이산요소 모형의 기하학적 구조는 수치 및 상사 모델링 기법으로 수행된 기존 연구들과 일치함을 확인하였다. 모형의 변형률을 분석하기 위해 연속체 역학(즉, 유한차분법 및 유한요소법)이 도입된 파이썬 기반 후처리 코드를 제작하였다. 이를 통해 습곡-충상단층대의 실제적인 취성-소성 변형을 모사하고 물질 대변형에 사용되는 유한변형률을 통해 단층의 운동 감각, 방향에 따른 변형의 규모, 체적변형률, 왜곡도를 구하는 등 물리적·정량적 분석을 성공적으로 수행하였다. 유한차분법을 사용한 변형률 계산 결과는 경계부에서는 부정확하나 전체적으로 균일한 분포를, 유한요소법에서는 경계부를 따라 해석해와 비교하였을 때 정확도가 높은 결과를 보였다. 이를 통해 본 연구는 유한요소법이 지질 구조의 불규칙한 경계 내 변형률을 계산하는데 효율적인 방법임을 제시하였다. 본 연구에서 적용한 변형률 계산법은 습곡-충상단층대와 같은 대변형이 일어난 지질 구조의 물리적 해석에 최적화되었다. 또한 퇴적물의 마찰성 입상 거동(예, 쪼개짐, 회전, 교란, 교결, 압밀 등)을 내재적으로 반영하는 이산요소법의 특성을 고려하면, 본 연구에서 개발한 변형률 계산 방법은 지표면 내 퇴적물 거동, 불연속 구조의 변형률, 시간에 따른 퇴적층의 부피 변화(예, 침하, 팽창) 등을 이해하는데 활용 가능하다.

Abstract

We develop a fold-and-thrust belts model with extensive brittle-plastic deformations using the discrete element method (DEM), widely adopted for simulating discontinuous granular materials. We confirm that our models are consistent with previous numerical and analog studies in terms of the surface angle, the spacing between two near-thrusts, and the number of thrusts. Furthermore, we build the Python post-processing code for analyzing the strain evolution within the model based on continuum mechanics, i.e., the finite difference method (FDM) and finite element method (FEM). Physical and quantitative analysis such as kinematic sense of fault, the amount of deformation, volumetric strain, and distortion is successfully performed by adopting the finite strain for large geological deformations. The strain calculated based on the FDM is inaccurate at the boundaries but shows a uniform distribution throughout the model, and the FEM shows high accuracy which is comparable with the analytic solution along the boundaries. This suggests that the FEM provides a more efficient tool for analyzing strain fields in geological domains with irregular boundaries such as fold-and-thrust belts. The strain calculation method applied in this study is optimized for the physical analysis of geological structures that have undergone large deviatoric deformations. Moreover, since the DEM inherently includes the frictional granular behaviors of the sediment (e.g., fracture, rotation, disturbance, cementation, and consolidation), our strain analysis method guarantees high feasibility in various cases, such as the observation of the sediments behavior in the surface, the strain analysis of discontinuous structures, and the volume change of sediment layers over time (e.g., settlement and expansion).

Keywords: fold-and-thrust belt, discrete element method, finite difference method, finite strain tensor, finite element method
키워드: 습곡-충상단층대, 이산요소법, 유한변형률텐서, 유한차분법, 유한요소법
1. 서 론

습곡-충상단층대(fold-and-thrust belts)는 조산대(Metcalf et al., 2009; Forte et al., 2013; Vaseghi et al., 2016; Roy et al., 2020)나 섭입대(Dahlen, 1984; Park et al., 2002; Borderie et al., 2019; Dominguez et al., 2000)와 같은 수렴형 경계에서 발달하는 대표적인 천부 취성 변형 구조이다(Buiter et al., 2016). 조산대에서는 전면지(foreland)에서 후방지(hinterland)로 향할수록, 섭입대에서는 해구에서 내륙으로 향할수록 단층 경사각이 증가하고(Davis et al., 1983; Park et al., 2002) 지각이 두꺼워져 습곡-충상단층대나 부가쐐기(accretionary wedge) 구조를 이룬다(그림 1; Ruh et al., 2018; Dominguez et al., 2000). 습곡-충상단층대 진화에 수반되는 부가 쐐기는 불도저가 흙더미를 밀었을 때 그 모양이 쐐기처럼 되는 현상과 같은 원리로 생성된다(Dalhen et al., 1984; Ruh et al., 2013). 시간에 따른 쐐기 구조 진화를 역학적으로 설명하는 임계 경사 이론(critical taper theory)에서는 쐐기 구조에 퇴적물이 부가될 때 구조적으로 안정된 상태를 유지하기 위해 시간에 따른 표면 경사각 α와 바닥 경사각 β의 합인 α + β가 보존된다(Dahlen, 1984). 습곡-충상단층대 및 부가쐐기 등 쐐기 형태로 진화하는 구조에 대한 이론 및 실험 모델은 임계 경사 이론에 쿨롱 파괴 기준(Coulomb failure criterion)을 도입(Chapple, 1978)하여 쐐기 구조의 이론적 진화 양상 뿐만 아니라 퇴적층의 취성변형까지 모사하였다(Lohrmann et al., 2003; McClay et al., 2004). Bonini (2007)는 상사 모델링 기법을 통해 습곡-충상단층대에서 횡압력을 받은 마찰성 미고결 퇴적물들은 분리단층(décollement)의 바닥면을 밀며 충상단층(thrust), 역방향 충상단층(back thrust), 돌출(pop-up), 융기(uplift), 습곡(fold), 비탈(ramp), 듀플렉스(duplex) 등 다양한 소규모 압축 구조를 복합적으로 발생시킴을 확인하였다.


Fig. 1. 
Schematic diagram of critical taper theory regarding the evolution of an orogenic wedge. The critical taper angle (α + β) for sustaining wedge stability is defined by the summation of the surface angle α and dip angle β. As compressional stress is applied to the foreland from the hinterland, faults splay out from the décollement, and a tectonic wedge appears in the foreland. The critical taper angle is conserved during the deformation.

유동학과 암석역학에서 습곡은 큰 규모의 단열없이 소성 변형이 영구적으로 누적된 상태의 연성전단대(ductile shear zone)이며, 충상단층은 응력이 모암의 항복강도를 초과하여 비가역적으로 파괴되는 취성전단대(brittle shear zone)에 속한다(Roy et al., 2020). 단열은 주응력에 수직이나 전단방향으로 갈라진 틈으로서(Liu, 2005), 퇴적물 입자의 전위나 회전을 일으킨다. 단열의 누적은 단층을 만들고 그 주변에 변형 집중화(strain localization) 현상을 야기한다. 미소 단열이 산발적으로 발생하나 구조 형태 파괴가 일어나지 않는 선에서 그치는 연성전단대와 달리, 취성전단대에서는 단열이 공극 유체의 통로(Stipp et al., 2013; Yamada et al., 2014)가 되어 윤활 작용이 일어나 단층대의 강도가 낮아지거나(Hori and Sakaguchi, 2011), 변형 집중 과정에서 전단열 증대(Regenauer-Lieb and Yuen, 1998; Jo and So, 2020) 및 마찰력 약화(Goldsby and Tullis, 2011; So and Capitanio, 2020)가 발생할 수 있다. 특히 일부 섭입대의 후방지 말단에서 시작된 물매 단층(splay fault)으로부터 전면지 말단까지 이어진 습곡-충상단층대는 여러 갈래로 갈라져 있으며 분리단층을 경계 삼아 상부 맨틀과 지각 사이 퇴적물과 공극유체 등 다상 물질의 이동 경로가 된다(Moore et al., 1990; Branquet et al., 1999). 횡압력과 정암압에 의해 과압 상태가 되어 상부지각으로 빠져나온 섭입 환경 내 슬랩의 유체는 습곡-충상단층대 내부에서 활발한 열역학 및 유체역학적 상호작용을 하면, 결과적으로 단층면에 작용하는 유효수직응력을 감소시켜 메가스러스트 지진의 발생원인이 되기도 한다(Sibson, 2020; Dorsey et al., 2021). 이러한 맥락에서, 습곡-충상단층대의 변형률 분포 계산은 단층 변위 규모와 운동 감각을 보여주며 더 나아가 단층대 강도의 열적 및 기계적 약화 기작을 추론하는 열쇠가 된다.

지질/지구물리학적 자료를 통해 변형률을 계산하는 방법은 측정 대상 규모에 따라 크게 상시 GPS (Global Positioning System) 관측과 변형타원체법으로 나뉜다. 조산 및 섭입 운동의 결과로 수반되는 습곡-충상단층대 구조는 일반적으로 전체 폭이 약 수십-수백 km 규모로 발달한다(Glasmacher et al., 2004; Kendall et al., 2020). 상시 GPS 관측을 통해 전체 구조에 대한 작용 힘의 방향, 속도, 변형률을 계산할 수 있다(Forte et al., 2013; McCaffrey et al., 2016). 반면 연구 지역 노두 상에서 관찰 가능한 소규모 변형 구조에서, 변형된 암석 또는 박편 상의 임의의 절단면은 가상의 응력변형타원을 가진다(Lisle, 1985; De Paor, 1990). 각 노두 상의 변형타원체를 통해 구한 변형률을 지질도 상에 표시했을 때 드러나는 변형률 분포로 전체 지질 구조를 해석할 수 있다(Ramsay et al., 1983; Burmeister et al., 2009). 앞선 두 방법으로는 연구 지역의 현재 지질학적 변형률 상태를 알 수 있으나 과거에 시작된 지질 구조 진화와 그에 따른 변형률 변화 이력을 추적하기 어렵고, 획득 가능한 관찰 자료의 개수에 한계가 있어 양질의 결과를 얻기 위해 자료의 보간이 불가피하다.

최근 암석/광물학 실험, 지진파 토모그래피 자료, 원격 탐사 자료 등 지질/측지학적 근거와 고속 컴퓨팅 기술을 융합한 수치해석법이 습곡-충상단층대를 포함한 여러 지구동역학적 현상 모사에 사용되고 있다(Reuber and Simons, 2020). 습곡-충상단층대의 경우, 수 미터 ~ 수백 킬로미터의 공간적 규모와 수십 년 ~ 수백만 년의 시간을 통해 진화하며(Cosgrove and Ameen, 2000) 지배방정식(예, 질량보존방정식, 운동보존방정식, 에너지보존방정식)을 이용한 수치해석법은 이를 동시에 다루기에 적합하다(Buiter et al., 2016; van Zelst et al., 2021). 이에 따라 연속체 역학을 기반으로 한 유한차분법과 유한요소법이 활발히 사용되고 있다. 그러나 Byerlee 법칙(Byerlee, 1978)을 따르는 마찰성 물질에 적용할 경우, (1) 취성 파괴 등 급격한 동적 상황에서 불안정한 계산의 수렴(Li and Gosh, 2006; Ma and Elbanna, 2019), (2) 물질 경계면의 변형에 따라 발생하는 과도한 격자 변형(Espinosa et al., 1998; Kaus et al., 2010), (3) 매질의 연속성 가정에 위배되는 단층의 본질적인 계산 한계(Gurnis et al., 2000; Lee and So, 2020)가 단점으로 언급되어왔다. 앞선 이유로, 유한차분법과 유한요소법은 추가적인 기법 없이 마찰성 물질이 취성거동을 하는 상부지각(깊이 15 km 이내)의 변형을 모사하기에 제한적이다(Feng et al., 2015; Zhao and Gao, 2014).

최근 이산요소법(discrete element method)을 도입한 습곡-충상단층대 수치 모사 연구가 활발히 진행 중이다(Buiter et al., 2006, 2016; Hardy and Finch, 2005; Egholm et al., 2007; Dean et al., 2013; Morgan, 2015; Meng and Hodgetts, 2019; Li et al., 2021). 이산요소법은 습곡-충상단층대를 구성하는 퇴적물 거동 해석에 적합한 수치모사법으로써, 운동 방정식(뉴턴 제2법칙)과 선형 탄성식에 따라 입자 간 상호작용을 계산한다. 연속된 격자가 아닌 독립된 입자 자체가 계산 요소이기 때문에 대규모 암반 파쇄와 미시적인 지형 변화 등 비선형 변형을 안정적으로 모사할 수 있다. 그러나 유한요소법이나 유한차분법을 사용한 습곡-충상단층대 모형에 비해, 이산요소법을 활용한 경우에서는 변형률 계산에 추가적인 단계가 요구된다. 각 입자의 변형률을 구하기 위해서는 변형구배텐서(deformation gradient tensor) 계산이 선행되어야 한다. 변형구배텐서는 입자의 변위를 공간증분에 대해서 편미분하여 얻어진다. 변위가 연속함수일 경우 해석적 방법으로 미분이 가능하지만, 이산요소모형의 경우 변위 분포가 일정하지 않고 불연속하여 미분을 하려면 수치 해석 기법이 필요하다. 이산요소모형의 변형률 계산을 위해 다양한 방법으로 연속체 역학을 도입한 연구 사례들을 소개한 O’Sullivan (2011)에 따르면, 연속된 변위장을 만들기 위해서는 불연속한 입자의 변위를 임의의 격자망으로 전사(mapping)하는 과정이 필요하다. 유한요소법은 이산요소모형에 사용된 입자 중 인접한 세개의 중심점을 서로 이어 만든 삼각형 요소로 구성된 격자를 사용한다. 따라서 격자의 모양이 실제 변형된 입자의 위치를 그대로 반영한다. 반대로 유한차분법에서는 이산요소모형의 입자 위에 균일한 간격으로 이루어진 가상의 격자망을 덧씌우기 때문에 격자의 모양이 변형된 입자 분포와 무관하다. 두 방법으로 조성된 격자 내 절점에 변위 정보를 부가하여 변형률 계산을 위한 미분을 수행하였다. 국내에서 이산요소법은 주로 절리암반 내 공동 주변의 전단영역 및 균열간극 추정 모형(Hong et al., 2021) 등 암반공학/안전성평가 등 공학적 목적으로 활용되었다. 그러나 국내에서 아직까지 이산요소법을 활용한 습곡-충상단층대 진화의 수치 모형 연구 및 변형률 분포 분석은 활발히 이어지지 않았다.

이에 따라 본 연구에서는 이산요소법을 이용한 습곡-충상단층 진화 모형을 제작하고, 연속체 역학 기반 수치계산법인 유한요소 및 유한차분법을 사용하여 물질 대변형의 변형률을 나타내는 유한변형률텐서(finite strain tensor; Zienkiewicz and Taylor, 2000)를 구하였다. 파이썬 프로그래밍 언어를 이용하여 변형률 해석 코드를 작성하였으며, 신뢰성 검증을 위해 해석해가 존재하는 모형 문제를 벤치마크 하였다. 벤치마크는 작성된 변형률 해석 파이썬 코드를 통해 모형 문제의 유한변형률텐서를 수치적으로 계산한 후 해석해와 비교하는 방식으로 수행되었다. 변형률 계산 코드는 격자 구성 방식에 따라 두개의 방법(유한차분법, 유한요소법)으로 나누어 개발되었다. 이후 타당성이 확인된 변형률 계산 코드를 습곡-충상단층대 이산요소 모형에 적용하고 압축률에 따른 운동학적 및 기하학적 구조 진화와 변형률 분포를 동시에 해석하였다.

2. 연구방법
2.1 이산요소법

입자법의 일종인 이산요소법은 암석의 역학적 거동을 모사하기위한 수치해석기법으로써 1971년에 Cundall에 의해 처음으로 제안되었다(Cundall, 1971). 현재는 분말 역학(Jerier et al., 2011), 분자 동역학(Depta et al., 2018), 기계 공학(Fleissner et al., 2007), 재료 공정(Mishra, 2003) 등 유동성 입자 군집을 대상으로 한 연구 분야에 다양하게 활용되고 있다. 이산요소법의 계산 영역은 연속된 격자가 아닌 서로 독립된 유한개의 입자 집합이다. 개개의 질량을 가진 입자들은 외부 힘(예, 중력, 충돌 힘)과 내부 힘(예, 관성, 회전 모멘트)의 합력에 따라 동시다발적으로 병진 및 회전 운동을 한다. 이산요소모형의 물질 내 입자는 다른 입자와 접촉되어 있거나 접촉이 파괴되어 있다. 접촉과 접촉 파괴를 설명하는 가장 간단한 역학 모형에는 질량-스프링 체계가 있다. 이는 질량을 가진 입자와 강성이 부여된 탄성 스프링이 직렬로 연결된 구조를 가지며, 외부 힘에 의해 탄성 한계 이상의 응력이 쌓인 스프링의 파괴 거동을 직관적으로 모사하기에 적합하다. 스프링의 높은 강성은 원상태를 유지하려는 복원력이 강함을 의미하며 이는 입자 간의 접촉력과 비례한다. 일반적으로 이산요소법을 활용한 지질구조모형에서는 시간에 따른 결합력의 감쇠 효과를 더하기 위하여 그림 2와 같이 질량-스프링 체계에 완충기를 병렬로 연결한 모형인 질량-스프링-완충기 체계를 사용하고 있다(Seyferth and Henk, 2004; Chiu et al., 2015; Furuichi et al., 2018).


Fig. 2. 
Mass-spring-damper system with the contact point between particles. Particles with individual masses are in contact; the red circle indicates the contact point. Each contact point comprises a pair consisting of spring and damper, representing restoring (green arrow) and damping forces (blue arrow), respectively. Particles attract each other by the restoring force of spring. In contrast, the damping force works in the opposite direction to buffer the effect of restoring force and stabilize the numerical system.

2.2 유한변형률텐서

물질의 변형은 크게 전위변형과 회전변형으로 나뉜다. 공학적인 목적으로 주로 도입하는 미소변형률은 무한소 단위의 전위변형만을 고려한다. 하지만 지질 구조 중 단층과 같이 대규모 변형을 겪는 물질의 변형률 계산에는 전위변형뿐만 아니라 회전변형까지 감안하는 유한변형률을 사용한다. 라그랑지안 관점의 유한변형률이론은 입자의 초기 위치를 기준공간좌표계로 사용한다. 이에 따라 시간이 t, 입자의 변형 전 좌표가 X, 변형 후 좌표가 x 일 때 변위장 u(X, t)는 x-X(x, t)로 정의된다. 유한변형률은 변형 전 좌표에 대해 변형 후 좌표를 공간 방향에 따라 미분한 변형구배텐서(F)로부터 유도된다. 2차원 공간에서 입자의 초기 위치와 시간에 대한 함수인 변형구배텐서는 F(X, t)로 표현된다. 이를 텐서표기법으로 나타내면 식 (1)과 같으며 i, j는 공간 방향을 표시하는 색인으로서 1은 수평방향, 2는 수직방향을 의미한다.

(1) 

이를 크로네커 델타(Kronecker delta)인 δi,j와 변위 u로 나타내면 식 (2)와 같다.

(2) 

라그랑지안 관점의 유한변형률텐서는 그린-라그랑지안 변형률텐서(Green-Lagrange strain tensor)라고도 불리며 식 (3)과 같다.

(3) 

이를 정리하면 식 (4)가 유도된다.

(4) 

습곡-충상단층모형 내 변형률의 물리적 해석을 위해 앞서 구한 변형구배텐서 및 유한변형률텐서로부터 체적변형률(volumetric strain)과 왜곡도(distortion)를 각각 계산하였다(Morgan, 2015). 변형 전 부피를 V, 변형 후 부피가 ν일 때 체적변형률 J는 자코비안으로 식 (5)와 같다.

(5) 

편차유한변형률(deviatoric finite strain, E')은 유한변형률에서 정적 변형률(hydrostatic strain,Ehyd)을 뺀 값으로서 텐서표기법으로는 식 (6)과 같다.

(6) 

편차유한변형률의 이차불변량(IIE')은 왜곡도를 의미하며 식 (7)에 따라 계산된다.

(7) 
2.3 유한차분법 및 유한요소법을 이용한 변형률 계산

본 연구에서는 파이썬 프로그래밍 언어를 이용하여 불규칙하게 배치된 이산요소모형 내 입자들의 변위 정보를 격자 내 절점에 전사한 후, 공간에 대한 미분을 수행하여 유한변형률(E)을 계산하는 코드를 작성하였다. 해당 과정에 대한 연속체 수치모사기법의 활용가능성을 파악하기 위해, 본 연구에서 사용한 알고리즘은 격자 생성 방식에 따라 유한차분법과 유한요소법으로 나뉜다.

유한차분법을 사용한 계산법에서는 이산요소모형의 불균일한 입자 배치가 격자로 투영된다. 격자는 차분이 용이하도록 일정한 간격으로 연결된 절점들로 구성된(그림 3a). 라그랑지안 관점의 유한변형률이론은 초기 위치를 기준좌표계로 사용하기 때문에, 변형 전 상태의 입자 변위가 절점으로 전사된다. 이 과정에서 발생하는 입자 값과 절점 값 사이의 오차를 줄이기 위해, 역거리 가중법을 사용하였다. 대상 절점과 주변 입자 간 거리에 따라 가중치가 산정된 변위를 절점에 내삽함으로써 절점에 가까운 입자일수록 그 변위가 더 크게 반영되게 하였다(그림 3a). 변위의 근사치가 부여된 임의의 한 절점을 수평 및 수직방향으로 중앙차분하여 변형구배텐서를 구하였다. 이어 계산한 유한변형률 및 물리적 해석을 위한 정량적 측도(체적변형률, 왜곡도)도 동일한 절점에 정의된다. 격자는 초기 입자의 위치를 토대로 생성되었기 때문에 절점 상에 놓인 값들을 그대로 도시할 경우 이산요소모형의 변형 전 구조에 나타나게 된다. 그러나 본 연구는 지질 구조 내 변형률 분포 가시화를 목적으로 하므로, 절점 값을 입자로 다시 전사하는 방법을 통해 변형된 모형 상에 변형률 및 기타 계산 값을 도시하였다.


Fig. 3. 
Mesh grids over particles in the discrete element method (DEM) used to obtain spatial differentials of the displacement field. Open and pink circles indicate DEM particles and nodes. At each node, strains are calculated. (a) In the finite difference method (FDM), particle displacements are mapped to the nodes of the uniform rectangular grid. According to the inverse distance weighting method, particle displacements around a node (e.g., a, b, c, and d) are interpolated regarding their distances from the node (r1, r2, r3, and r4). (b) In the finite element method (FEM), the centers of particles are overlapped with nodes. Three nodes connect to form a triangle-shaped element.

유한요소법을 이용한 변형률 계산법에서는 들로네 삼각화법(Delaunay triangulation method)를 통해 격자를 구축한다(그림 3b). 이는 인접한 입자 간 중심점을 서로 이어 삼각 요소를 구성하는 격자생성법으로 물질의 표면 기하를 명확하게 나타낼 수 있기 때문에 복잡한 형태를 가진 물질 계산에 빈번하게 사용된다(Wu and Lee, 1997; Fedorko et al., 2014). 물질경계면에서 격자 해상도가 높을수록 정확히 계산하는 유한차분법과 달리, 유한요소법에서는 입자 분포를 따라서 생성된 삼각 요소 내에서 원활히 계산을 수행한다. 임의의 요소 내 변위장은 삼각형 요소를 구성하는 3개의 절점 간의 형상함수로 내삽하였다. 요소에 형상함수의 적분값을 부여하는 경우 요소 전체에 걸쳐 평균 낸 결과값을 얻을 수 있지만, 이는 변형 전 입자 위치를 따라서 생성된 요소이므로 변형 후 모형에서 변형률 분포를 알기 어려운 문제점이 있다. 따라서 조성된 변위함수를 통해 유한변형률, 체적변형률, 왜곡도를 구하고 해당 요소 내 3개 절점에 지정하여 변형 후 입자 상에 나타내었다.

3. 결과 및 토의
3.1 후처리를 위한 임의의 연속체 모형의 개발 및 검증

유한차분법 및 유한요소법을 사용한 변형률계산코드를 이산요소모형에 적용하기 전, 본 연구에서 개발한 코드를 이용해 도출한 값과 해석해 간의 비교를 통해 정확성을 검증하는 단계가 필요하다. 그러나 습곡-충상단층대 내의 변형은 불균일하며 비선형성이 커서 해석적인 방법으로 변형률을 계산하는 것은 한계가 있다. 따라서 본 연구는 해석해가 존재하는 2차원 비선형 변형 모형 문제를 정의하고 유한요소법 및 유한차분법 코드를 사용하여 구한 변형률의 수치해를 해석해와 비교하였다. 모형 문제의 초기 상태는 X축과 Y축의 범위가 [0, 1]으로 정사각형을 이루며 초기 좌표가 (X, Y)일 때 변형 후 좌표 (x, y)는 식 (8)에 따라 정의된다.

(8) 

유한차분법 및 유한요소법을 사용하여 모형 문제의 유한변형률텐서, 체적변형률, 왜곡도를 계산하고 해석해와 비교하였다(그림 4). 유한변형률텐서의 대각성분인 E11과 E22은 각각 X방향과 Y방향의 변형률을 의미하며 이에 대한 해석해에서는 각 축과 평행한 방향으로 일정한 변형률을 보인다(그림 4a, 4b). 대칭행렬인 E의 비대각성분 E12과 E21는 서로 같으며 전단방향의 변형률을 나타낸다(그림 4c). 그림 4d와 4e는 각각 체적변형률과 왜곡도의 해석해이다. 유한차분법(그림 4f ~ 4j) 및 유한요소법(그림 4k ~ 4o)을 사용하여 계산한 변형률의 전반적인 분포는 해석해와 일치한다. 또한 해석해와 수치해의 범위를 비교해보면 E11, E22의 경우 0 ~ 20 사이 E21, E12의 경우 -4 ~ 4 사이, J의 경우 -40 ~ 40 사이, IIE'는 -100 ~ 0 사이에 동일하게 존재한다. 그러나 바깥 경계에 해당하는 부분 중 변형률이 높은곳에서는 두 방법의 특징에 따라 해석해와 다른 양상을 보였다. 보간을 사용하는 유한차분법을 통해 E11(그림 4f), E22(그림 4g), IIE'(그림 4j)를 계산할 경우 변형 전 형상에서 바깥 경계에 해당하는 부분(검은 점선 영역) 동일한 위치의 유한요소법을 사용한 결과(그림 4k, 4l, 4o)보다 평활화 되었다. 그에 반해 보간을 하지 않고 절점값 그대로를 사용하는 유한요소법으로 계산했을 경우 해당 위치의 변형률을 뚜렷하게 나타내지만 해석해 수준의 균일한 분포는 보이지 못했다. 두 방법의 정량적인 정확도를 평가하기 위해 해석해에 대한 수치해의 제곱평균제곱근오차(Root Mean Square error; 이하 RMS 오차)를 표 1에 제시하였다. IIE'를 제외한 모든 RMS 오차의 차수는 10-3로 매우 낮은 편에 속했으나 유한차분법의 평균 RMS 오차가 유한요소법보다 약 4배 더 높았다. 그림 5에서는 유한요소법, 유한차분법으로 계산한 E11, E22, IIE'에 대한 RMS 오차를 변형 전 좌표에 나타내고 그 분포를 해석해와 비교하며 분석하였다. 유한차분법을 사용했을 때 3개 지표에 대한 RMS 오차 분포의 공통점은 주로 경계에 집중 되어있다는 것이다(그림 5d, 5e, 5f). 이는 유한차분법의 격자에서 가장 마지막 행/열에 포함된 절점은 중앙차분을 수행하기에 한계가 있어 인접한 절점의 값으로 대신하기 때문으로 추정된다. 한편 유한요소법은 경계에 놓여있는 입자를 절점으로 삼아 그대로 계산한다. 이에 따라 E11(그림 5h), E22(그림 5i)에 대한 유한요소법의 오차 발생률은 각각 좌우 및 상하 경계에서 거의 0에 가깝다. 그러나 유한요소법을 이용한 IIE' RMS 오차 분포(그림 5g)와 해석해(그림 5a)를 비교해보면 변형 구배가 높은 곳에서 RMS 오차가 발생함을 확인하였다. 같은 지표에 대해서 유한차분법을 사용했을 때는 변형 구배가 높은 곳에서 오히려 균일하게 낮은 RMS 오차를 보이며, 이는 보간 과정되었다. 위와 마찬가지로 E11와 E22의 해석해(그림 5b, 5c)의 X축 및 Y축을 따라 점진적으로 변형률이 변화하는 곳에서 유한차분법을 이용했을 경우(그림 5e, 5f)는 평활화 된 값이 나타나며, 유한요소법의 경우(그림 5h, 5i)는 미량의 오차를 보였다. 각각의 경우에서 평균 RMS 오차는 유한차분법이 유한요소법보다 최소 2.5배에서 약 7배가 더 크며 유한차분법에서 경계의 변형률을 정확히 계산하지 못하는 것이 원인으로 고려된다. 예를 들어 변형 후 좌표에 유한차분법으로 구한 변형률을 도시한 그림 4f, 4g, 4j에서 경계선 중 변형률이 높은 부분을 보면 해석해(그림 4a, 4b, 4e)와 다소 상이한 양상을 보인다.


Fig. 4. 
Exact and numerical solutions of the finite strain tensor (E), volumetric strain (, and distortion (IIE'). The diagonal parts of E are E11 and E22. As E is symmetrical, the off-diagonal parts (E12 and E21) are identical. Volumetric strain and distortion are calculated for physical interpretations such as volume change and deviatoric deformation; (a) ~ (e) show exact solutions. (f) ~ (j) and (k) ~ (o) exhibit the corresponding FDM and FEM results, respectively. The black dashed curves in the FDM solutions refer to the region with large error compared to the exact and FEM solutions.

Table 1. 
Root mean square errors (RMSEs) of numerically calculated strains relative to exact solutions.



Fig. 5. 
Root-mean-square error (RMSE) distributions of E11, E22, and IIE' calculated using the FDM and FEM and compared to exact solutions. Panels (a) ~ (c) show exact solutions of strains, whereas (d) ~ (f) and (g) ~ (i) show RMSE for the FDM and FEM, respectively. Most errors of FDM results are on the boundary which showed high strains in the exact solutions. On the other hand, RMSEs of FEM results are concentrated on regions where the gradients of strain in the exact solution are large. The overall RMSEs of FDM results are about four times higher than those of FEM results in averages.

3.2 이산요소법을 이용한 습곡-충상단층대 모형 설정

다수의 지질구조모형 모사 연구에서는 습곡-충상단층대를 압축 하의 토조(sandbox)로 가정한다(Davis et al., 1983; Burbidge and Braun, 2002; Schreurs et al., 2006; Granado et al., 2017). 대표적인 압축 환경인 조산대와 섭입대에서는 각자 다른 압축 기작이 적용된다. 대륙판 사이의 충돌로 인한 조산대의 경우 지각의 연대가 높아 강한 강도를 지니며, 후방지에 해당하는 상대적으로 강한 지각이 전면지에 해당하는 약한 지각을 수평방향으로 수축시킨다. 이 때 전면지 하부층이 상부층과 함께 변형될 경우 단층의 경사각이 증가하며 두꺼운 습곡-충상단층대가 형성된다. 이에 반해 섭입대 환경에서는 밀도가 높은 해양판이 대륙판 하부로 섭입하면서 분리 단층이 활성화되고 상부층에만 변형이 집중된다. 따라서 습곡-충상단층대의 두께가 조산대 환경보다 상대적으로 얇고 평균적인 경사각이 낮다. 모의 실험에서 두 발생 기작에 따른 충상-습곡단층대의 진화를 파악하기 위해서는 모형에 각 경우의 작용 힘 방향에 맞는 경계조건을 부여해야한다. 예를 들어 조산대 내 습곡-충상단층대 모사는 토조의 좌우 벽면 중 한쪽을 후방지로 가정하여 변위를 가하는 방식으로 수행될 수 있다. 반면 섭입대의 경우 토조의 바닥면을 수평방향으로 이동시켜 모사할 수 있다. Davis et al. (1983)에 따르면 이 방법은 벽면에 변위를 가하는 경우보다 평균 10°이내로 낮은 경사각을 가진 단층을 유발한다. Schreurs et al. (2006) 와 Buiter et al. (2006) 는 각각 상사 모사 기법 및 수치 모사 기법을 활용하여 토조 압축 실험을 수행한 여러 연구들을 종합하고 결과를 비교하였다. 본 연구는 위 실험의 벤치마크를 수행하기 위해 이산요소법을 이용하여 토조의 한쪽 벽면에 1.07 × 10-6 m/time-step의 압축률을 가하면서 입자층의 형태 변화를 조사하였다. Schreurs et al. (2006) 및 Buiter et al. (2006) 의 토조는 실내 시험이 가능한 크기로 설정하였다. 이와 동일한 실험 조건을 갖추기 위하여, 그림 6과 같이 최대 높이 5 cm, 길이 35 cm의 2차원 토조 모형을 구축하였다. Schreurs et al. (2006)의 상사 모사 실험에서 토조의 주요 구성 물질은 쿨롱 파괴 법칙을 따르는 마찰성 규사이다. 분리단층 역할을 수행하는 0.5 cm 폭의 마이크로비드(microbead)층하나를 제외한 나머지 부분은 규사층으로서 변형의 가시화를 위해 0.5 cm 간격으로 색을 달리하였다. 또한 이미 존재하던 조산대 및 퇴적층을 나타내는 모래층이 오른쪽 최상부에 쐐기 모양으로 존재한다. 습곡-충상단층대의 진화 기작을 다루는 많은 선행 연구에서는 분리단층이 습곡-충상단층대의 발달을 촉진한다고 제시해왔다(Poblet and Lisle, 2011; Forte et al., 2013; Muñoz et al., 2013). 이는 분리단층이 주변부에 비했을 때 기계적인 약대로써 작용하기 때문이다. 수치모델링기법을 사용한 다수의 연구에서는 점성도, 마찰각, 밀도 등 해당 모델링 기법의 특정 변수를 낮게 조절함으로써 약대를 구현했다(Costa and Vendeville, 2002; Koyi and Vendeville, 2003; Li and Mitra, 2017). 따라서 본 연구는 Itasca사에서 개발한 이산요소법 상용 소프트웨어인 PFC2D (Particle Flow Code 2D)를 사용하였으며 각 층에 해당하는 물질의 특성(예, 밀도, 마찰각, 점착력 등)에 맞는 이산요소법의 입력 변수(예, 입자의 밀도, 접촉부 마찰각, 입자의 직경, 접촉 스프링의 강성, 점착력 등)를 적용하였다(표 2).


Fig. 6. 
Model setup for materials and boundary conditions of a two-dimensional (2D) sandbox for simulating fold-and-thrust belts evolutions using the DEM, which has conditions similar to the analog model developed by Schreurs et al. (2006). The dip angle is 0°, and the initial surface angle of the upper sand pile is 10°. The microbeads layer, which is embedded between layers of quartz sand, represents décollement which promotes significant deformation due to its distinct properties (e.g., density, frictional angle, and Young's modulus). The compressional rate applied on the right sidewall is 1.07 × 10-6 × m·time-step-1, which is sufficiently slow for maintaining a quasi-static state. The black box shows the distributions of contact bonds (see red bars) between adjacent particles.

Table 2. 
Discrete element model of fold-and-thrust belts model parameters.


3.3 습곡-충상단층대 모형 벤치마크 결과

그림 7은 압축에 따른 토조 내 모래층의 변화(그림 7a ~ 7e)와 입자 간 접촉 상태의 분포(그림 7f ~ 7j)를 보여준다. 토조의 초기 횡방향 길이 35 cm에서 25 cm가 될 때까지 지속적으로 1.07 × 10-6 m·time-step-1의 압축률을 가하였으며, 이에 따라 응력을 받은 입자 간 스프링이 강성 한계를 초과해 절단되어 인접한 입자와 분리되는 과정을 관찰하였다. 압축 변위 2 cm일 때(그림 7a) 이동 벽으로부터 약 8 cm 거리의 지점에서 전파단층(forward thrust)과 역방향단층(backthrust)이 동시에 발생하여 돌출구조가 생성되었다. 이 구조는 압축 변위 4 cm (그림 7b)까지 발달하였으며 압축 변위 6 cm (그림 7c)부터는 약 9 cm 거리의 지점에서 두번째 돌출구조가 발생하였지만 그 거리가 기존 돌출 구조와 매우 가까워 새로 발생한 역방향단층이 기존 역방향단층과 중첩되었다. 압축 변위 8 cm (그림 7d)에서는 3번째 돌출구조가 생성되었으며 앞선 구조와 약 8 cm 떨어진 지점에서 나타났기 때문에 역방향단층 및 돌출구조가 뚜렷하게 관찰되었다. 압축 변위 10 cm (그림 7e)에서는 약 4 cm 거리에서 4번째 돌출구조가 발생하였으며, 최종적으로 4개의 전파단층이 연이어 존재하는 인편상단층(imbricate fault)을 이루었다. 접촉 상태의 스프링은 그림 7f ~ 7j에서 붉은색으로, 파괴된 스프링은 푸른색으로 나타난다. 본 연구에서 단층은 가시적인 규모에서 파쇄된 무결암을 의미하며 이산요소모형을 이루는 입자 간 비 접촉 상태와 같다(그림 7f ~ 7j의 청색 영역). 따라서 그림 7a ~ 7e에서 관찰 가능한 단층 형태는 그림 7f ~ 7j에서 더욱 명확하게 확인할 수 있다. 또한 그림 7c ~ 7e에서 정확히 보이지 않았던 역방향 단층은 그림 7h ~ 7j에서 뚜렷하게 드러난다.


Fig. 7. 
Results of shortening experiment conducted using DEM. (a) ~ (e) Time evolution steps of the fold-and-thrust belts model, with increasing total displacement from 2 cm to 10 cm. A pop-up structure is shown in (a); (b) ~ (e) after 4 cm convergence, the successive forward thrusts develop in sequence, forming an imbricate thrust system. (f) ~ (j) States of the contact map at time steps equivalent to (a) ~ (e). The red and blue regions denote the areas where particles are in contact and separated, respectively. The area pervasive with separated particles represents faults, backthrusts comprising paired pop-up structures with forward thrusts clearly appear.

Buiter et al. (2006)은 7개의 수치 해석 코드를 사용하여 토조 압축 모형을 설정하고 압축에 따른 단층 개수, 표면 경사각, 단층 간격 등을 측정 후 Schreurs et al. (2006)의 상사 모형 결과와 비교하였다. 본 연구는 제작한 모형의 유효성을 검토하기 위해 이산요소모형의 표면 경사각, 단층 간격, 단층 개수를 총 압축 변위에 대해 도시하고 Buiter et al. (2006)가 정리한 수치 및 상사 모델링 결과와 대조하였다(그림 8). 본 연구에서 사용한 길이 30 cm 토조 모형의 총 압축 변위는 10 cm이며, 전체 압축량이 2, 4, 6, 8, 10 cm인 시점에서 표면 경사각, 단층 간격, 단층 개수를 측정하였다. 수치 및 상사 모델링을 통한 기존 연구 결과가 속한 영역은 채색하여 도시하였으며, 본 연구 결과는 붉은 점으로 표시하였다. 표면 경사각은 단층 상부 외곽선이 기울어진 각도로서 압축 변위가 2 cm 일 때 약 15°, 6 cm 일 때 21°, 10 cm 일 때 18°로 소폭 증가하였다가 다시 감소하는 경향을 보였다(그림 8a). 이는 기존에 수치 및 상사 모델링으로 수행한 결과 범위 내에 해당하였다. 압축이 진행됨에 따라 발생하는 단층 간 간격(그림 8b) 및 개수(그림 8c) 또한 기존 결과 내에 분포하였다. 단층 간격의 경우, 전파단층개수가 2개인 6 cm부터 측정하였으며 총 압축 변위인 10 cm에 도달할 때까지 단층의 간격이 평균 3 cm로 그 편차가 미미하였다. 본 토조모형을 압축할수록 전파단층과 역방향단층이 교대로 발생하였다. Buiter et al. (2006)의 단층 개수 측정 방법에 따라 정방향 단층만 고려하였으며 그림 8c에서 최대 4개까지 발생하였다. 이산요소법의 입자 간 강성은 입자 간 중심 거리에 반비례한 관계를 갖는다. 이는 입자의 반지름이 강성에 영향을 미치며, 입자의 크기가 달라질 시에는 입자를 새로 배치시키는 과정 및 강성 차이로 미시적인 접촉 힘 분포가 변화한다. 따라서 합리적인 크기의 입자를 사용하며 계산의 반복 수행을 통해 크기에 따른 민감도를 파악할 필요가 있다. 본 연구에서는 수행한 토조 압축 실험에 사용된 입자의 개수는 총 17,154개이고, 분리단층역할을 수행하는 마이크로비드층을 제외한 경우 입자의 크기는 3·10-4 ~ 6·10-4 m이다. Schreurs et al. (2006)에서 발췌한 Parma 대학에서 수행한 토조 압축 실험에서 사용된 입자의 크기는 6·10-5 ~ 2.5·10-4 m 이다. 따라서 본 연구는 실제 모래 수준과 비슷한 크기의 입자를 사용하였다. 이를 바탕으로 계산을 반복 수행하고 수렴된 결과를 사용하여 입자 크기에 따른 결과 변동성을 감소시켰다.


Fig. 8. 
Analysis of the DEM results compared to previous analog (Schreurs et al., 2006) and numerical (Buiter et al., 2006) studies. Geometric results (i.e., surface angle, spacing between two near-thrusts, and the number of thrusts) are measured at the time-steps with the same displacements from 2 to 10 cm (red circles). (a) Blue and purple bands denote the surface angles for the analog and numerical models, respectively. The top inset schematically describes the measurement of a surface slope. (b) Light green and green bands indicate the spacing between two near-thrusts in the analog and numerical models, respectively. The center inset represents the determination of the space. (c) The numbers of thrusts for the analog and numerical experiments are shown in the yellow and red bands, respectively.

3.4 습곡-충상단층대 모형의 변형률

이산요소법으로 수행한 습곡-충상단층대모형 결과의 유한변형률텐서, 체적변형률, 왜곡도를 계산하였다. 그림 9a ~ 9d에서 유한변형률텐서의 대각성분인 E11과 E22는 단층부에서 높은 변형률 값을 보여준다. 토조가 X축 방향으로 압축되었기 때문에 Y축방향으로 인장력이 작용하여 해당 방향의 변형률을 보여주는 E22가 X축 방향의 변형률인 E11 보다 평균적으로 높은 값을 보였다. E는 대칭행렬로서 E21과 E12가 서로 같기 때문에 E12만 도시하였다(그림 9e, 9f). E12는 절대값은 같지만 부호가 다른 결과 범위를 가지며 부호는 단층의 운동감각을 지시한다. 오른쪽 방향으로 경사진 단층은 양수이며 왼쪽 방향은 음수로 E12 값이 표시된다. 그림 9g, 9h에서 체적변형률(J)의 계산 범위는 E12와 같지만 방향에 따른 특징이 보이지 않고 양수 값과 음수 값 모두 단층에 혼재한다. 그림 9i, 9j의 왜곡도(IIE')는 모두 음수로 측정되었으며 E11과 E22과 반대로 단층에서 더 작은 값을 보였다. 유한요소법을 사용했을 때(그림 9b, 9d, 9f, 9h, 9j)가 유한차분법을 사용했을 때(그림 9a, 9c, 9e, 9g, 9i)보다 상부 경계의 변형률을 명확히 나타내지만 전체적으로는 유한차분법을 사용했을 때보다 불균일한 분포를 보인다. 이는 모형 문제의 변형률 계산 결과와 동일하게 보간 과정 및 경계 절점 계산의 유무를 원인으로 해석할 수 있다.


Fig. 9. 
Comparison between the strain distribution calculated using the FDM and FEM. Diagonal parts of the finite strain tensor are shown in (a) ~ (d). Panels (a) and (c) show the FDM results, whereas (b) and (d) show the FEM results. Panels (e) and (f) indicate the non-diagonal parts attained from the FDM and FEM, respectively. (g) represents volumetric strain of the FDM and (h) display that of the FEM. The amount of distortion, the second invariants of the deviatoric finite strain tensor, was also calculated by the FDM and FEM, which are displayed in (i) and (j), respectively.

4. 결 론

본 연구는 물질의 취성-소성 변형 모사에 적합한 이산요소법을 통해 습곡-충상단층대의 진화를 모사하였으며, 유한요소법 및 유한차분법을 도입하여 퇴적 입자 간 결합과 파쇄가 반영된 미소 단위 변형률을 계산하였다. 단층 연구에서 변형률은 일차적으로 단층의 변위 규모와 운동 감각을 보여주며, 이차적으로는 전단열 증대 현상(Jaquet et al., 2016), 마찰력 약화(Rathbun and Marone, 2010), 에너지 소산(Lieou et al., 2014) 등과 깊은 관련이 있다. 본 연구에서 계산한 입자 단위 변형률은 습곡-충상단층대 진화의 운동학적 특성을 반영하였으며 추후 위 현상들과 정량적으로 관계 지을 수 있는 가능성을 제시하였다. 지질 구조의 변형은 매우 오랜 시간에 걸쳐 미소하게 변형하여 해당 지대의 변형률과 응력 변화 이력을 추적하는 일은 매우 어렵다. 본 연구에서는 습곡-충상단층대 모형의 시간에 따른 변형률을 계산하였으며 교결, 압밀, 침식 등 퇴적물의 마찰성 입상 거동이 반영되었다는 점에서 다양한 활용가능성을 보여준다. 모형 모사의 결과로 획득한 정량적 자료들(예, 변위, 좌표 등)을 후처리 하여 모사 결과에 대한 정성적 의의(예, 운동감각, 체적 변형률, 왜곡도, 변형 규모 등)를 얻었다. 이에 불연속 물질을 기반한 이산요소법을 사용한 모형임을 고려하였을 때 변형률 계산을 위한 가장 효율적인 후처리 기법을 비교, 분석하였다.

본 연구에서는 습곡-충상단층대의 진화 및 변형률 분포를 파악하기 위해 이산요소법을 이용하여 습곡-충상단층대 모형을 제작하고 유한차분법 및 유한요소법을 통해 변형률을 계산하였다. Buiter et al. (2006)의 토조 압축 실험을 벤치마크 하였으며 압축률에 따른 세가지 지표(즉, 표면 경사각, 단층 간격, 단층 개수)의 비교를 통해 본 연구에서 개발한 습곡-충상단층대모형의 유효성을 검증하였다. 이산요소모형의 입자의 배열은 매우 불규칙하며 입자 간 간격이 일정하지 않아 추가적인 기법 없이 변형률을 구하기에 한계가 있다. 본 연구는 이산요소모형의 변형률 분석 방법으로써 유한차분법 및 유한요소법을 도입하여 변형률을 계산하였다. 두 계산법의 정확성 및 효율성을 정량적으로 평가하기 위하여 모형 문제의 변형률을 구하고 해석해와 비교하였다. 두 방법 모두 해석해와 전반적인 변형률 분포가 일치하였으나 변형률이 급격히 변화하는 지점에서는 두 결과에 다소 차이가 발생했다. 유한차분법에서는 격자 구성을 위해 변위 정보를 입자에서 절점으로 전사하는 과정에서 일부 변위가 손실되는 대신 평활화가 이루어지기 때문에 해석해보다 변형률 구배가 작게 나타났다. 유한요소법에서는 변위의 보간 과정이 불필요하며 입자 값이 곧 절점 값이 되어 자료의 손실을 피했지만 해석해에 완전히 상응하는 균일함은 구현하지 못했다. 그러나 두 방법 모두 RMS 오차의 단위가 10-3 ~ 10-2로 양질의 계산을 수행했다고 해석된다. 유한차분법의 RMS 오차는 유한요소법보다 약 4배가 높으며 특히 오차가 경계에 집중되어 있다. 이는 경계선에 존재하는 마지막 행/열은 중앙차분이 불가능하기에 가장 인접한 값을 대신 사용했기 때문이다. 유한차분법을 사용한 변형률 계산은 균일한 결과를 내어 결과 해석을 원활하게 할 수 있지만 보간 방법과 격자 해상도에 따라서 정확도 차이가 크게 발생하므로 신중한 처리가 요구된다. 특히 입자 개수 대비 격자 해상도가 매우 높을 경우 오차 발생률이 급격히 증가하였다. 이산요소법을 사용한 습곡-충상단층대모형은 입자 단위의 취성 변형을 본질적으로 계산하며, 특히 단층 내 암석 파쇄로 인한 불규칙한 경계를 사실적으로 나타낸다는 장점이 있다. 경계부는 침전, 침식, 교란 등과 같은 퇴적작용이 활발하게 일어나 변형의 가능성이 높은 부분이다. 유한요소법은 입자의 위치를 그대로 절점으로 사용한다는 점에 있어서 습곡-충상단층대의 불균일한 경계 내 변형률을 정확히 계산할 수 있다. 이는 유한요소법이 이산요소법을 사용한 모형의 후처리 수단으로써 높은 효율성을 가짐을 제시한다.

Acknowledgments

본 연구는 “행정안전부 방재안전분야 전문인력양성”사업과 기상청 기상·지진 See-At 기술개발사업(KMI2018-02110), 한국연구재단 신진연구지원 사업(NRF-2019R1C1C1010804)과 중점연구소지원 사업(No.2019R1A6A1A03033167)으로 수행되었습니다.

References
1. Bonini, M., 2007, Deformation patterns and structural vergence in brittle-ductile thrust wedges: An additional analogue modelling perspective. Journal of Structural Geology, 29, 141-158.
2. Borderie, S., Vendeville, B.C., Graveleau, F., Witt, C., Dubois, P., Baby, P. and Calderon, Y., 2019, Analogue modeling of large-transport thrust faults in evaporites-floored basins: Example of the Chazuta Thrust in the Huallaga Basin, Peru. Journal of Structural Geology, 123, 1-17.
3. Branquet, Y., Cheilletz, A., Giuliani, G., Laumonier, B. and Blanco, O., 1999, Fluidized hydrothermal breccia in dilatant faults during thrusting: the Columbian emerald deposits. Geological Society Special Publication, 155, 183-195.
4. Buiter, S.J.H., Babeyko, A.Y.U., Ellis, S., Gerya, T.V., Kaus, B.J.P., Kellner, A., Schreurs, G. and Yamada, Y., 2006, The numerical sandbox: comparison of model results for a shortening and an extension experiment. Geological Society, London, Special Publications, 253, 29-64.
5. Buiter, S.J.H., Schreurs, G., Albertz, M., Gerya, T.V., Kaus, B., Landry, W., le Pourhiet, L., Mishin, Y., Egholm, D.L., Cooke, M., Maillot, B., Thieulot, C., Crook, T., May, D., Souloumiac, P. and Beaumont, C., 2016, Benchmarking numerical models of brittle thrust wedges. Journal of Structural Geology, 92, 140-177.
6. Burbidge, D.R. and Braun, J., 2002, Numerical models of the evolution of accretionary wedges and fold-and-thrust belts using the distinct-element method. Geophysical Journal International, 148, 542-561.
7. Burmeister, K.C., Harrison, M.J., Marshak, S., Ferré, E.C., Bannister, R.A. and Kodama, K.P., 2009, Comparison of Fry strain ellipse and AMS ellipsoid trends to tectonic fabric trends in very low-strain sandstone of the Appalachian fold-thrust belt. Journal of Structural Geology, 31, 1028-1038.
8. Byerlee, J., 1978, Friction of rocks. In Rock friction and earthquake prediction. Birkhäuser, Basel., 615-626.
9. Chapple, W.M., 1978, Mechanics of thin-skinned fold-and-thrust belts. Geological Society of America Bulletin, 89, 1189-1198.
10. Chiu, C.C., Weng, M.C. and Huang, T.H., 2015, Biconcave bond model for cemented granular material. Journal of GeoEngineering, 10, 91-103.
11. Cosgrove, J.W. and Ameen, M.S., 2000, A comparison of the geometry, spatial organization and fracture patterns associated with forced folds and buckle folds. Geological Society, London, Special Publications, 169, 7-21
12. Costa, E. and Vendeville, B.C., 2002, Experimental insights on the geometry and kinematics of fold-and-thrust belts above weak, viscous evaporitic décollement. Journal of Structural Geology, 24, 1729-1739.
13. Cundall, P.A., 1971, A computer model for simulating progressive large scale movements in blocky rock systems. In: Proceedings of the Symposium of International Society of Rock Mechanics, vol. 1, Nancy: France, Paper No. II-8
14. Dahlen, F.A., 1984, Noncohesive critical Coulomb wedges: An exact solution. Journal of Geophysical Research: Solid Earth, 89(B12), 10125-10133.
15. Davis, D., Suppe, J. and Dahlen, F.A., 1983, Mechanics of fold-and- thrust belts and accretionary wedges. Journal of Geophysical Research, 88(B2), 1153-1172.
16. De Paor, D.G., 1990, Determination of the strain ellipsoid from sectional data. Journal of Structural Geology, 12, 131-137.
17. Dean, S.L., Morgan, J.K. and Fournier, T., 2013, Geometries of frontal fold and thrust belts: Insights from discrete element simulations. Journal of Structural Geology, 53, 43-53.
18. Depta, P.N., Jandt, U., Dosta, M., Zeng, A.P. and Heinrich, S., 2018, Toward Multiscale Modeling of Proteins and Bioagglomerates: An Orientation-Sensitive Diffusion Model for the Integration of Molecular Dynamics and the Discrete Element Method. Journal of chemical information and modeling, 59, 386-398.
19. Dominguez, S., Malavieille, J. and Lallemand, S.E., 2000, Deformation of accretionary wedges in response to seamount subduction: Insights from sandbox experiments. Tectonics, 19, 182-196.
20. Dorsey, M.T., Rockwell, T.K., Girty, G.H., Ostermeijer, G.A., Browning, J., Mitchell, T.M. and Fletcher, J.M., 2021, Evidence of hydrothermal fluid circulation driving elemental mass redistribution in an active fault zone. Journal of Structural Geology, 144, 104269.
21. Egholm, D.L., Sandiford, M., Clausen, O.R. and Nielsen, S.B., 2007, A new strategy for discrete element numerical models: 2. Sandbox applications. Journal of Geophysical Research: Solid Earth, 112, 1-12.
22. Espinosa, H.D., Zavattieri, P.D. and Emore, G.L., 1998, Adaptive FEM computation of geometric and material nonlinearities with application to brittle failure. Mechanics of Materials, 29, 275-305.
23. Fedorko, G., Stanova, E., Molnar, V., Husakova, N. and Kmet, S., 2014, Computer modelling and finite element analysis of spiral triangular strands. Advances in Engineering Software, 73, 11-21.
24. Feng, Y., Kowalsky, M.J. and Nau, J.M., 2015, Finite-element method to predict reinforcing bar buckling in RC structures. Journal of Structural Engineering, 141, 04014147.
25. Fleissner, F., Gaugele, T. and Eberhard, P., 2007, Applications of the discrete element method in mechanical engineering. Multibody system dynamics, 18, 81-94.
26. Forte, A.M., Cowgill, E., Murtuzayev, I., Kangarli, T. and Stoica, M., 2013, Structural geometries and magnitude of shortening in the eastern Kura fold-thrust belt, Azerbaijan: Implications for the development of the Greater Caucasus Mountains. Tectonics, 32, 688-717.
27. Furuichi, M., Nishiura, D., Kuwano, O., Bauville, A., Hori, T. and Sakaguchi, H., 2018, Arcuate stress state in accretionary prisms from real-scale numerical sandbox experiments. Scientific Reports, 8, 1-11.
28. Glasmacher, U.A., Matenaar, I., Bauer, W. and Puchkov, V.N., 2004, Diagenesis and incipient metamorphism in the western fold-and-thrust belt, SW Urals, Russia. International Journal of Earth Sciences, 93, 361-383.
29. Goldsby, D.L. and Tullis, T.E., 2011, Flash heating leads to low frictional strength of crustal rocks at earthquake slip rates. Science, 334, 216-218.
30. Granado, P., Ferrer, O., Muñoz, J.A., Thöny, W. and Strauss, P., 2017, Basin inversion in tectonic wedges: Insights from analogue modelling and the Alpine-Carpathian fold-and-thrust belt. Tectonophysics, 703-704, 50-68.
31. Gurnis, M., Zhong, S. and Toth, J., 2000, On the competing roles of fault reactivation and brittle failure in generating plate tectonics from mantle convection. In: The History and Dynamics of Global Plate Motions. Geophysical Monograph, No.121, American Geophysical Union, Washington, DC, 73-94.
32. Hardy, S. and Finch, E., 2005, Discrete‐element modelling of detachment folding. Basin Research, 17, 507-520.
33. Hong, S., Kwon, S., Min, K. and Ji, S., 2021, Effect of Excavation and Thermal Stress on Slip Zone and Aperture Change Around Disposal Hole and Tunnel in Fractured Rock. Tunnel & Underground Space, 31, 125-144 (in Korean with English abstract).
34. Hori, T. and Sakaguchi, H., 2011, Mechanism of décollement formation in subduction zones. Geophysical Journal International, 187, 1089-1100.
35. Itasca Consulting Group, 1999, PFC2D (Particle Flow Code in 2 Dimensions). Itasca Consulting Group, Minneapolis, MN, 1-124.
36. Jaquet, Y., Duretz, T. and Schmalholz, S.M., 2016, Dramatic effect of elasticity on thermal softening and strain localization during lithospheric shortening. Geophysical Journal International, 204, 780-784.
37. Jerier, J.F., Hathong, B., Richefeu, V., Chareyre, B., Imbault, D., Donze, F.V. and Doremus, P., 2011, Study of cold powder compaction by using the discrete element method. Powder Technology, 208, 537-541.
38. Jo, T. and So, B., 2020, Numerical Modeling of Shear Heating in 2D Elastoplastic Extensional Lithosphere using COMSOL Multiphysics®. Geophysics and Geophysical Exploration, 23, 1-12 (in Korean with English abstract).
39. Kaus, B.J., Mühlhaus, H. and May, D.A., 2010, A stabilization algorithm for geodynamic numerical simulations with a free surface. Physics of the Earth and Planetary Interiors, 181, 12-20.
40. Kendall, J., Vergés, J., Koshnaw, R. and Louterbach, M., 2020, Petroleum tectonic comparison of fold and thrust belts: the Zagros of Iraq and Iran, the Pyrenees of Spain, the Sevier of Western USA and the Beni Sub-Andean of Bolivia. Geological Society, London, Special Publications, 490, 79-103.
41. Koyi, H.A. and Vendeville, B.C., 2003, The effect of décollement dip on geometry and kinematics of model accretionary wedges. Journal of Structural Geology, 25, 1445-1450.
42. Lee, H. and So, B., 2020, Numerical simulation of coseismic and postseismic deformation using finite element modeling with weak elastic fault. Journal of the Geological Society of Korea, 56, 771-787 (in Korean with English abstract).
43. Li, C., Yin, H., Wu, C., Zhang, Y., Zhang, J., Wu, Z., Wang, W., Jia, D., Guan, S. and Ren, R., 2021, Calibration of the Discrete Element Method and Modeling of Shortening Experiments. Frontiers in Earth Science, 9, 394.
44. Li, J. and Mitra, S., 2017, Geometry and evolution of fold-thrust structures at the boundaries between frictional and ductile detachments. Marine and Petroleum Geology, 85, 16-34.
45. Li, S. and Ghosh, S., 2006, Extended Voronoi cell finite element model for multiple cohesive crack propagation in brittle materials. International Journal for Numerical Methods in Engineering, 65, 1028-1067.
46. Lieou, C.K.C., Elbanna, A.E. and Carlson, J.M., 2014, Grain fragmentation in sheared granular flow: Weakening effects, energy dissipation, and strain localization. Physical Review E - Statistical, Nonlinear, and Soft Matter Physics, 89, 1-14.
47. Lisle, R.J., 1985, The effect of composition and strain on quartz-fabric intensity in pebbles from a deformed conglomerate. Geologische Rundschau, 74, 657-663.
48. Liu, E., 2005, Effects of fracture aperture and roughness on hydraulic and mechanical properties of rocks: implication of seismic characterization of fractured reservoirs. Journal of Geophysics and Engineering, 2, 38-47.
49. Lohrmann, J., Kukowski, N., Adam, J. and Oncken, O., 2003, The impact of analogue material properties on the geometry, kinematics, and dynamics of convergent sand wedges. Journal of Structural Geology, 25, 1691-1711.
50. Ma, X. and Elbanna, A., 2019, Dynamic rupture propagation on fault planes with explicit representation of short branches. Earth and Planetary Science Letters, 523, 115702.
51. McCaffrey, R., King, R.W., Wells, R.E., Lancaster, M. and Miller, M.M., 2016, Contemporary deformation in the Yakima fold and thrust belt estimated with GPS. Geophysical Journal International, 207, 1-11.
52. McClay, K.R., Whitehouse, P.S., Dooley, T. and Richards, M., 2004, 3D evolution of fold and thrust belts formed by oblique convergence. Marine and Petroleum Geology, 21, 857-877.
53. Meng, Q. and Hodgetts, D., 2019, Combined control of décollement layer thickness and cover rock cohesion on structural styles and evolution of fold belts: A discrete element modelling study. Tectonophysics, 757, 58-67.
54. Metcalf, J.R., Fitzgerald, P.G., Baldwin, S.L. and Muñoz, J.A., 2009, Thermochronology of a convergent orogen: Constraints on the timing of thrust faulting and subsequent exhumation of the Maladeta Pluton in the Central Pyrenean Axial Zone. Earth and Planetary Science Letters, 287, 488-503.
55. Mishra, B.K., 2003, A review of computer simulation of tumbling mills by the discrete element method: part I—contact mechanics. International Journal of Mineral Processing, 71, 73-93.
56. Moore, G.F., Shipley, T.H., Stoffa, P.L., Karig, D.E., Taira, A., Kuramoto, S., Tokuyama, H. and Suyehiro, K., 1990, Structure of the Nankai Trough accretionary zone from multichannel seismic reflection data. Journal of Geophysical Research: Solid Earth, 95, 8753-8765.
57. Morgan, J.K., 2015, Effects of cohesion on the structural and mechanical evolution of fold and thrust belts and contractional wedges: Discrete element simulations. Journal of Geophysical Research: Solid Earth, 120, 3870-3896.
58. Muñoz, J., Beamud, E., Fernández, O., Arbués, P., Dinarès-turell, J. and Poblet, J., 2013, The Ainsa Fold and thrust oblique zone of the central Pyrenees: Kinematics of a curved contractional system from paleomagnetic and structural data. Tectonics, 32, 1142-1175.
59. O'Sullivan, C., 2011, Particulate discrete element modelling: a geomechanics perspective. CRC Press.
60. Park, J.O., Tsuru, T., Takahashi, N., Hori, T., Kodaira, S., Nakanishi, A., Miura, S. and Kaneda, Y., 2002, A deep strong reflector in the Nankai accretionary wedge from multichannel seismic data: Implications for underplating and interseismic shear stress release. Journal of Geophysical Research: Solid Earth, 107, ESE-3.
61. Poblet, J. and Lisle, R.J., 2011, Kinematic evolution and structural styles of fold-and-thrust belts. Geological Society, London, Special Publications, 349, 1-24.
62. Ramsay, J.G., Huber, M.I. and Lisle, R.J., 1983, The techniques of modern structural geology: Folds and fractures (Vol. 2). Academic press.
63. Rathbun, A.P. and Marone, C., 2010, Effect of strain localization on frictional behavior of sheared granular materials. Journal of Geophysical Research, 115, 1-16.
64. Regenauer-Lieb, K. and Yuen, D.A., 1998, Rapid conversion of elastic energy into plastic shear heating during incipient necking of the lithosphere. Geophysical Research Letters, 25, 2737-2740.
65. Reuber, G.S. and Simons, F.J., 2020, Multi-physics adjoint modeling of Earth structure: combining gravimetric, seismic, and geodynamic inversions. GEM-International Journal on Geomathematics, 11, 1-38.
66. Roy, S., Bose, S. and Saha, P., 2020, Spatial variations of ductile strain in fold-and-thrust belts: From model to nature. Journal of Structural Geology, 134, 104012.
67. Ruh, J.B., Gerya, T. and Burg, J.P., 2013, High-resolution 3D numerical modeling of thrust wedges: Influence of décollement strength on transfer zones. Geochemistry, Geophysics, Geosystems, 14, 1131-1155.
68. Ruh, J.B., Vergés, J. and Burg, J.P., 2018, Shale-related minibasins atop a massive olistostrome in an active accretionary wedge setting: Two-dimensional numerical modeling applied to the Iranian Makran. Geology, 46, 791-794.
69. Schreurs, G., Buiter, S.J.H., Boutelier, D., Corti, G., Costa, E., Cruden, A.R., Daniel, J.M., Hoth, S., Koyi, H.A., Kukowski, N., Lohrmann, J., Ravaglia, A., Schlische, R.W., Withjack, M.O., Yamada, Y., Cavozzi, C., Del Ventisette, C., Elder Brady, J.A., Hoffmann-Rothe, A., Mengus, J., Montanari, D. and Nilforoushan, F., 2006, Analogue benchmarks of shortening and extension experiments. Geological Society, London, Special Publications, 253, 1-27.
70. Seyferth, M. and Henk, A., 2004, Syn-convergent exhumation and lateral extrusion in continental collision zones-Insights from three-dimensional numerical models. Tectonophysics, 382, 1-29.
71. Sibson, R.H., 2020, Preparation zones for large crustal earthquakes consequent on fault-valve action. Earth, Planets and Space, 72, 31.
72. So, B.D. and Capitanio, F.A., 2020, Self-consistent stick-slip recurrent behaviour of elastoplastic faults in intraplate environment: a Lagrangian solid mechanics approach. Geophysical Journal International, 221, 151-162.
73. Stipp, M., Rolfs, M., Kitamura, Y., Behrmann, J.H., Schumann, K., Schulte-Kortnack, D. and Feeser, V., 2013, Strong sediments at the deformation front, and weak sediments at the rear of the Nankai accretionary prism, revealed by triaxial deformation experiments. Geochemistry, Geophysics, Geosystems, 14, 4791-4810.
74. van Zelst, I., Crameri, F., Pusok, A.E., Glerum, A., Dannberg, J. and Thieulot, C., 2021, 101 Geodynamic modelling: How to design, carry out, and interpret numerical studies. Solid Earth Discussions, 1-80.
75. Vaseghi, H., Maleki, Z. and Arian, M., 2016, Structural Style in the Zagros Fold-Thrust Belt: The Gavbast Anticline, Coastal Fars. Open Journal of Geology, 6, 109-116.
76. Wu, J.Y. and Lee, R., 1997, The advantages of triangular and tetrahedral edge elements for electromagnetic modeling with the finite-element method. IEEE Transactions on Antennas and Propagation, 45, 1431-1437.
77. Yamada, Y., Baba, K., Miyakawa, A. and Matsuoka, T., 2014, Granular experiments of thrust wedges: Insights relevant to methane hydrate exploration at the Nankai accretionary prism. Marine and Petroleum Geology, 51, 34-48.
78. Zhao, J. and Gao, Z., 2014, Modelling non-coaxiality and strain localisation in sand: The role of fabric and its evolution. In: 14th International Conference of the International Association for Computer Methods and Advances in Geomechanics (14IACMAG) (Abstracts), Kyoto, Japan, 22-25, 6 p.
79. Zienkiewicz, O.C. and Taylor, R.L., 2000, The finite element method: solid mechanics (Vol. 2). Butterworthheinemann, 445 p.