인장 분지의 열구조 진화에 대한 벤치마크 연구

1. 서 론

인장 분지(extensional basin)는 분지의 한 종류로서, 지각의 양 끝에 인장응력(extensional stress)이 작용하여 암석권(lithosphere)의 두께가 정단층(normal fault)을 비롯한 인장 구조에 의해 얇아져 형성된다. 예를 들어, 섭입대(subduction zone)에서 섭입하는 판의 후퇴(slab roll-back)에 의해 발생한 인장응력에 의하여 형성되는 배호분지(back-arc basin)가 있으며 대표적인 예로 우리나라의 동해가 있다( Uyeda and Kanamori, 1979). 따라서 분지는 동해의 형성과 같은 지구조사 연구에서 중요하다( Tamaki et al., 1992 ; Jolivet et al., 1994 ). 그리고 분지는 석유(oil), 천연가스(natural gas), 가스 하이드레이트(gas hydrate) 같은 화석 자원이 발견되는 곳으로서, 자원의 형성에 중요한 분지 퇴적물과 같은 기원암(source rock), 자원을 포함하고 있는 저류암(reservoir rock), 트랩 구조(trap) 등을 이해하기 위해서는 지구조(tectonic) 및 열적 과정(thermal processes)의 해석이 중요하다(Allen and Allen, 2013). 그러므로 이러한 분지의 지구조와 형성 과정은 지구물리학, 암석학, 퇴적학을 포함한 지질과학 여러 분야에서 꾸준히 연구되었다(e.g., Jolivet and Tamaki, 1992; Park et al., 2009 ; Sato et al., 2014 ; Yoon et al., 2014 ).

분지 형성 매커니즘(mechanism)으로 초기에는 암석권 하부 온도 구조의 변화에 관심이 집중되었다. 예를 들어, 암석권 하부에 열적 이상(thermal anomaly)이 발생하여 암석권 하부가 무거운 에클로자이트(eclogite)로 전환된 다음, 시간이 지남에 따라 지각이 냉각되면서 무거워진 지각 하부가 침강하면서 분지가 형성된다는 가설( Haxby et al., 1976 )과 에클로자이트의 형성 없이 지온의 증가로 암석권이 상승했다가 지온이 감소하면서 서서히 하강할 때, 침식과 열수축(thermal contraction)에 의해 분지가 형성된다는 가설(Sleep and Snell, 1976) 등이 있다. 하지만 이러한 가설들은 분지의 형성을 설명하는 데에 충분하지 않은데, 열적 이상에 의해 분지가 형성되려면 대륙 지각의 경우 해양지각에 비해 훨씬 지역적으로 집중된 큰 열적 이상이 존재해야 하고, 침식과 열수축에 의하여 분지가 형성되려면 상당한 침식이 발생해야 하는데 그러한 증거가 부족하기 때문이다.

이러한 가설 대신 McKenzie (1978)와 Jarvis and McKenzie (1980)는 분지가 암석권의 인장으로 인하여 형성되었다고 제시하였다. McKenzie (1978)는 분지가 순식간에 인장(instantaneous extension)되어 형성되었음을 가정하고, 인장이 완료된 이후부터 시간이 지남에 따라 깊이에 따른 지온(geotherm)과 표면에서의 열유동(heat flux)을 해석해(analytical solution)로 제시하였다. 그러나, McKenzie (1978)의 해석해는 분지의 인장이 순식간에 발생하고 이후 분지가 냉각된다는 가정에서 유도되었기 때문에 분지의 형성 기간이 매우 짧은 경우에만 유도된 해석해가 적용될 수 있는 한계를 지닌다. 이러한 문제를 보완하기 위해 Jarvis and McKenzie (1980)는 인장 기간 동안 분지 표면에서의 열손실을 고려하여 분지의 지온과 표면에서의 열유동 값을 수치해석해(numerical solution)로 계산하여 제시하였다.

하지만, 이 연구들은 여러 가지 가정들을 기반으로 수행되었기 때문에 보다 복잡한 분지 형성 및 진화사를 연구하는 데에는 한계를 지닌다. 예를 들어, 두 연구는 모두 분지의 인장이 분지 전체에 걸쳐 균일하게 일어났음을 가정하였다. 그러나, 실제로는 분지의 인장 과정에서 변형집중(strain localization)이 발생한다는 것이 지구물리 탐사 및 수치모델링 연구 등을 통하여 알려져 있다( Buiter et al., 2006 ; Sato et al., 2006 , 2014; So et al., 2012 ). 또한 정단층과 같은 인장 구조들이 고려되지 않아 인장 분지의 복잡한 지구조 작용이 분지 형성 및 진화에 미치는 영향이 고려되지 못하였다( Buiter et al., 2006 ; Choi et al., 2013 ). 그러므로, 궁극적으로 분지의 형성 및 진화사를 연구하기 위하여 정단층을 비롯한 지질 구조, 분지 형성 기간 및 이후에서 발생하는 퇴적 작용, 변형률(strain rate)과 온도에 의존하는 암석권 및 맨틀의 리올로지(rheology) 등을 고려하는 것이 필요하다. 최근 다양한 요소들을 고려하여 분지 형성 및 진화사에 대한 수치모델링 연구가 수행되어 과거 알지 못했던 분지 형성 및 진화사에 대한 이해가 높아지고 있다(e.g., Burov and Gerya, 2014; Göğüş, 2015).

그러므로, McKenzie (1978)와 Jarvis and McKenzie (1980)의 연구에서 고려되지 못하였던 분지 형성에서 발생하는 복잡한 현상들을 고려한 연구를 수행하는 것이 궁극적으로 요구된다. 이를 위한 초기 연구로서, 우리는 컴퓨터 수치모델링을 이용하여 McKenzie (1978)와 Jarvis and McKenzie (1980)가 제시한 해석해 및 수치해석해를 벤치마크 하였다. 또한 인장기간이 분지의 진화에 미치는 영향을 평가하고 벤치마크 결과를 울릉분지의 열구조 해석에 적용하였다. 이 연구를 통하여 분지 모델링의 초기 연구를 검증하고 울릉 분지를 포함한 인장 분지 해석에 활용될 수 있는 일차적인 해석 자료를 제시하고자 한다.

