콤솔 멀티피직스를 활용한 지자기장 모델링 사례 연구

1. 서 론

지구 상에서 관찰되는 지구 자기장은 지구 심부 구조의 영향을 반영하는 동시에 천부 지각이 가지는 자화 특성의 변화 및 공간적 분포 양상에 의해 다르게 관찰된다. 지구 자기장은 지구 표면을 따라 그 세기와 방향이 달라지는 특성을 보여주고 있으며, 또한 시간에 따라 작게는 일별로 크게는 연별로 주기적인 변화를 보인다(Turcotte and Schubert, 2014). 주기적인 변화와 더불어 자기폭풍과 같은 급작스러운 자기장의 변동도 함께 일어난다. 지구 자기장은 자성을 가진 지질 매질의 분포를 조사하거나, 기반암의 깊이 추정, 암석의 생성 시기 및 위치, 혹은 오염 환경의 분포 등을 조사하는 데 많이 활용되어 왔다(Telford et al., 1976). 이 중 가장 중요한 발견은 지구 자기장의 역전 현상을 활용하여 정의된 해저 확장(seafloor spreading)과 판구조(plate tectonics) 운동의 발견이라고 할 수 있을 것이다(Vine and Matthews, 1963). 중앙해령을 따라 새롭게 생성되는 해양 지각은 생성 당시의 지구 자기장 방향을 기록한다. 그러므로 지구 자기장의 역전으로 인해 해양 지각에 기록된 자기 이상도 함께 역전된다. 하지만 해양 지각은 지구 자기장 이외에도 단층, 열수 순환 등 다양한 지질학적 기작에 의해 변형된다(Fujii et al., 2015; Szitkar and Dyment, 2015). 이러한 다중 지구물리학적 기작과 해양 지각에 기록된 자기장 역전띠(geomagnetic reversals)와의 상관관계는 아직 뚜렷하게 밝혀진바 없다. 더욱이 다중 지구물리학적 연구는 유한요소 기반의 수치모델링 방식의 접근을 많이 택하는 경향을 보이고 있으나, 자기 이상 모델링 기법은 여전히 전통적인 접근 방식에 머물러 있는 경향이 있다(Butler and Sinha, 2012). 이 연구에서는 유한 요소 기반의 자기 이상 모델링 기법을 콤솔 멀티피직스 환경에서 구현하고 그 결과의 유효성을 검증하고자 한다.

콤솔 멀티피직스(COMSOL Multiphysics®)는 1998년 스웨덴의 콤솔사에 의해 개발된 소프트웨어로서 유한요소법을 바탕으로 하여 편미분 방정식으로 표현된 다중물리현상을 모델링 할 수 있는 상용 프로그램이다. 이 프로그램은 전자기, 열, 구조, 유체등의 복합적인 현상이 일어나는 다중물리현상을 지원할 수 있는 다양한 어플리케이션을 제공하고 있으며, 2차원 및 3차원의 모델링이 가능하다. 그리고 사용자가 해석하고자 하는 문제에 집중할 수 있도록, 기본적으로 지원되는 어플리케이션 이외에 직접 식을 만들어서 모델링에 적용할 수 있는 기반도 갖추고 있다. 또한 유한요소법은 모델을 구성하는 요소(element)의 크기와 갯수에 따라 계산의 유효성이 영향을 받는데, 콤솔 멀티피직스에서는 메쉬(mesh) 설정을 통해 사용자의 모델에 맞는 크기와 모양의 메쉬를 제공하고 있다. 콤솔 멀티피직스의 장점 중 하나는 기존의 모델링 소프트웨어들과는 달리 다중 물리 현상 모델링 방법이 비교적 간단하다는 점이다. 모든 물리현상은 그래픽 유저 인터페이스(Graphic User Interface)를 통해 쉽게 지정할 수 있으며, 간단한 마우스 클릭으로 모델의 완성 및 계산까지 가능하다. 그와 동시에 그래픽 유저 인터페이스 환경 하에서 제공되는 지배방정식을 원하는 대로 수정하거나 새로운 지배방정식을 설정 할 수도 있다.

이 연구에서는 콤솔 멀티피직스 환경에서 주어진 자기 이상체에 의한 자기장을 계산하고 그 유효성을 검증하고자 한다. 또한 중앙해령의 단층 작용으로 자기 이상체가 분리된 경우 나타나는 자기 이상띠의 변화를 수치모델링을 통하여 살펴보고자 한다.

