ASPECT를 이용한 열개 수치 모형의 자유 표면 모사 기법 간의 비교 연구

1. 서 론

열개 현상은 판 간의 상대 운동이나 플룸이 상승하는 지역의 상부에서 발생하는 인장력에 의해 지각 두께가 감소하거나 파괴되어 새로운 해양 지각을 생성하고 지구 물질을 순환시키는 지구동역학적 현상이다(Turcotte and Emerman, 1983; Wright et al., 2012; Frizon de lamotte et al., 2015). 열개는 섭입과 더불어 암석권 변형, 맨틀 대류와 플룸 상승이 복합적으로 작용하는 행성의 지구조적 진화 과정으로 판구조론 관점에서 지구조의 과거와 미래를 이해하는데 중요한 의미를 가진다. 판 운동 및 지구조 진화는 직접적인 관찰이 제한된 대규모의 시공간을 가지기 때문에, 이러한 제약을 극복할 수 있는 수치 모사나 상사 모사(Analogue modeling) 등을 사용하여 열개 현상을 재현하는 연구가 수행되어왔다(Corti, 2012; Liao and Gerya, 2015). 대륙 열개구조는 암석권의 대변형으로 인해 지표면에서 침강(van der beek et al., 1994)과 융기(Bohannon et al., 1989), 지구대(Graben; Rosendahl et al., 1986), 화성활동(Franke, 2013) 등의 형태로 나타나며 측지 자료와 같은 정량적 자료를 직접 획득할 수 있는 판 구조 진화 과정의 중요한 증거이다. 열개에 따른 지표변형을 모사를 통해 열개의 진화과정과 암석권 하부의 물성 및 구조를 파악하기 위해, 최근 열개 현상 수치 모사 연구는 실제 자연에서 관찰되는 지표면 변형을 모사하려는 노력을 기울이고 있다(예, Guyonnet-Benaize et al., 2015).