2. 실험 방법

2.1 McKenzie (1978)의 분지 형성 완료 이후의 열구조 해석

2.1.1 모델의 구성

McKenzie (1978)는 암석권이 순식간에 인장되어 인장 분지가 형성되었음을 가정하였다. 그는 인장이 완료된 이후부터 시간이 지남에 따라, 분지 하부에서 지표면으로의 열전달(heat transfer)이 열확산(heat diffusion)에 의해서만 이루어진다고 가정하여 분지하부의 지온과 표면에서의 열유동을 해석해로 제시하였다. McKenzie (1978)에서 암석권의 인장 정도는 인장인수(stretch factor, β)에 따라 결정된다. 예를 들어, 인장되기 전 암석권의 길이와 두께가 a이면(그림 1a), 인장이 완료된 이후 암석권의 너비는 aβ가 되고 두께는 a/ β가 된다(그림 1b). 그리고 인장이 일어나 암석권의 두께가 감소하면, 감소한 두께만큼 암석권 아래에 존재하는 일정한 온도를 지닌 연약권 맨틀(asthenospheric mantle)이 위로 상승하여 인장된 분지의 지온이 증가한다. 시간이 지남에 따라 분지 표면에서의 열유동에 의하여 분지 하부의 온도가 감소하기 때문에 무한대의 시간이 흐르면 다시 분지 인장 전의 지온으로 수렴한다(그림 1c).

Fig. 1.

(a)-(c) Schematic diagrams showing the process of the basin extension from McKenzie (1978). a and β indicate thickness of lithosphere and stretch factor, respectively. L and A mean lithosphere and asthenosphere, respectively. The basin extension instantaneously occurs from (a) to (b). After extension, the basin cools with time by heat loss through the surface. The right graphs show the geotherms of the basin and underlying mantle with depth. (d) The model domain formulated by using the COMSOL Multiphysics®. The mesh refinement is applied to the uppermost part of the domain for the accurate calculation of the heat flux at the surface. The right graph is the initial temperature condition relevant to the stretch factor (β) of 2.

2.1.2 콤솔 멀티피직스(COMSOL Multiphysics®)를 이용한 수치모델링

우리는 유한요소법 기반 상용 프로그램인 콤솔 멀티피직스(COMSOL Multiphysics®, www.comsol.com)를 사용하여 McKenzie (1978) 해석해에 대한 벤치마크를 수행하였다. 여기에서부터 표기로 인한 혼란을 막기 위해 McKenzie (1978)를 벤치마크한 우리의 수치모델은 ‘IBE (Instantaneous Basin Extension) 모델’로 표기한다.

IBE 모델에서 모델 범위(model domain)는 McKenzie (1978)가 제시한 것처럼 가로와 세로의 길이가 125 km인 2차원 정사각형으로 설정하였다( 그림 1d). 그림 1d에서 제시한 모델에서는 인장인수(β)를 2로 적용하여 모델 범위를 구성하였다. IBE 모델에서 사용한 상수들은 인장되기 전 암석권의 길이와 두께(a), 암석권 표면의 온도( T₀), 암석권-연약권 맨틀과의 경계와 연약권 맨틀전체의 온도( T₁), 암석권의 밀도( ρ_c ), 연약권 맨틀의 밀도( ρ_m )이다. 각각의 구체적인 값들은 McKenzie (1978)에서 사용된 값들과 같다(표 1). 그리고 모델에서 암석권과 연약권 맨틀은 강체(rigid body)로 가정되었으며 리올로지나 항복점(yielding point) 같은 물질의 고유 성질(material property)은 고려되지 않았다. 모델에서의 열전달은 다음의 열확산 방정식(heat diffusion equation)에 의해 계산된다.

Table 1.

Model Parameters.

$\[ \frac{\partial T}{\partial t}=\kappa \left(\frac{{\partial }^{2}T}{\partial x^2}+\frac{{\partial }^{2}T}{\partial y^2}\right) \]$

1) 열확산 방정식

κ는 열확산률(thermal diffusivity)이며 밀도와 열전도도(heat conductivity, k) 그리고 열용량(heat capacity, C _p )에 의해 $$ \kappa =\frac{k}{\rho C_p}$$ 로 정의된다.

IBE 모델의 경계조건으로서 상부 경계( T₀) 및 하부경계 온도( T₁)는 각각 0℃와 1,333℃로 고정되었으며 좌우 경계는 단열(insulated) 경계로 설정되었다. 그리고 분지 인장이 완료된 이후 뜨거운 연약권 맨틀의 상승에 의해 형성된 초기 지온 조건은 다음과 같다.

$\[ \begin{array}{cc}{T}_{initial}=\frac{\beta T_{1}depth}{a}\hfill & \left(depth\le \frac{\beta }{a}\right)\hfill \\ T_{initial}=T_1\hfill & \left(depth>\frac{\beta }{a}\right)\hfill \end{array} \]$

2) 초기 지온

깊이(depth)는 모델 상부 경계에서부터 측정된다. 그림 1b는 인장인수(β)가 2인 경우의 초기 지온을 나타낸 것이다.

