콤솔 멀티피직스를 활용한 2차원 수치 섭입모델링 벤치마크

1. 서 론

섭입(subduction)은 중앙해령(mid-ocean ridge)에서 생성된 해양판이 맨틀로 침강하여 소멸하는 작용으로 지구의 물질 및 에너지 순환에서 중요한 역할을 담당한다. 또한 섭입은 지진과 호화산(arc volcanism) 활동처럼 우리 인류의 삶과 밀접한 지질 현상을 발생시키므로 오랫동안 연구되었다. 그러나, 직접 관찰이 불가능한 맨틀 내부 섭입해양판의 거동과진화에 대해서는 여전히 잘 알려져 있지 않다.

호화산은 맨틀로 침강한 섭입해양판에서 탈수된 물이 맨틀 쐐기(mantle wedge)로 유입되어 발생하는 유체유입용융(flux melting)에 의한 것으로 여겨진다(Gaetani and Grove, 1998). 그러나, 섭입해양판의 탈수 및 유체유입용융을 이해하기 위해서는 섭입해양판에 포함된 함수 광물의 탈수(dehydration), 탈수된 물의 이동 그리고 그 물에 의한 맨틀의 부분용융(partial melting) 등을 정량적으로 이해하는 것이 필요하다. 이는 고체 및 유체의 거동 및 에너지의 정량적 계산을 요구하기 때문에 광물학 및 암석학 등의 지질과학과 유체 역학 및 열역학 등의 관련 이공학을 이용하여 연구되어 왔다. 유체의 거동과 에너지 교환을 해석하는 지배방정식(governing equations)은 편미분방정식(partial differential equations, PDE)의 형태로 표현되므로 일부의 경우를 제외하고는 해석해(analytic solution)가 존재하지 않는다. 그러므로 수치해석(numerical analysis)을 이용한 수치모델링(numerical modeling)이 과거 수십 년 동안 각광받아 왔으며 특히 컴퓨터 기술의 발달과 더불어 발전해 왔다.

이처럼 섭입 작용의 정량적 연구를 위하여 널리 활용되어 온 수치모델링 연구는 서로 다른 연구진들에 의해 다양한 공개(open) 및 상용(commercial) 소프트웨어를 이용하여 수행되어 왔다. 그러므로, 각기 다른 소프트웨어를 사용함으로써 얻어진 계산 결과의 신빙성을 평가하기 위한 벤치마크(benchmark)가 다수 수행되었다(e.g., Blankenbach et al., 1989; Travis et al., 1990; van Keken et al., 2008; King et al., 2010). van Keken et al. (2008)은 상수(constant) 및 포행(creep) 맨틀 점성(viscosity)을 사용하여 맨틀 쐐기(mantle wedge)에서의 거동을 해석학적(analytic) 및 수치해석적(numerical) 방법으로 계산하였다. 이 벤치마크 연구는 섭입 수치모델링을 수행하고자 하는 연구자들이 자신들의 연구 방법을 벤치마크하기 위하여 널리 활용되고 있다.

이 연구에서는 van Keken et al. (2008)에서 보고된 벤치마크 결과들을 유한요소법(finite element method) 기반 상용 소프트웨어인 콤솔 멀티피직스(COMSOL Multiphysics®)를 사용하여 재현하였다. 이미 van Keken et al. (2008)에서 우즈 홀 해양연구소(Woods Hole Oceanographic Institute, WHOI) 소속 Mark Behn 박사에 의해 콤솔 멀티피직스(ver. 3.2b)를 사용한 벤치마크 결과가 제시되었다. 그러나, Behn 박사에 의해 사용된 수치모델링의 구성에 대한 설명이 지나치게 간략하여 다른 연구진들이 이를 사용하여 벤치마크를 수행하는데 어려움이 있다. 그러므로 이 연구에서는 보다 자세한 설명을 부연하여 콤솔 멀티 피직스를 사용한 벤치마크를 보다 쉽게 수행할 수 있도록 하는데 목적을 두었다. 이 논문에서는 먼저 섭입 수치모델링의 설정 방법에 대해 서술하며, 그 다음 van Keken et al. (2008)에서 제시된 벤치마크 결과와 비교하였다.