2. 실험 방법

시불변장(static)에서의 자기장을 설명하기 위한 맥스웰 방정식(Maxwell’s equation)은 다음의 두 식으로 정의된다.

여기서 H는 자기장 세기(magnetic field intensity), B는 자속밀도(magnetic flux density)를 의미한다. 전류의 흐름이 없는 경우에는 식 1은 으로 표현되며, 자속밀도와 자기장 사이의 관계를 이용하면 다음의 관계식을 정의할 수 있다.

위 식에서 μ는 투자율(permeability)이며, B_r은잔류자속밀도(remnant magnetic flux density)로 자기장을 없앤 뒤에도 남아있는 자화에 의한 자속밀도를 의미한다.

전류의 흐름이 없는 정자계(magnetostatic) 문제에서는 자기 스칼라 포텐셜(magnetic scalar potential)V_m으로 다음과 같이 자기장 H를 표현할 수 있다.

위에서 주어진 식 2와 3을 활용하면 식 4는 다음의 미분방정식을 정의할 수 있다.

여기에서 전체 자기 스칼라 포텐셜을 외부에서 추가된 량과 내부에서 감소한 량의 합으로 나타낼 수 있고, 또한 식 4의 정의에 의하여 식 5는 다음과 같이 정의 할 수 있다.

여기서 V_red는 내부에서 감소한 자기 스칼라 포텐셜을, H_ext은 외부 자기장을 의미한다.

콤솔 멀티피직스의 AC/DC 모듈에서는 위의 미분방정식이 미리 내장되어 있어 사용자는 외부자기장과 내부 자기장의 세기 및 자화 방향 만을 입력하여 지표에서의 자기장의 변화를 계산할 수 있다. 식 6에 주어진 미분 방정식은 해는 콤솔에 내장되어 있는 정적 선형 솔버(Stationary Linear Solver)를 적용하여 획득하였다. 이 연구에서는 콤솔 기반의 자기장 계산 결과의 유효성 검증(사례 1)과 해양지각의 3차원 지질구조 변화에 따른 지자기 역전띠의 변화 양상(사례 2)에 대한 수치모델링을 수행하였다.

2.1 사례 1: 자기장 계산 유효성 검증

콤솔 멀티피직스를 활용한 자기장 수치 모델링 결과의 유효성 검증을 위하여 구면체 형태의 자성체를 활용하였다. 지표면 아래에 자성을 가진 구면체를 각각 다른 깊이에 위치시킨 후 지표면에서 측정 되는 자기 이상을 위에 열거된 실험방법으로 콤솔 멀티피직스에서 계산하였다. 모델링을 위한 전체 도메인은 x, y, z 축 방향으로 각각 100, 200, 30 m의 길이를 가진다. 구면체의 반지름은 5 m이며 지표면 아래에 구면체의 중심이 10 m 그리고 두 개의 구면체가 고려된 경우에는 동일 반지름을 가지며 구면체의 중심이 지하 10 m와 15 m에 위치하도록 구성하였다. 유한요소법을 위한 요소의 크기는 계산의 속도와 정확도를 고려하여 계산에 크게 영향을 주지 않는 도메인에는 비교적 성긴 메쉬를 설정하고, 지표면의 측선 주변부와 자성체 주위에는 메쉬 분포(mesh distribution) 항목을 이용하여 더 조밀하게 설정함으로써 모델 결과의 해상도를 조절하였다. 계산에는 자유 사면체 메쉬(free tetrahedral mesh)를 사용하였으며 요소의 최대 크기와 최소 크기는 각각 16.5m와 1.2 m이며 총 요소는 38,000개이다. 계산에 사용된 지자기장 정보는 IGRF (International Geomagnetic Reference Field) 모델을 통해 얻은 대전지역의 값을 사용하였고 표 1에 나타내었다(Thébault et al., 2015). 자기 이상체는 30 A/m의 자화 강도를 가지고 지자기장 방향과 평행하게 자화 되었다고 가정하였다. 수치모델링을 통해 계산된 자기 이상은 동일한 모델에 대한 해석해를 활용하여 계산된 자기 이상 결과와 비교하였다.

Table 1.

Parameters to compute magnetic anomalies due to spheres.

구형의 모양을 가지는 자기 이상체에 의한 측정점 P (x₀, y₀, z₀)에서의 자기 이상은 구형 중심에 쌍극자(dipole)로 근사화 할 수 있으며 다음의 식으로 표현된다(Blakely, 1995).