유한요소법을 이용한 수치모델링에서는 모델 범위를 구성하는 메쉬(mesh)를 어떻게 형성하느냐에 따라 계산의 정확도가 달라진다. 일반적으로 모델범위를 많은 개수의 요소(element)로 구성할수록 계산의 정확도가 높아지지만, 요소의 개수가 증가할수록 계산 시간과 비용이 증가하므로 지나치게 많은 개수의 요소로 모델 범위를 구성하는 것은 비효율적이다. IBE 모델에서는 콤솔 멀티피직스에서 기본적으로 지원하는 finer 크기의 삼각형 요소를 사용한 비정렬 격자(unstructured mesh)로 모델 범위를 구성하고 요소 내 온도 근사를 위하여 이차 형상 함수(quadratic shape function)를 사용하였다. 그리고 표면에서의 열유동을 보다 정확하게 계산하기 위해 메쉬 세밀화(mesh refinement)를 모델 범위의 상부에 적용하였다(그림 1b). 요소의 총 개수는 메쉬 세밀화 과정으로 인해 인장인수(β)의 값에 따라 조금씩 달라지며 인장인수(β)가 2인 경우, 총 요소의 개수는 4,032개이다.

수치해석은 시간에 따른 열전달을 계산하기 위하여 시간 의존성 솔버(time-dependent solver)를 사용하여 수행되었으며, 분지 인장이 완료된 시간( t)을 0으로 정의하여 시간이 지남에 따라 분지 하부의 지온과 표면에서의 열유동을 계산하였다. 이 때, 시간간격(time-step)은 1 Myr로 적용하였으며, 양해법(explicit method)을 사용하여 계산을 수행하였다. 그리고 인장인수(β)를 1.25, 1.5, 2, 4, 10으로 변화시켜 다양한 분지 인장 정도를 고려하였다.

2.2 Jarvis and McKenzie (1980)의 분지 형성 과정과 형성 이후의 열구조 해석

2.2.1 모델의 구성

McKenzie (1978)는 인장이 순식간에 일어났음을 가정하였기 때문에, 분지의 인장 기간 동안 열유동에 의한 분지 냉각을 고려할 수 없다. 그래서 Jarvis and McKenzie (1980)는 분지 인장 및 인장 완료 이후 분지 냉각 과정에서의 지온과 표면에서의 열유동을 해석하였다(그림 2a). 분지 인장 시기 동안 인장되는 암석권과 이에 따른 연약권 맨틀의 상승을 모사하기 위하여 암석권과 연약권 맨틀의 거동을 다음과 같은 공간에 대한 함수(u( x), v(y))로 표현하였다.

Fig. 2.

(a) Schematic diagrams showing the velocity trajectory of the lithosphere and asthenosphere from Jarvis and McKenzie (1980). The ‘V0’ is the vertical velocity at y=a. The blue line is the velocity trajectory of the lithosphere and asthenosphere. (b) The modeling domain formulated by using the COMSOL Multiphysics®. Except for the open boundary condition applied to the bottom and side walls, all the modeling conditions are the same as Figure 1b. The right graph is the initial temperature condition (geotherm).

3) 인장 속도

4) 상승 속도

Δt 는 암석권의 인장 기간(s)이며, β는 인장인수이다. 그림 2a의 유선(stream line)은 인장에 따른 암석권 및 연약권 맨틀의 거동을 보여준다. 인장이 완료된 이후에는 암석권과 연약권 맨틀의 거동이 발생하지 않기 때문에 열확산에 의해서만 열전달이 이루어진다. 모델에서의 열확산률(κ)과 열용량( C _p )을 제외한 각각의 상수들은 모두 McKenzie (1978)에서 제시된 것들과 동일하다(표 1).

2.2.2 콤솔멀티피직스(COMSOL Multiphysics®)를 이용한 수치모델링

IBE 모델의 표기법과 유사하게, Jarvis and McKenzie (1980)를 벤치마크한 우리의 수치모델은 ‘PBE (Protracted Basin Extension) 모델’로 표기한다. PBE 모델에서의 모델 범위는 2차원 정사각형으로 IBE 모델과 동일하게 설정하였으며, 암석권과 연약권 맨틀의 거동을 위에서 제시된 함수( u(x), v(y))로 모사하였다. 그리고 모델에서 암석권과 연약권 맨틀은 위에서 제시한 함수(식 3 및 4)에 의해서 거동이 결정되는 유체(fluid)로 가정되었기 때문에 점성도(viscosity)와 같은 유체의 성질은 고려되지 않았다. 그러므로, 암석권과 연약권 맨틀의 거동에 따른 열진화는 열대류와 열전도를 함께 고려한 다음의 에너지 방정식을 이용하여 계산된다.

$\[ \frac{\partial T}{\partial t}+u\frac{\partial T}{\partial x}+v\frac{\partial T}{\partial y}=\kappa \left(\frac{{\partial }^{2}T}{\partial x^2}+\frac{{\partial }^{2}T}{\partial y^2}\right) \]$

5) 에너지 방정식

왼쪽 항은 대류항(convection term)이며 오른쪽 항은 확산항(diffusion term)이다.

PBE 모델의 경계조건으로서 모델의 상부 경계( T₀) 및 하부 경계 온도( T₁)는 각각 0℃와 1,333℃로 고정되었으며, 연약권 맨틀은 열린(open) 경계인 하부 경계에서 유입된다. 암석권의 인장 및 인장에 의한 연약권 맨틀의 좌우 경계로의 유출을 고려하기 위하여 좌우 경계는 열린(open) 경계 조건으로 설정되었다. 그리고 초기 지온은 다음과 같이 정의된다.

$\[ T_{initial}=\frac{T_{1}depth}{a} \]$

6) 초기 지온

PBE 모델에 사용된 계수 및 상수 값들은 표 1에 제시되었다.

메쉬는 IBE 모델에 사용된 것과 유사하며, 총 요소의 개수는 3,972개이다. 해석은 시간 의존성 솔버를 사용하였으며, 시간 간격은 각 모델에 따른 인장기간의 1/100로 적용하고 양해법을 사용하여 계산을 수행하였다. 그리고 IBE 모델과 같이 인장인수(β)를 1.25, 1.5, 2, 4, 10으로 변화시켜 계산을 수행하였다.