열개 현상으로 인한 지표 변형은 지질학적 시간 단위(수백만 년 이상)에서 유체 거동을 보이므로, 수치 모사를 통해 지표면 변위를 정확히 모사하기 위해서는 자유 표면(Free surface) 경계조건을 도입해야만 한다. 구조적 대변형을 유발하는 열개 수치 모사는 자유표면 조건에 따라 수치 격자의 변형 및 물질의 큰 이동을 모사하기 때문에 유한요소법 바탕의 기존 지구동역학 수치 모형에 적용하기 난해한 것으로 알려져 있다(d'Acremont et al., 2003; May et al., 2014). 이를 극복하기 위해 Sticky-Air (이하, SA) 기법(Nikolaeva et al., 2010; Allken et al., 2011; Crameri et al., 2012)과 Arbitrary Lagrangian-Eulerian (이하, ALE) 기법(Fullsack, 1995; Sarrate et al., 2001; Wenker and Beaumont, 2018)이 개발되어 지구동역학 수치 모형에서 자유 표면을 모사하는 기법으로 주요하게 사용되고 있다(그림 1). ALE 기법은 구조 변형 및 응력 해석에 사용되는 라그랑지안(Lagrangian) 기법과 열-유체 흐름 해석에 사용되는 오일러리안(Eulerian) 기법을 모두 적용하여 격자 자체를 속도장에 따라 변형시켜 모형의 최상부를 자유 표면처럼 모사하는 기법이다(그림 1a). SA 기법과 비교하였을 때 지표면의 침강 및 융기를 실제와 유사하게 모사할 수 있지만 격자 변형으로 인해 강성행렬 재조립에 따른 연산 시간이 증가하고, 격자가 크게 뒤틀리면서 계산 안정성이 저해되는 단점이 있다. SA 기법은 모형의 최상부에 물 또는 대기와 같이 점성도가 낮고 열적 확산이 매우 빠른 가상의 층을 도입하여 하부의 지각층에 자유 표면 효과를 구현하는 기법이다(그림 1b). SA 기법은 비교적 도입이 쉽고 격자 변형이 없어 연산이 안정적이라는 장점이 있다. 그러나, SA층에 추가적인 연산 할당이 필요할 뿐만 아니라 수치 계산상 점성도를 낮추는 데 한계가 있어 SA층의 매우 낮은 점성도 조건을 완벽히 충족시킬 수 없다. 지각층과 물층/대기층 사이의 경계면에서 점성도는 10²⁸~10³⁰배 정도로 급격히 변하므로 실제 자연의 경계면을 수치 모형에 그대로 반영하게 되면 경계면 근처에 존재하는 격자의 강성행렬(Stiffness matrix) 요소 크기가 크게 차이 나서 역행렬을 구하는 수치 연산이 어려워진다.

Fig. 1.

Schematic diagrams for free-surface deformation using (a) ALE and (b) SA methods. The top and bottom rows show the models before and after deformation by a horizontal extension (grey arrows). The black lines and points indicate the mesh grid and nodes, respectively. The red region refers to the upper crust. The additional SA layers (blue region) are defined above the upper crust when using the SA method. The interface between the blue and red regions is assumed as the free surface of the SA method. (c) The meshes are deformed following the velocity field driven by mass flux. (d) Internal compositions are advected satisfying mass conservation.

지난 수십 년 간 열개 현상을 연구하기 위한 다수의 수치 모사 소프트웨어, 예를 들어 CITCOM (Moresi and Parsons, 1995), CONMAN (King et al., 1990), UNDERWORLD (Moresi et al., 2007)가 개발되었으며 최근 다수의 연구에서 ASPECT (Advanced Solver for Problems in Earth's ConvecTion; Kronbichler et al., 2012) 라는 지구동역학 소프트웨어를 활발히 이용하고 있다(예, Zhang and Li, 2018; Rajaonarison et al., 2020). ASPECT는 맨틀 대류를 모사하기 위해 열-대류 유한요소법을 기반으로 한 소프트웨어로서 다음과 같은 장점이 있다. i) 오픈 소스 소프트웨어로써 접근성이 용이하고 공용 라이브러리를 사용하기 때문에, 이용자와 개발자 간의 소통을 통한 즉각적인 기능 개발, 추가 및 적용으로 연산 기법의 최신화가 신속하다. ii) 취성 지각과 점성 맨틀을 동시에 포함하는 다상 유체 연산 기법, 자유 표면 변화를 구현하는 ALE 기법, 높은 해상도를 요구하는 국지적 영역에 대한 재격자화를 수행해주는 AMR (Adaptive Mesh Refinement) 기법(Gassmöller et al., 2018) 등을 지원하기 때문에 열개 현상을 모사하는데 용이하다.

자유 표면 관점에서 지표면 변형을 모사하기 위한 두 기법, 즉 ALE와 SA 기법의 정확성과 유효성은 현재까지 논쟁 중에 있으며(예, Crameri et al., 2012), ASPECT에서는 아직 열개 현상에 대해 두 기법을 비교한 연구가 부족하다. 열개 현상 등 시공간적으로 대규모인 지구조 운동은 긴 계산 시간이 필요하다. 수치 지구동역학 연구에 있어, 계산 시간과 모형 결과의 정확한 재현성은 연구 효율에 큰 영향을 미치기 때문에 두 요소를 모두 평가하고 적절한 자유 표면 모사 기법을 선택하는 것은 중요하다. 따라서 최근 열개 현상 수치 모형 제작에 널리 사용되고 있는 수치 모사 소프트웨어인 ASPECT를 사용해 동일한 계산 성능과 환경에서 두 기법을 사용한 2차원 열개 현상 수치 모형을 제작하고 각 모형 간의 정확성과 계산 시간의 차이를 비교하여 유효성과 효율성을 평가하고자 한다.

2. 연구방법

2.1 지배방정식

열개 현상을 지구동역학적으로 다룰 때 취성 지각과 지질학적 시간 단위(수백만 ~ 수천만 년)에서 점성 유체로 근사되는 특성을 지닌 맨틀을 동시에 포함하는 다상 유체 계산이 필요하다. 매우 느린 점성 유체의 거동을 기술하는(Kim and So, 2021) 스토크스(Stokes) 방정식(식 1)과 비압축성 유체를 가정한 질량 보존을 의미하는 연속 방정식(식 3)을 도입했다. 맨틀 물질 흐름이 온도차에 의한 부력에 의해서만 구동되는 부시네스크 근사(Boussinesq approximation)를 사용하였다.

$\[ -\text{∇}P+\text{∇}\cdot \tau =\rho g \text{where} \rho ={\rho }_0\left(1-\alpha \left(T-T_0\right.\right) \]$

(1)

$\[ \tau =2\eta \dot{\epsilon } \text{where} \dot{\epsilon }=\frac{1}{2}\left(\text{∇}u+{\left(\text{∇}u\right)}^T\right) \]$

(2)

$\[ \text{∇}\cdot u=0 \]$

(3)

P, τ, έ, u는 각각 압력, 편차응력 텐서, 변형률 속도 텐서, 속도장을 의미한다. α, T, g, η, ρ는 열팽창계수, 온도, 중력가속도, 점성도, 밀도를 표시한다. T₀는 기준 온도 293 K이며 ρ₀는 각 구성층의 T₀일 때의 밀도이다. 점성도 η는 연성 점성도 η_D (식 4)와 취성 점성도 η_B (식 5) 중 작은 값(식 6)으로 결정된다.

$\[ {\eta }_D=\frac{1}{2}{A}^{-\frac{1}{n}}{d}^{\frac{m}{n}}{\dot{\epsilon }}_{II}^{\frac{1-n}{n}}\exp \left(\frac{E+PV}{nRT}\right) \]$

(4)

$\[ {\eta }_B=\frac{{\sigma }_y}{2{\dot{\epsilon }}_{II}} \text{where} {\sigma }_y=C\text{cos}\varphi +P\text{sin}\varphi \]$

(5)

$\[ \eta =\min \left({\eta }_B,{\eta }_D\right) \]$

(6)

n, A, d, m, E, V는 각각 멱급수 지수, 전위 크리프 계수, 입자 크기, 입자 크기 계수, 활성화 에너지, 활성화 부피를 지시한다. R은 기체 상수이다. C는 점착도, ϕ는 마찰각이다. έ_II는 유효 변형률 속도로 $$ \sqrt{1/2\stackrel{\mathrm{\circ }}{\text{ε}}:\stackrel{\mathrm{\circ }}{\text{ε}}}$$로 표시되는 변형률 속도 텐서의 제2 불변량이다(So and Capitanio, 2016, 2017). 에너지 보존 방정식은 식 (7)과 같으며, 추가적인 열원으로 대륙 지각에 국한하여 동위원소 붕괴열만 고려한다.

$\[ \mathrm{\rho }{\mathrm{C}}_{\mathrm{p}}\left(\frac{\text{∂T}}{\text{∂t}}+\text{u}\cdot \text{∇T}\right)=\text{∇}\cdot \left(\text{k∇T}\right)+\mathrm{\rho }\mathrm{H} \]$

(7)

$\[ \frac{\text{∂}{\text{C}}_{\text{i}}}{\text{∂t}}+\text{u}\cdot \text{∇}{\text{C}}_{\text{i}}=0 \]$

(8)

C_p는 열 용량, k 및 H는 각각 열 전도도와 열 생산이다. 식 (8)은 조성의 이류 방정식으로, 첨자 i는 각 층, 즉 SA층, 상부 지각, 하부 지각, 암석권 맨틀을 지시한다. 자세한 물성은 표 1에 제시한다.

Table 1.

Values of parameters.a,b

2.2 수치 모형 구성

실험에서 사용된 2차원 열개 현상 모형은 일반적인 암석권 규모(Pascal et al., 2004)에서 발생하는 열개 현상을 기준으로 하여 수평과 수직 방향으로 600 km×200 km 크기를 가진다(그림 2a). 모형의 암석권 구성은 상부 지각과 하부 지각, 암석권 맨틀로 이루어져 있으며 각각의 두께는 20 km, 15 km, 165 km로 설정했다. SA 기법을 사용한 모형의 경우, 상부 지각 위에 다양한 두께(5 ~ 20 km)의 SA층을 추가하여 계산 속도와 열개 진화 과정에 SA층의 두께가 미치는 영향을 분석했다. ALE 기법을 사용한 경우에는 SA층을 정의하지 않았다. 스토크스 방정식에 대한 경계 조건으로써 양쪽 측벽에서 수평축 방향으로 각각 2 cm/yr의 속도로 일정하게 인장하였고 질량 보존 법칙을 만족시키기 위해 하부 경계에 양쪽 측벽에서 유출되는 질량만큼 유입시켰다. SA 기법을 사용한 모형의 경우 상부 SA 층의 질량을 보존하기 위해 상부 경계에 질량 유입 조건을 추가로 적용하였다(Wang et al., 2020). ALE 기법의 경우 상단 경계에 격자의 수직 변형만을 허용하는 자유표면 조건을 적용하였으며 SA층이 없기 때문에 추가적인 질량 유입 조건은 필요하지 않다.

Fig. 2.

(a) Model setup of lithospheric layer structure and boundary conditions. The blue, grey, beige, and brown regions indicate the SA layer, upper crust, lower crust, and lithospheric mantle. The red square located at the center of the lower crust refers to the mechanical weak zone defined by a positive temperature anomaly (i.e., +500 K). (b) Initial temperature profiles along yellow solid and green dashed lines in Fig. 2a. The temperature at a depth between 20 and 25 km in the green line displays a positive temperature anomaly. The temperature of the SA layer is fixed at 293 K. (c) Initial strength envelops with the reference strain-rate is 10-15 s-1. The strength reduction is exhibited at 20 to 25 km depth.

SA와 ALE 기법을 사용한 두 모형 모두 최상부의 온도는 293 K로 고정되어 있으며, SA층은 모든 부분에서 293 K로 구속하였다. 하부 지각의 중심부에 5 km×5 km 크기의 ~500 K의 양의 온도 이상을 적용하여 열개 발생의 씨앗(seed)으로 작동하게 하였다(Ammann et al., 2018; Pang et al., 2018). 초기 온도 조건으로는 정상 상태의 대륙 지각 온도 구배(Chapman, 1986)를 적용하였고 모형 하부 지각의 중심부에는 온도 이상을 위치시켜 지각의 기계적 강도를 주변부보다 약하게 만들어 열개의 시발점을 구축한다(그림 2b). 추가적인 열원은 방사성 동위원소 붕괴열만 고려하였다. SA층의 경우 열개 현상이 미치는 영향을 최소화하기 위해 매우 낮은 강도로 설정하였다. 초기 모형의 강도 계산을 위한 기준 변형률 속도는 10^-15 s^-1로 적용하였다(그림 2c). 모형의 총 진화 시간은 실험 모형 중에서 최초로 열개가 발생하고 해양지각의 폭이 100 km가 되는 시점인 5.3 Myrs로 설정했다. 모형에 적용된 자세한 기계적/열적 물성은 표 1에 제시하였다.

3. 연구 결과

3.1 기법 간 구조 및 진화 비교

열개 현상의 진화 초기에는 판의 상대 운동이나 암석권 하부의 플룸 상승 등이 발생시킨 인장력에 의해 지각이 얇아지고 지각 평형을 맞추기 위해 맨틀이 상승하는 과정을 거친다. 본 연구를 통해 제작된 열개 수치 모형에서 초기에 온도 이상이 위치한 지점의 지각이 주변에 비해 상대적으로 빠르게 인장하여 그 근방에 맨틀 상승이 시작된다. 그림 3은 지표의 변형을 모사하기 위해 ALE (그림 3a와 3c)와 SA (그림 3b, 3d) 기법을 도입한 모형의 각 층을 구성하는 조성 성분 변화를 열개 진화에 따라 2.7 Myrs와 5.3 Myrs 시점에서 보여준다. SA층의 두께는 20 km로 설정했다. 그림 3a, 3b는 각각 ALE와 SA 기법을 이용한 모형의 2.7 Myrs 시점으로 열개 진화 초기 단계로써 두 모형 모두 열개 중심부에서 지각이 얇아지는 현상을 보여준다. SA 기법을 도입한 모형의 경우 지각이 빠르게 얇아지고 하부 지각은 거의 분열된 반면에 ALE 기법을 도입한 모형은 지각 두께의 미세한 감소만이 관찰되었다. 수치 모형의 종료 시점인 5.3 Myrs에서 SA 기법을 사용한 모형(그림 3d)은 대륙붕괴 및 약 100 km의 대양저 확장을 보이는데 비해, ALE 기법을 사용한 모형(그림 3c)의 경우 대륙붕괴도 발생하지 않았음을 확인했다. ALE를 이용한 모형에서 열개 진화가 더 느린 현상은 SA층 상단에 적용된 질량 유입 조건의 영향으로 보이는데, 이 경계 조건은 경계 밖으로 향하는 운동을 상쇄하고, 수평 방향 흐름만을 허용하여 열개를 촉진시키기 때문이다. 모형의 상부로 향하는 운동을 완전히 제거하는 SA 기법과 달리 ALE 기법은 상승 운동을 허용함으로써 자유표면을 구현하기 때문에 진화 속도에 차이를 보이는 것으로 판단된다.

Fig. 3.

The distribution of compositions ALE (a and c) and SA (b and d) methods during rifting evolution in our numerical models. The top and bottom rows are the compositional distributions at time = 2.7 Myrs and 5.3 Myrs, respectively. At time = 5.3 Myrs, the model based on the SA method shows the 100 km wide fully-developed seafloor spreading (compare c and d).

SA와 ALE 기법을 사용한 모형은 진화 과정에서 보이는 암석권 맨틀 대류 양상에서도 상이한 결과를 보인다. 그림 4는 2.7 Myrs와 5.3 Myrs 시점의 온도 분포로 어두운색과 밝은색은 각각 저온(최저 293 K)과 고온(최대 1800 K)을 표현한다. 2.7 Myrs 시점에서, SA 기법을 사용한 모형은 지각이 빠르게 얇아지기 때문에 지각 평형을 만족하기 위한 빠른 하부 맨틀 상승류가 생성되지만(그림 4b), ALE 기법을 사용한 모형은 지각 두께의 변화가 느려서 맨틀 상승이 뚜렷하지 않다(그림 4a). 5.3 Myrs 시점에서는 두 모형은 공통적으로 열개 축 아래 맨틀 상승을 포함하는 명확한 맨틀 대류를 보여준다(그림 4c와 4d). SA 기법을 사용한 모형의 경우, 뜨거운 하부 맨틀이 열개 축 아래 기둥의 형태로 상승하여 지각에 도달한다. 그 후, 대양저 확장을 지시하는 지표의 수평 방향 확장이 발생하고 하부에서는 맨틀이 복잡하게 대류한다. ALE 기법을 사용한 모형에서도 열개 축 아래 중심 맨틀 상승 기둥과 수평 방향 맨틀 대류가 관찰되지만, 지표에서의 확장은 나타나지 않았다.

Fig. 4.

The thermal structure ALE (a and c) and SA (b and d) methods at time = 2.7 Myrs (top row) and 5.3 Myrs (bottom row). The faster mantle upwelling is developed in the SA model. At time = 5.3 Myrs, the hot mantle already reaches the SA layer, indicating the continental breakup.

그림 5는 SA와 ALE 기법을 사용한 두 모형에서 도출된 5.3 Myrs 시점의 점성도(그림 5a와 5b)와 변형률 속도(그림 5c와 5d) 분포를 도시함으로써 열개 진화 과정에서 형성되는 지질구조(대륙주변부 형태 및 전단대 등) 형태의 뚜렷한 차이를 보여준다. 5.3 Myrs 시점에서 SA 기법을 사용한 모형은 온도가 높아 낮은 점성도를 가진 물질이 지표에서 확장하며 열개 방향의 경사를 가진 성숙한 대륙주변부가 형성되는 데 비해(그림 5b), ALE 기법의 모형은 열개 축 주변에 낮은 점성도가 국지적으로 보이는 열개 진화의 초기 단계에 머무르고 있다(그림 5a). 정단층은 열개 현상이 진화하면서 발달하는 주요한 지질구조로 전단대의 형태로 추정할 수 있으며 지구대의 형태로 나타난다(그림 5c). 우리 모형에서는 낮은 점성도(그림 5a)와 높은 변형률 속도(그림 5c)를 가진 지역을 전단대로 간주한다(예, Brune et al., 2012; Lutz et al., 2021).

Fig. 5.

The distribution of viscosity (top row) and strain-rate (bottom row) in the numerical model using ALE (left column) and SA (right column) methods at time = 5.3 Myrs. In the ALE method, the crustal breakup did not occur (a). The SA method shows seafloor spreading (b) and high strain-rate (d) as the low viscosity mantle reaches the sticky-air.

3.2 Sticky-Air 층의 두께가 열개 형태와 연산 시간에 미치는 영향

SA 기법에서는 점성도가 낮은 층을 추가적으로 정의한다. 따라서, SA 기법을 사용한 모형은 같은 크기의 암석권에서 열개를 모사할 때 ALE 기법을 사용한 모형보다 더 큰 자유도(Degree of freedom)와 이에 따른 더 많은 연산 자원을 할당해야만 한다. DOF의 크기는 연산 시간에 직접적인 영향을 주고 연산 자원에 따라 한계 DOF가 제한되기 때문에, 모형에 따라 적절한 DOF를 설정하는 것은 연산 자원의 효율적 이용으로 이어진다. 두꺼운 SA층을 설정할수록 질량 유입 경계조건을 가진 모형 최상부와 상부 지각의 표면 사이의 거리가 멀어지고 그 영향을 적게 받기 때문에, 자유 표면에 가까워질 것으로 기대된다. 하지만, 더 두꺼운 SA층은 더 많은 연산 자원을 필요로 하기 때문에, SA층의 두께가 열개 형태 진화와 연산 시간에 미치는 영향을 정량적으로 파악해야 한다. 따라서, 우리는 SA에 다양한 두께(20 km, 15 km, 10 km, 5 km)를 적용하여 계산 결과를 비교하였다.

그림 6은 SA층의 두께를 20 km (1행), 15 km (2행), 10 km (3행), 5 km (4행)로 적용한 모형의 5.3 Myrs 시점에서의 결과로, 조성(1열)과 점성도(2열)의 분포를 보여준다. SA층의 두께가 감소할수록 열개 현상의 진행 속도가 느려지는 경향을 확인했다. 예를 들어, SA층이 20 km인 경우 대양저 확장이 100 km인데 비해(그림 6a), SA층이 5 km인 경우 30 km (그림 6g)에 불과하다. 뿐만 아니라, SA층의 두께 감소는 상부 지각의 변형에도 영향을 미친다. 이는, SA층의 두께 감소로 인해 열개 중심부가 상부의 경계에 가까워지면서 물질 유입 경계 조건에 더 지배적인 영향을 받은 것이 원인으로 추정된다.

Fig. 6.

The distribution of composition (left column) and viscosity (right column) with different values of SA layer thickness at time = 5.3 Myrs. The numerical result shows that the rate of rifting evolution becomes slower as the thickness of the SA layer decreases.

그림 7는 SA 기법을 사용한 모형에서 SA층 두께에 따른 DOF (가로축)가 전체 연산 시간(세로축)에 미치는 영향을 보여준다. 붉은색 동그라미는 SA 기법 기반의 열개 모형의 DOF를 표현한다. 정량적 비교를 위해 모든 모형은 동일한 CPU (Intel Xeon-E7-4850v4-2.1 Ghz, 코어 30개)와 RAM (SAMSUNG 1TB, 2133 Mhz)을 사용하여 MPI 병렬 처리 기반 환경에서 제작되었다. SA층의 두께가 20 km, 15 km이면, 총 연산 시간은 45,000 초 정도로 유사하다. 그러나, 두께가 10 km와 5 km로 감소하면 41,000 초 정도로 10%가량 감소하였다. 그에 비해 ALE 기법은 가장 낮은 DOF를 가짐에도 불구하고 약 86,000 초로 가장 긴 계산 시간을 사용했다. 이 결과는 ALE 기법에서 격자의 변형을 계산하기 위해 추가되는 계산 시간이 SA층에 의한 추가 계산 시간 보다 상대적으로 큰 것으로 해석된다.

Fig. 7.

A comparison of computing time depending on the thickness of the SA layer. Orange solid circles denote the total calculation time of a rifting model using SA methods. As the thickness of the SA layer increases from 5 km to 20 km, which leads to DOF increase from 6 million to 6.4 million, the calculation time increases by ~10%. The dashed line refers to the calculation time in the ALE method with DOF = 4.8 million.

4. 결론 및 토의

본 연구는 최근 열개 현상 수치 모사 연구에 주로 이용되고 있는 유한 요소 소프트웨어인 ASPECT를 사용하여 자유 표면 모사 기법인 SA와 ALE 기법을 도입한 열개 모형을 제작하고 SA의 유무 및 두께 차이에 따른 열개 구조 진화와 연산 시간 등을 비교하였다. 열개 현상 수치 모사에서 일반적으로 적용되는 점성도 제한 범위인 10¹⁸~10²⁴ (Candioti et al., 2020) 조건에서 수행된 수치 모사 실험의 5.3 Myrs 시점 결과에서, SA와 ALE 기법과 상관없이, 암석권 하부에서 열개 중심부를 향해 상승하는 맨틀 기둥과 열개 인장 방향으로 거동하는 맨틀 대류가 유사하게 관찰되었다. SA 기반 모형은 대륙 분열 및 해양 지각 형성까지 보였지만, ALE 기반 모형은 열개 중심부의 침하와 맨틀 상승만을 보였다. 이러한 열개 진화 속도의 차이는, SA 기법에서는 ALE 기법과 달리 상단에 적용된 질량 유입 조건이 지각 변형에 영향을 미쳤기 때문으로 해석된다.

SA기법에 있어서, 상단 경계 조건이 열개 진화 양상에 미치는 영향을 평가하기 위해 SA층의 두께를 20 km에서 5 km까지 다양하게 설정한 모형을 제작 및 계산하여 그 결과를 관찰했다. SA층이 얇아질수록, 상부 지각이 상단 경계와 가까워지면서 열개 진행 속도가 감소하였고 상부 지각의 지표 변형이 증가하였다. 이와 같은 결과는 SA층이 열개 현상과 분리되지 못하고 동역학적 상호작용하는 것을 지시하며, 열개 진화 과정에서 나타나는 지표면의 변형이 아닌 경계 조건에 의한 변형이기 때문에 열개 현상 수치 모사 결과 해석에 부정적인 영향을 미칠 수 있다. 따라서 SA층의 두께는 열개 현상에 영향을 주지 않을 만큼 충분히 크게 설정되어야 함과 동시에 연산 자원에 관한 고려도 필수적이다. 본 연구에 따르면, 200 km 두께의 암석권의 열개를 모사함에 있어서, 15 km와 20 km의 SA층 두께를 도입한 모형이 서로 유사한 결과를 도출했으므로, 해당 두께를 가진 SA층을 도입하는 것이 타당한 것으로 판단된다. 유사한 DOF를 가진 ALE 기법과 SA층의 두께가 5 km인 두 모형에서는 연산 시간이 각각 86,000 초와 45,000 초로 약 2배 차이가 발생함을 확인했다. ALE 기법에서는 격자 변형과 그에 따른 강성 행렬 재조립으로 인해, 많은 연산 시간이 소요되는 것으로 판단된다. 20 km의 SA층 두께에서는 계산 시간이 약 45,000 초, 5 km SA층 두께에서는 ~10% 감소하여 약 41,000 초가 소요되었으며, 이는 SA층 두께 감소에 따른 DOF 감소에 기인한 것으로 보인다.

현재 지구에는 태평양과 같은 대표적인 성숙한 열개부터 생성 초기의 홍해까지 다양한 열개 지역이 존재하며 그 진화와 생성 과정에 대한 연구는 아직까지 부족하다. 한반도의 동쪽에 위치한 동해 또한 태평양 판의 수렴 과정에서 발생한 해구 후퇴에 의해 인장하여 생성된 후열도 분지로서(Kim et al., 2013; Yoon et al., 2014) 그 생성 과정은 당겨 열림(Pull-Apart; Lallemand and Jolivet, 1986), 부채꼴 확장(Pan-shape; Otofuji and Matsuda, 1987), 두 단계 확장(Two-Stage Rifting; Hayashida et al., 1991; Lee et al., 1999) 등 다양한 가설들이 경쟁 중이다. 지구동역학적 컴퓨터 수치 모사 연구는 대규모 판 구조의 형성 및 진화를 연구하는데 효과적인 도구로써, 동해와 같은 열개 진화 및 생성 과정에 대한 가설들의 물리적 타당성을 검증할 수 있을 것으로 기대된다. 본 연구는 유한 요소 지구동역학 소프트웨어인 ASPECT를 사용하여 열개 현상 수치 모사에서 자유 표면 변형을 계산하기 위해 주로 사용되는 SA와 ALE 기법의 열개 진화 및 구조와 연산 시간을 비교하였다. SA 기법은 연산 시간이 비교적 적게 소요되지만, 열개 진화가 상단 경계 조건에 영향을 받는 것을 최소화하기 위한 적절한 SA층의 두께와 같은 조건이 필요하다. ALE 기법은 연산 시간은 느리지만 상부 지각의 표면에서 보이는 전단대와 같은 변형 구조 및 진화가 실제와 더 유사한 결과를 보였다. 본 연구 결과는 열개에 따른 지표 변형을 자유 표면 관점에서 모사하기 위한 두 기법인 SA와 ALE 간의 합리적 선택에 대한 방향성을 제시할 수 있다.

Acknowledgments

본 연구는 “행정안전부 방재안전분야 전문인력 양성”사업과 한국연구재단 중견연구지원사업(2022R1A2C1009742)과 중점연구소지원 사업(No.2019R1A6A1A03033167), 한국해양과학기술원의 지원을 받아 수행되었다(PEA0084).

References

• Allken, V., Huismans, R.S. and Thieulot, C., 2011, Three-dimensional numerical modeling of upper crustal extensional systems. Journal of Geophysical Research: Solid Earth, 116. [https://doi.org/10.1029/2011JB008319]
• Ammann, N., Liao, J., Gerya, T. and Ball, P., 2018, Oblique continental rifting and long transform fault formation based on 3D thermomechanical numerical modeling. Tectonophysics, 746, 106-120. [https://doi.org/10.1016/j.tecto.2017.08.015]
• Bohannon, R.G., Naeser, C.W., Schmidt, D.L. and Zimmermann, R.A., 1989, The timing of uplift, volcanism, and rifting peripheral to the Red Sea: A case for passive rifting?. Journal of Geophysical Research: Solid Earth, 94, 1683-1701. [https://doi.org/10.1029/JB094iB02p01683]
• Brune, S., Popov, A.A. and Sobolev, S.V., 2012, Modeling suggests that oblique extension facilitates rifting and continental break-up. Journal of Geophysical Research: Solid Earth, 117. [https://doi.org/10.1029/2011JB008860]
• Burov, E.B., 2011, Rheology and strength of the lithosphere. Marine and Petroleum Geology, 28, 1402-1443. [https://doi.org/10.1016/j.marpetgeo.2011.05.008]
• Candioti, L.G., Schmalholz, S.M. and Duretz, T., 2020, Impact of upper mantle convection on lithosphere hyperextension and subsequent horizontally forced subduction initiation. Solid Earth, 11, 2327-2357. [https://doi.org/10.5194/se-11-2327-2020]
• Chapman, D.S., 1986, Thermal gradients in the continental crust. Geological Society, London, Special Publications, 24, 63-70. [https://doi.org/10.1144/GSL.SP.1986.024.01.07]
• Corti, G., 2012, Evolution and characteristics of continental rifting: Analog modeling-inspired view and comparison with examples from the East African Rift System. Tectonophysics, 522-523, 1-33. [https://doi.org/10.1016/j.tecto.2011.06.010]
• Crameri, F., Schmeling, H., Golabek, G.J., Duretz, T., Orendt, R., Buiter, S.J.H., May, D.A., Kaus, B.J.P., Gerya, T.V. and Tackley, P.J., 2012, A comparison of numerical surface topography calculations in geodynamic modelling: an evaluation of the ‘sticky air’ method. Geophysical Journal International, 189, 38-54. [https://doi.org/10.1111/j.1365-246X.2012.05388.x]
• d’Acremont, E., Leroy, S. and Burov, E.B., 2003, Numerical modelling of a mantle plume: the plume head-lithosphere interaction in the formation of an oceanic large igneous province. Earth and Planetary Science Letters, 206, 379-396. [https://doi.org/10.1016/S0012-821X(02)01058-0]
• Franke, D., 2013, Rifting, lithosphere breakup and volcanism: Comparison of magma-poor and volcanic rifted margins. Marine and Petroleum Geology, 43, 63-87. [https://doi.org/10.1016/j.marpetgeo.2012.11.003]
• Frizon de Lamotte, D., Fourdan, B., Leleu, S., Leparmentier, F. and de Clarens, P., 2015, Style of rifting and the stages of Pangea breakup. Tectonics, 34, 1009-1029. [https://doi.org/10.1002/2014TC003760]
• Fullsack, P., 1995, An arbitrary Lagrangian-Eulerian formulation for creeping flows and its application in tectonic models. Geophysical Journal International, 120, 1-23. [https://doi.org/10.1111/j.1365-246X.1995.tb05908.x]
• Gassmöller, R., Lokavarapu, H., Heien, E., Puckett, E.G. and Bangerth, W., 2018, Flexible and scalable particle-in-cell methods with adaptive mesh refinement for geodynamic computations. Geochemistry, Geophysics, Geosystems, 19, 3596-3604. [https://doi.org/10.1029/2018GC007508]
• Guyonnet-Benaize, C., Lamarche, J., Hollender, F., Viseur, S., Münch, P. and Borgomano, J., 2015, Three-dimensional structural modeling of an active fault zone based on complex outcrop and subsurface data: The Middle Durance Fault Zone inherited from polyphase Meso-Cenozoic tectonics (southeastern France). Tectonics, 34, 265-289. [https://doi.org/10.1002/2014TC003749]
• Hayashida, A., Fukui, T. and Torii, M., 1991, Paleomagnetism of the Early Miocene Kani Group in southwest Japan and its implication for the opening of the Japan Sea. Geophysical Research Letters, 18, 1095-1098. [https://doi.org/10.1029/91GL01349]
• Huismans, R.S. and Beaumont, C., 2014, Rifted continental margins: The case for depth-dependent extension. Earth and Planetary Science Letters, 407, 148-162. [https://doi.org/10.1016/j.epsl.2014.09.032]
• Kim, D. and So, B., 2020, Effects of rheology and mantle temperature structure on edge-driven convection: Implications for partial melting and dynamic topography. Physics of the Earth and Planetary Interiors, 303, 106487. [https://doi.org/10.1016/j.pepi.2020.106487]
• Kim, G.B., Yoon, S.H., Sohn, Y.K. and Kwon, Y.K., 2013, Wave-planation surfaces in the mid-western East Sea (Sea of Japan): Indicators of subsidence history and paleogeographic evolution of back-arc basin. Marine Geology, 344, 65-81. [https://doi.org/10.1016/j.margeo.2013.07.008]
• King, S.D., Raefsky, A. and Hager, B.H., 1990, Conman: vectorizing a finite element code for incompressible two-dimensional convection in the Earth's mantle. Physics of the Earth and Planetary Interiors, 59, 195-207. [https://doi.org/10.1016/0031-9201(90)90225-M]
• Kronbichler, M., Heister, T. and Bangerth, W., 2012, High accuracy mantle convection simulation through modern numerical methods. Geophysical Journal International, 191, 12-29. [https://doi.org/10.1111/j.1365-246X.2012.05609.x]
• Lallemand, S. and Jolivet, L., 1986, Japan Sea: a pull-apart basin?. Earth and Planetary Science Letters, 76, 375-389. [https://doi.org/10.1016/0012-821X(86)90088-9]
• Lee, Y.S., Ishikawa, N. and Kim, W.K., 1999, Paleomagnetism of Tertiary rocks on the Korean Peninsula: tectonic implications for the opening of the East Sea (Sea of Japan). Tectonophysics, 304, 131-149. [https://doi.org/10.1016/S0040-1951(98)00270-4]
• Liao, J. and Gerya, T., 2015, From continental rifting to seafloor spreading: Insight from 3D thermo-mechanical modeling. Gondwana Research, 28, 1329-1343. [https://doi.org/10.1016/j.gr.2014.11.004]
• Lutz, B.M., Axen, G.J., van Wijk, J.W. and Phillips, F.M., 2021, Whole-lithosphere shear during oblique rifting. Geology, 50, 412-416. [https://doi.org/10.1130/G49603.1]
• May, D.A., Brown, J. and Le Pourhiet, L., 2014, pTatin3D: High-performance methods for long-term lithospheric dynamics. In SC'14: Proceedings of the International Conference for High Performance Computing, Networking, Storage and Analysis, p. 274-284. [https://doi.org/10.1109/SC.2014.28]
• Moresi, L. and Parsons, B., 1995, Interpreting gravity, geoid, and topography for convection with temperature dependent viscosity: Application to surface features on Venus. Journal of Geophysical Research: Planets, 100, 21155-21171. [https://doi.org/10.1029/95JE01622]
• Moresi, L., Quenette, S., Lemiale, V., Mériaux, C., Appelbe, B. and Mühlhaus, H.B., 2007, Computational approaches to studying non-linear dynamics of the crust and mantle. Physics of the Earth and Planetary Interiors, 163, 69-82. [https://doi.org/10.1016/j.pepi.2007.06.009]
• Nikolaeva, K., Gerya, T.V. and Marques, F.O., 2010, Subduction initiation at passive margins: Numerical modeling. Journal of Geophysical Research: Solid Earth, 115. [https://doi.org/10.1029/2009JB006549]
• Otofuji, Y. and Matsuda, T., 1987, Amount of clockwise rotation of Southwest Japan - fan shape opening of the southwestern part of the Japan Sea. Earth and Planetary Science Letters, 85, 289-301. [https://doi.org/10.1016/0012-821X(87)90039-2]
• Pang, Y., Zhang, H., Gerya, T.V., Liao, J., Cheng, H. and Shi, Y., 2018, The Mechanism and Dynamics of N-S Rifting in Southern Tibet: Insight From 3-D Thermomechanical Modeling. Journal of Geophysical Research: Solid Earth, 123, 859-877. [https://doi.org/10.1002/2017JB014011]
• Pascal, C., Cloetingh, S.A.P.L. and Davies, G.R., 2004, Asymmetric lithosphere as the cause of rifting and magmatism in the Permo-Carboniferous Oslo Graben. Geological Society, London, Special Publications, 223, 139-156. [https://doi.org/10.1144/GSL.SP.2004.223.01.06]
• Rajaonarison, T.A., Stamps, D.S., Fishwick, S., Brune, S., Glerum, A. and Hu, J., 2020, Numerical Modeling of Mantle Flow Beneath Madagascar to Constrain Upper Mantle Rheology Beneath Continental Regions. Journal of Geophysical Research: Solid Earth, 125, e2019JB018560. [https://doi.org/10.1029/2019JB018560]
• Rosendahl, B.R., Reynolds, D.J., Lorber, P.M., Burgess, C.F., McGill, J., Scott, D., Lambiase, J.J. and Derksen, S.J., 1986, Structural expressions of rifting: lessons from Lake Tanganyika, Africa. Geological Society, London, Special Publications, 25, 29-43. [https://doi.org/10.1144/GSL.SP.1986.025.01.04]
• Ruh, J.B., 2019, Effects of fault-weakening processes on oblique intracontinental rifting and subsequent tectonic inversion. American Journal of Science, 319, 315-338. [https://doi.org/10.2475/04.2019.03]
• Sarrate, J., Huerta, A. and Donea, J., 2001, Arbitrary Lagrangian-Eulerian formulation for fluid-rigid body interaction. Computer Methods in Applied Mechanics and Engineering, 190, 3171-3188. [https://doi.org/10.1016/S0045-7825(00)00387-X]
• So, B. and Capitanio, F.A., 2016, The emergence of seismic cycles from stress feedback between intra-plate faulting and far-field tectonic loading. Earth and Planetary Science Letters, 447, 112-118. [https://doi.org/10.1016/j.epsl.2016.05.002]
• So, B. and Capitanio, F.A., 2017, The effect of plate-scale rheology and plate interactions on intraplate seismicity. Earth and Planetary Science Letters, 478, 121-131. [https://doi.org/10.1016/j.epsl.2017.08.020]
• Turcotte D.L. and Emerman, S.H., 1983, Mechanisms of Active and Passive Rifting. Developments in Geotectonics, 19, 39-50. [https://doi.org/10.1016/B978-0-444-42198-2.50010-9]
• van der Beek, P., Cloetingh, S. and Andriessen, P., 1994, Mechanisms of extensional basin formation and vertical motions at rift flanks: Constraints from tectonic modelling and fission-track thermochronology. Earth and Planetary Science Letters, 121, 417-433. [https://doi.org/10.1016/0012-821X(94)90081-7]
• Wang, K., Chen, L., Xiong, X., Yan, Z. and Xie, R., 2020, The role of lithospheric heterogeneities in continental rifting: Implications for rift diversity in the North China Craton. Journal of Geodynamics, 139, 101765. [https://doi.org/10.1016/j.jog.2020.101765]
• Wenker, S. and Beaumont, C., 2018, Can metasomatic weakening result in the rifting of cratons?. Tectonophysics, 746, 3-21. [https://doi.org/10.1016/j.tecto.2017.06.013]
• Wright, T.J., Sigmundsson, F., Pagli, C., Belachew, M., Hamling, I.J., Brandsdóttir, B., Keir, D., Pedersen, R., Ayele, A., Ebinger, C., Einarsson, P., Lewi, E. and Calais, E., 2012, Geophysical constraints on the dynamics of spreading centres from rifting episodes on land. Nature Geoscience, 5, 242-250. [https://doi.org/10.1038/ngeo1428]
• Yoon, S.H., Sohn, Y.K. and Chough, S.K., 2014, Tectonic, sedimentary, and volcanic evolution of a back-arc basin in the East Sea (Sea of Japan). Marine Geology, 352, 70-88. [https://doi.org/10.1016/j.margeo.2014.03.004]
• Zhang, N. and Li, Z., 2018, Formation of mantle “lone plumes” in the global downwelling zone - A multiscale modelling of subduction-controlled plume generation beneath the South China Sea. Tectonophysics, 723, 1-13. [https://doi.org/10.1016/j.tecto.2017.11.038]