$\[ \mathbf{B}=\frac{{\mu }_0}{4\pi }\frac{m}{r^3}\left[3\left(\widehat{\mathbf{m}}\bullet \widehat{\mathit{r}}\right)\widehat{\mathbf{r}}-\widehat{\mathbf{m}}\right],r\ne 0, \]$

7) 자기 이상 계산식

여기서 μ₀는 진공상태에서의 투자율(magnetic permeability), $$ \mathbf{m}=m\widehat{\mathbf{m}}$$으로 자화강도(magnetization intensity) m을 가지는 쌍극자 모멘트이며, $$ \mathbf{r}=r\widehat{\mathbf{r}}$$로 쌍극자 혹은 구체의 중심에서 관측점으로 향하는 위치 벡터로 주어진다. 쌍극자 모멘트 m은 유도자기장의 방향코사인을 이용하여 다음과 같이 표현할 수 있다.

m_x = m cos(mi) cos(md)
m_y = m cos(mi) cos(md)
m_z = m sin(mi)

여기서 mi는 유도자기장의 복각, md는 편각을 의미한다.

2.2 사례 2: 해양지각 자기이상 모델링

해양지각은 생성되는 시기의 지구 자기장 방향을 기록하므로, 시간에 따른 지구 자기장 역전 현상으로 인해 중앙 해령을 중심으로 대칭적인 자기 이상 패턴이 나타난다. 이러한 해양 지자기 역전띠를 활용하면 해양지각의 생성 시기, 이동 속도 및 방향을 계산할 수 있다. 또한, 해령에서는 해양지각의 확장으로 인해 정단층과 열극이 생성되고, 대부분의 경우 단층면은 확장 방향에 수직하게 형성되며 경사는 해령 축을 향한다(Macdonald et al., 1991). 단층이 생긴 초기에는 45°-60°의 고각도를 유지하며 변위가 활발히 발생하다가 지각이 확장하는 장력에 의해 최종적으로는 25° 내외의 저각도를 가지는 분리 단층(detachment fault)이 발생하여 해양지각 하부에 위치하는 하반이 해저면에 노출되기도 한다(Tucholke et al., 1998; Canales et al., 2004; Reston and Ranero, 2011). 이 연구에서는 콤솔 멀티피직스 상에서 해양지각 단층 모델을 구성하고, 단층 형태와 각도, 그리고 해양 조사선의 측선 방향에 따른 자기 이상 패턴의 차이점을 살펴보았다. 전체 계산 도메인의 크기는 x, y, z 축 방향으로 각각 10,000 m, 8,000 m, 2,000 m 이며, 수심 1,000 m 깊이에 두께가 1,000 m인 상부 해양지각을 설정하였다. 해령 축 부근에서 발생하는 정단층의 변화를 반영하기 위하여 경사각이 각각 90°인 모델과 60°, 25°인 모델을 제작하였다(그림 2a-c). 해양지각이 생성되는 동안 반복적으로 발생한 지자기 역전을 가정하기 위해 해령 축을 기준으로 각각 4개의 블록에 정자기와 역자기 방향을 번갈아 가며 입력하였다. 정자기 방향은 IGRF 모델을 이용하여 태평양-남극 해령의 현재 지자기 방향 값을 사용하였고, 계산에 사용된 자화 정보는 표 2에 제시하였다.

Table 2.

Parameters to compute magnetic anomalies by faulted oceanic ridge crest.

3. 실험 결과

3.1 자기장 계산 유효성 검증

콤솔 멀티피직스를 통해 지표면에서 측정된 자기 이상과 분석적 방법을 통해 계산된 결과를 비교하여 도시하였다(그림 1). 자성체는 지구 자기장과 평행하게 자화된 것으로 가정하였기 때문에 지표면에서 음과 양의 자기 이상치가 관찰되고 자성체의 깊이가 깊어짐에 따라 자기 이상치가 감소함을 그림 1e에서 확인 할 수 있다. 지표면에서 측정된 자기 이상치를 분석적 방법과 비교한 결과는 그림 1c와 1f에 도시하였다. 두 결과 모두 ±30 nT의 오차 범위 내에서 일치하는 것을 확인 할 수 있다. 최대 오차가 구면체의 중심부에서 발생하는 것으로 바탕으로 고려해볼 때 보다 정확한 계산을 위해서는 충분히 조밀한 메쉬 설정이 필요함을 알 수 있다. 자기 이상치는 자화 강도에 의한 세기 변화에도 민감하지만 자화 방향에 따른 자기 이상 곡선의 모양도 변화하게 된다. 그림 1에 나타난 구형의 자기 이상체에 의한 자기 이상치의 공간상 형태는 수치모델링 및 해석해 사이에 큰 차이가 없음을 알 수 있다. 이것은 콤솔 멀티피직스를 활용한 자력 탐사 해석 및 정량적 모델링이 충분히 가능함을 지시한다.