2. 실험 방법

2.1 지배방정식 및 수치모델

벤치마크를 수행하기 위하여 van Keken et al. (2008)에서 사용된 지배방정식을 사용하였다. 사용된 지배방정식은 맨틀의 거동이 점성도가 매우 높은 비압축성 유체(incompressible fluid)의 유동으로 근사 될 수 있다는 가정에서 도입되었으며 맨틀 쐐기 내부에서의 온도 차에 의한 부력(buoyancy)은 무시되었다. 지배방정식은 다음의 연속(continuity), 모멘텀(momentum) 그리고 에너지 방정식으로 표현된다.

$\[ 0=\nabla \bullet \stackrel{⃑}{V} \]$

1) 연속 방정식

$\[ 0=-\nabla P+\nabla \bullet \left(2\eta \dot{ϵ}\right) \]$

2) 모멘텀 방정식

$\[ {\rho }_c{C}_p\left(\frac{{\partial }_T}{{\partial }_t}+\stackrel{⃑}{V}\bullet \nabla T\right)=\nabla \bullet \left(k\nabla T\right) \]$

3) 에너지 방정식

$$ \stackrel{⃑}{V}$$는 맨틀 속도(mantle velocity), P는 압력(pressure), η는 맨틀 점성(viscosity), $$ \dot{ϵ}$$는 변형율 텐서(strain rate tensor), ρ_c는 밀도(density), C_p는 등압열용량(specific heat at constant pressure), T는 온도(temperature), t는 시간(time), k는 열전도도(thermal conductivity)이다. 맨틀 내에서의 부력과 방사성 동위원소의 붕괴 및 광물의 상전이 등을 무시하였으며 지배방정식에 사용된 값들은 표 1에 제시되었다.

Table 1.

Parameters and reference values.

벤치마크에서 사용된 모델 구조(model geometry)는 2차원으로 구현되었으며 van Keken et al. (2008)에서 사용된 깊이 600 km, 너비 660 km로 설정되었다(그림 1a). 섭입해양판은 45° 각도로 5 cm/y의 속도로 섭입하며 섭입해양판 하부맨틀은 섭입해양판과 같은 속도와 방향으로 거동하는 것으로 가정되었다. 섭입해양판과 맨틀 쐐기와의 점성 결합(viscous coupling)에 의하여 맨틀 쐐기에서 구석 유동(corner flow)이 발생한다. 두께 50 km의 상부판(overriding plate)은 고정되어 맨틀 쐐기의 구석 유동에 의해 영향 받지 않는 것으로 설정되었다. 100 Ma의 연령에 해당하는 섭입해양판의 온도 구조는 모델 구조의 왼쪽 수직 경계에 반-공간 냉각 모델(half-space cooling model)을 이용하여 구현하였다(Turcotte and Schubert, 2002) (그림 1b). 섭입해양판과 모델 구조가 만나는 하부 경계는 유출 경계(outflow boundary)로 설정하여 열에너지가 맨틀과 함께 유출되도록 설정되었다. 상부판의 표면 온도는 0℃로 설정되었으며 상부판 하부 표면은 노-슬립(no slip) 조건을 부여하였다. 상부판과 모델 구조의 오른쪽 수직 경계에 부여한 온도 구조는 0에서 1,300℃까지 부여되었으며 맨틀 쐐기와 모델 구조의 경계는 200 km 깊이까지 1,300℃로 고정되었다. 경계 조건에 부여된 온도와 맨틀 쐐기 온도와의 간섭을 최소화시키기 위하여 200 km부터 모델 구조의 하부 경계까지 유출 경계를 설정하였다.

Fig. 1.

