콤솔 멀티피직스를 이용한 점탄성 자유 섭입 수치 모형에 관한 벤치마크 연구: 2차원 라그랑지안 유한요소법 틀 내에서

1. 서 론

해령에서 생성된 해양 암석권이 확장 과정과 그에따른 냉각을 거치면서 밀도가 증가하여 맨틀 심부로 침강하는 섭입은 판 구조 현상의 주된 구동력으로 작용한다는 점과 지구 내부 물질의 지구동역학적 순환이라는 관점에서 중요한 현상이다. 섭입대 하부를 시각화 해주는 다수의 기존 연구(Song and Simons, 2003; Huang and Zhao, 2006; Hsu et al., 2013; Tang et al., 2014; Chang et al., 2016; Barker et al., 2018)는 섭입대의 현시점 정보를 영상화하여 제공해주어 지구의 표면과 내부에서 발생하는 다양한 현상 사이의 지구동역학적 상관관계를 이해하는 데 결정적 역할을 했다. 하지만, 이러한 연구는 섭입대의 현재 상태만을 보여주기 때문에, 수치 모형(numerical model)을 통해 섭입 슬랩(subducting slab)의 변형과 이동 시간에 따른 진화 과정을 재현하는 것이 필요하다(Funiciello et al., 2003). 컴퓨터 계산 능력의 발전으로 인해 대용량 수치 계산이 가능해졌고, 고온/고압 환경에서의 암석실험(Ranalli and Murphy, 1987; Karato, 2010)과 후빙기 반동(postglacial rebound, Forte and Mitrovica, 1996), 지진 후 변형(postseismic deformation, Vergnolle et al., 2003; Wang et al., 2012)에 대한 측지학적 지구 물리 탐사, 슬랩의 하강 속도(Čížková et al., 2012) 등으로 유추한 맨틀 점성도를 구속하는 다양한 연구는 수치 모형에 있어 중요한 입력 변수인 유변학적 물성(rheology)을 제공하여 섭입에 관한 수치 모사는 계속 발전할 수 있었다.

섭입에 대한 수치 모사 기법에는 크게 두 가지 접근방법이 있는데 하나는 계산 영역(domain)을 변형하지 않는 고정된 격자(grid)로 나누고, 상판(overriding plate, 주로 대륙판)과 섭입하는 하판(subducting plate, 주로 해양판) 및 맨틀을 순수 점성 유체(purely viscous fluid)로 가정하는 오일러리안(Eulerian) 방식이다(Schmeling et al., 1999; Gerya and Yuen, 2003; Kneller et al., 2005; Čížková et al., 2007). 이러한 방법은 지진파 단층촬영 연구에서 유추한 다양한 형태의 슬랩을 동역학적으로 재현하는 데 성공하였고, 슬랩의 탈수화 작용, 맨틀 쐐기에서의 구석 유동(mantle corner flow)과 같은 슬랩과 맨틀 사이의 상호작용에 대한 해석에 기여하였다. 오일러리안 방법에서는 맨틀의 점성도 보다 특정 배율(100-10,000배 정도)로 높은 점성도를 갖는 부분을 슬랩으로 가정하고 섭입 현상을 해석한다(Billen and Hirth, 2007; Andrews and Billen, 2009). 점성도가 주로 온도에 의해 결정되고 열은 전도되는 특성이 있기 때문에 슬랩과 맨틀 사이의 경계가 모호한 것은 오일러리안 방법의 한계라 할 수 있다. 또한 여러 가지 물성(즉, 취성, 탄성, 소성 및 점성)이 복합적으로 작용하여 발생하는 섭입 현상을 점성이라는 물성 하나만으로 해석해야 한다는 단점이 있다.

다른 한 가지 접근법은 변형할 수 있는 격자를 사용하여 섭입 슬랩의 거동을 해석하는 라그랑지안(Lagrangian) 방법으로(Giunchi et al., 1996; Houseman and Gubbins, 1997; Toth and Gurnis, 1998; Buiter et al., 2001) 슬랩을 주로 점탄성 고체로 가정한다. 사실 지구의 암석은 시간 규모(time scale)에 따라 점성과 탄성 거동을 모두 보여준다(Farrington et al., 2014). 이 중 탄성은 섭입하는 해양판이 구부러지는 정도에 영향을 주고, 섭입하면서 슬랩 내에 저장된 탄성 변형에너지가 지진의 형태로 방출되기 때문에(Fourel et al., 2014) 섭입 수치 모형에서 점탄성 물질을 도입하는 것은 중요하다. 라그랑지안 접근법은 해석 도중에 초기 설정한 격자의 변형과 이동이 발생한다는 특징이 있다. 따라서 계산 시작 전에 설정한 슬랩과 맨틀의 경계를 시간에 따라서 추적할 수 있고, 이 경계는 계산 시간 내내 명확하게 유지될 것이기 때문에 섭입 슬랩 내부의 탄성에너지 저장과 방출에 따라 변하는 응력 분포를 높은 해상도로 파악하기에 적합하다고 판단된다.

섭입 모형에 대한 다양한 연구가 진행되는 가운데, 기존의 라그랑지안 유한요소법을 이용한 섭입 슬랩 모형 개발은 주로 ABAQUS라는 상용 유한요소 소프트웨어를 바탕으로 진행되었다. ABAQUS는 라그랑지안 틀 내에서도 섭입 슬랩의 거동을 성공적으로 모사한(Funiciello et al., 2003; Capitanio et al., 2007, 2009; Fourel et al., 2014) 유한요소 소프트웨어지만, 주로 공학을 다루는데 유용하고 다룰 수 있는 물성이 탄소성으로 제한된다. 그에 비해 본 연구에서 활용한 상용 유한요소 패키지인 콤솔 멀티피직스(COMSOL Mutiphysics)는 광범위한 물성(즉, 점성, 탄성과 소성)을 간편하게 다룰 수 있다. 이 소프트웨어는 편리한 GUI (Graphic User Interface)을 기반으로 사용자가 계산 영역의 형태, 물성과 지배방정식 및 요소의 크기, 형태(삼각형 혹은 사각형)와 형상 함수의 차수 등을 자유롭게 선택할 수 있고, 대규모 병렬 연산을 지원함으로써 큰 규모의 수치 해석이 가능하다는 장점이 있다(Lee, 2013). 뿐만 아니라, 자바(Java) 및 매트랩(MATLAB) 등과 같은 해석형 프로그래밍 언어와 연동하여 많은 수의 수치 모형 해석을 자동화할 수 있다. 지난 10여 년간 이 소프트웨어를 통하여 오일러리안 관점에서 슬랩을 점성 유체로 고려한 연구는 있었지만(Hall, 2012; Hebert and Montesi, 2013; Lee and Wada, 2017) 라그랑지안 접근법의 틀 안에서 섭입 슬랩을 점탄성 고체로 접근한 모형은 개발 및 평가되어 본 적이 없었다. 자유 섭입 모형의 신뢰성을 더욱 향상하기 위해서 기존 연구에서 이용한 소프트웨어와 다른 소프트웨어를 사용하여 벤치마크 연구를 수행하여야 한다. 이를 통하여 맨틀 안으로 섭입하는 슬랩의 하강 속도와 슬랩 내부의 응력 분포를 계산하여 이 값이 기존 연구 결과와 유사함을 확인하였으며, 전 세계 섭입대에서 흔하게 관찰되는 현상인 해구 후퇴 및 전진에 영향을 주는 변수(예, 해양 암석권과 상부 맨틀의 밀도 차, 해령 밀기 속도, 깊이 660 km 상변이대에서의 점성도 증가 정도)의 효과를 정량적으로 측정하였다.

본 연구에서는 슬랩 내부 온도의 진화에 대해서도 실험을 수행하였다. 우리의 수치 모형은 열적 초기조건 및 경계조건을 도입하여 뚜렷한 혀 구조(tongue structure)를 도출했으며, 이를 바탕으로 차가운 슬랩 코어와 고압상으로 상변이되지 못한 준안정 감람석의 영역에 관해서 논의하였다. 추가로 부시네스크 근사(Boussinesq approximation)를 통하여 열 팽창에 의한 섭입 슬랩의 밀도 변화를 적용하고 슬랩의 동적인 온도 구조 변화를 구현했다. 슬랩의 온도와 동역학적 거동 사이의 상관관계를 열적 매개변수(thermal parameter, 슬랩의 나이와 섭입 속도의 곱으로 표현)라는 슬랩 내부의 열적 구조를 지시하는 변수와 관련지어 논의하였다.

2. 실험 방법

2.1 지배방정식 및 응력-변형률 관계

2차원 평면 변형률(plane strain) 상태를 가정한 자유 섭입 모형을 모사하기 위하여 사용된 지배방정식은 운동량 보존 방정식(식 1)과 에너지 보존 방정식(식 2)이다.

$\[ \frac{\partial {\sigma }_{ij}}{\partial x_j}=-△\mathrm{\rho }g{\delta }_{i2} \]$

(1)

$\[ \frac{DT}{Dt}=\frac{k}{{\rho }_L{C}_p}\frac{\partial }{\partial x_i}\left(\frac{\partial T}{\partial x_i}\right) \]$

(2)