3. 실험 결과

3.1 McKenzie (1978) 벤치마크

McKenzie (1978)가 제시한 해석해를 IBE 모델을 이용하여 벤치마크한 결과로서 각각의 인장인수(β)에 해당하는 지온(그림 3a-e)과 상부 표면에서의 열유동(그림 3f)을 계산하였다. 인장인수( β)가 2인 경우, 인장 직후(t=0)의 초기 지온은 모델 표면에서부터 62.5 km 깊이(a/ β)까지 약 21℃/km의 지온 구배(geothermal gradient)를 보이며, 깊이 62.5 km부터 125 km까지는 1,333℃로 지온이 고정되었다(그림 3c). 이후 시간이 지나면서 인장된 분지가 냉각됨에 따라 지온 구배는 점차 완만해지고, 1억년(100 Myr)보다 더 오랜 시간이 지난 후에는 분지 인장 전의 지온인 정상상태(steady-state)에 도달한다. 모든 인장인수(β)에 해당하는 실험에서 깊이에 따른 지온은 인장 직후부터 시간이 지남에 따라 정상상태로 수렴하였다. 그리고 인장인수(β) 1.25와 10에 해당하는 계산결과를 비교하면 인장 직후부터 동일하게 1억년(100 Myr)이 지났을 때, 인장인수(β)가 1.25인 경우의 지온은 거의 정상상태의 지온에 근접하나 인장인수( β)가 10인 경우의 지온은 깊이 60 km 지점에서의 온도를 비교했을 때 정상상태의 지온과 약 150℃의 차이를 보인다. 이는 인장인수(β)의 값이 클수록 지온이 정상상태로 수렴하는데 걸리는 시간이 증가한다는 것을 지시한다.

다음으로 각각 인장인수(β)에 해당하는 실험에서 표면에서의 열유동의 시간에 따른 변화를 측정하였다( 그림 3f). 인장 직후의 열유동은 지온 구배에 비례하므로 각각의 인장인수(β)에 해당하는 열유동은 0.8 β HFU(Heat Flux Unit)가 된다. 예를 들어, 인장 인수(β)가 10인 경우 실험 초기 표면에서의 열유동은 8 HFU이며 인장인수(β)가 1.25인 경우 열유동은 1 HFU이다. 시간에 따라 열유동은 분지 하부가 냉각되어 지온이 완만해지면서 감소하게 된다. 특히 인장인수(β)가 클수록 냉각 초기의 지온의 감소가 크므로(e.g., 그림 3e), 열유동은 시간에 따라 급격하게 감소하게 된다. 그러나 약 1억년(100 Myr)의 시간이 지나면 모든 실험에서 열유동 값은 1.5 HFU 이하로 감소하여 점진적으로 0.8 HFU에 수렴한다. IBE모델을 이용한 모든 실험에서 측정된 지온 및 열유동은 McKenzie (1978)가 제시한 해석해와 1% 이내의 상대오차 이내에서 잘 맞는 것을 확인할 수 있었다.

Fig. 3.

(a)-(e) Geotherms calculated from the IBE (Instantaneous Basin Extension) models and analytic solutions from McKenzie (1978), corresponding to the varying stretch factors (β). Horizontal and vertical axes are temperature and depth, respectively. The green solid and red dotted lines correspond to our numerically calculated geotherms and the geotherms from McKenzie (1978), respectively. t is the cooling time after basin extension. (f) The heat fluxes at the surface with time that calculated using the IBE model with the heat fluxes from McKenzie (1978). The horizontal and vertical axes correspond to the cooling time and heat flux (HFU, heat flux unit), respectively. The black dotted line indicates the heat flux at the surface, corresponding to the stretch factor (β) of 1 or infinite cooling time.

3.2 Jarvis and McKenzie (1980) 벤치마크

암석권이 인장되어 분지가 형성되는 동안 분지표면에서 발생하는 열유동에 의하여 분지 하부의 지온은 변화한다. Jarvis and McKenzie (1980)는 인장 기간 동안 변화하는 지온을 벤치마크하기 위하여 무차원 수 G′를 다음과 같이 제시하였다.