Model geometry, boundary condition and model mesh used in the benchmark study modified from van Keken et al. (2008). a) Model geometry consists of kinematically subducted slab and underlying mantle at a rate of 5 cm/y, rigid overriding plate of which thickness is 50 km and mantle wedge of which behavior is analytically or dynamically calculated. For the viscosity of the mantle wedge, constant viscosity, diffusion and dislocation creep are used for the cases 1a-c, 2a and 2b, respectively. b) Boundary condition applied to the model domain. The temperature profile along the left vertical boundary is prescribed with half-space cooling model. A constant velocity is prescribed on the interface between the subducting slab and mantle wedge to solve the Stokes equation. The top and bottom boundaries of the overriding plate are fixed and prescribed with the no-slip boundary, respectively. The right vertical boundaries of the overriding plate and mantle wedge are prescribed with the linear temperature gradient and mantle potential temperature by a depth of 200 km. The other boundary of the mantle wedge is prescribed with the open boundary. c) Model mesh consisting of triangular elements of which sizes range from 1.8 to 18 km. Mesh refinement is applied to the subducting slab and the interfaces between the subducting slab/overriding plate and mantle wedge.

맨틀의 점성은 온도, 압력, 물의 유무 그리고 변형률 등에 크게 달라지며 실험을 통해 제안된 아레니우스 방정식(Arrhenius equation)으로 주로 표현된다(e.g., Karato and Wu, 1993; Hirth and Kohlstedt, 1995). 이 연구에서는 맨틀의 점성을 상수 점성과 아레니우스 방정식으로 표현되는 포행 점성으로 구분하여 사용하였다. 각 벤치마크에 다른 점성이 사용되었기 때문에 아래에서 필요에 따라 점성 계산식을 부연하였다.

수치계산을 위하여 유한요소법 기반 상용 소프트웨어인 콤솔 멀티피직스(COMSOL Multiphysics, ver. 5.3)가 사용되었다. 연속방정식은 이 연구에서 가정된 비압축성 유체 가정을 만족하므로 직접 계산되지 않는다. 모멘텀 방정식을 계산하기 위하여 우리는 유체거동 해석을 위해 사용되는 콤솔 멀티피직스의 내장모듈(module)인 포행 유동(creeping flow) 모듈을 사용하였으며, 식 2에서 제시된 것처럼 매우 느리게 거동하는 맨틀의 관성항(좌변항)을 무시한 스토크스 방정식(Stokes equation)을 사용하였다. 모멘텀 방정식의 이산화(discretization)는 P2+P1을 사용하였으며 이는 맨틀의 속도는 노드(node)에서 연속인 2차 다항 함수(second-order polynomial)로, 압력은 선형 함수로 근사된다는 것을 지시한다. 에너지 방정식을 계산하기 위하여 우리는 내장 모듈인 유체에서의 열전달(heat transfer in fluid) 모듈을 사용하였다. 에너지 방정식의 이산화는 선형(linear)을 사용하였으며 이는 노드(node)에서의 온도 값을 선형 근사하여 요소 내부의 온도 값을 계산한다는 것을 지시한다. 유한요소법에 기반한 수치해석을 수행하기 위하여 모델 구조를 삼각형 요소(triangular element)로 구성한 메쉬(mesh)를 생성하였다(그림 1c). 수치해석의 정확도를 높이기 위하여 온도 및 맨틀 거동의 변화가 크게 나타나는 지역인 섭입해양판과 맨틀 쐐기와의 경계는 메쉬 세밀화(mesh refinement)를 적용하여 최소 약 1.8 km의 삼각형 요소를 사용하였으며 섭입해양판과 거리가 먼 지역은 최대 약 18 km의 삼각형 요소를 사용하였다. 메쉬를 구성하는 전체 삼각형 요소의 갯수는 33,734개이다. 수치계산을 위하여 콤솔 멀티피직스에 내장된 시간의존성 솔버(time-dependent solver)를 내장된 일반화된 알파(generalizedalpha) 법을 적용하여 사용하였다. 실험 결과의 수렴성을 높이기 위하여 내장된 직접 솔버(direct solver)인 MUMPS (MUltifrontal Massively Parallel Sparse direct Solver)를 사용하였다. 거동 및 온도의 수렴을 보장하기 위하여 충분한 시간인 200 Ma의 시간을 부여하여 실험 결과가 정상 상태(steady-state)에 도달 하도록 하였다.

