The Geological Society of Korea
[ Article ]
Journal of the Geological Society of Korea - Vol. 49, No. 2, pp.245-265
ISSN: 0435-4036 (Print) 2288-7377 (Online)
Print publication date Apr 2013
Received 13 Dec 2012 Reviewed 20 Dec 2012 Accepted 14 Jan 2013
DOI: https://doi.org/10.14770/jgsk.2013.49.2.245

A Benchmark for 2-Dimensional Incompressible and Compressible Mantle Convection Using COMSOL Multiphysics®

LeeChangyeol
Faculty of Earth and Environmental Sciences, Chonnam National University, 333 Yongbong-no, Bukgu, Gwangju 500-757, Republic of Korea
콤솔 멀티피직스를 이용한 2차원 비압축성 및 압축성 맨틀 대류의 벤치마크

Correspondence to: ‡ +82-62-530-3451, E-mail: changyeol.lee@gmail.com

In this study, a benchmark of mantle convection is conducted using the COMSOL Multiphysics®, a commercial program based on the finite element method (FEM). For the benchmark, simplified mantle convection models are formulated by including incompressibility or compressibility of the mantle. The benchmarked governing equations include anelastic liquid approximation (ALA), truncated anelastic liquid approximation (TALA), extended Boussinesq approximation (EBA) and Boussinesq approximation (BA) by varying dissipation and Rayleigh numbers. Since the COMSOL Multiphysics® has user-friendly modeling interfaces and environments, non-specialists in computational geodynamics can use the program for their own geodynamic modeling. In addition, the program allows easy modifications of rheology and material properties as well as governing equations for the user’s specified modeling purposes. The results in this study are well consistent with the published results. The COMSOL Multiphysics® is promising for diverse modelings of computational geodynamics.

초록

본 연구는 최근 전산 지구동력학 연구에서 널리 활용되고 있는 유한요소법 기반의 상용 프로그램인 콤솔 멀티피직스를 이용한 벤치마크를 수행하였다. 벤치마크는 맨틀의 비압축성 및 압축성을 반영한 단순한 맨틀 대류를 모델링 함으로써 이루어졌다. 벤치마크에 사용된 지배방정식은 비탄성 유체 근사, 절단 비탄성 유체 근사, 확장 부시네스크 근사 그리고 부시네스크 근사이며 레일리히 수와 소산 수를 변화시켜 구해진 각각의 수렴해들을 비교하였다. 콤솔 멀티피직스는 전산 지구동력학의 비전문가도 쉽게 모델을 수립 및 사용할 수 있는 편리성을 내포하고 있으며 사용자가 자유롭게 물성과 지배방정식을 변형시킬 수 있기 때문에 이 연구에서 유용하게 활용되었다. 실험 결과는 기존의 벤치마크 결과와 전반적으로 잘 합치한다. 편리성과 정확성을 내포하는 콤솔 멀티피직스는 향후 전산 지구동력학 연구에 널리 활용될 수 있을 것으로 기대된다.

Keywords:

COMSOL Multiphysics®, compressible mantle, mantle convection, numerical model, benchmark, 콤솔 멀티피직스, 압축성 맨틀, 맨틀 대류, 수치 모델, 벤치마크

1. 서 론

20세기 이공학의 발전 역사에서 가장 중요한 발견을 꼽으라면 컴퓨터의 등장을 들 수 있다. 미육군에서 탄도 계산을 위한 최초의 범용 컴퓨터인 ENIAC (Electronic Numerical Integrator And Computer) (Goldstine, 1972)의 개발 이후 컴퓨터 기술은 크게 발달하였으며 현재는 거의 모든 이공학 분야에서 컴퓨터가 널리 사용되고 있다. 특히, 과거에는 불가능했던 막대한 양의 연산을 컴퓨터가 가능하게 해줌에 따라 수치해석 분야도 획기적인 발전을 이루었으며 이를 통하여 대부분 해석해(analytic solution)가 존재하지 않는 복잡한 편미분방정식의 수치해(numerical solution)를 찾는 것이 가능해졌다. 그리고 최근 더욱 강력한 연산 능력을 제공하기 위한 수퍼컴퓨터(supercomputer)의 등장은 불과 20여년 전에는 상상도 하지 못했던 대규모의 계산을 가능하게 해주었다.

이에 따라 20세기 후반 컴퓨터와 수치해석을 이용하여 이공학 분야의 복잡한 문제들을 해결하기 위한 노력들이 많이 이루어져 왔는데 지질과학(geological sciences)도 그 예외가 아니다. 컴퓨터 자체가 생소했던 1960년대에도 컴퓨터를 이용한 수치모델링(numerical modeling)이 선구적인 지질과학 연구에 일부 사용되었다(e.g., McKenzie and Sclater, 1968; Oxburgh and Turcotte, 1968). 컴퓨터가 보다 저렴해지고 널리 보급됨에 따라 컴퓨터와 수치해석을 이용한 지질과학 연구는 더욱 빠르게 확산되었는데 그 가운데 전산 지구동력학(computational geodynamics)은 컴퓨터와 수치해석을 이용하여 판의 이동과 변형, 맨틀의 대류를 연구하는 분야로써 미국, 일본 그리고 유럽 선진국에서 중요한 지질과학의 연구분야로 자리잡았다. 현재에도 맨틀 대류 및 섭입 작용 등과 같은 다양한 전산 지구동력학 연구들이 수행되고 있다(Sleep, 2006 and references therein; Billen, 2008).

그러나 전산 지구동력학 연구는 다양한 분야의 관련 지식들이 요구되는 융합학문이다. 우선 Fortran이나 C와 같은 컴퓨터 프로그래밍에 대한 지식이 필요하며, 선형대수학 및 미분방정식 등 수치해석에 필수적으로 요구되는 수학 지식이 필요하다. 더 나아가서는 판과 맨틀의 거동을 지배하는 물성(rheology)이 온도, 압력 그리고 화학 조성에 의존하기 때문에 지질과학 및 재료과학적 지식도 함께 필요하다. 이러한 특성들 때문에 전산 지구동력학 연구를 위하여 필요한 프로그램을 개발 및 이용하는 것조차 큰 노력과 시간이 소요되며, 결과적으로 지질과학 전공자들이 전산 지구동력학 연구를 수행하며 다양한 지식을 습득하는 데 걸림돌로 작용하고 있다. 이를 극복하기 위하여 선진국 지질과학계에서는 전문 컴퓨터 프로그래머, 응용수학자 및 지질과학자 등으로 특성화된 연구 그룹을 조직하여 전산 지구동력학 프로그램 개발과 이를 이용한 연구에 힘쓰고 있다(e.g., Computational Infrastructure for Geodynamics, www.geodynamic.org and Earthbyte Group, www. earthbyte.org).

다른 나라들과 달리, 현재 우리나라의 전산 지구동력학 연구는 걸음마 수준에 머무르고 있어 선진국처럼 전문화된 연구 그룹을 구성하여 프로그램을 개발 및 이용하는 종합적인 지구동력학 연구는 전무한 실정이다. 그러므로 전산 지구동력학을 이용한 지질과학 연구를 활발하게 국내 지질과학계에 소개하여 저변확대를 꾀하는 것이 먼저 필요하다. 그러나 전산 지구동력학 연구자들의 수가 절대적으로 부족한 국내환경을 고려하면 여러 전문가의 많은 노력과 시간이 필요한 전산 지구동력학 연구를 계획하는 단계조차도 어렵다. 그러므로 전산 지구동력학을 전문적으로 전공한 사람뿐만이 아니라 비전공자들까지도 보다 편리하게 연구를 계획하고 수행할 수 있는 환경을 조성하여 인식의 확대를 꾀하는 것이 무엇보다도 시급하다.

이러한 현실에서 필자는 유한요소법(Finite Element Method, FEM)에 기반을 둔 상용 수치해석 프로그램인 콤솔사(COMSOL, www.comsol.com)사의 콤솔 멀티피직스(COMSOL Multiphysics®)를 이용한 전산 지구동력학 연구 수행을 하나의 대안으로 제시한다. 콤솔 멀티피직스는 이미 이공학 분야에서 널리 활용되고 있는 상용 프로그램으로써 1, 2 및 3차원의 정교한 기체, 유체 및 고체의 거동을 정량적으로 계산할 수 있다. 특히 그래픽 유저 인터페이스(Graphic User Interface, GUI)를 기본적으로 지원함으로써 사용자의 데스크탑 모니터에서 마우스와 키보드를 이용하여 쉽게 모델을 수립 및 해석할 수 있게 해준다. 이뿐만이 아니라 막대한 계산을 요구하는 모델링의 효율적 계산을 위하여 수퍼컴퓨터를 이용한 병렬연산(parallel computation)을 지원함으로써 보다 정교하고 복잡한 수치모델링을 수행하는 것도 가능하다. 실제 콤솔 멀티피직스는 섭입 작용 및 중앙해령에서 일어나는 판구조 운동 연구에 다수 활용되어 전산 지구동력학 연구에서의 활용 가능성은 매우 높다고 할 수 있다(e.g., Carminati and Petricca, 2010; Montési et al., 2011).