$\[ G\text{'}=\frac{a^2\ln \beta }{\kappa \mathrm{\Delta }t} \]$

7) 무차원수 G′

a는 인장되기 전 암석권의 두께, β는 인장인수, κ는 열확산률, Δ t는 암석권의 인장 기간이다. 그리고 G′를 이용하여 암석권 표면에서 인장 속도를 다음과 같이 계산할 수 있다.

$\[ V=\frac{\kappa G\text{'}}{a} \]$

8) 인장 속도

McKenzie (1978)에서 제시된 암석권의 두께와 열확산율을 적용하면 G′가 5, 10, 20, 50, 100일 때 인장 속도는 각각 0.1, 0.2, 0.41, 1.01, 2.03 cm/y이다. 그리고 인장인수(β)가 같은 경우 인장 기간은 인장 속도가 빠를수록 작은 값을 가진다. 예를 들어, 인장 인수( β)가 10이고 인장 속도가 2.03 cm/y일 때의 인장 기간은 약 14.2 Myr이지만, 인장 속도가 0.1 cm/y일 때의 인장 기간은 약 283.6 Myr로 증가한다. 또한, 같은 인장 속도를 가질 때 인장인수(β)가 클수록 인장 기간은 증가한다. 인장 속도가 2.03 cm/y일 때, 인장인수(β)가 10인 경우의 인장 기간은 약 14.2 Myr이지만, 인장인수( β)가 1.25인 경우 인장 기간은 약 1.4 Myr이다. G′=∞인 경우는 인장 속도가 무한대이며, 인장 기간이 무한소라는 것을 의미하므로 McKenzie (1978)의 초기 지온과 같으며, G′=0인 경우는 정상상태의 지온과 같다.

각각의 인장인수(β)에 해당하는 깊이에 따른 지온을 측정하였다(그림 4). PBE 모델을 이용하여 계산된 지온은 인장이 완료된 직후의 지온이며, IBE 모델의 초기 지온(G′=∞)과 비교하였다. 인장인수(β)가 2인 경우에서 G′값이 클수록 즉, 인장 속도가 빠를수록 인장 기간은 짧아지기 때문에(e.g., Δt= 4.3 Myr)전도에 의한 열손실이 상대적으로 적어 지온이 IBE모델의 초기 지온과 큰 차이를 보이지 않는다. 그러나, 인장 속도가 느릴수록 인장 기간 동안(e.g., Δt=85.4 Myr) 전도에 의한 열손실이 발생하여 지온이 감소하므로 인장 완료 직후의 지온은 IBE 모델의 초기 지온에 비해 훨씬 정상상태에 가까운 값을 보인다(그림 4c).

Fig. 4.

Geotherms calculated from the PBE (Protracted Basin Extension) models, corresponding to the varying stretch factors (β) after basin extension. The axes are the same as Figure 3. The black, blue, green, pink and red solid lines correspond to G′=5, 10, 20, 50 and 100, respectively. The black dashed and dotted lines indicate the infinity and 0 of G′, respectively. The calculated velocities of the basin extension equivalent to G′s are described in the legend. Δ t is duration of the basin extension.

다음으로 각각의 인장인수(β)에 해당하는 표면에서의 열유동을 측정하였다(그림 5). 그래프의 가로축은 인장이 완료된 시간을 0으로 기준하여 음수(-)에 해당하는 부분은 인장 기간에 해당하며, 양수(+)에 해당하는 부분은 인장이 완료된 이후 분지가 냉각되는 기간에 해당된다. 인장인수(β)가 2에 해당하는 그래프를 보면, 인장 속도가 2.03 cm/y일 때 인장이 완료된 직후의 표면에서의 열유동은 약 1.6 HFU이지만 인장 속도가 0.1 cm/y일 때는 약 1.3 HFU로 인장 속도가 느릴수록 인장 직후의 열유동이 감소하는 것을 확인할 수 있다(그림 5c). 그리고 이러한 인장 속도에 따른 각각의 열유동 값의 차이는 인장인수(β)가 커질수록 증가한다. 인장인수(β)가 1.25인 실험에서는 인장 속도가 가장 빠를 때와 작을 때의 열유동의 차이가 약 0.0006 HFU정도로 매우 작지만(그림 5a), 인장인수(β)가 10인 경우에는 약 4.03 HFU로 인장인수(β)가 1.25일때보다 훨씬 크다(그림 5e). 이처럼 인장인수( β)와 인장 속도에 따라서 인장 완료 직후의 표면에서의 열유동은 차이를 보인다. 그러나 시간이 지남에 따라 분지가 냉각되고, 그에 따라 표면에서의 열유동도 감소한다. 그리고 시간이 지남에 따라 모든 경우에서 표면에서의 열유동은 정상상태의 열유동인 0.8 HFU에 수렴한다. 이 때 0.8 HFU에 수렴하기까지 걸리는 시간은 인장인수(β)가 클수록, 인장 속도가 빠를수록 증가한다.

Jarvis and McKenzie (1980)가 제시한 수치해석해는 복잡한 계산을 요구하므로, 우리의 수치모델인 PBE 모델의 결과값과 직접 비교를 수행하지는 못하였다. 그러므로, 우리는 Jarvis and McKenzie (1980)의 그림 3과 그림4를 각각 우리의 그림 4 및 그림 5와 직접 비교하였으며 대체로 잘 일치하는 것을 확인하였다.

Fig. 5.

Surface heat fluxes corresponding to the varying stretch factors ( β) calculated from the PBE (Protracted Basin Extension) model. Horizontal and vertical axes correspond to the time from the basin extension to cooling and heat flux at the surface expressed as heat flux unit (HFU). The other notations are the same as in Figure 4.

3.3 메쉬 해상도(mesh resolution) 테스트

벤치마크 실험과 함께 메쉬 해상도(mesh resolution) 테스트를 수행하였다. 이를 통해 수치모델링을 위하여 어느 정도의 해상도를 가진 메쉬를 사용하는 것이 효율적인지 검증하였다. 메쉬 해상도 테스트를 위하여 사용된 실험 조건과 상수들은 이미 제시된 실험에서 쓰인 값들과 같으며, 단지 콤솔 멀티피직스에서 제공하는 메쉬 해상도 옵션인 coarse, finer, extremely fine를 사용하여 메쉬 해상도를 변화시킴으로써 실험을 수행하였다. 여기에서부터 coarse, finer, extremely fine의 메쉬 해상도를 각각 low, intermediate, high로 표기한다.

IBE 모델은 인장인수(β)가 2인 모델을 사용하여 메쉬 해상도 테스트를 진행하였다. 그리고 PBE 모델에서는 인장인수( β)가 2이고, 인장 속도가 0.41 cm/y(G′=20)에 해당하는 모델을 이용하여 테스트를 수행하였다. 그리고 각 테스트의 메쉬 해상도는 low, intermediate, high로 나누어 구성하였다. 각 해상도에 따른 요소의 개수와 요소의 최소/최대 면적은 다음과 같으며, 간결한 표기를 위해 요소의 개수, 가장 작은 요소의 면적, 가장 큰 요소의 면적의 순서로 표기하였다. IBE 모델에서 low는 568개, 4.6 km², 74.2 km ²로 구성하고, intermediate는 4,032개, 0.6 km², 10.2 km² 그리고 high는 54,494개, 0.038 km², 0.717 km²로 구성하였다. 또한 PBE 모델은 low는 578개, 4.6 km², 75.1 km ²로 구성하고, intermediate는 3,972개, 0.6 km², 9.07 km² 그리고 high는 54,446개, 0.037 km², 0.711 km²로 구성하였다.

메쉬 해상도 테스트에 의해 계산된 분지의 지온과 표면에서의 열유동 결과값을 그림 6에 도시하였다. McKenzie (1978)의 해석해에 대해 IBE 모델의 해상도 low, intermediate, high에 따른 결과의 상대오차를 각각 계산해보았다. 그 결과 지온의 경우, low일 때는 0.038%, intermediate일 때는 0.032% 그리고 high 일 때는 0.031%의 상대오차를 보였다(그림 6a). 그리고 표면에서의 열유동은 low일 때는 0.042%, intermediate일 때는 0.041%, high 일 때는 0.037%의 상대오차를 보였으며 지온과 표면에서의 열유동 모두 높은 해상도일수록 McKenzie (1978)의 해석해와의 상대오차가 감소하였다(그림 6b). 그리고 모든 해상도 조건에서 상대오차는 0.1% 이내로 IBE 모델의 결과와 McKenzie (1978) 해석해가 잘 일치함을 확인하였다.

Jarvis and McKenzie (1980)의 수치해석해는 계산을 수행할 수 없어 그래프의 육안 비교만을 수행하였고, 수치해석해와 PBE 모델의 결과 간의 직접비교를 수행하지 못하였다. 따라서 가장 높은 해상도인 high의 결과가 Jarvis and McKenzie (1980)의 수치해석해에 가장 가까운 값임을 가정하고 이에 대해 low와 intermediate에 따른 결과의 상대오차를 각각 계산해 보았다. 그 결과 지온의 경우, low일 때는 0.156%, intermediate일 때는 0.025%의 상대오차를 보였다(그림 6c). 또한 표면에서의 열유동은 low일 때는 0.086%, intermediate일 때는 0.014%의 상대오차를 보이며 모두 0.1~0.2%이내로 PBE 모델의 결과와 Jarvis and McKenzie (1980)의 수치해석해가 잘 일치함을 확인하였다(그림 6d). 이러한 결과는 우리가 앞에서 수행한 IBE 모델 및 PBE 모델의 벤치마크 실험에서 사용한 메쉬 해상도가 지온 및 열유동을 계산하기에 충분히 높다는 것을 지시한다.

Fig. 6.

Mesh resolution tests of the IBE (Instantaneous Basin Extension) and PBE (Protracted Basin Extension) models. The stretch factor (β) is 2. The green, red and blue lines indicate the mesh resolutions of high, intermediate and low, respectively. (a) The geotherms with time after basin extension calculated from the IBE model. (b) The heat fluxes at the surface with time calculated from the IBE model. (c) Same as (a) but the PBE model. Duration of the basin extension and extension velocity for the PBE model are 21.3 Myr and 0.41 cm/y, respectively. (d) Same as (b) but the PBE model.

4. 토의 및 결론

우리는 분지가 순식간에 인장되었음을 가정한 McKenzie (1978)의 해석해와 분지의 인장 기간을 고려한 Jarvis and McKenzie (1980)의 수치해석해를 콤솔 멀티피직스를 이용하여 벤치마크하였다. 인장된 분지의 지온과 표면에서의 열유동을 측정하여 비교한 결과, 우리의 IBE 모델의 결과는 McKenzie (1978)와 0.1%이내의 상대오차를 보이며 잘 일치하였다. 비록 PBE 모델의 결과와 Jarvis and McKenzie (1980)의 수치해석해와의 직접적인 비교는 이루어지지 못하였지만 제시된 결과와의 육안 비교를 통하여 잘 일치함을 확인하였다. 이는 우리가 McKenzie (1978)와 Jarvis and McKenzie (1980)가 제시한 해석해 및 수치해석해를 컴퓨터 수치모델링을 이용하여 성공적으로 재현하였음을 의미한다.

McKenzie (1978)와 Jarvis and McKenzie (1980)가 제시한 해석들 간의 가장 큰 차이점은 Jarvis and McKenzie (1980)는 분지의 인장 기간 동안 표면에서의 열손실을 고려할 수 있다는 것이다. McKenzie (1978)가 가정한 것처럼 순식간에 분지가 인장되는 경우는 존재하지 않는다. 특히 분지의 형성 기간이 충분히 길다면( 그림 5), 분지 형성 기간 중의 열유동에 의한 열손실을 고려해야 한다. 인장 기간 동안의 열손실에 의한 분지의 열구조 차이를 검증하기 위하여 McKenzie (1978)와 Jarvis and McKenzie (1980)를 각각 벤치마크한 IBE 및 PBE 모델을 이용하여 열유동의 진화를 비교해보았다. 우선 인장인수( β)가 1.25로 작은 경우에서 인장 기간 동안 발생한 열손실이 인장 완료 이후 열유동에 미치는 영향을 검증하였다( 그림 5a). 비교를 위해 인장 속도가 0.1, 0.2, 0.41, 1.01, 2.03 cm/y일 때 계산된 PBE 모델과 IBE 모델의 열유동과의 상대오차를 계산하였다. 상대오차는 고려된 인장속도에 따라 각각 1.96, 0.77, 0.3, 0.09, 0.03%로 측정되었다. 인장 속도가 느릴수록 IBE 및 PBE 모델에서 계산된 열유동과의 차이가 증가한다는 것은 느린 인장 속도에 의해 형성된 분지의 경우 인장 기간에서의 열손실이 인장 완료 이후 열유동에 큰 영향을 미친다는 것을 보여준다. 이와 반대로, 상대적으로 빠른 인장 속도에 의해 형성된 분지의 경우 IBE 모델로 분지 형성 이후의 열유동을 해석하는 것이 유용함을 의미한다.

다음으로 인장인수(β)가 10으로 큰 경우, 앞의 경우와 같은 인장속도들을 사용하여 인장 기간 동안 발생한 열손실로 인한 인장 완료 후 열유동의 변화를 계산하였다(그림 5e). 이를 위하여 IBE 및 PBE 모델에서 계산한 열유동간의 상대오차를 계산하였다. 상대오차는 고려된 인장속도에 따라 각각 47.97, 35.97, 23.95, 12.13, 6.74%의 값을 보인다. 인장속도가 2.03 cm/y로 가장 빠른 경우, 상대오차는 인장인수(β)가 1.25로 작을 때는 0.03%를 나타냈지만 인장인수(β)가 10으로 클 때는 6.74%로 상대적으로 큰 값을 보인다. 인장인수(β)가 클 때 IBE 및 PBE 모델간의 열유동의 차이가 크다는 것은 인장인수(β)가 큰 분지형성의 경우 인장 속도가 빠르더라도 인장 기간 동안의 열손실이 인장 완료 이후 열유동에 큰 영향을 미친다는 것을 지시한다. 따라서 인장인수(β)가 큰 분지일수록 빠른 인장 속도에 의해 형성된 분지와 마찬가지로 PBE 모델이 분지 형성 기간 및 완료 이후의 열진화 해석에 적합하다. 그러나, 인장인수가(β)가 상대적으로 작은 분지 형성에서는 IBE 모델이 인장완료 이후 열유동을 계산하는데 유용함을 의미한다.

벤치마크를 위하여 개발된 IBE 및 PBE 모델을 동해 울릉 분지의 열유동 계산에 적용해보았다. 동해는 28 Ma에서 12 Ma까지 태평양 판(Pacific plate)과 필리핀 판(Philippine plate)의 후퇴로 인해 발생한 인장응력에 의해 형성되었다( Tamaki et al., 1992 ; Jolivet et al., 1994 ). 울릉 분지는 동해의 인장 기간 동안 형성된 분지 중 하나로 정확한 인장 기간이 알려지지 않았다. 그러므로 울릉 분지의 화산 활동이 시작된 20 Ma을 분지 형성 시작 시기로 가정하고 동해의 형성이 완료되었을 때 울릉 분지의 인장이 완료되었음을 가정하여 울릉 분지의 인장 기간을 20 Ma부터 12 Ma까지 8 Myr로 가정하였다( Chough and Lee, 1992). 울릉 분지 인장 전의 암석권의 두께가 알려지지 않았으므로 일반적인 대륙암석권(continental lithosphere)의 두께인 80 km와 100 km로 암석권의 두께를 가정하였다. 현재 울릉 분지의 너비의 근사값인 200 km를 사용하여 인장인수(β)를 인장 전의 암석권의 두께가 80 km일 때는 2.5, 100 km일 때 2로 가정하였다. 그리고 아직까지 알려지지 않은 암석권 아래 연약권 맨틀의 온도는 1,350℃와 1,450℃ 두 가지의 경우로 가정하였다. 이외의 나머지 조건은 모두 앞서 제시된 두 모델과 동일하게 적용하였다. 각각의 인장인수(β)에 따라 맨틀의 온도를 다르게 적용하여 실험을 수행하고 표면에서의 열유동을 계산하였다(그림 7a와 7b). 두 모델 모두에서 인장인수(β)가 2일 때 보다 2.5일 때, 맨틀의 온도가 1,350℃일 때보다 1450℃일 때 표면에서의 열유동이 더 높게 나타났다. 그리고 이 결과를 울릉 분지의 열유동과 비교해보았다. Kim et al. (2010) 에 따르면 울릉 분지 주변부의 열유동은 약 2.3 HFU이고, 중심부는 약 2.8 HFU이다(그림 7c). 이는 IBE 및 PBE 모델 모두에서 맨틀의 온도가 1450℃이고, 인장인수(β)가 2.5일 때의 결과와 대체로 일치한다.

Fig. 7.

(a) and (b) The calculated heat fluxes at the surface corresponding to the estimated stretch factors (β) of 2 and 2.5 for the Ulleung basin using the IBE (Instantaneous Basin Extension) and PBE (Protracted Basin Extension) models, respectively. The horizontal and vertical axes correspond to the time from the basin extension to cooling and heat flux at the surface expressed as heat flux unit (HFU). The red and blue lines correspond to the calculated heat fluxes using mantle potential temperatures of 1450℃ and 1350℃, respectively. The dashed and solid lines correspond to the IBE and PBE models, respectively. (c) Measured heat fluxes (mW/m2) at the surface in the Ulleung basin and nearby ( Kim et al., 2010 ).

기존의 McKenzie (1978)와 Jarvis and McKenzie (1980)에서는 실제 분지 형성 과정을 단순화 하여 열구조를 해석하였기 때문에 실제 분지 형성 과정에서 일어나는 암석권과 맨틀의 리올로지에 의한 정단층과 같은 인장 구조, 응력 또는 변형 집중 그리고 퇴적작용 등의 지구조 작용을 고려할 수 없다. 그러나 분지 형성을 수치모델링을 통해 연구하기 위하여 이러한 요소들은 필수적으로 고려되어야 한다. 구체적으로, 암석권에서의 응력 또는 변형 집중을 일으킬 수 있는 취약대(weak zone) 혹은 단층대(fault zone), 퇴적물에 의해 가해지는 퇴적 하중(sediments loading) 그리고 전단 탄성 계수(shear modulus)의 차이 등에 의한 응력집중 효과 등을 고려해야 한다(e.g., So et al., 2012 ). 또한, 동해의 경우 두 개의 판이 섭입하여 형성되기 때문에 3차원 모델링의 수행이 필수적으로 요구된다. 우리가 개발한 IBE 및 PBE 모델은 이와 같은 요소들을 향후에 개발 및 고려함으로써, 보다 사실적인 분지 형성 및 진화사 연구에 활용될 수 있을 것이다.

Acknowledgments

세심하게 본 논문을 검토해주신 김승섭, 소병달 교수님께 감사드립니다. 그리고 편집위원이신 최은서 교수님께 감사드립니다. 이 연구는 2011년 정부(교육과학기술부)의 재원으로 한국연구재단의 지원을 받아 수행되었습니다(NRF-35B-2011-1-C00043).

References

• Allen, P.A., Allen, J.R., (2013), Basin analysis: Principles and application to petroleum play assessment, John Wiley & Sons.
• Buiter, S.J.H., Babeyko, A.Y., Ellis, S., Gerya, T.V., Kaus, B.J.P., Kellner, A., Schreurs, G., and Yamada, Y., (2006), The numerical sandbox: comparison of model results for a shortening and an extension experiment , Geological Society, London, Special Publications, 253, p29-64. [https://doi.org/10.1144/GSL.SP.2006.253.01.02]
• Burov, E., and Gerya, T., (2014), Asymmetric three-dimensional topography over mantle plumes, Nature, 513, p85-89. [https://doi.org/10.1038/nature13703]
• Choi, E., Tan, E., Lavier, L.L., and Calo, V.M., (2013), DynEarthSol2D: An efficient unstructured finite element method to study long-term tectonic deformation , Journal of Geophysical Research: Solid Earth, 118, p2429-2444. [https://doi.org/10.1002/jgrb.50148]
• Chough, S.K., and Lee, K.E., (1992), Multi-stage volcanism in the Ulleung Back-arc Basin, East Sea (Sea of Japan) , Island Arc, 1, p32-39. [https://doi.org/10.1111/j.1440-1738.1992.tb00055.x]
• Göğüş, O.H., (2015), Rifting and subsidence following lithospheric removal in continental back arcs , Geology, 43, p3-6. [https://doi.org/10.1130/G36305.1]
• Haxby, W., Turcotte, D., and Bird, J., (1976), Thermal and mechanical evolution of the Michigan Basin, Tectonophysics, 36, p57-75. [https://doi.org/10.1016/0040-1951(76)90006-8]
• Jarvis, G.T., and McKenzie, D.P., (1980), Sedimentary basin formation with finite extension rates, Earth and Planetary Science Letters, 48, p42-52. [https://doi.org/10.1016/0012-821X(80)90168-5]
• Jolivet, L., and Tamaki, K., (1992), Neogene kinematics in the Japan Sea region and volcanic activity of the northeast Japan arc , Proc. ODP, Sci Results, 127/128(Part 2), Ocean Drilling Program, College Station, TX, p1311-1331.
• Jolivet, L., Tamaki, K., and Fournier, M., (1994), Japan Sea, opening history and mechanism: A synthesis, Journal of Geophysical Research: Solid Earth, 99, p22237-22259. [https://doi.org/10.1029/93JB03463]
• Kim, Y.-G., Lee, S.-M., and Matsubayashi, O., (2010), New heat flow measurements in the Ulleung Basin, East Sea (Sea of Japan): relationship to local BSR depth, and implications for regional heat flow distribution , Geo-Marine Letters, 30, p595-603. [https://doi.org/10.1007/s00367-010-0207-x]
• McKenzie, D., (1978), Some remarks on the development of sedimentary basins, Earth and Planetary Science Letters, 40, p25-32. [https://doi.org/10.1016/0012-821X(78)90071-7]
• Park, G.S., Park, J.S., Kwon, B.D., Kim, C.H., and Park, C.H., (2009), Tectonic structure modeling around the Ulleung Basin and Dokdo using potential data , J. Korean Earth Sci. Soc., 30, p165-175, (in Korean with English abstract). [https://doi.org/10.5467/JKESS.2009.30.2.165]
• Sato, T., No, T., Kodaira, S., Takahashi, N., and Kaneda, Y., (2014), Seismic constraints of the formation process on the back-arc basin in the southeastern Japan Sea , Journal of Geophysical Research: Solid Earth, 119, p1563-1579. [https://doi.org/10.1002/2013JB010643]
• Sato, T., Sato, T., Shinohara, M., Hino, R., Nishino, M., and Kanazawa, T., (2006), P-wave velocity structure of the margin of the southeastern Tsushima Basin in the Japan Sea using ocean bottom seismometers and airguns , Tectonophysics, 412, p159-171. [https://doi.org/10.1016/j.tecto.2005.09.001]
• Sleep, N.H., and Snell, N.S., (1976), Thermal contraction and flexure of mid-continent and Atlantic marginal basins , Geophysical Journal International, 45, p125-154. [https://doi.org/10.1111/j.1365-246X.1976.tb00317.x]
• So, B.-D., Yuen, D.A., Regenauer-Lieb, K., and Lee, S.-M., (2012), Asymmetric lithospheric instability facilitated by shear modulus contrast: implications for shear zones , Geophysical Journal International, 190, p23-36. [https://doi.org/10.1111/j.1365-246X.2012.05473.x]
• Tamaki, K., Suyehiro, K., Allan, J., Ingle Jr, J.C., and Pisciotto, K.A., (1992), Tectonic synthesis and implications of Japan Sea ODP drilling, in Tamaki, K., Suyehiro, J., Allan, M., al., M.e. (Eds.), Proc. ODP, Sci Results, 127/128(Part 2), Ocean Drilling Program, College Station, TX, p1333-1348.
• Uyeda, S., and Kanamori, H., (1979), Back-arc opening and the mode of subduction, Journal of Geophysical Research: Solid Earth, 84, p1049-1061. [https://doi.org/10.1029/JB084iB03p01049]
• Yoon, S., Sohn, Y., and Chough, S., (2014), Tectonic, sedimentary, and volcanic evolution of a back-arc basin in the East Sea (Sea of Japan) , Marine Geology, 352, p70-88. [https://doi.org/10.1016/j.margeo.2014.03.004]