2.2 벤치마크 사례

벤치마크를 위해 수행된 모델 사례(case)는 van Keken et al. (2008)에서 수행된 것과 같으며 사례에 따라 필요한 설명을 부연하였다. 사례 1은 상수 맨틀 점성을 사용한 실험들로 구성되었으며 사례 2는 포행 점성을 사용한 실험들로 구성되었다.

2.2.1 사례 1a

일부 초기 섭입 수치모델링 연구(e.g., McKenzie and Sclater, 1968; Oxburgh and Turcotte, 1968)에서 맨틀 쐐기의 구석 유동을 해석하기 위하여 해석 구석유동식(analytic cornerflow solution, 부록 참조)을 사용하였다. 이 해석 구석 유동식은 임의의 섭입 각도와 속도가 주어지면 맨틀 쐐기 내의 유체의 속도 분포를 계산할 수 있으므로 초기 섭입대 연구에 다수 활용되었으나, 맨틀 점성이 상수이여야 하며 섭입해양판이 모델 구조에서 제시된 것처럼 단순한 직선 구조로 정의되어야 한다는 한계를 지니고 있다. 그럼에도 불구하고 벤치마크를 위한 목적에 유용하게 활용될 수 있으므로 이 사례에서는 해석 구석 유동식을 이용하여 맨틀 쐐기(mantle wedge)와 모델 경계에서의 유동이 계산되었다(그림 2). 그러므로 에너지 방정식만을 수치해석적으로 계산한다.

Fig. 2.

Schematic diagram used to calculate the cornerflow solution of the mantle wedge. The subducting rate is U and the x- and y-component of the velocity in the mantle wedge are expressed as u and v, respectively. The θ is the angle between the subducting slab and fixed overriding plate.

2.2.2 사례 1b

이 사례에서는 맨틀 쐐기의 유동을 사례 1a에서 사용된 해석 구석 유동식을 사용하지 않고 모멘텀 방정식을 수치해석을 통하여 계산하였다. 그러나 맨틀 쐐기와 모델 경계에서의 맨틀의 유입 및 유출은 해석 구석 유동식을 이용하여 계산되었다. van Keken et al. (2008)에서 토의된 것처럼 일정한 속도로 유입되는 섭입해양판과 고정된 상부판 그리고 맨틀 쐐기가 만나는 점에서 압력 특이점(pressure singularity)이 발생한다. 압력 특이점이 맨틀 거동 및 온도 계산에 미치는 영향을 감소시키기 위하여 van Keken et al. (2008)에서 사용된 것과 같이 압력 특이점에서 섭입해양판 표면을 따라 5 km에 걸친 구간에 섭입 속도를 0에서 5 cm/y로 선형적으로 증가시켜 부여하였다.

2.2.3 사례 1c

이 사례에서는 사례 1b와 동일한 조건이 사용되었으며 다만 맨틀 쐐기와 모델 경계에서의 맨틀의 유입 및 유출이 맨틀 쐐기의 거동에 따라 자유롭게 이루어지도록 압력을 0으로 부여하였다.

2.2.4 사례 2a

과거 연구들은 맨틀의 거동이 크게 확산(diffusion) 포행과 전위(dislocation) 포행으로 표현될 수 있다는 것을 보였다(e.g., Karato and Wu, 1993; Hirth and Kohlstedt, 1995). 이 사례에서는 나머지 모델 조건을 사례 1c에서 쓰인 것과 같은 것을 사용하고 다만 맨틀의 점성을 다음의 확산 포행식을 사용하여 계산하였다.

$\[ {\eta }_{dif}\left(\mathrm{T}\right)A_{dif}\mathrm{e}\mathrm{x}\mathrm{p}\left(\frac{E_{dif}}{RT}\right) \]$

4) 확산 포행식