이처럼 높은 활용 가능성을 가진 콤솔 멀티피직스를 이용한 본격적인 전산 지구동력학 연구에 앞서 필자는 가장 기본적으로 요구되는 벤치마크(benchmark)를 수행하고자 한다. 섭입대 모델링에 적용한 콤솔 멀티피직스의 벤치마크는 이미 수행되었다(van Keken et al., 2008). 그러나 맨틀의 압축성(compressibility)과 점성소산(viscous dissipation)을 고려한 콤솔 멀티피직스의 벤치마크는 아직까지 보고되지 않았다. 그러므로 필자는 본 연구에서 2차원 비압축성 및 압축성 맨틀 대류 모델의 선행 연구 결과(King et al., 2010)와 콤솔 멀티피직스를 사용하여 얻은 실험 결과들을 비교하고자 하였다.


2. 실험 방법

2.1 콤솔 멀티피직스(COMSOL Multiphysics®)의 간략한 소개

콤솔 멀티피직스는 1996년 스웨덴에서 설립된 콤솔(COMSOL, www.comsol.com)사에서 1998년 소개한 유한요소법에 기반한 상용 소프트웨어로써 2012년 12월 현재 4.3a 버전(version)까지 출시되었다. 콤솔 멀티피직스는 1, 2 및 3차원의 다양한 수치 모델링을 그래픽 유저 인터페이스를 이용하여 구현할 수 있도록 지원하며 전기, 기계, 유체, 연소, 구조해석 및 화학 반응에까지 이공학에서 필요한 다양한 수치 해석을 손쉽게 구현할 수 있게 한다. 이미 프로그램 내부에 다양한 모델링 예제들이 탑재되어 있어 사용자가 목적에 맞게 예제를 변형하여 사용할 수 있다. 그리고 상수나 변수의 형태로 원하는 재료의 물성을 입력할 수 있기 때문에 다양한 물성을 가진 수치모델 수립이 가능하다.

콤솔 멀티피직스가 탑재하고 있는 모델링 예제도 방대하며 유용하지만 또 하나의 중요한 장점은 사용자가 목적에 맞게 상미분 및 편미분 방정식을 자유롭게 모델에 삽입할 수 있다는 것이다. 상미분 및 편미분 방정식의 자유로운 삽입은 사용자의 목적에 특성화된 수치 모델 수립을 가능하게 한다. 그리고 이미 콤솔 멀티피직스에 내장된 상미분 및 편미분 방정식도 간단한 수식 기호의 조합으로 변경이 가능하기 때문에 상미분 및 편미분 방정식에 대한 기본적인 수학적 지식이 있으면 누구나 수식을 삽입 및 수정이 가능하다.

콤솔 멀티피직스는 수치 모델링을 통하여 도출된 결과들에 대한 후처리(post-processing) 기능을 제공하여 다른 소프트웨어의 도움 없이도 해석이 가능하다. 예를 들어 유체 거동 모델의 경우, 거동을 해석하는데 중요한 자료인 온도, 압력, 속도 등을 그래픽을 통하여 도시하는 것이 가능하며 각 수치 계산 값들을 외부 파일로 출력할 수 있다. 이 뿐만이 아니라 선적분, 면적분 그리고 부피적분을 통해 유체 및 열의 유입 및 유출량, 엔탈피(enthalpy), 총에너지 등 다양한 이공학적으로 유용한 값들을 손쉽게 계산할 수 있다. 실제로 이러한 기능들은 이 연구에서 수행된 결과와 기존 벤치마크에서 도출된 결과와의 상호비교를 위하여 유용하게 활용되었다. 콤솔 멀티피직스에 대한 보다 자세한 소개는 콤솔 멀티피직스의 메뉴얼이나 출판된 참고서적(Pryor, 2012)들을 통하여 접할 수 있다.

2.2 비압축성 및 압축성 맨틀 대류 지배방정식

2.2.1 맨틀 거동 지배방정식

판구조론의 성립 이후 섭입 작용과 중앙해령의 확장은 전지구적으로 일어나는 상부 해양판과 하부 맨틀간의 열화학적 순환(thermochemical circulation)으로 설명될 수 있다. 이러한 열화학적 순환의 가장 큰 부분을 담당하는 맨틀 거동은 온도 및 압력에 의존하는 점성도가 매우 높은 유체의 거동으로 근사 될 수 있다(Schubert et al., 2001). 점성도가 매우 높은 유체의 거동은 압축성 층류 유동(compressible laminar flow)으로 표현될 수 있으며 이 유동의 지배방정식(governing equations)은 질량, 모멘텀 그리고 에너지의 보존을 지시하는 편미분 방정식인 연속(continuity), 모멘텀(momentum) 그리고 에너지(energy) 방정식으로 다음과 같이 표현된다(Schubert et al., 2001).

ρ는 밀도(kg/m3), t는 시간(s), 는 속도 벡터(m/s), P는 압력(Pa), 는 중력가속도 벡터(m/s2), τ는 차응력 텐서(deviatoric stress tensor) (Pa), Cp는 등압열용량 (J/kg.K), T는 온도(K), α는 열팽창률(/K), k는 열전도도(J/m.K.s), φ는 단위 부피당 점성소산(viscous dissipation) (J/m3.s), H는 단위 질량당 내부 열에너지 발생률(J/kg.s)이다. 점성소산은 유체의 거동으로 발생하는 전단응력에 의하여 유체의 모멘텀이 열에너지로 전환되는 양을 지시하며 다음과 같이 정의된다.

만약 맨틀이 거동이 없으며 화학적으로 균질하다면(homogeneous) 맨틀의 밀도는 온도와 압력에만 의존하는 함수로써 나타낼 수 있으며 이를 기준 상태(reference state)로 정의한다. 맨틀 거동(대류)은 위치에 따른 밀도 차이에 의하여 발생하며 밀도 차이는 기준 상태에서 벗어난 온도와 압력에 의하여 발생한다. 그러므로 실제 맨틀의 밀도는 기준 상태를 지시하는 함수와 기준 상태에서 벗어난 온도와 압력에 의하여 발생한 밀도의 차이를 선형적으로 합하여 나타낼 수 있다. 이를 수식으로 나타내면 아래와 같이 서술된다.

위 수식에서 , 그리고 는 기준 상태에 해당하는 밀도, 온도 그리고 압력을 각각 지시하며 ρ', T ' 그리고 P'는 기준 상태에서 벗어난 밀도, 온도 그리고 압력을 각각 지시한다. 맨틀의 거동이 발생하지 않을 때 기준 상태 밀도는 Adams-Williamson 상태방정식(Birch, 1952)으로 나타낼 수 있다.

ρ0는 표면에서의 맨틀의 밀도이며 HT는 압축성 유체의 온도 규모 높이(temperature scale height of compressible liquid)이며 Γ는 그루네이센(Grüneisen) 변수(parameter)로써 다음과 같이 각각 정의된다.

Ks는 체적압축성 등엔트로피 탄성률(isentropic modulus of bulk compressibility)이다. 위에서 제시된 지배방정식을 변수인 맨틀의 밀도, 온도 및 압력을 대입하여 풀면 맨틀의 거동을 계산할 수 있으며 다음은 맨틀 거동 지배방정식에서 유도된 여러 가지 맨틀 거동 근사식들에 대하여 설명하고자 한다.

2.2.2 비탄성 유체 근사

(Anelastic Liquid Approximation, ALA)

맨틀 거동 지배방정식은 맨틀이 다음에 제시되는 가정들을 만족하므로 보다 간략하게 서술될 수 있다. (Schubert et al., 2001) 첫째, 맨틀의 속도에 비하여 맨틀 밀도의 변화량이 매우 작으므로 연속방정식의 좌변항은 0으로 근사 될 수 있다. 둘째, 맨틀의 가속도는 매우 작으므로 맨틀은 마하 수(Mach number)가 사실상 0이며 프란틀 수(Prandtl number)가 1022 이상으로 매우 큰 유체에 해당한다. 이 때문에 모멘텀 방정식의 좌변항도 역시 0으로 근사 될 수 있다. 이러한 두 가지 가정을 지배방정식에 적용하면 비탄성 유체 근사(anelastic liquid approximation, ALA)로 일컬어지는 다음의 방정식이 도출된다.

지배방정식을 효율적으로 계산하기 위하여 무차원화(non-dimensionalization)를 다음의 관련식들을 이용하여 수행한다.

ΔT는 맨틀의 상하표면 온도차, d는 맨틀의 두께(m), κ는 열확산률(heat diffusivity, m2/s), μ는 역학점성도(dynamic viscosity, Pa.s), ρ0는 맨틀 표면의 밀도(kg/m3)이다. 무차원화를 수행한 후 표기의 편의를 위하여 '를 제외하고 지배방정식을 서술하면 다음과 같다.

무차원화는 소산 수(dissipation number) Di와 레일리히 수(Rayleigh number) Ra를 도출시키는데 각각 다음과 같이 정의된다.

무차원화는 밀도에 해당하는 기준 상태식에도 적용되므로 기준 상태식은 다음과 같이 재정의된다.

y'는 맨틀의 깊이를 맨틀 두께로 무차원화시킨 값으로써 맨틀 표면에서는 0, 맨틀 바닥에서는 1이다. 이 기준 상태식은 맨틀의 밀도는 깊이에 선형적으로 비례하여 증가함을 지시한다. 맨틀의 밀도 증가는 단열압축(adiabatic compression)에 의한 것으로 근사 될 수 있는데 단열압축에 의하여 생성되는 온도의 기준 상태식은 다음과 같이 유도된다.

T0는 맨틀 표면의 온도이다. 이 식은 밀도와 마찬가지로 맨틀의 온도는 깊이가 증가함에 따라 단열압축에 의하여 증가함을 지시한다.

2.2.3 절단 비탄성 유체 근사