σ_ij는 코시 응력 텐서(Cauchy’s stress tensor)이며 i와 j는 공간적 방향을 표시하는 색인으로 1일 때는 수평 방향, 2일 때는 수직 방향을 가리킨다. x_i는 i방향 좌표를 의미하고, ∂ / ∂x_i는 i방향으로의 편미분을 지시한다. δ_ij는 크로네커 델타(Kronecker delta) 함수이며 δ_i2는 수평 방향 (i = 1)일 때는 0, 수직 방향(i = 2)일 때는 1이다. 즉 식 (1)의 우항은 i=2일 때만 작동하며 이는 중력의 작용 방향을 반영하고 있다. Δρ는 해양 암석권과 상부 맨틀 사이의 밀도차, g는 중력가속도, T는 온도, t는 시간, k는 열전도도, ρ_L은 해양 암석권의 밀도, C_p는 열용량, D/Dt는 전미분(material derivative)이다. 방사성 동위원소 붕괴열 및 전단열, 단열 팽창 및 압축에 의한 열 감소와 증가는 무시하였다. 수치 모형에서 맨틀은 실재하지 않으므로 맨틀 내 대류에 의한 열전달은 고려하지 않고 오직 맨틀과 슬랩 경계면을 통한 전도만을 고려하였다. 지배방정식에 사용된 변수의 기호와 상세한 값은 표 1에 포함하였다.

Table 1.

Variable descriptions and values.

맨틀과 암석권은 거시적 관점에서 연속체로 가정할 수 있는데, 연속체 역학에서 물체의 변위 u_i(i방향의 변위)와 온도 T는 각각 식 (1)과 (2)를 유한요소법을 이용해서 요소마다 이산화한 후 전체 요소를 대상으로 조립한 거대한 선형계를 풀어서 구한다. 식 (2)의 경우 T에 관한 편미분 방정식이기 때문에 원래 형태대로 이산화하여 선형계를 풀이하면 T의 공간적 분포를 간단하게 구할 수 있다. 하지만, 식 (1)의 경우 응력 텐서에 관한 편미분으로 표현되어 있으므로 이 식을 이산화하여 u_i를 구하기 위해서는 u_i와 응력 텐서 σ_ij사이의 관계식이 필요하다. 따라서, 물체에 응력이 가해질 때 변형률 속도(strain rate)를 결정하는 조성 관계(constitutive relation)를 도입한다(Beuchert and Podladchikov, 2010; Schmalholz et al., 2014; So and Capitanio, 2017). 먼저 변형률(strain) 텐서 δ_ij는 식 (3)으로 정의한다.

$\[ {\epsilon }_{ij}=\frac{1}{2}\left(\frac{\partial u_j}{\partial x_i}+\frac{\partial u_i}{\partial x_j}\right) \]$

(3)