η_dif는 확산 포행 점성(viscosity of diffusion creep, Pa·s), A_dif는 확산 포행 선지수(pre-exponent constant, Pa·s), E_dif는 확산 포행 활성화 에너지(activation energy, J·mol^-1), R은 기체 상수(gas constant, J·mol^-1·K^-1) 그리고 T는 온도(temperature, K)이다. 실제 모델에서 사용된 점성은 다음의 식을 이용하여 지나치게 높은 점성이 수치계산에서 사용되는 것을 방지하였다.

$\[ {\eta }_{dif,effective}={\left(\frac{1}{{\eta }_{diff}}+\frac{1}{{\eta }_{max}}\right)}^{-1} \]$

5) 실질 확산 포행식

η_{dif,effective}는 실질 확산 포행 점성(effective viscosity of diffusion creep, Pa·s) 그리고 η_max는 최대 점성(maximum viscosity, Pa·s)이다. 확산 점성 포행 계산을 위한 수치 값들은 표 1에 제시되었다.

2.2.5 사례 2b

이 사례에서는 나머지 모델 조건을 사례 1c에서 쓰인 것과 같은 것을 사용하고 다만 맨틀의 점성을 다음의 전위 포행식을 사용하여 계산하였다.

$\[ {\eta }_{dis}\left(T,\dot{ϵ}\right)=A_{dis}\mathrm{e}\mathrm{x}\mathrm{p}\left(\frac{E_{dis}}{nRT}\right){\dot{ϵ}}^{\left(1-n\right)/n} \]$

6) 전위 포행식

η_dis는 전위 포행 점성(viscosity of dislocation creep, Pa·s), A_dis는 전위 포행 선지수(pre-exponent constant, Pa·s^1/n), E_dis는 전위 포행 활성화 에너지(activation energy, J·mol^-1), n은 전위 포행 멱법칙 지수(powerlaw exponent for dislocation creep), 그리고 $$ \dot{ϵ}$$은 변형률(strain rate, s^-1)이다. 사례 2a에 적용된 것과 마찬가지로 실제 모델에 사용된 점성은 다음의 식을 이용하여 지나치게 높은 점성이 사용되는 것을 방지하였다.

$\[ {\eta }_{dis,effective}={\left(\frac{1}{{\eta }_{dis}}+\frac{1}{{\eta }_{max}}\right)}^{-1} \]$

7) 실질 전위 포행식

η_{dis,effective}는 실질 전위 포행 점성(effective viscosity of dislocation creep, Pa·s)이다. 전위 점성 포행 계산을 위한 수치 값들은 표 1에 제시되었다.

3. 실험 결과

3.1 속도 및 열구조: 상수 맨틀 점성

사례 1의 경우는 상수 맨틀 점성을 사용하였으며 맨틀 쐐기의 거동과 경계 조건을 변화시켜 그 결과 들을 벤치마크 하였다. 모든 사례 1의 실험에서 섭입해양판과 맨틀 쐐기와의 점성 결합에 의한 구석 유동이 발생하는 것을 관찰할 수 있으며 온도 경계조건이 부여된 맨틀 쐐기와 모델 경계에서 유입된 뜨거운 맨틀이 맨틀 쐐기의 구석으로 유입되어 차가운 섭입해양판을 가열시키는 현상을 관찰할 수 있다(그림 3). 전반적인 맨틀의 거동 및 온도 분포는 사례 1에서 변화된 조건과는 크게 상관 없이 유사한 분포를 보이는 것을 확인할 수 있다. 다만 맨틀 쐐기와 모델 경계에서 맨틀 거동과 온도 분포가 달라지는 것을 확인 할 수 있다.

Fig. 3.

Velocity and temperature distributions obtained from the cases 1a, b and c. a) Velocity and temperature distributions from the case 1a. Left panel: velocity distribution described with red arrows of which sizes are proportional to velocity. Right panel: temperature distribution described with colors from blue (cold) to red (hot). b) The same with a except for the case 1b. c) The same with a except for the case 1c.