(Truncated Anelastic Liquid Approximation, TALA)

일부 수치 모델링 패키지의 경우 기술적 문제 등의 이유로 역학압력 요소를 모멘텀 방정식에 고려하는데 어려움이 존재한다(King et al., 2010). 이 경우 비탄성 유체 근사의 모멘텀 방정식에서 역학압력 요소( )를 생략하여 절단 비탄성 유체 근사(truncated anelastic liquid approximation, TALA)를 이용한다. 나머지 방정식은 비탄성 유체 근사와 동일하다.

2.2.4 확장 부시네스크 근사

(Extended Boussinesq Approximation, EBA)

확장 부시네스크 근사(extended Boussinesq approximation, EBA)는 절단 비탄성 유체 근사에서 기준 상태의 밀도를 상수로 가정하여 다시 정의한 근사식으로써 지배방정식은 다음과 같이 정의된다.

νy는 y방향으로의 속도 벡터이다. 이 지배방정식은 깊이에 따른 밀도의 증가 효과가 크지 않은 상대적으로 얕은 깊이의 맨틀 거동을 근사 하는데 유용하다.

2.2.5 부시네스크 근사

(Boussinesq Approximation, BA)

부시네스크 근사(Boussinesq Approximation, BA)은 비압축성 유체의 거동을 설명하는 가장 간결화된 지배방정식이다. 이 근사는 확장 부시네스크 근사의 에너지 방정식에서 소산 수를 0으로 가정함으로써 얻어진다. 그러므로 이 가정은 점성 소산에 의한 열의 발생이 크지 않을 것으로 기대되는 맨틀 거동에 적합하며 다음과 같이 서술된다.

위에서 설명된 지배방정식에 대한 보다 자세한 내용은 과거의 연구 문헌을 참고함으로써 얻을 수 있다(Jarvis and McKenzie, 1980; Schubert et al., 2001; King et al., 2010). 전술하였듯이 맨틀의 거동에 가장 가까운 지배방정식은 비탄성 유체 근사이다. 그러나 기술적 한계 혹은 모델의 간결화 등의 이유로 여전히 부시네스크 근사와 같은 간결화된 지배방정식을 이용한 모델이 널리 쓰이고 있다.

2.3 수치 모델

유도된 지배방정식과 기준 상태 방정식을 콤솔 멀티피직스에 삽입하고 2차원 유한요소 모델링을 이용하여 벤치마크를 수행하였다. 벤치마크를 위한 모델은 King et al. (2010)에 제시된 것과 같으며 간단하게 요약하자면 다음과 같다. 모델 범위(domain)는 2차원 정사각형으로써 상부와 하부 벽(wall)은 273K 및 3273K의 온도로 고정되었으며 좌우 벽은 단열(insulated) 환경으로 설정되었다. 그리고 모든 벽의 응력 환경은 자유슬립(free slip) 환경으로 가정되었으며 이는 어떠한 전단 응력도 유체의 거동에 작용하지 않음을 뜻한다. 맨틀 대류는 상부와 하부 벽의 온도차(3000K)에 의하여 발생한다. 이 연구에서 사용된 지배방정식은 무차원화된 수를 사용하기 때문에 실제로 모델링에서 사용되는 값들은 모두 무차원화된 값들이다. 맨틀의 점성은 온도, 압력 그리고 성분에 따라 크게 달라지나(Karato and Wu, 1993; Ranalli, 2001), 벤치마크의 편의를 위하여 점성도는 상수로 가정되었다. 그리고 방사능 동위원소의 붕괴에 의하여 내부 열에너지가 맨틀 내부에서 발생하나 벤치마크의 편의를 위하여 무시되었다(H = 0). 모델링에 사용된 계수 및 상수 값들은 표 1에 서술되었다.

Model parameters.

유한요소법을 이용한 수치해석에서는 모델링 범위를 구성하는 메쉬(mesh)의 효율적 형성 정도에 따라 계산의 정확도가 크게 영향 받는다. 일반적으로 많은 요소(element)를 사용할수록 정확도가 높지만 요소의 양이 증가함에 따라 모델링에 요구되는 계산 시간이 크게 증가하므로 비효율적이다. 콤솔 멀티피직스는 삼각형 및 사각형 등 다양한 요소로 구성된 메쉬의 형성을 기본적으로 지원하고 있다. 이 연구에서는 사용자-조정(user-controlled) 메쉬를 사용하였으며 메쉬의 크기는 콤솔 멀티피직스에서 제공하는 유체 역학 모델링에 적합한 삼각형 형태의 요소들로 전체 도메인을 균질하게 구성하였다. 메쉬 세밀화는 도메인 상부 벽에서 측정되는 대류에 의한 열 에너지 유출량과 전도에 의한 열 에너지 유출량의 비인 넛셀(Nusselt) 수와 상부 벽에서의 대류 속도의 보다 정확한 측정을 위하여 필요하다. 이를 위해 도메인의 상부는 3단계의 메쉬 세밀화(refinement)를 적용하여 보다 작은 메쉬들로 구성되었다(그림 1). 메쉬 세밀화를 통하여 형성된 총 요소의 수는 16,065개이다.

Fig. 1.

The modeling domain formulated in this study using the COMSOL Multiphysics®. The non-dimensionalized length of all the walls of the square box is 1. The uppermost region of the modeling domain consists of very finer elements for the accurate calculation of the Nusselt number at the top wall.

수치 모델링을 위한 변수로써 과거 연구에서 변수로 사용된 레일리히 수를 104, 2×104, 5×104, 105, 2×105, 5×105 그리고 106으로 변화시켜 사용하였다. 그리고 각 레일리히 수에 해당하는 실험에 각각 소산 수를 0.25, 0.50, 1.0, 1.5 그리고 2.0을 사용하여 별도의 실험을 수행하였다. 소산 수가 존재하는 경우 압축단열에 의한 맨틀 온도의 기준 상태식이 각각의 소산 수에 대응하여 각각 존재한다. 벤치마크를 위하여 이 연구에서는 맨틀 온도의 기준 상태식을 맨틀 온도에 선형적으로 더한 것이 맨틀 하부에서 3273K (무차원값: 1)이 되도록 하였다. 즉, 맨틀 거동을 일으키는 알짜 온도차는 총 맨틀 온도에서 맨틀 온도의 기준 상태 값을 제외한 값이 된다.

부시네스크 근사 지배방정식을 사용한 실험의 경우 소산 수와 관련된 항들이 존재하지 않으므로 레일리히 수만 변화시켜 실험을 수행하였다. 실험은 시간의존성 솔버(time-dependent solver)를 사용하여 무차원화된 시간 1.0까지 계산하여 수행하여 실험 결과가 실질적인 수렴에 이른 정상상태(steady-state)에 도달하였다. 일반적으로 맨틀 대류를 발생시키기 위하여 콤솔 멀티피직스는 자체적으로 온도 교란(perturbation)을 초기값으로 부여하며 이는 대부분 단일 대류 세포(convection cell)를 모델링 범위에 형성시키게 한다. 그러나, 레일리히 수와 소산 수가 높은 실험의 경우 단일 대류 세포보다는 두 개 이상의 대류 세포로 수렴하는 경우가 있다. 가급적 단일 대류 세포로 수렴하도록 실험 결과값을 유도하기 위하여 일부 실험들은 낮은 레일리히 수를 이용하여 계산된 온도 분포를 초기 온도 분포로 사용하였으며 이는 기존 벤치마크 연구(King et al., 2010)에서도 사용되었다.


3. 실험 결과

3.1 벤치마크

3.1.1 비탄성 유체 근사

(Anelastic Liquid Approximation, ALA)