Fig. 1.

Magnetic anomalies from the COMSOL-based approach and analytic solution for uniformly magnetized spheres. (a) Model geometry for one sphere located at 10 m depth. Along profile A-A’, the computed magnetic anomalies are compared. (b) Magnetic anomalies along the profile A-A’. (c) The difference between the analytic and numerical results along profile A-A’. (d) Model geometry for two spheres placed at 10 and 15 m depths, respectively. Along profile B-B’, the computed magnetic anomalies are compared. (e) Magnetic anomalies along profile B-B’. (f) The difference between the analytic and numerical results along profile B-B’.

3.2 해양지각 자기이상 모델링

정단층의 경사각을 다르게 설정한 해령 자력탐사 모델의 해수면에서의 자기 이상 계산 값을 그림 2에 나타내었다. 세 모델 모두 해령을 중심으로 대칭적인 자기 이상 분포를 보였고 정단층의 경사각이 작아짐에 따라 해령 축 부근의 자기 이상치가 감소함을 확인할 수 있다(그림 2d-f). 일반적으로 중앙해령의 지자기장 역전 현상의 해석 시에는 동일한 자화 방향을 가지는 지질 구조가 수직 단면을 가지는 것으로 가정한다(그림 2a). 그러나 중앙해령 확장 속도에 따라 해저 지질구조는 다른 양상을 가질 수 있다. 단순한 수직 구조가 아닌 각기 다른 자화 방향을 가진 자성체가 서로 중첩되어 존재할 수도 있기 때문이다(그림 2b-c). 이 연구에서는 중앙 해령 주변부의 정단층 각도가 저각 일수록, 혹은 확장 속도가 느릴수록, 지자기 역전에 의한 세기는 감소할 수 있음을 간접적으로 확인하였다.

Fig. 2.

Synthetic faulted oceanic crusts used in COMSOL Multiphysics. (a-c) Geometries of the synthetic oceanic crust models with dip of 90°, 60° and 25°, respectively. (d-f) Magnetic anomalies at sea surface predicted by the corresponding models. (g) The difference between magnetic anomaly profiles along the vertical and diagonal direction (i.e., ship tracks) with respect to the ridge axis.

또한 해양 지자기 측정 방향(즉, 해양 조사선의 측선 진행 방향)에 따라 관찰되는 자기 이상치의 변화를 알아보기 위해 해령을 기준으로 수직(perpendicular) 방향과 대각선(diagonal) 방향으로 측정된 값을 비교하였다(그림 2g). 대각선 방향 측선의 경우에는 측정거리가 길어짐에 따라 자기 이상의 정점이 해령 축을 중심으로 점차 멀어졌지만 자기 이상의 크기는 수직방향 측선의 결과와 유사함을 확인 할 수 있었다.

4. 토의 및 결론

