Aims & Scope

2차원 평면 변형률(plane strain) 상태를 가정한 자유 섭입 모형을 모사하기 위하여 사용된 지배방정식은 운동량 보존 방정식(식 1)과 에너지 보존 방정식(식 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에 포함하였다.

맨틀과 암석권은 거시적 관점에서 연속체로 가정할 수 있는데, 연속체 역학에서 물체의 변위 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)으로 정의한다.

u_i는 i방향의 변위를 의미한다. 변형률 텐서를 시간에 대해 미분하면 변형률 속도(strain-rate) 텐서 ε˙ij가 도출된다. 식 (4)에서 ε˙ij와 ε˙ij'는 각각 변형률 속도 텐서와 편향 변형률 속도 텐서이다.

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

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

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

대부분의 기존 섭입 모사 연구는 섭입 각도 및 속도를 부가하는 운동학적 경계조건(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)를 사용하여 구한 수치해를 바탕으로 섭입 과정을 모사할 수 있었다.

그림 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 (V_p = 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 V_p are fixed at 80 kg/m³ and , 4 cm/yr repectively. (c) Time evolution of trench position in the cases with varying values of V_p and a fixed value of Δρ = 80 kg/m³.

그림 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., V_p = 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 V_p show larger stress (e.g., ~3.5 GPa with V_p = 4 cm/yr and η660 = 100 × η_M).

수치 모형에 적용한 열적 경계조건 및 초기조건을 통해 섭입 과정 동안 진화하는 슬랩 내부의 온도 구조를 얻을 수 있었다. 그림 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.