사례 1에서 계산된 맨틀 거동과 온도 분포를 van Keken et al. (2008)에서 보고된 벤치마크 결과들과 비교하기 위하여 첫째, (60, 60 km) 위치에서의 온도값을 측정하였으며, 둘째, 0 km에서 210 km 깊이의 섭입해양판과 맨틀 쐐기의 경계면에 6 km 간격으로 분포하는 36개 지점에서 온도의 L2 놈(norm)을 다음의 방정식을 이용하여 구하였으며,

$\[ ‖T_{\mathrm{s}\mathrm{l}\mathrm{a}\mathrm{b}}‖=\sqrt{\frac{{\sum }_{i=1}^{36}{T}_{ii}^2}{36}} \]$

8) 섭입해양판 표면에서의 온도의 L2 놈

셋째, 맨틀 쐐기의 구석에 너비 66 km인 대칭 직각삼각형 형태의 공간을 설정하고 그 공간에 6 km간격으로 균등하게 분포된 78개 지점에서 측정된 온도의 L2 놈을 다음의 방정식을 이용하여 구하였다.

$\[ ‖T_{\mathrm{w}\mathrm{e}\mathrm{d}\mathrm{g}\mathrm{e}}‖=\sqrt{\frac{{\sum }_{\mathrm{i}=10}^{21}\sum _{\mathrm{j}=10}^{\mathrm{i}}{T}_{\mathrm{i}\mathrm{j}}^2}{78}} \]$

9) 맨틀 쐐기에서의 온도의 L2 놈

van Keken et al. (2008)에서는 참여 연구 그룹들이 계산한 온도 값들을 가로 111, 세로 101 개의 구조화된 메쉬(structured mesh)에 투영한 값들을 비교에 활용하였으나, 이 연구에서는 주어진 좌표값에서 직접 측정된 온도 값들을 비교에 활용하였다. 표 2는 van Keken et al. (2008)에서 보고된 온도 값들과 함께 이 연구에서 측정된 값들을 나타내었다. 사례 1 실험들에서 측정된 온도 값들은 보고된 온도값들과 유사한 값들을 보여주며 이전 버전의 콤솔멀티피직스(ver. 3.2b)를 사용하여 계산한 온도 값(WHOI)과 유한요소법에 기반한 공개 소프트웨어인 Sepran 소프트웨어(Cuvelier et al., 1986)를 사용한 온도 값(UM)들과 전반적으로 잘 일치하는 것을 확인할 수 있었다. 이는 콤솔 멀티피직스를 사용하여 재현한 섭입대 벤치마크가 성공적으로 수행되었음을 지시한다. van Keken et al. (2008)에서 보고된 맨틀 쐐기 구석에서의 메쉬 해상도(mesh resolution)에 따른 온도 값들을 도시한 그림에 우리가 측정한 온도 값들을 도시하였다(그림 4). 콤솔을 사용한 WHOI의 실험 결과는 메쉬 해상도에 거의 영향을 받지 않았다. 그러므로 이 연구에서는 하나의 메쉬 해상도를 사용하였고 그림의 패널 각각에 하나의 측정 값들을 도시하였다.

Table 2.

The chart of temperature distributions for cases 1a, 1b and 1c.

Fig. 4.

Slab surface temperature at a depth of 60 km, L2 norm of the slab temperatures along the slab surface and L2 norm of the temperatures in the triangle part of the mantle wedge varying mesh resolution near the corner of the mantle wedge from the cases 1a, b and c in van Keken et al. (2008). In this study, the benchmark experiments are conducted using one mesh resolution and the calculated temperatures are depicted with the purple diamonds.

3.2 속도 및 열구조: 포행 맨틀 점성

상수 맨틀 점성을 사용한 실험과 함께 포행 맨틀 점성식을 사용한 사례 2 실험에서 계산된 맨틀 거동과 온도 분포를 나타내었다(그림 5). 확산과 전위 포행식을 사용한 실험 모두 유사한 맨틀 거동과 온도 분포를 보여준다. 상수 맨틀 점성을 사용한 실험 결과와 대조적으로 상부판의 배호에 두꺼운 열암석권(thermal lithosphere)의 발달이 관찰된다. 그리고 섭입해양판과 맨틀 쐐기 사이의 강한 점성 결합에 의하여 섭입해양판 표면에 발달하는 열경계층(thermal boundary layer)의 두께가 얇으며, 이로 인하여 섭입해양판의 온도는 상수 맨틀 점성을 사용한 실험과 대조하여 높은 온도에 빠르게 도달하게 된다.