이 연구에서는 자력 탐사를 통해 얻게 되는 지구 자기장의 변화 특성을 유한요소 기반의 수치모델링을 통하여 정의할 수 있음을 살펴보았다. 구형의 자성체에 의한 자기 이상을 유한요소 기반의 수치 모델링을 통해 계산하였고, 그 결과를 기존에 알려진 해석해를 통해 계산된 결과와 상호 비교하여 수치 모델링 결과를 검증하였다. 또한 중앙해령 확장 환경에서 확장 속도에 따라 달라질 수 있는 단층의 경사각 차이를 고려한 자기 이상 역전띠를 모사하였다. 경사각이 낮을 수록 각기 다른 자화 방향을 가지는 해양 지각 구조가 중첩되어 자기 이상치를 감소시키는 결과를 유도함을 확인하였다. 또한 관측 측선의 방향에 따라 지자기 역전띠의 위치가 다르게 정의될 수 있으며, 이는 해양지각 확장 속도를 정의하는 데 오차의 원인으로 작용할 수 있음을 확인하였다. 그러므로 해양지각의 자기 이상 정밀하게 파악을 위해서는 지표에 드러난 구조 뿐 만 아니라 해양지각 내부의 단층면의 경사 방향, 경사각 또한 고려해야 함을 알 수 있다. 이러한 해양 지각의 내부 구조는 지자기장 뿐 아니라 다른 지질학적 기작 및 자료가 함께 사용되어야만 정량적으로 제한될 수 있다. 그리고 다중 지질 요소 사이의 상호작용은 콤솔 멀티피직스에서 수치모델로서 구현되어 다중지질 작용에 대한 모사 가능성과 재현성을 확보할 수 있을 것이라고 기대한다. 이 연구에서는 콤솔 기반의 자기 이상 수치 모델링이 기존의 자기 이상 모델링 기법보다 우월하다는 것을 보이는 것이 목적이 아니다. 콤솔 기반의 자기 이상 수치모델링이 가능할 경우 다른 물리적 현상과 함께 복합 연산이 가능하게 되며 이를 통하여 다른 물리 인자들 사이의 상관관계를 보다 정량화 할 수 있으므로 그 적용성을 이 연구에서는 검증하여 본 것이다.

Acknowledgments

이 논문에 대하여 세심한 수정과 의견을 제공해주신 이창열, 소병달 심사위원님들께 감사드립니다. 이 연구는 2018년도 충남대학교 학술연구비에 의해 지원되었습니다.

References

• Blakely, R., (1995), Potential Thoery in Gravity and Magnetic Applications, Cambridge University Press, Cambridge, UK.
• Butler, S.L., and Sinha, G., (2012), Forward modeling of applied geophysics methods using Comsol and comparison with analytical and laboratory analog models, Comput. Geosci, 42, p168-176. [https://doi.org/10.1016/j.cageo.2011.08.022]
• Canales, J.P., Tucholke, B.E., and Collins, J.A., (2004), Seismic reflection imaging of an oceanic detachment fault: Atlantis megamullion (Mid-Atlantic Ridge, 30°10′N), Earth Planet, Sci. Lett, 222(2), p543-560. [https://doi.org/10.1016/j.epsl.2004.02.023]
• Fujii, M., Okino, K., Honsho, C., Dyment, J., Szitkar, F., Mochizuki, N., and Asada, M., (2015), High-resolution magnetic signature of active hydrothermal systems in the back-arc spreading region of the southern Mariana Trough, J. Geophys. Res. Solid Earth, 120(5), p2821-2837. [https://doi.org/10.1002/2014JB011714]
• Macdonald, K.C., Scheirer, D.S., and Carbotte, S.M., (1991), Mid-Ocean Ridges: Discontinuities, Segments and Giant Cracks, Science, 253(5023), p986-994. [https://doi.org/10.1126/science.253.5023.986]
• Reston, T.J., and Ranero, C.R., (2011), The 3-D geometry of detachment faulting at mid-ocean ridges, Geochemistry, Geophys., Geosystems, 12(7), pQ0AG05. [https://doi.org/10.1029/2011GC003666]
• Szitkar, F., and Dyment, J., (2015), Near-seafloor magnetics reveal tectonic rotation and deep structure at the TAG (Trans-Atlantic Geotraverse) hydrothermal site (Mid-Atlantic Ridge, 26°N), Geology, 43(1), p87-90. [https://doi.org/10.1130/G36086.1]
• Telford, W.M., Geldart, L.P., Sheriff, R.E., and Keys, D.A., (1976), Applied Geophysics, Cambridge University Press, Cambridge, UK.
• Thébault, E. et al., (2015), International Geomagnetic Reference Field: the 12th generation, Earth, Planets Sp, 67(1), p79. [https://doi.org/10.1186/s40623-015-0228-9]
• Tucholke, B.E., Lin, J., and Kleinrock, M.C., (1998), Megamullions and mullion structure defining oceanic metamorphic core complexes on the Mid-Atlantic Ridge, J. Geophys. Res. Solid Earth, 103(B5). [https://doi.org/10.1029/98JB00167]
• Turcotte, D.L., and Schubert, G., (2014), Geodynamics, 3rd ed, Cambridge University Press, Cambridge, UK.
• Vine, F.J., and Matthews, D.H., (1963), Magnetic Anomalies Over Oceanic Ridges, Nature, 199(4897), p947-949. [https://doi.org/10.1038/199947a0]