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

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.2 비압축성 및 압축성 맨틀 대류 지배방정식

2.2.1 맨틀 거동 지배방정식

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

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

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

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

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

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

2.2.2 비탄성 유체 근사

(Anelastic Liquid Approximation, ALA)

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

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

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

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

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

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

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

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에 서술되었다.

Table 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.

수치 모델링을 위한 변수로써 과거 연구에서 변수로 사용된 레일리히 수를 10⁴, 2×10⁴, 5×10⁴, 10⁵, 2×10⁵, 5×10⁵ 그리고 10⁶으로 변화시켜 사용하였다. 그리고 각 레일리히 수에 해당하는 실험에 각각 소산 수를 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 #), 대류하는 맨틀의 상부 표면에서의 최대 속도(V_surf-max), 대류하는 맨틀의 상부 표면에서의 선적분 값(V_surf), 대류 세포 전체의 평균 온도(<T>), 대류 세포 전체에서 점성소산에 의하여 발생한 열에너지(<φ>), 마지막으로 대류 세포 전체에서 맨틀이 중력에 대하여 한 일(<W>)을 보여준다. 물론 여기서 제시된 값들은 모두 무차원화된 값이다. 각 값들의 계산 방법은 본 연구와 비교 대상이 되는 기존 벤치마크 연구를 참고하길 바란다(King et al., 2010).

Table 2.

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으로 고정되어 있으며 레일리히 수가 10⁴, 10⁵ 그리고 10⁶으로 각각 10배씩 증가하였을 때의 대류 세포의 온도 분포를 나타내고 있다. 이미 알려진 것과 같이 레일리히 수가 증가할수록 대류 세포의 상부 및 하부에 존재하는 열경계층(thermal boundary layer)의 두께가 감소하는 것이 확인되며 이는 넛셀 수의 증가로 잘 표현된다(넛셀 수 4.410, 9.229 그리고 14.867이 레일리히 수 10⁴, 10⁵, 그리고 10⁶을 사용한 실험 결과에 각각 대응한다). 그리고 대류의 강도가 강해짐에 따라 맨틀 대류 속도가 증가하는데 이는 레일리히 수의 증가와 함께 상부 표면에서의 대류의 최고 속도 및 선적분 값이 증가하는 것으로 잘 표현된다(V_surf-max의 값인 58.065, 252.509 그리고 749.419이 레일리히 수 10⁴, 10⁵, 그리고 10⁶을 사용한 실험 결과에 각각 대응하며, V_surf의 값인 38.826, 184.406 그리고 540.730이 레일리히 수 10⁴, 10⁵, 그리고 10⁶을 사용한 실험 결과에 각각 대응한다). 점성소산에 의한 열에너지 (<φ>)는 레일리히 수의 증가와 함께 선형적인 증가 관계를 보인다. 그리고 맨틀이 중력에 대하여 한 일(<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인 경우, 대류의 강도는 급격하게 약화되어 레일리히 수가 10⁴인 경우 매우 약한 대류만이 관찰된다. 넛셀 수가 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×10⁵ 이상인 경우 상대오차는 2% 이내로써 매우 작아진다. 지구 맨틀의 레일리히 수가 약 10⁷ 정도로 추정되므로(Schubert et al., 2001), 이 실험은 절단 비탄성 유체 근사 지배방정식이 비탄성 유체 근사 지배방정식으로 근사 될 수 있음을 지시한다. 실제로 표 3에서 관찰되듯이 점성소산에 의한 열에너지는 맨틀이 중력에 대하여 한 일에 비하여 상대오차가 0.1~6% 내외이다. 위와 같이 기존 실험 결과(King et al., 2010)와 비교하였을 때 비록 소산 수가 큰 실험에서는 상대오차가 증가하는 경향을 보였으나 대부분 1~4% 정도로 작으며 소산 수 2.0와 작은 레일리히 수를 사용한 실험(예, 레일리히 수: 10⁴)에서만 상대오차가 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).

Table 3.

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

Table 4.

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% 이내의 상대오차를 보였다.

Table 5.

Results from the experiments using BA.

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

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

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

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

Table 6.

Resolution test.ALA, Di = 0.5

TALA, Di = 0.5

EBA, Di = 0.5

BA

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

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

4. 토의 및 결론

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

다양한 지배방정식을 사용하여 수행한 실험들은 맨틀의 압축성이 맨틀의 온도 및 대류 구조에 큰 영향을 미친다는 것을 보여주었다. 맨틀의 소산 수로 추정되는 값은 약 0.5 정도이며 레일리히 수는 10⁷ 정도이다. 그림 5는 소산 수를 0.5, 레일리히 수 5×10⁵로 가정하여 계산한 각 지배방정식에 대응하는 맨틀 대류의 수렴 결과이다. 온도 분포에서 잘 드러나듯이 압축성을 고려한 실험(그림 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]