Fig. 5.

Velocity and temperature distributions obtained from the cases 2a and 2b. The notations are the same with those in Figure 2.

van Keken et al. (2008)에서 보고된 온도 값들과 비교하기 위하여 사례 1에서 벤치마크를 위하여 사용한 방법으로 온도 값들을 계산하였다(표 3). 사례 1에서와 마찬가지로 사례 2의 경우에도 과거 보고된 온도 값들과 잘 일치하는 것을 확인할 수 있었다. 메쉬 해상도에 의한 영향을 나타낸 van Keken et al. (2008)에서 보고된 결과와 이 연구에서 측정된 온도 값들을 같이 도시하였으며 사례 1과 마찬가지로 메쉬 해상도가 온도에 미치는 영향이 작은 것을 확인 할 수 있었다(그림 6).

Table 3.

The chart of temperature distributions for cases 2a and 2b.

Fig. 6.

Slab surface temperature at a depth of 60 km, L2 norm of the slab temperatures along the slab surface and L2 norm of the temperatures in the triangle part of the mantle wedge varying mesh resolution near the corner of the mantle wedge from the cases 2a and b in van Keken et al. (2008). The notations are the same with those in Figure 3.

4. 결론 및 제언

이 연구에서는 유한요소법 기반 상용 소프트웨어인 콤솔 멀티피직스를 이용하여 과거 보고된 섭입대 수치모델링 벤치마크를 재현하였다. 실험 결과는 상수 맨틀 점성과 포행 맨틀 점성을 사용한 실험 모두에서 과거 연구에서 보고된 결과들과 잘 일치하는 것을 확인할 수 있었다. 특히 콤솔 멀티피직스는 메쉬 해상도의 차이에 의한 온도 값의 차이가 크지 않다는 것을 잘 보여주었으며 이는 콤솔 멀티피직스를 이용한 실험이 상대적으로 낮은 메쉬 해상도를 사용하여도 높은 계산 정확도를 보장할 수 있다는 것을 잘 지시한다. 이는 메쉬 해상도의 증가에 따라 기하 급수적으로 증가하는 메모리 및 계산 시간을 절약할 수 있다는 것을 뜻하므로 콤솔 멀티피직스의 활용 가능성을 잘 보여준다. 또한 콤솔 멀티피직스가 지원하는 사용자 친화적인 그래픽 사용자 인터페이스(graphic user interface)와 유용한 모듈의 탑재는 수치모델링 연구를 위하여 필요한 시간과 노력을 크게 절약할 수 있게 해준다는 점에서 수치모델링을 보다 쉽게 지질과학 연구에 활용할 수 있는 가능성을 잘보여준다.

5. 부 록

본문에서 언급한대로 맨틀 쐐기의 유동을 설명하기 위하여 사용한 해석 구석 유동식($$ \stackrel{⃑}{U}$$)은 맨틀의 점성이 상수이며, 섭입해양판이 직선 형태로 가정(그림 2)되었을 때 다음과 같이 표현된다(Turcotte and Schubert, 2002).

u와 v 는 해석 구석 유동식의 x및 y 방향 성분이며, A, B, C, D는 경계 조건에 의해 결정되는 상수들이다. 이 때 맨틀 쐐기 상부에 존재하는 상부판이 거동하지 않는다고 가정하면 표면(y=0)에서의 해석 구석 유동식의 경계 조건을 다음과 같이 정의할 수 있다.

위의 경계 조건을 식 A2와 A3에 대입하여 풀면 아래와 같이 상수가 정의된다.

A = 0, B = -C

맨틀 쐐기 내부(y>0)에서의 속도 성분 u와 v가 섭입 각도(θ)에 관한 식으로 표현된다.