먼저 비탄성 유체 근사 지배방정식을 사용하여 수행한 벤치마크에 대하여 레일리히 수와 소산 수를 변화시킨 결과에 대하여 살펴보자. 표 2는 비탄성 유체 근사 지배방정식을 사용하여 수행한 벤치마크 결과를 넛셀 수(Nu #), 대류하는 맨틀의 상부 표면에서의 최대 속도(Vsurf-max), 대류하는 맨틀의 상부 표면에서의 선적분 값(Vsurf), 대류 세포 전체의 평균 온도(<T>), 대류 세포 전체에서 점성소산에 의하여 발생한 열에너지(<φ>), 마지막으로 대류 세포 전체에서 맨틀이 중력에 대하여 한 일(<W>)을 보여준다. 물론 여기서 제시된 값들은 모두 무차원화된 값이다. 각 값들의 계산 방법은 본 연구와 비교 대상이 되는 기존 벤치마크 연구를 참고하길 바란다(King et al., 2010).

Results from the experiments using ALA.Di = 0.25

Di = 0.50

Di = 1.00

Di = 1.50

Di = 2.00

소산 수와 상관없이 모든 실험에서 레일리히 수가 증가하면 대류의 강도가 증가함이 확인되었다(그림 2; 표 2). 레일리히 수가 증가함에 따라 대류의 강도가 증가한다는 것은 이미 널리 알려진 사실이다(e.g., Blankenbach et al., 1989; King et al., 2010). 그림 2는 소산 수가 0.25으로 고정되어 있으며 레일리히 수가 104, 105 그리고 106으로 각각 10배씩 증가하였을 때의 대류 세포의 온도 분포를 나타내고 있다. 이미 알려진 것과 같이 레일리히 수가 증가할수록 대류 세포의 상부 및 하부에 존재하는 열경계층(thermal boundary layer)의 두께가 감소하는 것이 확인되며 이는 넛셀 수의 증가로 잘 표현된다(넛셀 수 4.410, 9.229 그리고 14.867이 레일리히 수 104, 105, 그리고 106을 사용한 실험 결과에 각각 대응한다). 그리고 대류의 강도가 강해짐에 따라 맨틀 대류 속도가 증가하는데 이는 레일리히 수의 증가와 함께 상부 표면에서의 대류의 최고 속도 및 선적분 값이 증가하는 것으로 잘 표현된다(Vsurf-max의 값인 58.065, 252.509 그리고 749.419이 레일리히 수 104, 105, 그리고 106을 사용한 실험 결과에 각각 대응하며, Vsurf의 값인 38.826, 184.406 그리고 540.730이 레일리히 수 104, 105, 그리고 106을 사용한 실험 결과에 각각 대응한다). 점성소산에 의한 열에너지 (<φ>)는 레일리히 수의 증가와 함께 선형적인 증가 관계를 보인다. 그리고 맨틀이 중력에 대하여 한 일(<W>)도 레일리히 수의 증가와 함께 선형적인 증가 관계를 보이는데 점성소산에 의한 열에너지와 거의 비슷한 값을 보인다. Jarvis and McKenzie (1980)는 점성소산에 의한 열에너지와 맨틀이 중력에 대하여 한 일이 정확하게 서로 상쇄한다는 것을 보였는데 본 연구에서 두 일이 서로 잘 상쇄된다. 또한 이 논문 후반부(단락 3.2)에서 제시될 메쉬의 해상도(resolution) 테스트는 본 실험에서 나타난 오차가 메쉬 해상도에 의한 것일 뿐 지배방정식의 오류에 의한 것이 아님을 보인다.

Fig. 2.

Temperature distributions in the modeling domain using the anelastic liquid approximation (ALA). Dissipation number is fixed as 0.25 but Rayleigh numbers are varied as 104, 105 and 106. The single clockwise convection cell is obtained in the experiments using Rayleigh numbers of 104 and 105. The temperature distribution of the experiment using a Rayleigh number of 106 shows count-clockwise convection cell but the temperature distribution is vertically inverted for a comparison among the temperature distributions. The black arrows indicate direction and vigor of convecting mantle. For a better view, the scale factors for black arrows are differently used as 0.0020, 0.0004 and 0.0002 for 2a, 2b and 2c, respectively. As Rayleigh number increases, upwelling and downwelling as well as the thermal boundary layers become thinner. The unit of temperature is K without non-dimensionalization.

레일리히 수가 증가함에 따라 대류의 강도가 증가하는 것과 달리, 소산 수의 증가는 대류의 강도를 약화시킨다(그림 3). 이는 각 벤치마크 결과 값들이 소산 수가 증가함에 따라 감소되는 것으로 잘 나타난다. 예를 들어, 소산 수가 2.0인 경우, 대류의 강도는 급격하게 약화되어 레일리히 수가 104인 경우 매우 약한 대류만이 관찰된다. 넛셀 수가 1에 근접하는 경우 맨틀 내부의 열전달(heat transfer)이 거의 대부분 대류가 아닌 전도에 의해서만 이루어짐을 뜻한다(넛셀 수가 1인 경우 맨틀 대류는 발생하지 않으며, 전도만이 열전달을 가능하게 한다). 이는 대류의 강도와 양의 관계를 가지는 레일리히 수와 달리 소산 수는 음의 관계를 가지는 것을 지시한다.

Fig. 3.

Temperature distributions in the modeling domain using anelastic liquid approximation (ALA). All the Rayleigh numbers are fixed as 2×105 but dissipation number is varied as 0.50, 1.00, 1.50 and 2.00. With increases in dissipation number, vigor of mantle convection is substantially weakened, expressed sluggish mantle convection and layered temperature distribution with depth. The experiments using higher dissipation number tends to develop multiple convection cells such as two or three convection cells. The black arrows indicate direction and vigor of convecting mantle. For a better view, the scale factors for black arrows are differently used as 0.0002, 0.0004, 0.0008 and 0.0016 for 3a, 3b, 3c and 3d, respectively. The black arrows in 3c and 3d indicate two convection cells; the experiments developing three convection cells are not shown here. The experiment using a dissipation number of 2.00 apparently develops strong layered temperature distributions with depth. The unit of temperature is K without non-dimensionalization.

본 실험에서 획득한 결과를 기존의 벤치마크 결과와 비교하였다. 일부의 경우를 제외하고는 기존 벤치마크 결과와 잘 합치하며 특히 King et al. (2010) 연구에서 ConMan 프로그램(King et al., 1990)을 사용한 필자의 VT (Virginia Tech)의 연구 결과와 1% 이내의 상대오차를 보일 정도로 매우 잘 맞는 것을 확인하였다. 그리고 두 개의 대류 세포를 발달시키는 실험의 경우 Citcom 프로그램(Moresi et al., 1996)을 사용한 CU (University of Colorado)의 연구 결과와도 잘 맞는 것을 확인할 수 있었다. 소산 수가 큰 실험에서 상대오차가 다소 커지는 것을 확인할 수 있으나 대부분 1~4% 정도로 작다.

3.1.2 절단 비탄성 유체 근사

(Truncated Anelastic Liquid Approximation, TALA)

여기서는 절단 비탄성 유체 근사 지배방정식을 사용한 실험 결과에 대하여 요약하였다. 절단 비탄성 유체 근사 지배방정식은 비탄성 유체 근사 지배방정식의 모멘텀 방정식에서 역학압력 요소만을 제거한 근사이다. 일반적으로 역학압력 요소는 다른 요소들에 비하여 상대적으로 유체의 거동에 미치는 영향이 작다고 알려져 있다(King et al., 2010). 그 영향을 정량적으로 평가하기 위하여 절단 비탄성 유체 근사 지배방정식을 사용한 실험과 비탄성 유체 근사 지배방정식을 사용한 실험 결과를 비교하였다. 비교가 가능한 결과 값들의 상대오차를 도시한 그림 4에서 알 수 있듯이, 레일리히 수가 작은 경우에 역학압력 요소에 의한 차이가 상대적으로 큼을 알 수 있다. 그러나, 레일리히 수가 증가함에 따라 상대오차는 감소하는 경향을 보인다. 넷셀 수, 점성소산에 의한 열에너지 그리고 맨틀이 중력에 대하여 한 일간의 상대오차 감소는 레일리히 수가 증가할 경우에 특히 두드려지며 레일리히 수가 5×105 이상인 경우 상대오차는 2% 이내로써 매우 작아진다. 지구 맨틀의 레일리히 수가 약 107 정도로 추정되므로(Schubert et al., 2001), 이 실험은 절단 비탄성 유체 근사 지배방정식이 비탄성 유체 근사 지배방정식으로 근사 될 수 있음을 지시한다. 실제로 표 3에서 관찰되듯이 점성소산에 의한 열에너지는 맨틀이 중력에 대하여 한 일에 비하여 상대오차가 0.1~6% 내외이다. 위와 같이 기존 실험 결과(King et al., 2010)와 비교하였을 때 비록 소산 수가 큰 실험에서는 상대오차가 증가하는 경향을 보였으나 대부분 1~4% 정도로 작으며 소산 수 2.0와 작은 레일리히 수를 사용한 실험(예, 레일리히 수: 104)에서만 상대오차가 8~13% 정도로 크게 나타남을 확인할 수 있었다. 그러나 전반적으로 본 연구의 실험 결과는 기존 벤치마크의 결과와 서로 유사함을 확인할 수 있었다.

Fig. 4.

Relative errors between the experiments using anelastic liquid approximation (ALA) and truncated anelastic liquid approximation (TALA) with varying Rayleigh and dissipation numbers. All the parameters except for the mean temperatures of the modeling domain show a clear tendency to decreasing relative errors between ALA and TALA with increases in Rayleigh number. Contrary to the effect of Rayleigh number, increase in dissipation number contributes to large relative errors between ALA and TALA. Missed data corresponding to smaller Rayleigh numbers are due to much higher relative errors between ALA and TALA such as 37.19% for Rayleigh number of 104 in 4c. The other missed data corresponding to higher Rayleigh number are due to time-dependent convection or different style of convection (e.g., one cell vs. two cells).

Results from the experiments using TALA.Di = 0.25

Di = 0.50

Di = 1.00

Di = 1.50

Di = 2.00

3.1.3 확장 부시네스크 근사

(Extended Boussinesq Approximation, EBA)

확장 부시네스크 근사는 절단 비탄성 유체 근사에서 맨틀이 비압축성이라는 가정을 도입한 것으로 맨틀의 단열압축에 의한 단열 온도 곡선은 존재하지 않는다. 이러한 가정 때문에 실제 단열 압축이 존재하는 맨틀의 온도 구조를 근사하기 위하여 맨틀의 포텐셜(potential) 온도만을 고려한 확장 부시네스크 근사를 적용한 실험을 수행한 다음 그 결과에 단열압축에 의해 발생하는 단열 온도 곡선을 선형적으로 더하는 방법이 널리 사용되었다(e.g., Billen and Hirth, 2007; Lee and King, 2011). 그러나 이 연구는 콤솔 멀티피직스 프로그램의 벤치마크에 목적을 두고 있기 때문에 이러한 차이를 직접 비교하는 연구는 수행하지 않았다. 그러므로 여기서는 확장 부시네스크 지배방정식을 사용한 실험 결과와 앞서 설명된 실험 결과와 직접 비교하지 않았다.

그러나 레일리히 수와 소산 수가 맨틀 대류에 미치는 효과는 위에서 언급한 실험결과와 유사하다. 레일리히 수가 증가함에 따라 맨틀 대류의 강도는 증가하며 소산 수가 증가함에 따라 맨틀 대류의 강도는 감소한다(표 4). 특히 소산 수가 2.0인 경우, 레일리히 수 104와 2×104를 사용한 실험에서는 맨틀 대류가 발생하지 않는데 이는 점성소산이 유체의 거동을 방해하는 효과가 극대화된 경우에 해당한다. 주목할 점은 확장 부시네스크 근사는 점성소산에 의하여 발생한 열에너지와 맨틀이 중력에 대하여 한 일간의 오차가 매우 작다는 것인데 이는 레일리히 수와 소산 수가 큰 값을 사용한 실험에서도 잘 관찰된다. 실험 결과는 기존 벤치마크 결과와 비교되었으며 대부분의 경우 1% 이내의 상대오차를 보인다는 것을 확인할 수 있었다(소산 수 1.5 및 2.0을 사용한 결과는 기존 벤치마크 결과에서 보고되지 않았으므로 비교하지 않았다).

Results from the experiments using EBA.Di = 0.25

Di = 0.50

Di = 1.00

Di = 1.50

Di = 2.00

3.1.4 부시네스크 근사

(Boussinesq Approximation, BA)

마지막으로 가장 일반적으로 쓰이는 유체 거동의 근사된 지배방정식인 부시네스크 근사 지배방정식을 이용하여 벤치마크를 수행하였다. 표 5에서 알 수 있듯이, 부시네스크 근사는 다른 근사 지배방정식을 이용한 실험과 비교하여 가장 강한 맨틀 대류를 발생시킨다. 그리고 다른 근사 지배방정식은 레일리히 수가 커질수록 대류 세포가 두 개 혹은 세 개로 수렴하거나 아니면 시간의존성(time-dependent) 대류를 보이는 것과 대조적으로 부시네스크 근사는 큰 레일리히 수를 사용한 실험에서도 한 개의 대류 세포로 안정적으로 수렴하는 것을 관찰할 수 있었다. 이 실험 결과도 기존 벤치마크 결과와 비교되었으며 모두 1% 이내의 상대오차를 보였다.

Results from the experiments using BA.

3.2 메쉬 해상도(resolution) 테스트

위에서 소개한 다양한 지배방정식에 대한 벤치마크와 더불어 메쉬 해상도(resolution) 테스트를 수행하였다. 특히 해상도 테스트를 통하여 비탄성 유체 근사 지배방정식에서 점성 소산에 의한 열에너지와 맨틀이 중력에 대하여 한 일이 과연 정확하게 일치하는지 검증해 볼 수 있다. 상기 실험들은 상부 메쉬 부분을 세밀화함으로써 넛셀 수를 보다 정확하게 측정하였지만 메쉬 해상도 실험에서는 상부 메쉬 부분을 세밀화하지 않고 단지 요소의 크기를 세 단계로 세분화하여 해상도 테스트를 수행하였다. 가장 적은 개수의 요소를 사용한 실험은 콤솔 멀티피직스가 제공하는 메쉬 해상도 중 ‘coarse’ 옵션을 사용하였으며 총 요소의 개수는 578개이다. 이보다 많은 요소들로 이루어진 ‘fine’ 옵션을 사용한 경우 총 요소의 개수는 2124개이며 가장 많은 요소를 사용한 실험인 ‘extra fine’ 옵션은 총 15244개의 요소를 형성시킨다. 메쉬 해상도 테스트의 편의를 위하여 소산 수는 0.5로 고정시켰으며 단지 레일리히 수만을 변화시켜 실험을 수행하였다.

표 6은 각각의 지배방정식에 해당하는 해상도 테스트의 결과를 표로 나타낸 것이다. 비탄성 유체 근사 지배방정식을 이용한 실험의 경우 넛셀 수가 메쉬 해상도가 증가함에 따라 증가하는 것을 관찰할 수 있다. 이는 넛셀 수가 상부 표면에 가장 가까운 요소들을 이용하여 계산되므로 요소의 크기가 작을 수록 열경계층에서의 온도 변화를 정확하게 측정할 수 있기 때문이다. 즉, 넛셀 수를 계산하는 부분에서는 요소의 크기가 작을수록 보다 정확한 넛셀 수가 계산된다.

넛셀 수가 메쉬 해상도의 영향을 크게 받음에 반하여 다른 결과 값들이 차이는 상대적으로 작음이 관찰되었다. 이는 상대적으로 적은 개수의 요소들로 모델을 구성하여도 맨틀 거동 모델링 결과는 큰 차이가 없음을 지시한다. 그러나 레일리히 수가 증가함에 따라 오차 값이 커지는 경향이 관찰되므로 레일리히 수가 107 정도인 맨틀의 거동을 계산하기 위해서는 어느 정도의 메쉬 해상도는 꼭 확보되어야 한다.

Resolution test.ALA, Di = 0.5

TALA, Di = 0.5

EBA, Di = 0.5

BA

비탄성 유체 근사 지배방정식을 사용한 실험의 경우 레일리히 수와 관계 없이 메쉬 해상도가 증가할수록 점성소산에 의한 열에너지와 맨틀이 중력에 대하여 한일 간의 오차가 작아짐을 확인할 수 있다. 특히 ‘extra fine’의 메쉬 해상도를 사용할 경우 모든 실험에서 두 계산 값간의 상대오차가 0.001% 이내로써 비탄성 유체 근사 지배방정식이 콤솔 멀티피직스에서 정확하게 구현되었음을 지시한다. 이와 달리 모멘텀 방정식에서 역학압력 요소가 제거된 절단 비탄성 유체 근사의 경우에는 메쉬 해상도와 상관없이 거의 일정한 점성소산에 의한 열에너지와 맨틀이 중력에 대하여 한 일 간의 상대오차가 관찰된다

확장 부시네스크 근사 지배방정식을 사용한 실험의 경우에도 메쉬 해상도의 증가와 점성소산에 의한 열에너지와 맨틀이 중력에 대하여 한 일 간의 상대오차가 감소하는 것이 비례한다는 점은 메쉬 해상도의 증가가 실험의 정확도 향상에 기여함을 잘 지시한다. 부시네스크 근사 지배방정식을 사용한 실험의 경우에는 점성소산이 존재하지 않으며 메쉬 해상도의 증가는 다른 실험에서와 같이 넷셀 수의 정확도 향상에 크게 기여한다.


4. 토의 및 결론

유한요소 기반 상용 프로그램인 콤솔 멀티피직스를 이용하여 비압축성 및 압축성 맨틀 대류의 벤치마크를 수행한 결과는 기존 벤치마크 결과들과 아주 잘 합치하였다. 이는 콤솔 멀티피직스가 지니고 있는 장점(그래픽 유저 인터페이스, 손쉬운 물성 및 지배방정식의 수정 및 삽입 등)을 활용하면서도 맨틀 대류 모델링을 수행하기 위한 정확성과 신뢰성을 확보할 수 있음을 뜻한다.

다양한 지배방정식을 사용하여 수행한 실험들은 맨틀의 압축성이 맨틀의 온도 및 대류 구조에 큰 영향을 미친다는 것을 보여주었다. 맨틀의 소산 수로 추정되는 값은 약 0.5 정도이며 레일리히 수는 107 정도이다. 그림 5는 소산 수를 0.5, 레일리히 수 5×105로 가정하여 계산한 각 지배방정식에 대응하는 맨틀 대류의 수렴 결과이다. 온도 분포에서 잘 드러나듯이 압축성을 고려한 실험(그림 5a, b)은 서로 비슷한 결과를 보여주나 비압축성 가정에서 수행된 실험(그림 5d)과 대조해볼 때 큰 온도 분포의 차이를 보인다. 이뿐만 아니라 비탄성 유체 근사 지배방정식에 비압축성 가정을 고려한 확장 부시네스크 근사 지배방정식을 이용한 실험(그림 5c)조차 비압축성 가정에서 수행된 실험(그림 5d)과 비교하여 큰 온도 분포의 차이를 보인다. 이는 맨틀의 압축성 이외에도 점성소산이 맨틀의 온도 및 유동 구조에 큰 영향을 미친다는 것을 잘 지시한다.

Fig. 5.

Temperature distributions in the modeling domain using anelastic liquid approximation (ALA), truncated anelastic liquid approximation (TALA), extended Boussinesq approximation (EBA) and Boussinesq approximation (BA). Dissipation and Rayleigh numbers are 0.5 and 2×105, respectively. Mantle compressibility (ALA and TALA) contributes to higher temperatures to the deeper region of the modeling domain (5a and 5b), contrast to the experiment without mantle compressibility (5d). The experiment (EBA) including viscous dissipation only also contributes to higher mantle temperature to deeper region of the modeling domain but weaker than those of the experiments considering mantle compressibility (ALA and TALA).

이 연구는 맨틀 대류 연구에 쓰이는 각 프로그램 간의 벤치마크를 목적으로 수행되었으므로 실제 맨틀의 물성이나 맨틀 내부에서 발생하는 상전이(phase transformation) 현상은 고려하지 않았다. 실제 맨틀의 거동에 큰 영향을 미치는 점성도는 확산(diffusion) 및 전위(dislocation) 포행(creep)으로 근사 될 수 있으며 온도, 압력, 입자 크기(grain size) 그리고 변형률(strain rate)에 크게 영향을 받는다(Karato and Wu, 1993; Ranalli, 2001). 그러므로 점성을 상수로 가정한 이 실험은 실제 맨틀의 거동을 단순화시킨 실험임을 밝혀둔다. 또한 맨틀에서 온도와 압력의 증가에 의하여 발생하는 감람석(olivine)-와드슬레이아이트(wadsleyite) 등의 상전이와 링우다이트(ringwoodite)-페로브스카이드(perovskite)+마그네시오우스타이트(Magnesiowüstite)등의 상해리(phase dissociation)는 섭입과 맨틀 풀룸(plume)으로 대표되는 맨틀 대류를 강화시키거나 반대로 약화시킨다(King and Ita, 1995; Ita and King, 1998). 그러므로 이러한 상전이 및 상해리 효과와 함께 맨틀의 압축성과 점성소산에 의한 효과들이 복합적으로 작용하여 맨틀 온도 및 거동에 영향을 미칠 것이다.

단순화 및 이상화된 실험에도 불구하고 본 연구는 콤솔 멀티피직스가 향후 맨틀 대류 및 섭입을 모델링하는데 편리성과 정확성을 내포하고 있으며 앞으로 관련 연구들에 큰 기여를 할 가능성이 있음을 보여준다. 대부분의 기존 연구들이 맨틀의 비압축성을 가정한 지배방정식(부시네스크 근사 등)을 이용하여 맨틀 대류 및 섭입 작용을 연구하였으며 현재까지도 널리 활용되고 있다. 다만 압축성에 의한 맨틀 온도 및 거동의 차이가 일부 연구에서만 제한적으로 비교 및 대조하여 평가되었지만 (Lee and King, 2009), 장기적으로 맨틀 및 섭입 작용 연구는 맨틀의 압축성이 중요한 모델링의 요소로써 고려되어야 한다. 향후 전산 지구동력학 연구들의 시작점으로써 본 연구의 벤치마크 결과는 고무적이다.

Acknowledgments

본 논문의 수준 향상에 크게 도움을 주신 김승섭, 김영희 교수님께 감사드립니다. 또한 편집위원장이신 손영관 교수님께 감사드립니다. 이 연구는 2011년 정부(교육과학기술부)의 재원으로 한국연구재단의 지원을 받아 수행되었습니다(NRF-2011-35B-1-C00043).

References

  • M.I Billen,, Modeling the Dynamics of Subducting Slabs, Annual Review of Earth and Planetary Sciences, (2008), 36,, p325-356. [https://doi.org/10.1146/annurev.earth.36.031207.124129]
  • M.I. Billen,, G Hirth,, Rheologic controls on slab dynamics, Geochemistry Geophysics Geosystems, 8, Q08012. [https://doi.org/10.1029/2007GC001597]
  • F Birch,, Elasticity and constitution of the Earth's interior, Journal of Geophysical Research, (1952), 57, p227-286. [https://doi.org/10.1029/JZ057i002p00227]
  • B. Blankenbach,, F. Busse,, U. Christensen,, L. Cserepes,, D. Gunkel,, U. Hansen,, H. Harder,, G. Jarvis,, M. Koch,, G. Marquart,, D. Moore,, P. Olson,, H. Schmeling,, T Schnaubelt,, A BENCHMARK COMPARISON FOR MANTLE CONVECTION CODES, Geophysical Journal International, (1989), 98, p23-38. [https://doi.org/10.1111/j.1365-246X.1989.tb05511.x]
  • E. Carminati,, P Petricca,, State of stress in slabs as a function of large-scale plate kinematics, Geochem. Geophys. Geosyst., 11, Q04006, (2010.). [https://doi.org/10.1029/2009GC003003]
  • H Goldstine,, he Computer: from Pascal to von Neumann, Princeton University Press, (1972.).
  • J. Ita,, S.D King,, The influence of thermodynamic formulation on simulations of subduction zone geometry and history, Geophysical Research Letters, (1998), 25, p1463-1466. [https://doi.org/10.1029/98GL51033]
  • G. Jarvis,, D McKenzie,, Convection in a compressible fluid with infinite Prandtl number, Journal of Fluid Mechanics, (1980), 96.. [https://doi.org/10.1017/S002211208000225X]
  • S.-I. Karato,, P Wu,, Rheology of the upper mantle: a synthesis, Science, (1993), 260, p771-778. [https://doi.org/10.1126/science.260.5109.771]
  • S.D. King,, J Ita,, Effect of slab rheology on mass-transport across a phase-transition boundary, Journal of Geophysical Research, (1995), 100, p20211-20222. [https://doi.org/10.1029/95JB01964]
  • S.D. King,, C. Lee,, P.E.v. Keken,, W. Leng,, S. Zhong,, E. Tan,, N. Tosi,, M.C Kameyama,, A community benchmark for 2-D Cartesian compressible convection in the Earth's mantle, Geophysical Journal International, (2010), 180, p73-87. [https://doi.org/10.1111/j.1365-246X.2009.04413.x]
  • S.D. King,, A. Raefsky,, B.H Hager,, Conman: vectorizing a finite element code for incompressible two-dimensional convection in the Earth's mantle, Physics of the Earth and Planetary Interiors, (1990), 59, p195-207. [https://doi.org/10.1016/0031-9201(90)90225-M]
  • C. Lee,, S.D King,, Effect of mantle compressibility on the thermal and flow structures of the subduction zones, Geochem. Geophys. Geosyst., 10, Q01006, (2009.). [https://doi.org/10.1029/2008GC002151]
  • C. Lee,, S.D King,, Dynamic buckling of subducting slabs reconciles geological and geophysical observations, Earth and Planetary Science Letters, (2011), 312, p360-370. [https://doi.org/10.1016/j.epsl.2011.10.033]
  • D.P. McKenzie,, J.C Sclater,, Heat-flow inside the island arcs of the northwest Pacific, Journal of Geophysical Research, (1968), 73, p3173-3179.
  • L.G.J. Montési,, M.D. Behn,, L.B. Hebert,, J. Lin,, J.L Barry,, Controls on melt migration and extraction at the ultraslow Southwest Indian Ridge 10°-16°E.J. Geophys, Res., 116, B10102, (2011.). [https://doi.org/10.1029/2011JB008259]
  • L. Moresi,, S.J. Zhong,, M Gurnis,, The accuracy of finite element solutions of Stokes' flow with strongly varying viscosity, Physics of the Earth and Planetary Interiors, (1996), 97, p83-94. [https://doi.org/10.1016/0031-9201(96)03163-9]
  • E.R. Oxburgh,, D.L Turcotte,, Problem of high heat flow and volcanism associated with zones of descending mantle convective flow, Nature, 218, 1041-1043, (1968.). [https://doi.org/10.1038/2181041a0]
  • R Pryor,, Multiphysics Modeling Using COMSOL v.4: A First Principles Approach, Mercury Learning & Information, (2012.).
  • G Ranalli,, Mantle rheology: radial and lateral viscosity variations inferred from microphysical creep laws, Journal of Geodynamics, (2001), 32, p425-444. [https://doi.org/10.1016/S0264-3707(01)00042-4]
  • G. Schubert,, D. Turcotte,, P Olson,, Mantle Convection in the Earth and Planets, Cambridge Univ. Press, Cambridge, 1st Ed, (2001.). [https://doi.org/10.1017/CBO9780511612879]
  • N.H Sleep,, Mantle plumes from top to bottom, Earth-Science Reviews, (2006), 77, p231-271. [https://doi.org/10.1016/j.earscirev.2006.03.007]
  • P.E. van Keken,, C. Currie,, S.D. King,, M.D. Behn,, A. Cagnioncle,, J. He,, R.F. Katz,, S.-C. Lin,, E.M. Parmentier,, M Spiegelman,, K Wang,, A community benchmark for subduction zone modeling, Physics of The Earth and Planetary Interiors, (2008), 171, p187-197. [https://doi.org/10.1016/j.pepi.2008.04.015]

Fig. 1.

Fig. 1.
The modeling domain formulated in this study using the COMSOL Multiphysics®. The non-dimensionalized length of all the walls of the square box is 1. The uppermost region of the modeling domain consists of very finer elements for the accurate calculation of the Nusselt number at the top wall.

Fig. 2.

Fig. 2.
Temperature distributions in the modeling domain using the anelastic liquid approximation (ALA). Dissipation number is fixed as 0.25 but Rayleigh numbers are varied as 104, 105 and 106. The single clockwise convection cell is obtained in the experiments using Rayleigh numbers of 104 and 105. The temperature distribution of the experiment using a Rayleigh number of 106 shows count-clockwise convection cell but the temperature distribution is vertically inverted for a comparison among the temperature distributions. The black arrows indicate direction and vigor of convecting mantle. For a better view, the scale factors for black arrows are differently used as 0.0020, 0.0004 and 0.0002 for 2a, 2b and 2c, respectively. As Rayleigh number increases, upwelling and downwelling as well as the thermal boundary layers become thinner. The unit of temperature is K without non-dimensionalization.

Fig. 3.

Fig. 3.
Temperature distributions in the modeling domain using anelastic liquid approximation (ALA). All the Rayleigh numbers are fixed as 2×105 but dissipation number is varied as 0.50, 1.00, 1.50 and 2.00. With increases in dissipation number, vigor of mantle convection is substantially weakened, expressed sluggish mantle convection and layered temperature distribution with depth. The experiments using higher dissipation number tends to develop multiple convection cells such as two or three convection cells. The black arrows indicate direction and vigor of convecting mantle. For a better view, the scale factors for black arrows are differently used as 0.0002, 0.0004, 0.0008 and 0.0016 for 3a, 3b, 3c and 3d, respectively. The black arrows in 3c and 3d indicate two convection cells; the experiments developing three convection cells are not shown here. The experiment using a dissipation number of 2.00 apparently develops strong layered temperature distributions with depth. The unit of temperature is K without non-dimensionalization.

Fig. 4.

Fig. 4.
Relative errors between the experiments using anelastic liquid approximation (ALA) and truncated anelastic liquid approximation (TALA) with varying Rayleigh and dissipation numbers. All the parameters except for the mean temperatures of the modeling domain show a clear tendency to decreasing relative errors between ALA and TALA with increases in Rayleigh number. Contrary to the effect of Rayleigh number, increase in dissipation number contributes to large relative errors between ALA and TALA. Missed data corresponding to smaller Rayleigh numbers are due to much higher relative errors between ALA and TALA such as 37.19% for Rayleigh number of 104 in 4c. The other missed data corresponding to higher Rayleigh number are due to time-dependent convection or different style of convection (e.g., one cell vs. two cells).

Fig. 5.

Fig. 5.
Temperature distributions in the modeling domain using anelastic liquid approximation (ALA), truncated anelastic liquid approximation (TALA), extended Boussinesq approximation (EBA) and Boussinesq approximation (BA). Dissipation and Rayleigh numbers are 0.5 and 2×105, respectively. Mantle compressibility (ALA and TALA) contributes to higher temperatures to the deeper region of the modeling domain (5a and 5b), contrast to the experiment without mantle compressibility (5d). The experiment (EBA) including viscous dissipation only also contributes to higher mantle temperature to deeper region of the modeling domain but weaker than those of the experiments considering mantle compressibility (ALA and TALA).

Table 1.

Model parameters.

Symbol Explanation and Si Unit Value
d Depth (m) 106
g Gravity (m/s2) 10
Cp Heat capacity (J/kg.K) 1250 (BA)
1250 (Di = 0.25)
625 (Di = 0.50)
312.5 (Di = 1.00)
208.3 (Di = 1.50)
156.3 (Di = 2.00)
ΔT Temperature difference (K) 3000
κ Heat diffusivity (m2/s) 10-6
ρ0 Density at the top wall (kg/m3) 4000
α Thermal expansivity (/K) 3.125×10-5
Γ Grüneisen’s parameter (.) 1
μ Dynamic viscosity (Pa.s) 3.750×10-23 (Ra = 104)
1.875×10-23 (Ra = 2×104)
7.500×10-22 (Ra = 5×104)
3.750×10-22 (Ra = 1×105)
1.875×10-22 (Ra = 2×105)
7.500×10-21 (Ra = 5×105)
3.750×10-21 (Ra = 106)

Table 2.

Results from the experiments using ALA.Di = 0.25

Ra Nu # Vsurf-max Vsurf <T> <φ> <W>
1.0E+04 4.410 58.065 38.826 0.515 0.848 0.848
2.0E+04 5.536 91.640 62.908 0.519 1.129 1.129
5.0E+04 7.424 164.002 116.752 0.526 1.599 1.599
1.0E+05 9.229 252.509 184.406 0.532 2.050 2.050
2.0E+05 11.413 386.336 288.060 0.538 2.596 2.596
5.0E+05 14.712 657.238 494.947 0.545 3.425 3.425
1.0E+06 14.867 749.419 540.730 0.539 3.469 3.468

Di = 0.50

Ra Nu # Vsurf-max Vsurf <T> <φ> <W>
1.0E+04 3.816 52.789 35.048 0.522 1.380 1.380
2.0E+04 4.717 83.314 56.553 0.529 1.824 1.824
5.0E+04 6.200 148.868 104.068 0.539 2.557 2.556
1.0E+05 7.559 227.216 161.682 0.547 3.231 3.231
2.0E+05 8.919 331.026 236.411 0.555 3.914 3.914
5.0E+05 9.279 397.198 274.344 0.563 4.116 4.115+
1.0E+06 no convergence to steady state

Di = 1.00

Ra Nu # Vsurf-max Vsurf <T> <φ> <W>
1.0E+04 2.466 37.584 24.391 0.510 1.360 1.360
2.0E+04 2.909 59.083 38.775 0.515 1.781 1.781
5.0E+04 3.541 101.569 67.095 0.523 2.399 2.398
1.0E+05 3.876 137.531 89.126 0.529 2.759 2.758
2.0E+05 4.078 171.944 106.126 0.534 2.994 2.992
5.0E+05* 6.144 226.274 143.416 0.525 5.035 5.031+
1.0E+06 no convergence to steady state

Di = 1.50

Ra Nu # Vsurf-max Vsurf <T> <φ> <W>
1.0E+04 1.309 17.216 10.292 0.478 0.415 0.415
2.0E+04 1.419 26.703 15.493 0.476 0.570 0.570
5.0E+04* 1.710 29.485 18.548 0.462 0.984 0.984+
1.0E+05* 1.992 48.287 30.175 0.461 1.388 1.387+
2.0E+05* 2.324 76.860 47.637 0.461 1.877+ 1.874+
5.0E+05* 2.801 133.614 80.672 0.460 2.602+ 2.593+
1.0E+06* no convergence to steady state

Di = 2.00

Ra Nu # Vsurf-max Vsurf <T> <φ> <W>
1.0E+04- 1.027 4.643 2.612 0.501 0.023 0.023
2.0E+04 1.100+ 14.189+ 7.213 0.487 0.152+ 0.151+
5.0E+04* 1.186 19.319 12.136 0.463 0.322 0.321
1.0E+05* 1.245 30.403 18.856 0.451 0.440 0.438
2.0E+05* 1.308 46.060 28.029 0.441 0.565 0.561+
5.0E+05* 1.398 76.790 44.900 0.429 0.749 0.737+
1.0E+06* 1.468+ 109.888 60.950 0.420 0.893 0.869+

Table 3.

Results from the experiments using TALA.Di = 0.25

Ra Nu # Vsurf-max Vsurf <T> <φ> <W>
1.0E+04 4.421 58.722 39.290 0.513 0.853 0.850
2.0E+04 5.550 92.568 63.589 0.518 1.134 1.131
5.0E+04 7.441 165.499 117.878 0.525 1.606 1.602
1.0E+05 9.249 254.713 186.059 0.531 2.058 2.053
2.0E+05 11.436 389.604 290.464 0.537 2.605 2.600
5.0E+05 14.731 662.195 498.296 0.543 3.433 3.427
1.0E+06 no convergence to steady state

Di = 0.50

Ra Nu # Vsurf-max Vsurf <T> <φ> <W>
1.0E+04 3.858 54.290 36.094 0.519 1.408 1.391
2.0E+04 4.766 85.445 58.083 0.526 1.857 1.837
5.0E+04 6.258 152.280 106.566 0.537 2.597 2.572
1.0E+05 7.621 232.028 165.171 0.545 3.276 3.246
2.0E+05 8.966 336.533 240.096 0.554 3.951 3.922
5.0E+05 9.295 399.578 275.571 0.563 4.131 4.117
1.0E+06 no convergence to steady state

Di = 1.00

Ra Nu # Vsurf-max Vsurf <T> <φ> <W>
1.0E+04 2.562 40.596 26.423 0.509 1.463 1.398
2.0E+04 3.011 63.224 41.617 0.515 1.894 1.817
5.0E+04 3.634 107.046 70.794 0.523 2.506 2.422
1.0E+05 3.922 141.621 91.506 0.530 2.820 2.754
2.0E+05 4.106 175.265 107.672 0.534 3.032 2.986+
5.0E+05 no convergence to steady state
1.0E+06 no convergence to steady state

Di = 1.50

Ra Nu # Vsurf-max Vsurf <T> <φ> <W>
1.0E+04 1.360 19.551 11.705 0.478 0.479 0.447
2.0E+04 1.469 29.651 17.221 0.476 0.632 0.596
5.0E+04* 1.762 31.715 19.955 0.462 1.054 1.012+
1.0E+05* 2.047 51.344 32.099 0.461 1.463 1.412+
2.0E+05* 2.378 80.849 50.117 0.461 1.953+ 1.897+
5.0E+05* 2.847 138.735 83.653 0.460 2.668+ 2.615+
1.0E+06 no convergence to steady state

Di = 2.00

Ra Nu # Vsurf-max Vsurf <T> <φ> <W>
1.0E+04 1.053+ 7.585+ 4.159+ 0.497 0.062++ 0.060++
2.0E+04 1.122+ 16.421+ 8.413 0.484 0.182+ 0.174+
5.0E+04* 1.213+ 21.384 13.435 0.461 0.360 0.345
1.0E+05* 1.274 33.138 20.559 0.450 0.479 0.459
2.0E+05* 1.335 49.459 30.094 0.440 0.601 0.580
5.0E+05* 1.421 81.026 47.276 0.428 0.779 0.756+
1.0E+06* 1.482 114.680 63.425 0.419 0.916+ 0.891+

Table 4.

Results from the experiments using EBA.Di = 0.25

Ra Nu # Vsurf-max Vsurf <T> <φ> <W>
1.0E+04 4.092 54.870 36.590 0.491 0.773 0.773
2.0E+04 5.156 86.928 59.605 0.494 1.039 1.039
5.0E+04 6.937 155.673 110.964 0.500 1.484 1.484
1.0E+05 8.644 239.712 175.592 0.505 1.911 1.911
2.0E+05 10.721 367.238 275.175 0.510 2.432 2.432
5.0E+05 14.016 635.358 483.195 0.518 3.260 3.260
1.0E+06 15.563 841.575 621.144 0.519 3.650 3.650

Di = 0.50

Ra Nu # Vsurf-max Vsurf <T> <φ> <W>
1.0E+04 4.092 54.870 36.590 0.491 0.773 0.773
2.0E+04 5.156 86.928 59.605 0.494 1.039 1.039
5.0E+04 6.937 155.673 110.964 0.500 1.484 1.484
1.0E+05 8.644 239.712 175.592 0.505 1.911 1.911
2.0E+05 10.721 367.238 275.175 0.510 2.432 2.432
5.0E+05 14.016 635.358 483.195 0.518 3.260 3.260
1.0E+06 15.563 841.575 621.144 0.519 3.650 3.650

Di = 1.00

Ra Nu # Vsurf-max Vsurf <T> <φ> <W>
1.0E+04 2.190 34.243 22.236 0.467 1.187 1.187
2.0E+04 2.646 55.373 36.676 0.468 1.638 1.638
5.0E+04 3.354 99.492 67.650 0.474 2.341 2.341
1.0E+05 3.957 150.240 103.289 0.482 2.940 2.940
2.0E+05 4.434 209.343 141.282 0.494 3.420 3.419
5.0E+05** 5.675 219.238 139.203 0.488 4.667 4.667
1.0E+06* 6.686 320.165 204.078 0.492 5.684 5.683+

Di = 1.50

Ra Nu # Vsurf-max Vsurf <T> <φ> <W>
1.0E+04 1.315 18.886 11.908 0.466 0.482 0.482
2.0E+04 1.510 32.509 20.585 0.458 0.771 0.771
5.0E+04 1.755 58.539 37.120 0.456 1.131 1.131
1.0E+05 1.921 84.854 52.570 0.458 1.375 1.375
2.0E+05* 2.238 75.124 47.044 0.443 1.861 1.861
5.0E+05* 2.717 132.376 81.562 0.445 2.574 2.574
1.0E+06* 3.079 192.452 114.548 0.447 3.112 3.111

Di = 2.00

Ra Nu # Vsurf-max Vsurf <T> <φ> <W>
1.0E+04 1.000 0.000 0.000 0.500 0.000 0.000
2.0E+04 1.000 0.000 0.000 0.500 0.000 0.000
5.0E+04 1.118 28.281 15.154 0.466 0.262 0.262
1.0E+05 1.168 44.151 22.470 0.457 0.367 0.367
2.0E+05* 1.291 47.817 29.943 0.433 0.609 0.609
5.0E+05* 1.386 79.895 48.541 0.422 0.800 0.800
1.0E+06* 1.455 112.477 64.821 0.416 0.936 0.936

Table 5.

Results from the experiments using BA.

Ra Nu # Vsurf-max Vsurf <T>
1.0E+04 4.884 61.632 41.422 0.500
2.0E+04 6.195 97.154 67.289 0.500
5.0E+04 8.401 173.316 125.194 0.500
1.0E+05 10.529 266.450 198.334 0.500
2.0E+05 13.157 408.690 312.246 0.501
5.0E+05 17.603 719.580 565.064 0.501
1.0E+06 21.900 1105.788 882.340 0.502

Table 6.

Resolution test.ALA, Di = 0.5

Ra Mesh Nu # Vsurf-max Vsurf <T> <φ> <W> Error (%)
1.0E+04 coarse 3.566 52.726 34.983 0.521 1.374 1.374 0.014
fine 3.730 52.767 35.031 0.521 1.379 1.379 0.001
extra fine 3.798 52.780 35.044 0.521 1.380 1.380 0.000
2.0E+04 coarse 4.287 83.273 56.472 0.528 1.817 1.817 0.025
fine 4.573 83.284 56.526 0.528 1.823 1.823 0.002
extra fine 4.689 83.283 56.536 0.528 1.824 1.824 0.000
5.0E+04 coarse 5.331 149.060 104.051 0.537 2.552 2.550 0.050
fine 5.908 148.860 104.036 0.538 2.556 2.556 0.005
extra fine 6.146 148.769 104.001 0.539 2.557 2.557 0.000
1.0E+05 coarse 6.135 228.010 161.925 0.545 3.233 3.230 0.080
fine 7.063 227.300 161.680 0.547 3.233 3.232 0.009
extra fine 7.470 226.991 161.501 0.547 3.231 3.231 0.000
2.0E+05 coarse 6.822 332.476 236.940 0.552 3.926 3.921 0.127
fine 8.153 331.148 236.349 0.555 3.918 3.918 0.016
extra fine 8.782 330.471 235.862 0.556 3.912 3.912 0.000
5.0E+05 coarse 7.180 408.889 281.777 0.559 4.190 4.178 0.297
fine 8.456 397.386 274.146 0.562 4.119 4.118 0.035
extra fine 9.135 396.931 273.622 0.563 4.117 4.117 0.001

TALA, Di = 0.5

Ra Mesh Nu # Vsurf-max Vsurf <T> <φ> <W> Error (%)
1.0E+04 coarse 3.596 54.232 36.030 0.518 1.403 1.386 1.195
fine 3.767 54.272 36.078 0.519 1.407 1.391 1.183
extra fine 3.838 54.281 36.090 0.519 1.409 1.392 1.182
2.0E+04 coarse 4.318 85.417 58.009 0.525 1.851 1.830 1.131
fine 4.614 85.420 58.057 0.526 1.857 1.836 1.109
extra fine 4.736 85.414 58.066 0.526 1.858 1.838 1.106
5.0E+04 coarse 5.358 152.507 106.571 0.535 2.593 2.566 1.045
fine 5.954 152.281 106.539 0.536 2.597 2.571 1.000
extra fine 6.202 152.180 106.500 0.537 2.598 2.573 0.994
1.0E+05 coarse 6.156 232.899 165.459 0.543 3.279 3.247 0.973
fine 7.109 232.124 165.178 0.545 3.278 3.248 0.902
extra fine 7.529 231.800 164.988 0.545 3.276 3.247 0.892
2.0E+05 coarse 6.831 337.983 240.621 0.550 3.963 3.930 0.829
fine 8.183 336.651 240.026 0.553 3.955 3.926 0.723
extra fine 8.826 335.960 239.529 0.554 3.949 3.921 0.707
5.0E+05 coarse 7.183 411.506 283.302 0.559 4.208 4.181 0.622
fine 8.466 399.771 275.414 0.562 4.135 4.120 0.367
extra fine 9.150 399.308 274.832 0.563 4.132 4.118 0.329

EBA, Di = 0.5

Ra Mesh Nu # Vsurf-max Vsurf <T> <φ> <W> Error (%)
1.0E+04 coarse 3.166 47.986 31.717 0.482 1.183 1.183 0.002
fine 3.306 48.046 31.774 0.482 1.187 1.187 0.000
extra fine 3.363 48.063 31.790 0.482 1.188 1.188 0.000
2.0E+04 coarse 3.841 76.446 51.818 0.485 1.599 1.598 0.004
fine 4.089 76.502 51.900 0.486 1.604 1.604 0.000
extra fine 4.189 76.514 51.921 0.486 1.605 1.605 0.000
5.0E+04 coarse 4.822 137.344 96.366 0.493 2.286 2.285 0.008
fine 5.333 137.286 96.465 0.494 2.292 2.292 0.001
extra fine 5.545 137.259 96.479 0.494 2.294 2.294 0.000
1.0E+05 coarse 5.595 211.397 151.686 0.500 2.935 2.935 0.014
fine 6.437 211.046 151.743 0.501 2.940 2.939 0.002
extra fine 6.807 210.894 151.699 0.502 2.940 2.940 0.000
2.0E+05 coarse 6.335 320.710 233.111 0.508 3.692 3.691 0.020
fine 7.628 319.826 233.130 0.510 3.692 3.692 0.003
extra fine 8.250 319.370 232.917 0.511 3.690 3.690 0.000
5.0E+05 coarse 6.828 446.725 312.185 0.519 4.208 4.207 0.033
fine 8.454 446.726 313.152 0.522 4.213 4.213 0.006
extra fine 9.258 447.287 313.721 0.523 4.218 4.218 0.000

BA

Ra Mesh Nu # Vsurf-max Vsurf <T>
1.0E+04 coarse 4.572 61.533 41.312 0.500
fine 4.796 61.602 41.396 0.500
extra fine 4.872 61.623 41.420 0.500
2.0E+04 coarse 5.595 97.053 67.113 0.500
fine 6.019 97.108 67.242 0.500
extra fine 6.170 97.121 67.276 0.500
5.0E+04 coarse 7.043 173.374 124.892 0.500
fine 7.969 173.265 125.107 0.500
extra fine 8.338 173.210 125.140 0.500
1.0E+05 coarse 8.135 266.841 197.876 0.500
fine 9.698 266.461 198.230 0.500
extra fine 10.405 266.235 198.212 0.500
2.0E+05 coarse 9.148 409.507 311.458 0.500
fine 11.601 408.908 312.198 0.500
extra fine 12.910 408.302 312.012 0.500
5.0E+05 coarse 10.554 723.565 565.773 0.500
fine 14.377 720.646 565.493 0.500
extra fine 17.027 718.923 564.690 0.500