u_i는 i방향의 변위를 의미한다. 변형률 텐서를 시간에 대해 미분하면 변형률 속도(strain-rate) 텐서 $$ {\dot{\epsilon }}_{ij}$$가 도출된다. 식 (4)에서 $$ {\dot{\epsilon }}_{ij}$$와 $$ {\dot{\epsilon }}_{ij}^{\text{'}}$$는 각각 변형률 속도 텐서와 편향 변형률 속도 텐서이다.

$\[ {\dot{\epsilon }}_{ij}^{\text{'}}={\dot{\epsilon }}_{ij}-\frac{1}{2}\left({\dot{\epsilon }}_{11}+{\dot{\epsilon }}_{22}\right){\delta }_{ij} \]$

(4)

본 연구에서 도입한 섭입 슬랩은 탄성과 점성을 모두 가지고 있는 점탄성체이기 때문에, 대표적인 점탄성 거동 모형인 멕스웰 점탄성체(Maxwell viscoelastic body) 가정을 이용한다. 맥스웰 점탄성체의 경우, $$ {\dot{\epsilon }}_{ij}^{\text{'}}$$는 탄성 변형률 속도 텐서($$ {\dot{\epsilon }}_{ij}^{ela}$$)와 점성 변형률 속도 텐서($$ {\dot{\epsilon }}_{ij}^{vis}$$)의 단순 합(예, So and Yuen, 2015)으로 표현된다(식 5).

$\[ \begin{array}{l}{\dot{\epsilon }}_{ij}^{\text{'}}={\dot{\epsilon }}_{ij}^{ela}+{\dot{\epsilon }}_{ij}^{vis}\\ \text{where} {\dot{\epsilon }}_{ij}^{ela}=\frac{1}{2G}\frac{D{\tau }_{ij}}{Dt} \text{and} {\dot{\epsilon }}_{ij}^{vis}=\frac{{\tau }_{ij}}{2\eta }\end{array} \]$

(5)

τ_ij는 편향 응력 텐서(deviatoric stress tensor), G는 전단 계수, η는 동적 점성도(dynamic viscosity)다. 코시 응력 텐서와 편향 응력 텐서의 관계는 식 (6)와 같으며 p는 압력으로 식 (7)과 같이 체적 변형률(volumetric strain, ε₁₁ + ε₂₂)과 체적 탄성 계수(bulk modulus, K)의 음의 곱으로 표현된다.

$\[ {\sigma }_{ij}={\tau }_{ij}-p{\delta }_{ij} \]$

(6)

$\[ p=-K\left({\epsilon }_{11}+{\epsilon }_{22}\right) \]$

(7)

식 (1)에 식 (3)-(7)을 대입하면 σ_ij에 관한 식 (1)은 u_i에 대한 편미분 방정식으로 변환되며 이를 이산화하고 선형계를 조립하여 풀면 u_i의 공간적 분포를 도출할 수 있다.

2.2 수치 모형의 구성

대부분의 기존 섭입 모사 연구는 섭입 각도 및 속도를 부가하는 운동학적 경계조건(kinematic boundary condition)을 택하였다(Gerya and Yuen, 2003; Kneller et al., 2005; Běhounková and Čížková, 2008; Lee and King, 2009). 이러한 방법은 경계조건을 통해서 슬랩의 움직임을 미리 정해 주기 때문에 섭입하는 슬랩의 거동이 지배방정식 내에서 자체적으로 결정되지 않는다. 반면에, 본 연구에서는 중력에 의한 판 자체의 무게 때문에 발생하는 슬랩 당김 힘(slab pull force) 즉, 음의 부력(negative buoyancy) 효과에 의해 자발적으로 슬랩의 침강과 변형이 발생하는 자유 섭입 모형(free subduction model)을 구성하였다.

먼저, 역학적 경계조건으로 150 km 깊이까지 슬랩의 왼쪽 끝부분에 중력 방향의 힘을 부가하여(그림 1a의 F_imposed참조) 초기 섭입대의 형성에 중요한 슬랩의 굽힘 모멘트를 강제적으로 극복하였고, 그 이후부터는 중력만 작용하게 설정했다(Capitanio et al., 2007). 섭입하는 슬랩에 작용하는 주된 힘은 슬랩 당김힘(slab pull force)에 의한 음의 부력으로 이 힘은 해령 밀기 힘(ridge push force), 해구에서의 흡입력(trench suction force) 등 여러 가지 판 구조론적 구동력과 함께 주요한 역할을 하고 있다(Elsasser, 1971; Forsyth and Uyeda, 1975; Schellart, 2004). 맨틀 안으로 가라앉는 슬랩은 맨틀의 점성 저항을 크게 받고 이에 의해서 맨틀 쐐기에서의 구석 유동, 스테그넌트 슬랩(stagnant slab) 등 다양한 지구동역학적 현상이 발생한다. 자유 섭입 모형의 라그랑지안적 접근법에서는 맨틀이 실재하지 않으므로 슬랩의 경계에 10²¹ Pa·s의 점성도(Mitrovica, 1996)를 가지는 뉴토니안 대쉬팟(Newtonian dash-pot)을 설치하여 섭입하는 슬랩의 움직임에 대해 반대로 작용하는 점성 맨틀의 저항력을 근사하였다. 해양 암석권의 점성도는 10²³ Pa·s으로 설정하였다(Royden and Husson, 2006; Capitanio et al., 2007). 지구에서 발생하는 섭입은 판의 끝부분에서부터 해구를 따라서 순차적으로 가라앉는 형태로 발생하는데, 이를 라그랑지안 접근법으로 재현하기 위해서는 윈클러 파운데이션(Winkler foundation)을 적용해야 한다(Funiciello et al., 2003; Royden and Husson, 2006; Capitanio et al., 2007; van Dinther et al., 2010; Fourel et al., 2014; Keum, 2019). 윈클러 파운데이션은 퇴적물의 두께(3 km, 그림 1a)와 밀도(2,000 kg/m³, 그림 1a) 효과를 재현해주는데, 판의 침하를 억제하는 복원력으로써 작용한다. 이를 통해, 중력 방향으로 연속적인 섭입이 발생하게 유도하고, 해구에서의 해양 암석권의 휨을 촉진하는 역할을 한다.

Fig. 1.

(a) Mechanical model set-up. The model is assumed to 2D plain-strain, dominated by gravity, based on Linear Maxwell viscoelasticity. Isostatic restoring force (i.e. Winkler foundation) supports oceanic lithosphere. Viscous mantle resistance against subduction is reproduced by viscous dash-pot elements attached on surface of oceanic plate. (b) Thermal model set-up. In each time-step, depth-temperature boundary conditions T(y) enveloping the slab induces thermal conduction through the interface between slab and mantle. (c) Viscosity jump at the 660 km phase boundary with different magnitude.

열적 초기 및 경계조건으로는 반-공간 냉각 모델(half-space cooling model)을 사용하여 해양판의 길이는 2,500 km, 두께는 80 km로 가정하고 양쪽 끝에서 가장 어린 부분을 0 Ma, 가장 나이가 많은 부분을 150 Ma로 설정하여 해양판의 확장속도를 고려한 온도 분포를 적용하였다(Turcotte and Schubert, 2002). 깊이에 따른 맨틀 지온구배를 도입하여(Turcotte and Schubert, 2002) 해석 시간 단계(time step)마다 온도가 실시간으로 반영되도록 하여 슬랩과 맨틀의 경계에서부터 맨틀 열이 슬랩 내부로 전도되어 슬랩의 가운데 부분에서 열적 혀 구조(tongue structure)가 나타남을 확인하고자 한다(그림 1b).

한편, 지오이드, 후빙기 반동, 프리에어 중력 이상(free-air gravity anomaly), 동역학적 고도(dynamic topography) 등의 관측 결과는 깊이 660 km 상변이대를 경계로 하는 상부 맨틀과 하부 맨틀의 점성도 차이가 10 - 100배임을 제시하였고(Chen and King, 1998; Panasyuk and Hager, 2000; Mitrovica and Forte, 2004), 이를 우리 연구의 수치 모형에 반영하여 많은 지진파 단층 촬영 연구에서 제시된 슬랩의 형태를 재현하였다(그림 1c).

모든 모형 영역에 10 km × 10 km 크기를 갖는 동일한 사각형 요소를 사용하였으며, 총 요소 개수는 2,000 개로 계산을 수행하였다. 해석 시간 간격은 최소 32 kyrs, 최대 320 kyrs로 설정하였다. 수치 해석을 위하여 콤솔 멀티피직스에서 지원하는 내삽 방식을 이용한 시간 의존성 솔버(implicit time-dependent solver)와 유한요소법을 통한 이산화의 결과인 희소 행렬(sparse matrix)을 효율적으로 풀기 위한 직접 솔버(direct solver)인 MUMPS (MUltifrontial Massively Parallel Sparse)를 사용하여 구한 수치해를 바탕으로 섭입 과정을 모사할 수 있었다.

3. 결 과

3.1 섭입 슬랩의 역학적 거동

그림 2a는 시간에 따른 섭입 슬랩의 변형 과정과 슬랩 내부의 2차 편향 응력 불변량(second invariant deviatoric stress, τ_ⅠⅠ)의 공간적 분포를 보여 준다. 해양 암석권의 탄성 영률(E)과 동적 점성도(η_L)는 각각 2 × 10¹⁰ Pa, 10²³ Pa·s이며 상부 맨틀의 동적 점성도(η_M)는 10²¹ Pa·s, 해양 암석권과 상부 맨틀의 밀도 차(∆ρ)는 80 kg/m³이고 밀도는 온도에 의존하지 않도록 설정하였다. 660 km 깊이의 상변이대 즉, 상부와 하부 맨틀의 경계에서 맨틀의 점성도가 약 100배 증가하기(η₆₆₀ = 100 × η_M) 때문에 스테그넌트 슬랩이 형성되었으며, 이는 지진파 단층 촬영 연구에서 제시된 일부 지역에서의 섭입대 형상(예, 동아시아 및 보닌 지역 밑으로 섭입하는 태평양 슬랩, Zhao et al., 2009; Obayashi et al., 2017)과 유사함을 확인하였다. 한편, 전세계 섭입대 중 통가-커마덱(Tonga-Kermadec), 뉴 헤브리디스(New Hebrides) 섭입대 지역을 포함한 25% 이상의 섭입대에서 해구 전진이 관찰되기는 하지만(Cížková and Bina, 2015), 대부분의 섭입대에서 해구 후퇴 현상이 나타난다고 알려져 있다(Schellart et al., 2007, 2008). 본 연구에서 구현한 대부분의 자유 섭입 모형에서도 섭입이 진행됨에 따라 안정적인 해구의 후퇴를 보여주고 있으며 슬랩의 상부와 하부 힌지(hinge)에 각각 약 1.2 GPa와 0.9 GPa의 응력 집중이 나타났다(그림 2a). 이 응력 값은 Capitanio et al. (2009)과 van Dinther et al. (2010) 등의 기존 자유 섭입 모형에서 제시한 슬랩의 힌지에서 집중된 응력 크기의 차수와 일치한다. 또한 맨틀의 점성도에 의한 저항력을 극복하여 섭입하는 슬랩의 하강 속도는 최대 4.65 cm/yr였다(그림 2b). 슬랩의 하강 속도는 시간에 따라 증가하는데, 이는 섭입 과정 동안 맨틀 안으로 가라앉은 섭입 슬랩의 양이 증가함에 따라 슬랩 당김 힘도 증가하여 나타나는 현상이다.

Fig. 2.

Free subduction model. The subduction of oceanic lithosphere is self-consistently driven by slab pull force only. (a) Time evolution of second invariant of deviatoric stress in subducting slab. The model shows obvious trench retreat motion and the stress concentration along lower slab hinges. (b) Distribution of sinking velocity of slab with time. The sinking velocity is accelerated by slab pull force that increases as the slab descends into the upper mantle.

맨틀 내 슬랩의 하강 속도는 슬랩 당김 힘을 대변하는 상부 맨틀과 섭입판의 밀도차(∆ρ)와 가라앉는 슬랩에 저항하는 맨틀 점성도(η_M)에 의존한다(Capitanio et al., 2007). 슬랩의 하강 속도를 정량적으로 파악하는 것이 중요한 이유는 맨틀과 슬랩의 경계를 통해서 맨틀의 높은 온도가 슬랩 내부로 전달되는 시간에 따라 차갑고 기계적 강도가 높은 슬랩의 코어(core)의 존재 깊이와 폭이 결정되고, 나아가 슬랩 전체의 기계적 강도에 큰 영향을 주기 때문이다. 하강 속도가 느린 슬랩의 경우 전도에 의한 열 전달의 상대적 시간 단위가 길어지므로, 하강 속도가 빠른 슬랩보다 온도 구조가 더 뜨거울 것이다. η_M = 1 × 10²¹ Pa·s, ∆ρ = 50 kg/m³이 모형에서 슬랩 하강 속도(V_sink)는 2.54 cm/yr이며, 같은 η_M값과 ∆ρ = 80 kg/m³ 조건에서는 V_sink = 4.7 cm/yr로 슬랩의 하강 속도가 약 87% 증가하였다(그림 3a). η_M = 0.61 × 10²¹ Pa·s, ∆ρ = 50 kg/m³일 때 V_sink = 4.01 cm/yr, ∆ρ = 80 kg/m³ 조건에서는 V_sink = 7.41 cm/yr로 약 85% 증가하였으며, 맨틀의 점성도가 다른 경우에서도 밀도 차에 따라 비슷한 V_sink증가량을 보였다(그림 3a). 본 연구의 자유 섭입 모형에서 도출된 V_sink값의 범위는 Capitanio et al. (2007)의 2차원 유한요소 수치모형 연구에서 제시한 ∆ρ = 80 kg/m³일 때, 4.5 - 5 cm/yr의 범위인 것과 유사함을 확인하였다. 즉, 해양 암석권과 상부 맨틀의 밀도 차가 클수록 섭입하는 해양판 자체의 무게 때문에 슬랩 당김 힘이 증가하여 슬랩의 하강 속도가 상승하고, 맨틀의 점성도가 낮을수록 섭입하는 슬랩에 대한 점성 저항력이 감소하여 슬랩의 하강 속도가 상승함을 알 수 있으며, 이는 기존 자유 섭입 모형(Stegman et al., 2010; Goes et al., 2013) 뿐만 아니라 섭입 각도와 속도를 직접 부가하는 운동학적 경계조건(kinematic boundary condition)을 이용하는 순수 점성 유체역학을 기반으로 한 수치모형(Běhounková and Cížková, 2008; Lee and King, 2009)과도 유사한 결과이다. 추가로, 슬랩 당김 힘의 크기와 방향에 중요한 역할을 하는 해양 암석권과 상부 맨틀의 밀도 차에 따른 해구의 위치 및 후퇴 속도(V_r)를 측정하였다(그림 3b). 먼저, ∆ρ = 50 kg/m³일 때 V_r = 3.2 cm/yr (그림 3b의 초록색 선)이고, ∆ρ = 80 kg/m³일 때, V_r = 5.4 cm/yr (그림 3b의 노란색 선)로 계산되었다. 이를 그림 3a와 연관 지어 생각하면 음의 부력 효과에 의한 슬랩 당김 힘이 증가 할수록 슬랩의 하강 속도와 해구 후퇴 속도가 증가하는 것을 확인하였다. 본 연구의 모형에서 도출한 V_r값은 전 세계 섭입대에서 관찰된 해구 후퇴 속도 범위(0-15 cm/yr) 안에 포함된다(Schellart et al., 2007, 2008).

Fig. 3.

(a) Sinking velocity of slab with different values of Δρ and ηM. (b) Change of trench position during subduction progress in the cases with different values of Δρ.

섭입대 역학을 지배하는 주된 힘은 슬랩 당김 힘이라고 알려져 있지만, 해령 밀기 힘 역시 섭입 해양판의 거동에 중요한 역할을 한다(Elsasser, 1971; Forsyth and Uyeda, 1975). 특히 해령 밀기 힘은 해구의 운동방향에 영향을 미친다고 알려져 있다(Cížková and Bina, 2015). 본 연구에서는 해령 밀기 힘을 운동학적 경계조건에 반영하여 수치 해석을 진행하였다. 그림 4에 도시된 모형에는 해양 암석권의 영률과 점성도는 각각 E = 2 × 10¹⁰ Pa, η_L = 10²³ Pa·s이며 ∆ρ = 80 kg/m³으로 설정하였다. 이 모형의 슬랩의 오른쪽 끝에 해령 밀기 속도 V_p = 4 cm/yr를 부가하였을 때 섭입 슬랩의 시계열 진화과정을 보여준다(그림 4). 시간이 지남에 따라 섭입 슬랩 전체가 왼쪽(즉, 섭입 방향)으로 이동하는 것을 확인할 수 있다. V_p를 부가한 경우 슬랩이 더 고각으로 휘며 660 km 상변이대에 닿는 시간이 약 10 Myrs로 나타났고, 슬랩의 하부 힌지에 최대 3.5 GPa의 강한 응력 집중을 보여주었다. 3.5 GPa의 응력집중은 V_p = 0 cm/yr으로 설정한 모형에서 구한 0.9 GPa (그림 2 참조) 보다 강한 응력 집중으로 섭입하는 슬랩의 지표면에서의 움직임이 슬랩 심부까지 영향을 미침을 암시한다.

Fig. 4.

The spatial distribution of τⅠⅠ in different time- steps of free subduction model with an imposed kinematic boundary condition (i.e., ridge push).

우리는 육지방향으로 해령 밀기 속도(V_p)를 부가한 자유 섭입 모형에서 해구 전진 현상을 관찰하였으며, 해구 전진에 영향을 미치는 여러가지 변수에 따른 해구의 위치 변화 양상과 슬랩의 응력 분포에 관해 수치모사 하였다(그림 5). 먼저, V_p = 4 cm/yr로 고정하고, 해양 암석권과 상부 맨틀의 밀도 차, ∆ρ에 따른 해구의 시계열 위치 변화를 계산하였다(그림 5a). ∆ρ = 65 kg/m³일 때(그림 5a의 파란색 선과 빨간 점참조) 해구 후퇴 속도(V_r)은 0.7 cm/yr이고 깊이 660 km 상변이대에 닿는 시간은 12.5 Myrs인 반면, ∆ρ = 80 kg/m³일 때(그림 5a의 노란색 선과 빨간 점 참조), V_r = 2.0 cm/yr이고 9.6 Myrs에 슬랩 끝이 660 km 상변이대에 도달하는 것으로 나타났다. 이는 중력에 의한 슬랩 당김 힘이 증가할수록 해구의 후퇴 속도가 빠르고, 깊이 660 km 상변이대에 빠르게 도달한다는 것을 의미한다. 해양 암석권의 밀도가 높을수록, 즉 슬랩 당김 힘이 클수록 해구 전진의 평균 속도가 느리게 나타났다. 예를 들어, ∆ρ = 65 kg/m³일 때 해구 전진의 평균 속도 V_ad = 4.57 cm/yr이며, ∆ρ = 80 kg/m³일 때는 V_ad = 2.52 cm/yr로 계산되었다. 그림 5a에서 보여주는 결과로부터 ∆ρ의 증가에 따른 슬랩 당김 힘의 증가는 해구 후퇴를 촉진하였으며 반대로 해구 전진은 억제함을 추측할 수 있었다. 그림 5b는 ∆ρ = 80 kg/m³, V_p = 4 cm/yr일 때, 깊이 660 km 상변이대(phase transition zone)에서의 점성도 점프(viscosity jump)의 정도(10-100배)가 다른 4개의 모형에서 계산한 해구의 시계열 위치 변화를 보여준다. 모든 모형에서 ∆ρ는 동일하므로 깊이 660 km에 닿기 전까지는 해구의 후퇴 속도는 일정하다(V_r = 2.0 cm/yr). 그러나 슬랩이 660 km에 닿은 이후(그림 5b의 빨간 점 참조) 해구 전진은 상부 맨틀과 하부 맨틀의 점성도 차이가 10배 일 때, 해구의 전진 속도가 느리고(V_ad = 1.56 cm/yr, 그림 5b의 파란색 선), 100배 일 때 빠른 속도(V_ad = 2.52 cm/yr, 그림 5b의 보라색 선)로 계산되었다. 상부와 하부 맨틀의 점성도 차이가 100배 일 때 의 변화에 따른 해구의 위치를 시간에 따라 측정하였다(그림 5c). V_p = 0, 1, 2 cm/yr인 경우(각각 그림 5c의 초록, 파란, 주황색 선), 660 km 상변이대에 닿은 후에는 후퇴 속도에 변화가 있기는 하지만 전체적으로 안정적인 해구 후퇴 현상이 나타난다. V_p = 3 cm/yr인 경우(그림 5c의 노란색 선), 660 km 상변이대에 닿은 이후로 해구의 위치가 거의 고정되어 있다. V_p = 4 cm/yr인 경우(그림 5c의 보라색 선), 660 km 상변이대에 슬랩이 닿은 이후에 해구 속도 방향에 역전이 발생하여 해구 전진이 나타난다(V_ad = 2.52 cm/yr). 이는 상변이대에 닿기 전에는 슬랩의 하강 속도가 육지 방향으로 부가한 속도보다 빠르기 때문에 해구의 후퇴가 발생하고, 닿은 후부터는 슬랩 앵커링(slab anchoring) 현상이 발생하여 슬랩의 왼쪽 끝 부분이 상변이대에 고정됨과 동시에 해령 밀기 속도에 의한 효과 때문에 해구의 전진이 나타나는 것으로 해석된다.

Fig. 5.

(a) Change of trench position in the models with an imposed ridge push velocity (Vp = 4 cm/yr). Red dots indicate the time of touching 660 km phase boundary. Different colors indicate the different values of Δρ. (b) Trench position with different magnitude of viscosity jump (i.e., × 10 - 100) at 660 km phase boundary. The values of Δρ and Vp are fixed at 80 kg/m3 and , 4 cm/yr repectively. (c) Time evolution of trench position in the cases with varying values of Vp and a fixed value of Δρ = 80 kg/m3.

그림 6a는 V_p값을 달리하여 해석한 수치 모형 결과를 보여주며, AA' 곡선을 따라(그림 6a의 좌상부 참조) 편차 응력 텐서의 제 2 불변량 τ_ⅠⅠ를 도시하였다. V_p = 1 cm/yr일 때(그림 6a의 초록색 선), 슬랩의 하부 힌지에서 최대 0.97 GPa의 응력이 집중되고, V_p = 4 cm/yr일 때(그림 7a의 빨간색 선), 최대 3.5 GPa의 응력이 집중되는 것을 확인하였다. 즉, 해령 밀기 속도가 증가 할수록 더 강한 응력 집중이 발생함을 확인하였다. 그림 6b-c에서는 깊이 660 km 상변이대에서의 점성도 점프에 따른 슬랩의 하부 힌지의 응력 집중 양상을 계산하였다. 그림 6b는 해령 밀기 속도를 부가하지 않은 자유 섭입 모형에서의 결과이며, 상부와 하부 맨틀의 점성도 차이가 10배 일 때 0.5 GPa (그림 6b의 파란색 선), 100배 일 때 0.71 GPa (그림 6b의 초록색 선)의 응력이 계산되었다. 그림 6c는 V_p = 4 cm/yr의 해령 밀기 속도를 부가한 자유 섭입 모형의 계산 결과이며, 슬랩의 하부 힌지에 나타나는 응력 집중은 10배와 100배의 점성도 점프 경우 각각 1.9 GPa (그림 6c의 초록색 선)와 3.5 GPa(그림 6c의 빨간색 선)로 계산되었다. 일반적으로 660 km 깊이에서의 점성도 증가가 높을수록 슬랩 하부 힌지에 높은 크기의 응력이 집중되는 현상이 나타났다.

Fig. 6.

Second invariant deviatoric stress profiles along slab hinge line AA' at time = 19 Myrs. (a) The values of Vp controls the magnitude of stress concentration. (b) When the edge of subducting slab remains free (i.e., Vp = 0 cm/yr), the larger viscosity jump at 660 km phase boundary causes the larger stress along the hinge. (c) The models with non-zero Vp show larger stress (e.g., ~3.5 GPa with Vp = 4 cm/yr and η660 = 100 × ηM).

3.2 섭입 슬랩 내 온도 분포의 시계열 진화

수치 모형에 적용한 열적 경계조건 및 초기조건을 통해 섭입 과정 동안 진화하는 슬랩 내부의 온도 구조를 얻을 수 있었다. 그림 7a는 해양판의 초기 온도 구조를 150 Ma의 나이에 해당하도록 설정하였을 때, 섭입하는 슬랩의 시계열 온도 구조 변화이다. 맨틀의 높은 온도가 슬랩의 외부로부터 내부로 전도되는 것을 확인할 수 있다. 이를 통해 온도가 낮은 슬랩의 중심부가 맨틀에 가까워질수록 온도가 상승하여 섭입 슬랩의 특징적인 열적 구조로 알려진 혀 구조(tongue structure)와 그에 따른 차가운 슬랩 코어(cold slab core)가 생성되었다(그림 7a의 점선 영역). 슬랩 내부의 온도 구조에서 관찰된 혀 구조는 타 섭입대 수치 모사 논문에서 제시된 형태와 유사함을 확인했다(Funiciello et al., 2003; Běhounková and Cížková, 2008). AA' - FF' 각 측선을 따라 온도를 도시하면 혀 구조가 명확하게 존재함을 알 수 있다(그림 7b). 가령, 깊이가 가장 깊은 AA' 측선의 경우(그림 7b의 파란색 선) 다른 측선에 비해 섭입 시간이 길었기 때문에 맨틀로부터 열이 전도되는 시간이 길어 슬랩의 중심에서 1,606 K의 높은 온도를 보인다. 반면에, 상대적으로 깊이가 얕아 섭입 시간이 짧은 FF' 측선에서는 (그림 7b의 초록색 선) 슬랩의 중심 부분의 온도가 1,228 K로 뚜렷한 혀 구조가 나타났다.

Fig. 7.

(a) Temperature distribution of subducting slab calculated from free subduction model. (b) Temperature profiles along the lines AA'-FF'. (c) Temperature profiles along the line GG' with different values of thermal parameters. The model with larger thermal parameter shows colder core of subducting slabs (e.g., ~950 K and ~1000 K for thermal parameter of 825 km and 420 km, respectively).

본 연구에서 계산한 슬랩 내부의 온도 구조의 현실성을 평가하기 위해서 수치 모델 결과와 실제 섭입대를 비교하여 Kirby et al. (1996)가 제안한 열적 매개변수(thermal parameter)를 도입하였다. 열적 매개변수는 Φ = V_sub × age = V_sink × sinθ × age와 같이 정의된다. V_sub는 섭입 속도, θ는 해양 암석권의 나이이며 는 평균 섭입 각도이다. 해양 암석권의 나이가 많을수록 지상에서 오랜 시간 냉각되어 밀도가 높고, 이로 인해 슬랩 당김 힘의 크기가 증가하여 섭입 속도가 빠르다. 따라서 나이가 많은 해양 암석권은 섭입 시간이 짧기 때문에 맨틀 온도가 전도되는 시간 역시 짧아 차가운 슬랩의 코어가 크고 길게 발달한다. 그림 7c에서 해양 암석권의 나이는 150 Ma로 고정하고, ∆ρ = 50 - 95 kg/m³로 설정한 다수의 모형에서, 슬랩의 섭입 속도에 의해 결정된 열적 매개변수를 바탕으로 슬랩 내부의 GG' 측선을 따라 온도 분포를 도시하였다. 섭입 속도를 임의로 지정할 수 없는 자유 섭입 모형의 경우에는 섭입 속도와 각도가 물리적 틀 내에서 자발적으로 계산되기 때문에 일반적인 열적 매개변수의 개념을 적용하기에 어려움이 있지만, 최대 섭입 속도 및 평균 섭입 각도를 활용하였다. 그 결과, 슬랩의 하강 속도가 느린 일 때 열적 매개변수는 420 km이며 측선 ∆ρ = 50 kg/m³를 따라 가장 낮은 부분의 온도가 1,068 K로 나타났다(그림 7c의 자주색 선). 반면, ∆ρ = 80 kg/m³의 무거운 슬랩일 때 열적 매개변수는 825 km으로 계산되었으며, 같은 측선을 따라 984 K의 최저 온도를 보였다(그림 7c의 파란색 선). 이렇게 우리는 열적 매개변수와 슬랩 내부 온도 사이에 상관관계가 있음을 확인했다. 수치 모형에서 얻어진 열적 매개변수의 값은 56개의 섭입대에 대하여 지구 물리탐사 결과를 종합한 Syracuse et al. (2010)의 연구에서 제시한 캐스캐디아(Cascadia) 섭입대의 100 km와 통가-커마덱 섭입대의 14,000 km의 열적 매개변수 범위와도 부합된다는 것을 확인하였다.

수치 모형에 적용한 열적 초기 및 경계 조건으로부터 얻은 슬랩 내부의 온도 구조를 바탕으로 밀도가 온도와 열팽창 t 계수에 의존하는 부시네스크 근사(Boussinesq approximation)를 적용할 수 있다 (So and Yuen, 2015). 즉, 슬랩이 섭입하면서 맨틀과 인접한 경계에서 전도의 형태로 열이 전달되고, 슬랩 내부의 온도에 따라 열팽창하는 정도가 달라져 밀도에 변화가 발생한다. 본 연구에서는 부시네스크 근사에 따라 온도에 의존하는 상부 맨틀과 슬랩의 밀도차를 ∆ρ = ∆ρ₀(1 - α(T - T₀))로 표현했다. ∆ρ₀는 지표 온도를 기준으로한 해양 암석권과 상부 맨틀 사이의 기준 밀도차이며, α는 열팽창 계수, 그리고 T₀는 지표에서의 온도이다. 부시네스크 근사를 위한 물성값은 표 1에 제시하였다. 그림 8에서 확인할 수 있듯 온도가 낮은 슬랩의 중심부에서는 밀도가 높고, 온도가 높은 맨틀과의 경계 부분에서는 밀도가 낮기 때문에 균일한 밀도를 가지는 자유 섭입 모형보다 평균적으로 낮은 밀도를 보인다. 부시네스크 근사를 적용하지 않은 경우보다 적용한 경우에 중력 방향의 음의 부력 효과가 더 약하다는 것을 의미하며, 이는 슬랩의 하강 속도와 더불어 섭입대의 형상 및 해구의 위치 차이로 나타난다. 예를 들어, 부시네스크 근사를 적용한 경우, 21.2 Myrs에서의 섭입 슬랩의 왼쪽 끝과 해구의 수평 위치는 각각 111 km과 1,660 km 인 반면 그림 2a의 부시네스크 근사를 미적용한 모형의 경우, 같은 시간에서 각각 114 km와 1,730 km이다. 부시네스크 근사를 사용하지 않은 경우, 슬랩의 왼쪽 끝은 3 km, 해구는 70 km가 더 물러났다. 특히 부시네스크 근사의 사용 여부에 따라 나타나는 70 km 정도의 해구의 위치 변화는 무시할 수 없을 정도로 크다.

Fig. 8.

Free subduction model with Boussinesq approximation. The density within subducting slab changes dynamically during subduction process by thermal conduction.

그림 9는 균일한 밀도를 가정한 해양 암석권(점선)과 부시네스크 근사를 적용한 밀도 구조(실선)를 적용했을 때 섭입하는 슬랩의 하강 속도(V_sink)를 계산한 결과이다. 전체적으로 부시네스크 근사를 적용한 수치 모형에서 하강 속도가 느리다. 예를 들어, η_M = 10²¹ Pa·s이고 ∆ρ = 80 kg/m³일 때 V_sink = 4.7 cm/yr로 계산되었으며(그림 9의 하늘색 실선), 같은 상부 맨틀의 점성도이고, 부시네스크 근사를 적용한 밀도 차에서는 V_sink = 4.6 cm/yr으로(그림 9의 하늘색 점선), 2.1%의 미세한 섭입 속도 감소가 발생하였다. 이는 부시네스크 근사에 의한 밀도의 온도 의존성으로 인해 슬랩 당김 힘이 약해져 발생하는 것으로 판단된다.

Fig. 9.

Sinking velocity of subducting slabs with Boussinesq approximation (solid lines) for the cases with a wide range of ηM. The dotted lines mean the same cases without Boussinesq approximation.

4. 토 의

우리가 개발한 자유 섭입 모형은 전 지구상의 섭입대에서 흔하게 발견되는 해구 후퇴와 전진 현상을 성공적으로 재현하였으며, 그 속도를 정량적으로 측정할 수 있었다. 해구의 후퇴와 전진은 상판에 가해지는 인장력 및 압축력에 중요한 역할을 한다. 태평양판의 슬랩 롤백(slab rollback)에 의해 발생한 후열도 분지(back-arc basin)인 서 필리핀 분지(West Philippine Basin)와 한반도와 일본 열도 사이에 위치한 동해의 울릉 분지(Ulleung Basin)에 관한 기존 연구에서는 상판에 인장력이 작용하게 하는 해구 후퇴가 주요 분지 생성 기작이라고 주장하고 있다(Faccenna et al., 2009; Kim et al., 2015). 해구 전진은 해령 밀기 힘 또는 판 사이의 상호작용에 의해 발생하는 압축력에 의해 나타난다고 알려져 있다(Cížková and Bina, 2015). 한편, Choi et al. (2012)의 진원 기구 해석 연구와 Kim et al. (2018)의 울릉 분지의 좌굴 구조 연구에서는 동해에 동서 방향의 압축력이 가해지고 있다고 주장하였다. 이는 과거 인장력에 의해 확장되었던 동해가 압축력에 의해 닫히고 있으며 인장과 압축을 반복적으로 경험하였음을 암시한다. 따라서 동해를 포함한 후열도 분지의 시계열적 진화를 이해하기 위해서는 섭입하는 슬랩과 해구의 움직임(해구 전진 및 후진)을 해석하는 것이 필수적이다. 이와 관련하여 라그랑지안 섭입대 수치 모형을 통하여 운동학적 경계조건을 부가하는 것이 아니라 물리적 틀 안에서 자발적으로 슬랩의 변형이 계산되는 자유 섭입 모형에서 해구 후퇴와 전진 속도를 제시하는 것은 동해 등의 후배호 분지를 포함한 섭입대 일원의 판구조적 진화에 대한 이해를 증진할 것이다.

맨틀 대부분을 차지하는 규산염 광물인 감람석(Mg₂SiO₄)에 대한 410-660 km 상변이대에서의 암석 변형 실험은 14 GPa 정도의 고압 환경에서 수행해야 하는 기술적 어려움이 따르기 때문에, 주기율표상 규산염과 같은 족에 속하는 게르마늄(germanium)을 포함한 감람석과 동일한 구조의 Mg₂SiO₄를 이용하여 고온/고압 광물 실험을 수행해 왔다(Green and Burnley, 1989; Tingle et al., 1993). Schubnel et al. (2013)은 Mg₂SiO₄을 사용하여 차가운 슬랩 코어에 해당하는 1,000-1,250 K의 온도와 2-5 GPa의 압력 범위에서 암석 변형 실험을 수행하여, 400-700 km 깊이에서 2.5-3 GPa의 높은 편차 응력이 주어진다면 암석의 파괴 가능성이 거의 없는 고온/고압 환경일지라도 지진이 발생할 수 있다고 주장하였다. 우리의 수치 모형에서 슬랩 하부 힌지의 높은 응력 집중은 해령 밀기 속도의 부가 여부와 관계없이 발생하였고, Schubnel et al. (2013)의 고압 광물학 실험 결과를 종합하여 보면, 슬랩 힌지에서의 대 변형에 따른 힌지 근방의 암석 파괴와 심발지진 발생 가능성을 유추할 수 있다. 실제로, 이주-보닌-마리아나(Izu-Bonin-Mariana) 섭입계에 위치하는 보닌(Bonin) 섬 아래에서 2015년 5월 30일에 발생한 규모 7.9의 오가사와라(Ogasawara) 지진은 진원의 깊이가 680 km로 관측 이래 최고 심도이며, 심발지진 중 최대 규모로 알려져 있다(Ye et al., 2016). 이 지진은 지진파 단층 촬영 및 진원 기구 연구를 통해 보닌 슬랩의 하부 힌지에서 발생한 것으로 제시되었다(Obayashi et al., 2017). 오가사와라 지진의 발생 위치와 우리의 수치 모형에서 계산된 하부 힌지에서의 강한 응력 집중(0.9-3.5 GPa) 결과를 바탕으로 하부 힌지의 응력 집중이 오가사와라 심발 지진에 영향을 주었다고 판단할 수 있다.

본 연구에서 모사한 수치 모형은 슬랩의 열적 구조를 성공적으로 재현하였으므로 함수광물의 탈수화 조건, 온도 의존성 점성도, 준안정 감람석(metastable olivine)의 영역 파악 및 410 km와 660 km 상변이대의 고도(topography) 파악 등 다양한 지구동역학적 현상을 파악할 수 있을 것이다. 사실, 해양 암석권과 맨틀의 점성도는 온도와 응력에 의존하는 멱급수(power-law) 형태의 전위 크리프(dislocation creep) 또는 확산 크리프(diffusion creep)로 표현되어야 현실적인 점성도 구조를 반영할 수 있다. 그러나 기존 연구에서도 계산의 편의성을 위하여 슬랩의 점성도를 하나의 상수 또는 2-3층 구조로 나누어 섭입 현상을 모사하였다(Funiciello et al., 2003; Capitanio et al., 2007; OzBench et al., 2008). 추후에 우리의 모형에서 얻어진 슬랩의 응력 분포와 온도 구조를 바탕으로 다양한 크리프 모형을 도입하면 현실적인 점성도 구조와 강도 곡선(strength envelop)을 제시할 수 있을 것이다.

맨틀 지온구배를 도입한 열적 경계조건에 의해 계산된 슬랩 내부의 혀 구조는 차가운 슬랩 코어의 준안정 감락석의 영역을 파악할 수 있게 하였다. 준안정 감락석 쐐기의 존재는 온도가 낮고 단단한 슬랩 코어의 높은 기계적 강도를 의미하며 이는 슬랩 내부의 높은 응력이 집중시킬 수 있다. 따라서 상변이대 깊이의 섭입대에서 발생하는 심발지진에 대한 가능성을 높일 수 있다.

섭입대 수치 모사에서 탄성은 해양 암석권의 섭입이 시작할 때 휨 모멘트 극복과 슬랩 내부의 응력 분포의 결정, 저장된 탄성에너지의 방출이 지진으로 나타남에 있어 고려해야 할 필수 요소이다. 우리의 모형은 영 계수를 고정하여 실험을 진행하였지만, 이를 실제 암석의 영 계수 범위안에서 변화시키면서 실험을 수행하면 탄성이 섭입대 역학에 어떤 역할로써 작용하는지 구속할 수 있을 것이다. 상판 위에 생성되는 후열도 분지의 형성과 섭입하면서 상판과의 마찰력에 의해 발생하는 거대 지진(megathrust earthquake), 섭입하는 슬랩과 상판 사이의 결합도(interplate coupling) 등을 고려해보면 섭입대 역학에서 상판의역할은 중요하다(van Dinther et al., 2010). 마찰을 직접 도입하는 것은 상판과 섭입하는 슬랩을 유체로 가정하는 오일러리안 접근법으로는 구현하기에 부적합하여 라그랑지안 틀 내에서 접근해야 한다. 본 연구의 수치 모형을 바탕으로 상판을 도입하면 마찰이 섭입에 미치는 영향에 대하여 연구할 수 있을 것이다. 모든 기하 구조를 고체로 가정하는 수치 모사 기법에서는 점성 유체의 맨틀은 실재할 수 없다. 이러한 약점을 극복하기 위해 섭입하는 슬랩의 경계에 점성 대쉬팟을 설치하여 하강하는 슬랩에 반대 방향으로 작용하는 맨틀의 항력(drag force)을 근사하여 경계조건으로 도입함으로써 섭입 과정 중에 슬랩과 맨틀 사이의 상호작용을 재현하였다. 그러나, 슬랩과 맨틀의 경계로부터 먼 거리에서 발생하는 구석 유동(Lee and King, 2009), 플룸과 슬랩의 상호작용(Druken et al., 2014; Chang et al., 2016)과 같은 슬랩에 의해 유도되는 맨틀의 대류 현상을 포함하지 못하였다. 본 연구에서 개발한 자유 섭입 모형을 바탕으로 해양판을 고체, 맨틀을 유체로 가정하는 유체-구조 상호작용(Morra and Regenauer-Lieb, 2006)방법을 도입하면 거대 맨틀 쐐기(Big Mantle Wedge, Zhao et al., 2009)와 더불어 플룸이 슬랩을 뚫고 올라오는 현상(Obrebski et al., 2010; Tang et al., 2014), 플룸에 의한 섭입 각도 변화 및 슬랩의 기울어짐(Mériaux et al., 2016), 플룸에 의한 섭입 시작(Ueda et al., 2008; Burov and Cloetingh, 2010; Gerya et al., 2015)과 같은 플룸과 슬랩 사이의 상호작용을 파악할 수 있을 것으로 기대된다.

5. 결 론

우리는 유한요소법에 기반한 상용 수치해석 소프트웨어인 콤솔 멀티피직스를 활용하여 라그랑지안 접근법으로 점탄성 자유 섭입 모형에 대한 벤치마크 연구를 수행하였다. 이를 통하여 기존 연구와 유사한 슬랩의 거동과 해양 암석권의 밀도에 따른 슬랩의 하강 속도를 도출했다. 우리의 자유 섭입 모형은 지구상 섭입대에서 흔하게 관찰되는 해구 후퇴 현상을 3.2-6.0 cm/yr의 후퇴 속도로 보여주었으며, 이는 현실적인 관측 결과와 일치한다는 것을 확인하였다. 뿐만 아니라, 660 km 상변이대에서 점성도 증가(10-100배)로 인해 슬랩이 가라앉지 못하여 발생하는 스테그넌트 슬랩 현상을 모사 하였고, 이는 일본 해구를 통해 섭입해 들어오는 태평양 슬랩의 형태와 유사하다.

해령 밀기 힘을 모사하기 위하여 운동학적 등속도 경계조건을 부가하여 수치 모사를 진행하였다. 속도를 부가하지 않았을 때와 비교하여 마리아나 해구로 섭입하는 태평양 슬랩의 형상과 같이 고각으로 섭입하는 것을 보여주었으며, 슬랩이 660 km 상변이대에 최초로 닿을 때까지 걸리는 시간이 단축됐다(예, 해령 밀기 속도를 4 cm/yr로 준 경우 2.7 Myrs 감소). 해령 밀기 속도를 충분히 부가 하였을 때 해구 전진 현상이 나타났으며 해양 암석권과 상부 맨틀의 밀도차, 660 km 상변이대에서 점성도 증가 정도가 해구 전진에 미치는 영향에 대해 파악하였다. 밀도차가 작을 수록 (즉, 슬랩 당김 힘이 작을 수록), 660 km 상변이대에서 점성도 증가 정도가 클수록 해구 전진 속도가 빠르게 나타났다.

점성도 증가의 정도, 해령 밀기 속도 및 슬랩의 밀도 조건에 따라 다양한 값이 나타나긴 하지만 슬랩의 하부 힌지에서 응력 집중이 조건에 따라 0.5-3.5 GPa의 강도로 발생했다. 암석 실험에서 제시한 바에 따르면, 우리의 수치 모형에서 나타난 강한 응력 값은 지구 심부의 고온/고압 환경일지라도 슬랩이 크게 변형하는 경우 슬랩 물질의 파괴 및 심발지진이 발생할 수 있음을 지시한다.

맨틀 지온구배에 기반한 열적 경계조건과 나이에 따른 해양판의 온도 분포를 열적 초기조건으로 적용하여 얻은 슬랩 내부의 시계열적 온도 구조는 차가운 슬랩 코어인 혀 구조를 재현하였다. 이를 바탕으로 슬랩 내부의 온도와 역학적 거동 사이의 상관 관계를 파악하기 위하여 열적 매개변수를 측정하였고, 해양 암석권의 밀도가 높을수록 슬랩의 하강 속도가 빠르기 때문에 열적 매개변수 값이 증가하고, 맨틀로부터의 열이 슬랩 내부로 전도되는 시간이 짧아 더 차가운 온도 구조를 얻을 수 있었다. 부시네스크 근사를 적용하여 온도에 의존하는 슬랩의 밀도 구조를 얻었고, 해양 암석권의 밀도를 상수로 적용했을 때보다 슬랩의 하강 속도가 감소하는 것을 확인하였다.

복잡한 물성(탄성, 점성, 소성)과 두 가지 이상의 지배 방정식(예, 열과 유체 혹은 열과 구조 해석)을 동시에 반영할 수 있을 뿐 아니라 병렬 연산의 지원과 해석형 프로그래밍 언어와의 연동을 통하여 수치 해석의 자동화가 가능한 콤솔 멀티피직스를 활용하면 향후 보다 복잡한 수치 계산이 가능할 것이며, 우리가 재현한 자유 섭입 모형은 섭입대 역학에 기반한 다양한 현상(슬랩의 하강속도, 해구 전진 및 후퇴, 응력 분포, 열적 구조 등)을 파악했다는 점에서 범용성있게 활용될 것으로 판단된다.

Acknowledgments

본 논문의 심사에 참여 해주신 익명의 두 분 위원님과 편집위원님께 깊이 감사드립니다. 이 연구는 한국연구재단의 대통령POST_DOC펠로우쉽(NRF-2014R1A6A3A04055841), 2017년도 강원대학교 대학회계 학술연구조성비, 행정안전부 ‘지진방재분야 전문인력 양성’ 사업과 기상청 기상·지진See-At기술개발사업(KMI2018-02110)의 지원으로 수행되었습니다.

References

• Andrews, E.R., and Billen, M.I., (2009), Rheologic controls on the dynamics of slab detachment, Tectonophysics, 464, p60-69. [https://doi.org/10.1016/j.tecto.2007.09.004]
• Barker, D.H. N., Henrys, S., Tontini, F.C., Barnes, P.M., Bassett, D., Todd, E., and Wallace, L., (2018), Geophysical constraints on the relationship between seamount subduction, slow slip and tremor at the north Hikurangi subduction zone, New Zealand, Geophysical Research Letters, 45, p12804-12813. [https://doi.org/10.1029/2018gl080259]
• Běhounková, M., and Cížková, H., (2008), Long-wavelength character of subducted slabs in the lower mantle, Earth and Planetary Science Letters, 275, p43-53. [https://doi.org/10.1016/j.epsl.2008.07.059]
• Beuchert, M.J., and Podladchikov, Y.Y., (2010), Viscoelastic mantle convection and lithospheric stresses, Geophysical Journal International, 183, p35-63. [https://doi.org/10.1111/j.1365-246x.2010.04708.x]
• Billen, M.I., and Hirth, G., (2007), Rheologic controls on slab dynamics, Geochemistry, Geophysics, Geosystems, 8. [https://doi.org/10.1029/2007gc001597]
• Buiter, S.J.H., Govers, R., and Wortel, M.J.R., (2001), A modelling study of vertical surface displacements at convergent plate margins, Geophysical Journal International, 147, p415-427. [https://doi.org/10.1046/j.1365-246x.2001.00545.x]
• Burov, E., and Cloetingh, S., (2010), Plume-like upper mantle instabilities drive subduction initiation, Geophysical Research Letters, 37, p1-6. [https://doi.org/10.1029/2009gl041535]
• Capitanio, F.A., Morra, G., and Goes, S., (2007), Dynamic models of downgoing plate-buoyancy driven subduction: Subduction motions and energy dissipation, Earth and Planetary Science Letters, 262, p284-297. [https://doi.org/10.1016/j.epsl.2007.07.039]
• Capitanio, F.A., Morra, G., and Goes, S., (2009), Dynamics of plate bending at the trench and slab-plate coupling, Geochemistry, Geophysics, Geosystems, 10, pQ04002. [https://doi.org/10.1029/2008gc002348]
• Chang, S.J., Ferreira, A.M.G., and Faccenda, M., (2016), Upper- and mid-mantle interaction between the Samoan plume and the Tonga-Kermadec slabs, Nature Communications, 7, p10799. [https://doi.org/10.1038/ncomms10799]
• Chen, J., and King, S.D., (1998), The influence of temperature and depth dependent viscosity on geoid and topography profiles from models of mantle convection, Physics of the Earth and Planetary Interiors, 106, p75-92. [https://doi.org/10.1016/s0031-9201(97)00110-6]
• Choi, H., Hong, T., He, X., and Baag, C., (2012), Seismic evidence for reverse activation of a paleo-rifting system in the East Sea (Sea of Japan), Tectonophysics, p572-573, p123-133.
• Cížková, H., and Bina, C.R., (2015), Geodynamics of trench advance: Insights from a Philippine-Sea-style geometry, Earth and Planetary Science Letters, 430, p408-415.
• Cížková, H., van den Berg, A.P., Spakman, M., and Matyska, C., (2012), The viscosity of Earth’s lower mantle inferred from sinking speed of subducted lithosphere, Physics of the Earth and Planetary Interiors, p200-201, p56-52.
• Čížková, H., van Hunen, J., and van den Berg, A., (2007), Stress distribution within subducting slabs and their deformation in the transition zone, Physics of the Earth and Planetary Interiors, 161, p202-214. [https://doi.org/10.1016/j.pepi.2007.02.002]
• Druken, K.A., Kincaid, C., Griffiths, R.W., Stegman, D.R., and Hart, S.R., (2014), Plume-slab interaction : The Samoa-Tonga system, Physics of the Earth and Planetary Interiors, 232, p1-14. [https://doi.org/10.1016/j.pepi.2014.03.003]
• Elsasser, W.M., (1971), Sea-floor spreading as thermal convection, Journal of Geophysical Research, 76, p1101-1112. [https://doi.org/10.1029/jb076i005p01101]
• Faccenna, C., Di Giuseppe, E., Funiciello, F., Lallemand, S., and Van Hunen, J., (2009), Control of seafloor aging on the migration of the Izu-Bonin-Mariana trench, Earth and Planetary Science Letters, 288, p386-398. [https://doi.org/10.1016/j.epsl.2009.09.042]
• Farrington, R.J., Moresi, L.-N., and Capitanio, F.A., (2014), The role of viscoelasticity in subducting plates, Geochemistry, Geophysics, Geosystems, 15, p4291-4304. [https://doi.org/10.1002/2014gc005507]
• Forsyth, D.W., and Uyeda, S., (1975), On the relative importance of the driving forces of plate motion, Geophysical Journal International, 43, p163-200. [https://doi.org/10.1111/j.1365-246x.1975.tb00631.x]
• Forte, A.M., and Mitrovica, J.X., (1996), New inferences of mantle viscosity from joint inversion of long-wavelength mantle convection and post-glacial rebound data, Geophysical Research Letters, 23, p1147-1150. [https://doi.org/10.1029/96gl00964]
• Fourel, L., Goes, S., and Morra, G., (2014), The role of elasticity in slab bending, Geochemistry, Geophysics, Geosystems, 15, p4507-4525. [https://doi.org/10.1002/2014gc005535]
• Funiciello, F., Morra, G., Regenauer-Lieb, K., and Giardini, D., (2003), Dynamics of retreating slabs: 1. Insights from two-dimensional numerical experiments, Journal of Geophysical Research: Solid Earth, 108, p1-17. [https://doi.org/10.1029/2001jb000898]
• Gerya, T.V., Stern, R.J., Baes, M., Sobolev, S.V., and Whattam, S.A., (2015), Plate tectonics on the Earth triggered by plume-induced subduction initiation, Nature, 527, p221-225. [https://doi.org/10.1038/nature15752]
• Gerya, T.V., and Yuen, D.A., (2003), Rayleigh-Taylor instabilities from hydration and melting propel “cold plumes” at subduction zones, Earth and Planetary Science Letters, 212, p47-62. [https://doi.org/10.1016/s0012-821x(03)00265-6]
• Giunchi, C., Sabadini, R., Boschi, E., and Gasperini, P., (1996), Dynamic models of subduction: geophysical and geological evidence in the Tyrrhenian Sea, Geophysical Journal International, 126, p555-578. [https://doi.org/10.1111/j.1365-246x.1996.tb05310.x]
• Goes, S., Capitanio, F.A., Morra, G., Seton, M., and Giardini, D., (2013), Signatures of downgoing plate-buoyancy driven subduction in Cenozoic plate motions, Physics of the Earth and Planetary Interiors, 184, p1-13.
• Green Ii, H.W., and Burnley, P.C., (1989), A new self-organizing mechanism for deep focus earthquakes, Nature, 341, p733-737. [https://doi.org/10.1038/341733a0]
• Hall, P.S., (2012), On the thermal evolution of the mantle wedge at subduction zones, Physics of the Earth and Planetary Interiors, p198-199, p9-27.
• Hebert, L.B., and Montesi, L.G.J., (2013), Hydration adjacent to a deeply subducting slab: The roles of nominally anhydrous minerals and migrating fluids, Journal of Geophysical Research: Solid Earth, 118, p1-18. [https://doi.org/10.1002/2013jb010497]
• Houseman, G.A., and Gubbins, D., (1997), Deformation of subducted oceanic lithosphere, Geophysical Journal International, 131, p535-551. [https://doi.org/10.1111/j.1365-246x.1997.tb06598.x]
• Hsu, S.K., Yeh, Y.C., Sibuet, J.C., Doo, W.B., and Tsai, C.H., (2013), A mega-splay fault system and tsunami hazard in the southern Ryukyu subduction zone, Earth and Planetary Science Letters, 362, p99-107. [https://doi.org/10.1016/j.epsl.2012.11.053]
• Huang, J., and Zhao, D., (2006), High-resolution mantle tomography of China and surrounding regions, Journal of Geophysical Research: Solid Earth, 111, pB09305. [https://doi.org/10.1029/2005jb004066]
• Karato, S., (2010), Tectonophysics Rheology of the deep upper mantle and its implications for the preservation of the continental roots: A review, Tectonophysics, 481, p82-98.
• Keum, J.Y., (2019), Evoultion process of stress distribution and thermal structure in oceanic subducting plate using 2D finite element method: in a framework of viscoelastic structural mechanics, Master D. thesis, Kangwon National University, Chuncheon, p23, (in Korean with English abstract).
• Kim, G.B., Yoon, S.H., Kim, S.S., and So, B.D., (2018), Transition from buckling to subduction on strike-slip continental margins: Evidence from the East Sea Japan Sea, Geology, 467, p603-606.
• Kim, H.-J., Lee, G.H., Choi, D.-L., Jou, H.-T., Li, Z., Zheng, Y., Kim, G.-Y., Lee, S.-H., and Kwon, Y.K., (2015), Back-arc rifting in the Korea Plateau in the East Sea (Japan Sea) and the separation of the southwestern Japan Arc from the Korean margin, Tectonophysics, 638, p147-157. [https://doi.org/10.1016/j.tecto.2014.11.003]
• Kirby, S.H., Stein, S., Okal, E.A., and Rubie, D.C., (1996), Metastable mantle phase transformations and deep earthquakes in subducting oceanic lithosphere, Reviews of Geophysics, 342, p261-306.
• Kneller, E.A., van Keken, P.E., Karato, S., and Park, J., (2005), B-type olivine fabric in the mantle wedge: Insights from high-resolution non-Newtonian subduction zone models, Earth and Planetary Science Letters, 237, p781-797. [https://doi.org/10.1016/j.epsl.2005.06.049]
• Lee, C., (2013), A Benchmark for 2-Dimensional Incompressible and Compressible Mantle Convection Using COMSOL Multiphysics^®, Journal of the Geological Society of Korea, 49, p245-265, (in Korean with English abstract).
• Lee, C., and King, S.D., (2009), Effect of mantle compressibility on the thermal and flow structures of the subduction zones, Geochemistry, Geophysics, Geosystems, 10, pQ01006. [https://doi.org/10.1029/2008GC002151]
• Lee, C., and Wada, I., (2017), Clustering of arc volcanoes caused by temperature perturbations in the back-arc mantle, Nature Communications, 8, p15753. [https://doi.org/10.1038/ncomms15753]
• Mériaux, C.A., Mériaux, A.S., Schellart, W.P., Duarte, J.C., Duarte, S.S., and Chen, Z., (2016), Mantle plumes in the vicinity of subduction zones, Earth and Planetary Science Letters, 454, p166-177. [https://doi.org/10.1016/j.epsl.2016.09.001]
• Mitrovica, J.X., (1996), Haskell [1935] revisited, Journal of Geophysical Research, 101, p555-569. [https://doi.org/10.1029/95jb03208]
• Mitrovica, J.X., and Forte, A.M., (2004), A new inference of mantle viscosity based upon joint inversion of convection and glacial isostatic adjustment data, Earth and Planetary Science Letters, 225, p177-189.
• Morra, G., and Regenauer-Lieb, K., (2006), A coupled solid-fluid method for modelling subduction, Philosophical Magazine, 86, p3307-3323. [https://doi.org/10.1080/14786430500256359]
• Obayashi, M., Fukao, Y., and Yoshimitsu, J., (2017), Unusually deep Bonin earthquake of 30 May 2015: A precursory signal to slab penetration?, Earth and Planetary Science Letters, 459, p221-226. [https://doi.org/10.1016/j.epsl.2016.11.019]
• Obrebski, M., Allen, R.M., Xue, M., and Hung, S.H., (2010), Slab-plume interaction beneath the Pacific Northwest, Geophysical Research Letters, 37, p1-6. [https://doi.org/10.1029/2010gl043489]
• OzBench, M., Regenauer-lieb, K., Stegman, D.R., Morra, G., Farrington, R., Hale, A., May, A., Freeman, J., Bourgouin, R., Hale, A., Mühlhaus, H., and Moresi, L., (2008), A model comparison study of large-scale mantle - lithosphere dynamics driven by subduction, Physics of the Earth and Planetary Interiors, 171, p224-234. [https://doi.org/10.1016/j.pepi.2008.08.011]
• Panasyuk, S.V., and Hager, B.H., (2000), Inversion for mantle viscosity profiles constrained by dynamic topography and the geoid, and their estimated errors, Geophysical Journal International, 143, p821-836. [https://doi.org/10.1046/j.0956-540x.2000.01286.x]
• Ranalli, G., and Murphy, C., (1987), Rheological stratification of the lithosphere, Tectonophysics, 132, p281-295. [https://doi.org/10.1016/0040-1951(87)90348-9]
• Royden, L.H., and Husson, L., (2006), Trench motion, slab geometry and viscous stresses in subduction systems, Geophysical Journal International, 167, p881-905. [https://doi.org/10.1111/j.1365-246x.2006.03079.x]
• Schellart, W.P., (2004), Quantifying the net slab pull force as a driving mechanism for plate tectonics, Geophysical Research Letters, 31, p10-14. [https://doi.org/10.1029/2004gl019528]
• Schellart, W.P., Freeman, J., Stegman, D.R., Moresi, L., and May, D., (2007), Evolution and diversity of subduction zones controlled by slab width, Nature, 446, p1-4. [https://doi.org/10.1038/nature05615]
• Schellart, W.P., Stegman, D.R., and Freeman, J., (2008), Global trench migration velocities and slab migration induced upper mantle volume fluxes: Constraints to find an Earth reference frame based on minimizing viscous dissipation, Earth-Science Reviews, 88, p118-144. [https://doi.org/10.1016/j.earscirev.2008.01.005]
• Schmalholz, S.M., Duretz, T., Schenker, F.L., and Podladchikov, Y.Y., (2014), Kinematics and dynamics of tectonic nappes: 2-D numerical modelling and implications for high and ultra-high pressure tectonism in the Western Alps, Tectonophysics, 631, p160-175. [https://doi.org/10.1016/j.tecto.2014.05.018]
• Schmeling, H., Monz, R., and Rubie, D.C., (1999), The influence of olivine metastability on the dynamics of subduction, Earth and Planetary Science Letters, 165, p55-66. [https://doi.org/10.1016/s0012-821x(98)00249-0]
• Schubnel, A., Brunet, F., Hilairet, N., Gasc, J., Wang, Y., and Green, H.W., (2013), Deep-Focus Earthquake Analogs Recorded at High Pressure and Temperature in the Laboratory, Science, 341, p1377-1380. [https://doi.org/10.1126/science.1240206]
• So, B.D., and Capitanio, F.A., (2017), The effect of plate-scale rheology and plate interactions on intraplate seismicity, Earth and Planetary Science Letters, 478, p121-131. [https://doi.org/10.1016/j.epsl.2017.08.020]
• So, B.D., and Yuen, D.A., (2015), Generation of tectonic over-pressure inside subducting oceanic lithosphere involving phase-loop of olivine-wadsleyite transition, Earth and Planetary Science Letters, 413, p59-69. [https://doi.org/10.1016/j.epsl.2014.12.048]
• Song, T.R.A., and Simons, M., (2003), Large trench-parallel gravity variations predict seismogenic behavior in subduction zones, Science, 301, p630-633. [https://doi.org/10.1126/science.1085557]
• Stegman, D.R., Farrington, R., Capitanio, F.A., and Schellart, W.P., (2010), A regime diagram for subduction styles from 3-D numerical models of free subduction, Tectonophysics, 483, p29-45. [https://doi.org/10.1016/j.tecto.2009.08.041]
• Syracuse, E.M., Keken, P.E. Van, Abers, G.A., Suetsugu, D., Editor, G., and Bina, C., (2010), The global range of subduction zone thermal models, Physics of the Earth and Planetary Interiors, 183, p73-90. [https://doi.org/10.1016/j.pepi.2010.02.004]
• Tang, Y., Obayashi, M., Niu, F., Grand, S.P., Chen, Y.J., Kawakatsu, H., Tanaka, S., and Ning, J., (2014), Changbaishan volcanism in northeast China linked to subduction- induced mantle upwelling, Nature Geoscience, 76, p470-475.
• Tingle, T.N., Green II, H.W., Scholtz, C.H., and Koczynski, T.A., (1993), The rheology of faults triggered by the olivine-spinel transformation in Mg₂GeO₄ and its implications for the mechanism of deep focus earthquakes, Journal of Structural Geology, 15, p1249-1256. [https://doi.org/10.1016/0191-8141(93)90167-9]
• Toth, J., and Gurnis, M., (1998), Dynamics of subduction initiation at preexisting fault zones, Journal of Geophysical Research: Solid Earth, 103, p18053-18067. [https://doi.org/10.1029/98jb01076]
• Turcotte, D., and Schubert, G., (2002), Geodynamics, Cambridge Univ. Press, Cambridge, 2nd Ed.
• Ueda, K., Gerya, T., and Sobolev, S.V., (2008), Subduction initiation by thermal-chemical plumes: Numerical studies, Physics of the Earth and Planetary Interiors, 171, p296-312. [https://doi.org/10.1016/j.pepi.2008.06.032]
• van Dinther, Y., Morra, G., Funiciello, F., and Faccenna, C., (2010), Role of the overriding plate in the subduction process: Insights from numerical models, Tectonophysics, 484, p74-86. [https://doi.org/10.1016/j.tecto.2009.08.038]
• Vergnolle, M., Pollitz, F., and Calais, E., (2003), Constraints on the viscosity of the continental crust and mantle from GPS measurements and postseismic deformation models in western Mongolia, Journal of Geophysical Research: Solid Earth, 108, p1-14. [https://doi.org/10.1029/2002jb002374]
• Wang, K., Hu, Y., and He, J., (2012), Deformation cycles of subduction earthquakes in a viscoelastic Earth, Nature, 484, p327-332. [https://doi.org/10.1038/nature11032]
• Ye, L., Lay, T., Zhan, Z., Kanamori, H., and Hao, J.L., (2016), The isolated ~680 km deep 30 May 2015 MW 7.9 Ogasawara (Bonin) Islands earthquake, Earth and Planetary Science Letters, 433, p169-179. [https://doi.org/10.1016/j.epsl.2015.10.049]
• Zhao, D., Tian, Y., Lei, J., Liu, L., and Zheng, S., (2009), Seismic image and origin of the Changbai intraplate volcano in East Asia: Role of big mantle wedge above the stagnant Pacific slab, Physics of the Earth and Planetary Interiors, 173, p197-206. [https://doi.org/10.1016/j.pepi.2008.11.009]