U는 섭입해양판의 섭입 속력(subducting rate)이며 해석 구석 유동식의 절대값(U=||$$ \stackrel{⃑}{U}$$||)이다. 표면에서의 경계 조건을 해석 구석 유동식의 x, y 방향 성분식 A2와 A3에 대입하면 다음의 임의의 섭입 속도와 각도에 대한 일반식으로 정의된다.

u = Csin²θ - D(θ + sinθcosθ)
v = C(θ - sinθcosθ) - Dsin²θ
C = (θsinθ/(θ² - sin²θ))U
D = ((sinθ - θcosθ)/(θ² - sin²θ))U

Acknowledgments

본 논문을 심사해주시고 수준 향상에 큰 도움을 주신 두 분의 심사위원님과 편집위원장님께 감사드립니다. 본 연구는 한국연구재단의 지원으로 수행되었습니다(과제 번호: 2016R1D1A1B03930906).

References

• Blankenbach, B., Busse, F., Christensen, U., Cserepes, L., Gunkel, D., Hansen, U., Harder, H., Jarvis, G., Koch, M., Marquart, G., Moore, D., Olson, P., Schmeling, H., and Schnaubelt, T., (1989), A benchmark comparison for mantle convection codes, Geophysical Journal International, 98, p23-38. [https://doi.org/10.1111/j.1365-246x.1989.tb05511.x]
• Cuvelier, C., Segal, A., and Van Steenhoven, A.A., (1986), Finite Element Methods and the Navier–Stokes Equations, 22, Reidel, Norwell, Mass.
• Gaetani, G.A., and Grove, T.L., (1998), The influence of water on melting of mantle peridotite, Contribution to Mineralogy and Petrology, 131, p323-346. [https://doi.org/10.1007/s004100050396]
• Hirth, G., and Kohlstedt, D., (1995), Experimental constraints on the dynamics of the partially molten upper mantle 2. Deformation in the dislocation creep regime, Journal of Geophysical Research, 100, p15441-15449. [https://doi.org/10.1029/95jb01292]
• Karato, S.-I., and Wu, P., (1993), Rheology of the upper mantle: a synthesis, Science, 260, p771-778. [https://doi.org/10.1126/science.260.5109.771]
• King, S.D., Lee, C., van Keken, P.E., Leng, W., Zhong, S., Tan, E., Tosi, N., and Kameyama, M.C., (2010), A community benchmark for 2-D Cartesian compressible convection in the Earth's mantle, Geophysical Journal International, 180, p73-87. [https://doi.org/10.1111/j.1365-246x.2009.04413.x]
• McKenzie, D.P., and Sclater, J.C., (1968), Heat-flow inside the island arcs of the northwest Pacific, Journal of Geophysical Research, 73, p3173-3179.
• Oxburgh, E.R., and Turcotte, D.L., (1968), Problem of high heat flow and volcanism associated with zones of descending mantle convective flow, Nature, 218, p1041-1043. [https://doi.org/10.1038/2181041a0]
• Travis, B.J., Anderson, C., Baumgardner, J., Gable, C.W., Hager, B.H., O'Connell, R.J., Olson, P., Raefsky, A., and Schubert, G., (1990), A benchmark comparison of numerical-methods for infinite Prandtl number thermal-convection in two-dimensional Cartesian geometry, Geophysical and Astrophysical Fluid Dynamics, 55, p137-160. [https://doi.org/10.1080/03091929008204111]
• Turcotte, D., and Schubert, G., (2002), Geodynamics, Cambridge Univ. Press, Cambridge, 2nd Ed.
• van Keken, P.E., Currie, C., King, S.D., Behn, M.D., Cagnioncle, A., He, J., Katz, R.F., Lin, S.-C., Parmentier, E.M., Spiegelman, M., and Wang, K., (2008), A community benchmark for subduction zone modeling, Physics of The Earth and Planetary Interiors, 171, p187-197. [https://doi.org/10.1016/j.pepi.2008.04.015]