1. 서 론

최근 심부 지하 활용에 대한 수요가 전 세계적으로 증가하면서 심부 암반에서 발생하는 열-수리-역학적 연계거동에 대한 연구 필요성이 증대되고 있다. 심부 지층은 천부 지층에 비해 고온, 고압의 극한환경에 노출되어 있을 뿐만 아니라, 제한된 접근성으로 인한 불확실성이 높기 때문에 그 거동 특성을 적절히 이해하고 예측할 수 있는 신뢰성 높은 수치해석 모델의 개발은 매우 중요하다. 특히, 국내와 같이 저투수성 또는 불투수성 결정질 암석이 많이 분포하는 환경에서는 단층이나 절리와 같은 불연속면이 암반의 열-수리-역학적 거동 특성에 지배적인 영향을 미치게 되므로, 불연속면의 특성을 적절히 반영할 수 있는 수치모델의 개발은 우선적으로 해결되어야 할 문제 중 하나이다.

암반 불연속면에서 발생하는 열-수리-역학적 연계거동은 주로 지층 내로의 유체 주입이 동반되는 이산화탄소 지층처분(Carbon Capture and Storage, CCS), 지열발전(Enhancd Geothermal Energy, EGS), 원유의 회수증진(Enhnaced Oil Recovery, EOR) 분야에서 주변 단층의 재활성화와 유발 지진 문제와 관련해 많은 연구가 이루어져 왔다. 단층 재활성이란 유체 주입에 의한 압력 증가와 유효수직응력 감소로 인해 단층의 전단강도가 감소하여 역학적 미끄러짐이 발생하는 현상으로 관련 시설의 안정성을 저해할 뿐만 아니라 주변 지역에 인공지진을 유발할 수 있으므로 프로젝트의 성패를 가르는 중요한 설계인자라고 할 수 있다(Rutqvist et al., 2007, Cladoushos et al., 2009, Rinaldi and Rutqvist, 2019). 방사성폐기물 처분장과 관련하여서도, 단층 재활성은 수리적 요인뿐만 아니라 열적 요인에 의한 처분장 규모의 열팽창으로 인해 발생할 수 있는 것으로 알려져 있으며(Ghabezloo et al., 2009), 핵종의 이동 통로와 저류장의 역할을 할 수 있어 부지선정과 운영 단계에서 단층과 균열의 거동을 적절히 이해하고 예측하는 일은 처분장의 안정성과 안전성 확보를 위해 무엇보다 중요한 연구 테마이다.

본 연구는 유체 주입에 따른 단층의 미끄러짐과 수리역학적 거동 특성을 모사할 수 있는 해석기법을 개발하기 위한 연구로, DECOVALEX(DEvelopment of COupled models and their VALidation against EXperiments)-2019 프로젝트의 일환으로 수행되었다. DECOVALEX 프로젝트는 고준위 방사성폐기물의 지층처분을 고려하는 세계 각국의 원자력 관련 기관들이 참여하는 국제공동연구로서 1992년에 착수되어 4년을 주기로 새로운 주제에 대한 연구들이 수행되고 있다. 현재 7가지 연구 주제(Task)에 대해 DECOVALEX-2023이 진행중에 있으며, 방사성 폐기물 처분장에서의 공학적 방벽과 지하 암반의 열-수리-역학-화학적(thermal-hydrological-mechanical-chemical, THMC) 연계거동 해석기법 개발을 목표로 한다. 본 논문에서는 2016년에서 2019년까지 DECOVALEX-2019 Task B를 통해 수행된 국내 참여기관(한국지질자원연구원)의 연구내용을 소개하였다. 이와 관련된 연구내용은 Park et al.(2018a), Park et al.(2018b), Park et al.(2019), Park et al.(2020)의 논문을 통해 보고된 바 있으며, 여기에서는 Task B의 전반적인 연구결과와 함께 일부 업데이트된 해석결과들을 소개하고, 개발된 해석기법의 의의와 적용성에 대해 고찰하였다.

2. DECOVALEX-2019 Task B 개요

DECOVALEX-2019 Task B(Fault Slip Test: Modelling the induced slip of a fault in argillaceous rock)는 스위스의 방사성폐기물 처분장 부지로 고려되고 있는 불투수성(또는 저투수성) 점토질 암반에서 단층의 수리역학적인 거동을 모사할 수 있는 수치해석기법을 개발하는 데에 목표가 있다. Task B는 Swiss Federal Nuclear Safety Inspectorate(ENSI)이 제안하였으며, 총 7개국(독일, 캐나다, 스위스, 대만, 한국, 미국, 스웨덴)에서 7개의 연구그룹이 참여하였다: Federal Institute for Geosciences and Natural Resources/Helmholtz Centre for Environmental Research(BGR/UFZ, 독일), Canadian Nuclear Safety Commission (CNSC, 캐나다), ENSI(스위스), Sinotech Engineering Consultants/Institute of Nuclear Energy Research(SINOTECH/INER, 대만), 한국지질자원연구원/한국원자력연구원(KIGAM/KAERI, 한국), Lawrence Berkeley National Laboratory/Department of Energy(LBNL/DOE, 미국), DynaFrax/Swedish Radiation Safety Authority(DynaFrax/SSM, 스웨덴).

Task B의 추진일정은 크게 세 단계의 연구로 구성된다(Fig. 1). 1단계는 벤치마크 해석 연구로, 각 연구팀은 독자적인 수치해석기법을 통해 동일한 문제를 모델링하게 되며 이 과정에서 해석결과의 상호 비교를 통해 지속적으로 코드를 개선하게 된다. 2단계와 3단계는 각각 스위스 Mont Terri 지하연구시설의 단층 손상대(fault damage zone)와 단층핵(fault core)에서 수행된 물 주입시험을 모델링하는 심화연구로, 1단계에서 개발된 해석코드를 통해 현장에서 취득한 수리역학적 거동을 모델링하게 된다. 배포된 현장시험 자료는 주입공에서의 주입율 및 주입압, 모니터링 지점에서의 압력변화, 변위 곡선 등이며 이를 적절히 재현함으로써 단층 재활성과 관련된 수리 역학적 메커니즘을 평가할 수 있는 수치모델을 개발하게 된다.

Table 1은 Task B 참가팀들의 수치해석코드와 단층 모델링 방법을 정리한 것이다. 수치해석코드의 특성에 따라 단층을 일정 부피(또는 두께)를 갖는 연속체모델(solid element)이나 불연속면(interface element)으로 가정하는 방법이 사용되었고, 입자결합모델(particle flow code, PFC)에서는 약한 결합력을 갖는 입자의 그룹을 통해 단층이 모델링되었다. 한국지질자원연구원의 경우, 연속체모델과 불연속면으로 모델링하는 두 방식을 모두 사용하였으나, 본 논문에서는 불연속면으로 모델링하는 경우에 대해서만 논의하였다. Fig. 2는 2단계 연구에 사용된 각 팀의 해석 메쉬를 보여주는 것이다.

Fig. 1.

DECOVALEX-2019 Task B modeling steps

Fig. 2.

Mesh discretization applied by the different teams for Step 2 modeling (after Graupner et al., 2020)

3. 벤치마크 해석을 통한 해석모델 개발

3.1 해석 개요

벤치마크 해석에서는 단일 불연속면(단층)을 포함하는 정육면체(암반) 해석모델에 단계적 압력의 물을 주입하는 과정을 모델링한다. 여기에서 암반은 불투수성으로 간주되고, 주입수의 유동은 단층을 통해서만 이루어지는 것으로 가정된다. 단층 내 유체의 흐름은 Darcy의 법칙을 따르며, 유량과 수리간극의 관계가 평행한 두 평판사이의 유체 흐름 모델인 삼승법칙(Witherspoon, 1980)을 따른다. Fig. 3은 해석모델의 형상과 적용된 주입압을 보여주는 것으로 주입압은 Mont Terri 지하연구시설 현장시험에 적용된 조건과 동일하다. 암반은 20×20×20 m³의 정육면체 형태로 가정되고 중심에 65° 경사의 단층이 관통하며, 단층의 주향은 X축과 동일하다. 물은 단층 중심(모델 중심)을 교차하는 반경 0.07 m의 시추공을 통해 주입되며, 주입수의 압력은 9단계에 걸쳐 변화한다. 각 주입단계에서 수십 초 동안 일정한 압력이 유지되며, 여덟 번째 단계까지 압력이 최대 6.302 MPa까지 증가하다가 마지막 단계에서는 354초 동안 3.382 MPa의 낮은 주입압이 적용된다.

Fig. 3.

Geometric information of benchmark simulations and injection pressure steps applied at the injection point

단층 내에서 발생하는 수리역학적 커플링을 구현하기 위한 관계식으로서 두 가지 수리간극모델(FM1, FM2)이 제안되었다. FM1과 FM2의 가장 큰 차이는, 전자의 경우 초기의 단층을 수리적으로 닫힌 균열(불투수성)로 간주하는 반면, FM2에서는 이미 일정한 간극을 가지고 있는 열린 균열(투수성)로 모델링한다는 점이다. FM1 수리간극모델은 Guglielmi et al.(2015)이 프랑스 Tournemire 지하연구시설에서 수행한 현장시험 결과를 바탕으로 제안한 것으로, 전단파괴(hydroshearing)나 인장파괴(tensile opening)가 발생하기 전까지 균열면이 닫혀 있는 상태로 가정한다. 즉, 일정 압력에 도달하여 파괴가 발생하기 이전에는 단층의 수리간극은 0이며, 파괴가 발생하면 일시에 일정량의 수직팽창이 발생하여 수리간극으로 작용하는 모델이다. 본 논문에서는 이를 ‘creation aperture’로 표기하였다. 수직팽창량은 경험적으로 산출하게 되며, 벤치마크 해석에서는 28 μm가 적용되었다. 벤치마크 모델에서는 주입공 주변으로 반경 0.5 m 거리에서 시추공 굴착 및 교란으로 인해 파괴 상태에 있는 것으로 가정하며, 열린 균열로서 28 μm의 수리간극을 갖는다. 이를 수식으로 나타내면 다음과 같다.

여기서 b_h는 수리간극, Δb_he는 탄성적 수직변위이며, Δb_hc는 파괴로 인한 수직변위(creation aperture), r은 주입지점에서의 거리, r_f는 초기 파괴구간의 반경이다.

탄성변형인 Δb_he는 단층의 수직강성 K_n과 단층에 작용하는 유효 수직응력 변화 Δσ_n′를 이용해 다음과 같이 계산되며, FM1의 경우 열린 균열에서만 탄성변형이 계산된다.

FM2는 암석역학 및 암반공학 분야에서 보다 널리 사용된 모델로서, 수리간극의 크기는 다음과 같이 정의된다. 단, 수리간극과 역학적 간극의 크기는 동일한 것으로 가정한다.

여기서 b_hi는 초기수리간극, Δb_he는 탄성변형, Δb_hs는 전단파괴 이후 전단변위와 함께 단층면의 거칠기에 의해 수직방향으로 발생하는 수직팽창을 의미한다.

Δb_hs는 전단변위 Δu_s와 팽창각 ψ를 이용하면 다음과 같이 수식화할 수 있다.

해석모델에 적용된 현지응력 조건은 σ_x = 3.3 MPa, σ_y = 6.0 MPa, σ_z = 7.0 MPa이며, 초기압력은 0.5 MPa이다. 해석모델의 경계에서는 각 면의 수직방향으로 변위가 구속되고 0.5 MPa의 압력이 일정하게 유지된다. Jaeger et al.(2007)의 이론식에 따르면 단층에 작용하는 초기 수직응력과 전단응력은 각각 6.179 MPa, 0.389 MPa로 계산된다. Table 2는 벤치마크 해석에 사용된 입력파라미터를 정리한 것으로 단층, 모암, 유체의 물성이 포함되어 있다.

3.2 TOUGH-FLAC 인터페이스 연동해석 모델

본 연구에서는 단층의 수리역학적 상호작용을 모델링하기 위하여 TOUGH-FLAC 연동해석 기법을 이용하였다. TOUGH-FLAC 연동해석은 Rutqvist et al.(2002)이 제안한 해석기법으로 현재까지 다양한 암반공학 문제에 적용되어 왔다. 위 해석기법은 다상 다성분 열수리유동 해석 코드인 TOUGH2와 암반의 역학해석 프로그램인 FLAC3D를 연동하여 동일한 해석모델(메쉬)를 대상으로 TOUGH2의 비선형 반복계산과 FLAC3D의 준정적해석을 순차적으로 수행함으로써 열-수리-역학적 연계거동 특성을 해석하는 기법이다. 이 과정에서 FLAC3D는 TOUGH2에서 계산된 열-수리적 해석 결과(온도, 압력 등)를 입력받아 역학적 거동(응력, 변위 등) 해석을 수행하고, TOUGH2는 FLAC3D의 결과를 바탕으로 입력물성이나 해석 조건을 업데이트하여 다음 time step 또는 Newton iteration에 대한 열수리 유동해석을 실시하게 된다.

TOUGH-FLAC 시뮬레이터를 통해 단층의 수리역학적 거동을 모사한 선행연구들에서는 대부분 단층을 연속체 모델로 가정하여 왔다(Rinaldi et al., 2014, Rutqvist et al., 2015). 일반적인 암반공학 스케일에서 수리간극의 크기는 단층의 폭과 연장에 비해 매우 작기 때문에 불량한 종횡비 문제를 해결하기 위해 수리간극에 비해 상당히 큰 해석요소 두께를 가정한다. 단층 내 수리유동은 TOUGH2에서 다공성매질 내의 흐름으로 간주되고, 수리간극의 크기와 해석모델의 두께를 고려해 이론적 수리물성이 부여된다. FLAC3D의 역학해석 모델에서도 불연속면의 물성(강성, 전단강도 등)을 연속체 모델의 물성(탄성계수, 내부마찰각, 점착력 등)으로 변환하여 사용한다. 이러한 접근법은 단층을 두께가 없는 불연속면으로 모델링하는 경우에 비해 모델링이 간편하다는 장점이 있으나(Cappa and Rutqvist, 2011), 입력물성 산정과 해석결과 처리 과정에서 상당히 많은 가정이 중첩될 수 있고, 본 연구와 같이 물 주입으로 인해 급격한 수리간극의 변화가 예상되는 문제에서 쌍방향 수리역학 커플링(two-way hydro-mechanical coupling)의 효과를 적절히 예측하는 데에 한계가 따른다.

따라서 본 연구에서는 가상의 단층 두께를 가정하지 않고, 단층의 거동을 직접적으로 재현하기 위한 새로운 접근법을 제시하였다. 적분유한차분법(integral finite difference method)을 이용하는 TOUGH2의 장점을 활용하여 실제 수리간극에 상응하는 매우 얇은 해석 요소를 단층의 유체 유동경로로 모델링하는 한편, FLAC3D에서는 인터페이스 요소(interface element)를 이용하여 단층의 전단파괴(slip)와 인장파괴(separation)를 직접적으로 재현하고자 하였다. FLAC3D의 인터페이스 요소는 연속체모델의 단점을 보완하기 위한 기능으로, 이를 통해 미끄러짐이나 개폐가 발생할 수 있는 불연속면(절리, 단층)이나 상이한 재료의 경계면의 거동을 모사할 수 있다. 연속체 요소(zone 요소)들 사이의 경계면(target surface)에 인터페이스를 할당하면, 세 개의 절점(interface node)을 갖는 삼각형 인터페이스가 정의되며, 경계면과 인터페이스 절점간의 접촉(contact)에서 불연속면의 거동(힘과 변형)을 결정하는 구성방정식들을 정의할 수 있다. 또한, 전단 또는 인장방향의 강도정수나 물성을 할당함으로써 미끄러짐(slip)과 분리(separation)의 발생을 직접적으로 구현할 수 있다.

Fig. 4는 역학해석과 수리유동해석 사용된 해석메쉬를 보여준다. FLAC3D에서 암반과 단층은 각각 탄성모델(zone 요소), 탄소성모델(interface 요소)로 모델링되었다. 역학적 거동해석에서는 암반의 탄성 거동이 단층의 응력과 변위에 복합적 영향을 미치므로 암반과 단층의 형상을 모두 고려하였으며, 수리해석에서는 암반을 불수투성으로 가정하므로 단층 내에서의 수리유동만을 해석하였다. TOUGH2에서 단층은 공극률이 1.0이고 수리간극에 상응하는 두께를 갖는 부피 요소로 모델링된다.

Fig. 4.

Numerical models used in FLAC3D and TOUGH2

본 연구에서는 TOUGH2의 수리유동해석을 위해 단층 내 유체의 흐름을 두께가 b_h이고 공극률인 1.0인 대수층 내 2차원 평면유동으로 단순화하고, 대수층 유동이론을 적용하여 단층의 투수량계수(transmissivity)와 저류계수(storativity)를 다음과 같이 수식화하였다.

여기서 T는 투수량계수, ρ_f는 유체의 밀도, g는 중력가속도이며, µ는 점성계수(dynamic viscosity)이다.

여기서 K_n은 단층의 수직강성, β는 물의 압축률이다.

식 (5)과 (6)에 의하면, 단층의 투수계수(k_f)와 압축률(α_f)은 다음과 같이 표현된다(Park et al., 2018a)

단층의 수리특성을 결정하는 주요인자인 저류계수와 투수량계수는 모두 수리간극의 함수로 표현되며, TOUGH2 코드에서는 위 인자들의 영향을 압축률과 투수계수 변화를 통해 반영할 수 있다. 수리간극모델(FM1, FM2)의 특징에 따라 단층 해석 요소의 두께를 모델링하였고, 입력물성인 투수계수와 압축률은 식 (7)과 식 (8)에 따라 결정하였다.

Fig. 5는 TOUGH2와 FLAC3D 프로그램 간의 연동해석 과정을 보여주는 것으로, 두 프로그램의 호출(실행)과 데이터 교환은 FLAC3D의 인터페이스 요소와 TOUGH2의 전체 요소 사이에서 매 time step마다 순차적, 반복적으로 이루어진다. 먼저 TOUGH2에서 수리해석을 통해 각 요소 부피 중심에서의 압력을 계산하여 FLAC3D로 전송하면, 역거리가중법(inverse-distance weighting)에 따라 모든 인터페이스 절점의 압력이 환산된다(Fig. 6). FLAC3D는 유효응력해석을 통해 각 절점에서의 수직/전단방향의 응력과 변위를 계산하고 인터페이스의 파괴 여부를 결정한다. 각 절점에서 수리간극 크기가 계산되면, 다시 역거리가중법에 의해 TOUGH2 요소 중심에서의 간극 크기가 결정되어 다음 time step의 TOUGH2 해석에 반영된다. 이 과정에서 식 (7)과 (8)을 통해 해석요소의 투수계수와 압축률도 업데이트된다.

Fig. 5.

Hydro-mechanical coupling process and data transfer between volume elements of TOUGH2 and interface elements of FLAC3D (after Park et al., 2020)

Fig. 6.

Inverse distance weighting method for data transfer between volume elements of TOUGH2 and interface elements of FLAC3D

본 연구에서는 TOUGH2 계산 시 단층 요소의 두께를 실제 수리간극과 동일하게 모델링하였으므로 해석 메쉬도 지속적으로 변화하게 된다. 매 time step마다 단층면과 인터페이스 절점의 변위 정보를 이용해 3차원 해석메쉬를 재구성하였으며, 해석메쉬의 부피 정보와 연결성(connectivity)을 업데이트하였다. 이를 위해 FLAC3D의 해석 결과를 토대로 TOUGH2 요소들의 메쉬 정보와 투수계수, 압축률을 동시에 업데이트할 수 있는 별도의 subroutine을 개발하였다.

3.3 해석 결과 및 논의

본 절에서는 두 가지 수리간극모델에 대한 벤치마크 해석결과와 함께 수리역학 커플링 조건에 따른 해석결과(FM2 모델)를 살펴봄으로써 수리간극모델과 커플링 인자의 영향을 고찰하였다. 벤치마크 해석결과와 관련하여 Park et al.(2018b)의 논문을 통해 상세한 연구내용이 소개된 바 있으나, 그들의 연구에서는 수리역학 커플링 시 단층의 압축률 변화(식 (8))가 반영되지 않았으며, FM1 수리간극모델에서 파괴 이전에 탄성적 수직변위를 허용하였으므로 FM1 모델과 FM2 모델이 거의 유사한 거동을 보이는 것으로 해석되었다. 여기에서는 압축률과 수정된 FM1 수리간극모델을 반영하여 업데이트된 해석결과를 기술하고, 커플링 인자들의 영향을 함께 살펴보았다.

3.3.1 수리간극모델에 따른 해석결과

Fig. 7은 두 수리간극모델의 해석결과로 주입지점(P1)의 압력 증가에 따른 유효수직응력, 전단응력, 전단강도의 변화를 보여준다. 여기서 전단강도는 유효수직응력 해석 결과와 단층의 강도(마찰각, 점착력)를 이용해 이론적으로 계산된 Coulomb sliding 모델의 전단강도이다. 압력의 증가에 따라 유효수직응력과 전단강도가 감소하였고, 두 수리간극모델 모두 6.302 MPa의 최대 압력이 적용된 420초에 주입공 주변에서 전단파괴가 개시되었다. FM1 모델의 경우, 초기 파쇄 영역(주입공으로부터 반경 0.5 m) 이내에서만 전단파괴가 관찰될 뿐 파쇄 영역이 확장되지는 않았으며, 두 모델에서 모두 전단파괴를 유발하는 이론적인 최소압력(critical injection pressure, P_crit)에 비해 높은 주입압 조건에서 전단파괴가 유도되었다. Park et al.(2018b)은 높은 전단강도가 발현되는 원인을, 단층의 탄성거동과 수직변위의 국부적 구속에 따라 유효 수직응력의 증가분이 발생하기 때문인 것으로 분석하였다. 즉, 압력의 증가량에 비하여 유효수직응력의 감소량이 적고, 단층의 전단강도는 이론적인 수치보다 높게 해석되었다.

Fig. 7.

Variations of injection pressure, total normal stress, effective normal stress, shear stress and shear strength of fault

Fig. 8은 최대 주입압의 적용이 종료된 453초에 예측된 전단파괴 범위를 보여주는 그림이며, Fig. 9는 420, 453, 807초에 예측된 주입공과 암반의 변위를 보여준다. FM1 모델에서는 소규모 영역에서 작은 전단변위가 예측되며, 주입공 주변의 국부적인 영역에서만 암반의 변위가 발생하는 것으로 계산되었다. 반면, FM2의 경우에는 압력의 발달과 암반의 변위가 모두 단층면의 넓은 범위에 걸쳐 유도되는 것으로 나타났다.

Fig. 8.

Shear failure zone estimated at 453 s of water injection under injection pressure of 6.302 MPa

Fig. 9.

Contour of displacement; the scaled arrow denotes the direction and magnitude of fault displacement at P1

Fig. 10은 단층의 주향방향을 따라 주입공을 교차하는 프로파일 상에 시간에 따른 수리간극의 크기 변화를 나타낸 것이다. FM1 모델에서는 주입압의 증가에 따라 수리간극의 크기도 비례하여 증가하였다. 그러나 초기 파쇄구간 이외에 새로운 균열이 확장되지는 않았고, 탄성적 변형에 의한 수리간극 변화만이 관찰되었다. FM2 모델의 경우, 단층 전체를 초기수리간극을 갖는 열린 균열로 가정하므로 시간의 흐름에 따라 전 단층면에 걸쳐 압력이 발달하게 되며, 압력과 수리간극 크기분포는 주입공을 중심으로 하는 동심원 형태를 갖는다.

Fig. 10.

Hydraulic aperture evolutions of FM1 and FM2 (σ_x = 3.3 MPa, σ_y = 6.0 MPa, σ_z = 7.0 MPa)

Fig. 11은 모니터링 지점의 압력 변화를 보여주는 것이며, 모니터링 지점인 P2와 P3은 각각 단층의 주향방향과 경사방향으로 주입공과 1.5 m 이격거리에 위치한다. Fig. 12는 주입공에서 구한 주입율 변화 곡선이다. Fig. 11에는 FM2 수리간극모델의 결과만을 도시하였는데 이는 FM1 모델에서 해석이 종료될 때까지 새로운 균열이 확장되지 않아 모니터링 지점에서 수리유동이나 압력 발달이 이루어지지 않았기 때문이다. 주입율 곡선에서도 확인할 수 있듯이, FM1 모델에서는 초기파쇄구간의 제한된 부피로 인해 각 주입의 초기단계에서만 일시적으로 주입율이 증가할 뿐 수초 이내에 주입율이 0으로 수렴하는 결과를 보였다. FM2 모델에서는 초기수리간극과 탄성적 변형으로 인해 압력구배가 지속적으로 발달하며 각 주입단계에서 시간이 경과하더라도 주입율이 0으로 수렴하지 않고 일정 수준을 유지하게 된다.

Fig. 11.

Variations of pressures monitored at monitoring points, P2 and P3 (FM2) (σ_x = 3.3 MPa, σ_y = 6.0 MPa, σ_z = 7.0 MPa)

Fig. 12.

Variations of injection flow rates (σ_x = 3.3 MPa, σ_y = 6.0 MPa, σ_z = 7.0 MPa)

상기한 해석결과에서는 FM1 모델의 단층이 파괴조건에 도달하지 않아 균열의 확장에 따른 수리역학적 거동을 면밀히 살펴보는데 제약이 따른다. Park et al. (2020)의 논문에서는 벤치마크 모델의 초기응력 조건을 조정하여(σ_x = 3.3 MPa, σ_y = 5.3 MPa, σ_z = 7.0 MPa) 낮은 주입압 조건에서 전단파괴가 개시되도록 유도하였다. Fig. 13, Fig. 14, Fig. 15는 각각 Fig. 10, Fig 11, Fig 12에 대응하는 해석결과로서 두 수리간극모델의 거동 특성의 차이점을 보다 효과적으로 제시한다.

Fig. 13에서 FM1 모델의 수리간극 변화를 Fig. 10과 비교하여 보면, 주입 초반에는 초기파쇄영역 구간의 탄성적 변형에 의해 반경 0.5 m 이내에서만 수리간극이 증가한다. Fig. 13에서는 주입압의 증가로 인해 단층이 파괴조건에 도달하면서 새로운 균열이 생성되기 시작하고, 453초에는 파쇄구간이 주입공으로부터 약 1.8 m까지 확장되는 것을 확인할 수 있다. 이에 따라 420초를 전후하여 모니터링 지점 P2와 P3의 압력이 순간적으로 급격히 증가하며, 주입율 역시 균열의 성장과 함께 높은 수준을 유지하게 된다(Fig. 14, Fig. 15).

Fig. 13.

Hydraulic aperture evolutions of FM1 and FM2 (σ_x = 3.3 MPa, σ_y = 5.3 MPa, σ_z = 7.0 MPa) (Park et al., 2020)

Fig. 14.

Variations of pressures monitored at monitoring points, P2 and P3 (σ_x = 3.3 MPa, σ_y = 5.3 MPa, σ_z = 7.0 MPa) (Park et al., 2020)

Fig. 15.

Variations of injection flow rates (σ_x = 3.3 MPa, σ_y = 5.3 MPa, σ_z = 7.0 MPa) (Park et al., 2020)

3.3.2 수리역학 커플링 조건에 따른 해석

본 연구에서는 수리역학 커플링 모듈을 통해 수리간극의 변화에 따른 투수계수, 압축률, 해석메쉬의 변화를 업데이트하였다. 개발된 해석모델이 수리역학적 커플링 프로세스를 적절히 재현하는지 검토하기 위하여 FM2 수리간극모델에 대하여 Table 3과 같은 해석 조건을 설정하여 커플링 조건에 따른 해석결과를 비교, 검토하였다. 여기에서는 1.0 MPa의 일정한 압력의 물을 800초 동안 주입되는 것으로 가정하였으며, 낮은 수준의 주입압 적용으로 인해 전단파괴가 유도되지는 않는다.

Case 1의 경우 단층의 변형으로 인한 수리간극의 변화를 고려하지 않는 경우로 역학적 요인에 의해 수리특성이 변화하지는 않지만, 수리해석을 통해 계산된 압력을 이용하여 유효응력해석을 수행하는 일방향 연동해석(one-way hydro-mechanical coupling)이다. Case 2와 Case 3은 단층의 변형을 고려하는 쌍방향 연동해석(two-way hydro-mechanical coupling)으로 Case 2에서는 저류계수(storativity)의 변화만을, Case 3에서는 저류계수와 투수량계수의 변화를 모두 업데이트한다. 따라서 Case 1과 Case 2의 결과 비교를 통해 단층의 저류계수 변화가 해석결과에 미치는 영향을 검토할 수 있으며, Case 2와 Case 3의 결과 비교를 통해 투수량계수의 변화가 해석결과에 미치는 영향을 검토할 수 있다. 본 연구에서는 수리해석 시 단층 요소의 두께를 수리간극과 동일하게 모델링하였으므로 저류계수의 변화를 고려하는 쌍방향 연동해석 시 해석메쉬도 함께 업데이트된다. 기타 해석모델의 형상과 물성은 벤치마크 모델과 동일하다.

Fig. 16.

Results of Cases 1, 2 and 3 in Table 3; Pressure distributions estimated at 200s and 800 s of water injection

Fig. 17.

Results of Cases 1, 2 and 3 in Table 3; Variations of injection flow rate

Fig. 18.

Results of Cases 1, 2 and 3; Variations of effective normal stress estimated at injection point P1

Fig. 16은 주입 후 후 200초와 800초 경과 시 단층면의 압력분포를 보여주며, Fig. 17과 Fig. 18은 각각 주입지점 P1에서 예측된 주입율과 유효응력 변화를 보여준다. 위 결과를 종합적으로 살펴보면, 수리간극의 변화가 고려되지 않는 Case 1의 경우, 주입 후 매우 짧은 시간 이내에 정상류 상태(steady state flow)에 수렴하는 것으로 나타났다. Case 1의 경우 200초가 경과하기 전 이미 평형상태에 도달하여 800초에 도달할 때까지 더 이상 압력분포의 변화가 발생하지 않으며, 물 주입율은 수 초 이내에 일정한 값으로 수렴하였다. 반면, 단층의 변형으로 인한 압축률(또는 저류계수)의 변화를 고려한 Case 2의 경우 Case 1에 비하여 압력 전달 속도가 현저히 느리고 주입공 주변에서 압력의 소산이 이루어지지 않고 있음을 확인할 수 있다. 이에 따라 압력수두 경사가 크게 발달하고 주입율은 높게 유지된다. 이는 개발된 수치해석모델이 단층 압축률의 영향을 적절히 반영하고 있음을 시사한다. Case 2에서는 높은 압축률로 인해 주입공 주변 압력이 증가하면 수리간극의 부피가 팽창하고, 이에 따라 주입수가 원거리로 전달되기보다는 주입공 주변에 저장되는 현상이 나타났다. 식 (6)에 따르면 Case 1과 같이 단층의 변형을 무시하는 경우(K_n = ∞), 저류계수는 약 4.4×10^-11이며, Case 2의 경우(K_n = 20 GPa/m)의 초기 조건의 저류계수는 약 5.0×10^-7로 10,000배 이상의 저류계수를 갖는다. 다시 말해 Case 2는 이론적으로 Case 1에 비해 극도로 느린 압력 전달 속도를 갖게 된다. 하지만 투수량계수가 일정하므로 오랜 시간이 흐른 이후에는 Case 1과 동일한 압력 분포를 갖게 되며, 주입율 역시 Case 1의 수렴치에 도달하게 될 것이다.

한편, 압축률, 투수계수, 해석메쉬가 지속적으로 업데이트되는 Case 3은 Case 1이나 Case 2와는 다른 압력 발달 양상을 보여준다. Case 1과 같이 즉각적으로 압력이 전달되지는 않지만, 투수계수의 증가로 인해 Case 2에 비해 현저히 빠른 유동을 보여준다. 주목할만한 점은 Case 1과는 달리 주입공 주변의 압력이 소산되지 않고 높은 수준을 유지하고 있다는 점이다. 이는 주입공 주변에서 단층의 탄성변형으로 인한 수리간극의 부피 증가와 투수계수 증가가 동시에 일어나기 때문이다.

각 해석 시나리오에서 주입공 주변의 유효수직응력(Fig. 18)을 살펴보면, Case 1에서는 압력이 0.5 MPa(초기 압력)에서 1.0 MPa(주입압)로 증가하였을 때 정확히 0.5 MPa(5.679 MPa → 5.179 MPa)의 유효수직응력 감소가 발생한다. 그러나 Case 2와 Case 3에서는 주입공 주변의 유효응력 감소량이 압력 증가량(초기압력 0.5 MPa → 주입압 1 MPa)보다 낮은 수치를 보였다. 앞서 언급한 바와 같이 이러한 현상은 단층과 주변 암반의 탄성거동과 단층 수직변위의 구속에 따른 전응력의 증가에 기인한다. 이론적, 해석적 접근법에서는 주입 이외의 외력이 작용하지 않을 때, Biot 상수를 1.0으로 가정하여 전응력이 일정하다고 가정하지만, 단층의 변형과 주변암반의 구속에 의하여 단층 주변의 전응력은 상승하게 되며, 이에 따라 유효응력의 증가분이 유도되는 것을 모든 연동해석 케이스에서 관찰할 수 있었다.

상기한 해석결과는 단층을 대상으로 하는 수치해석 연구에서 수리간극모델 및 수리역학 커플링관계식의 적절한 선택이 얼마나 중요한지를 시사한다. 단층의 경우, 압력 수준에 따라 초기수리간극의 크기가 수십 ~ 수백 배 이상 증가할 수 있으며, 수리유동경로의 압축성과 투수성 변화에 따라 압력 발달의 속도와 분포는 전혀 다른 형태로 발달할 수 있으므로, 단층 모델링 시 이에 대한 이해와 주의가 반드시 수반되어야 한다.

4. 스위스 Mont Terri 지하연구시설 단층재활성 시험 모델링

Task B의 2단계, 3단계 연구는 각각 스위스 Mont Terri 지하연구시설의 단층 손상대(fault damage zone)와 단층핵(fault core)에서 수행된 물 주입시험 결과를 모델링하는 심화 단계의 연구로서 각각 1단계에서 개발한 해석코드를 이용해 현장에서 얻어진 단층 재활성 시험 결과를 모델링하는 연구이다. 주어진 실험자료를 수치적으로 재현해 냄으로써 단층 내 수리역학적인 메커니즘을 이해하고 예측할 수 있는 수치모델을 개발하는 것이 주요 연구목표이다. Fig. 19는 단층 재활성 시험이 실시된 시추공 위치와 시험장치를 보여준다. 단층핵을 교차하는 BFS1 시추공(심도 47.2 m)과 BFS2 시추공(심도 37.2, 40.6, 44.65 m)을 통해 총 4회의 물 주입시험이 수행되었으며, ‘Step-rate Injection Method for Fracture In-situ Properties(SIMFIP)’ 시험이 시행되었다(Guglielmi et al., 2017). SIMFIP은 국제암반역학회(International Society for Rock Mechanics, ISRM)의 표준시험법 중 하나로 단층이나 절리와 같은 암반 균열에 다양한 압력의 물을 주입하여 탄성, 비탄성 변형을 유도함으로써 단층의 역학적․수리적 물성을 측정하는 현장시험법이다(Guglielmi et al., 2014).

Fig. 19.

Mont Terri ‘Main Fault’ reactivation experiment (after Guglielmi et al., 2014, Guglielmin et al., 2017): (a) Mont Terri Main Fault plane with the URL; (b) cross section of the Main Fault; (c) SIMFIP test equipment; (d) three-dimensional deformation unit

본 절에서는 Task B 2단계의 연구결과(Park et al., 2020)를 간략히 소개하고, 그 의의에 대하여 논의하고자 한다. 2단계 연구의 해석대상은 BFS2 시추공의 37.2 m 지점에서 이루어진 주입실험이며, 수평방향으로 약 1.5 m 이격된 BFS4 시추공에서 압력 측정이 이루어졌다. Fig. 20은 현장시험 결과(Guglielmi et al., 2017)를 보여주는 것으로 각 그림에서는 시간의 경과에 따른 주입압(injection chamber pressure)과 주입율(injection flow rate), 주입지점에서 일정 거리에 있는 모니터링 지점의 압력(monitoring pressure) 변화 곡선, 단층의 상하반에 각각 설치된 앵커들 간 상대변위 변화 곡선이 제시되어 있다. 여기서 앵커는 주입공에서 수직방향으로 0.5 m 거리에 위치하며 상반과 하반에 각각 설치되었다.

Fig. 20.

Field experimental results selected for Task B step 2 modelling: Injection and monitoring pressures, injection flow rate and relative displacement of upper anchor to lower anchor (Guglielmi et al., 2017)

Fig. 20의 현장시험 결과를 벤치마크 모델의 해석결과(Fig. 14, Fig. 15)와 비교하여 살펴보면, Mont Terri 지하연구시설 단층의 거동이 FM2 모델보다 FM1 모델의 거동과 유사한 특징을 보임을 확인할 수 있다. FM1 모델의 거동과 같이, 현장시험에서도 낮은 주입압 조건에서 주입율과 모니터링 지점 압력의 변화가 거의 발생하지 않다가 주입압이 증가하면서 급격한 변화가 관찰된다. 본 연구에서는 FM1 수리간극모델이 현장시험 결과(Fig. 20)를 모델링하기에 더욱 적합한 것으로 판단하였으며, 다음의 같은 세 가지 현상을 재현하는 데에 주안점을 두고 연구를 진행하였다: 1) 420 ~ 453초 사이의 주입율과 주입부피, 2) 420초에 발생하는 급격한 모니터링 압력의 증가, 3) 상하반 앵커의 변위의 방향과 규모.

2단계 연구에서는 현장시험의 기초적인 정보와 실험결과만이 주어지고, 나머지 입력물성과 해석조건은 각 참여팀이 독자적으로 결정하는 방식으로 연구가 진행되었다. 본 연구에서는 문헌연구 등을 통해 획득한 Mont Terri 현장의 단층과 암반 물성(단층 방향, 단층마찰각, 탄성계수 등), 현지응력 조건 등을 바탕으로 여러 해석 케이스들을 선정하고 시행착오법(trial and error)을 통해 파라미터를 개선해가면서 현장시험 결과를 재현하고자 하였다. Table 4는 해석모델에 적용한 입력 변수들의 범위와 함께 현장자료를 가장 잘 반영한 해석케이스의 파라미터 조합을 제시한 것이다. 현지응력의 방향과 크기 범위는 Yong et al.(2010), Martin and Lanyon(2003)의 연구를 참고하였다.

Fig. 21, Fig. 22, Fig. 23은 해석케이스 중 현장시험 결과를 가장 적절히 재현한 경우(Table 4의 Best matching case)의 해석 결과로, 각각 주입압(주입 챔버 내 압력)과 모니터링 지점 압력, 주입율, 상하반 앵커의 상대변위를 현장 취득 자료와 비교하여 나타낸 그림이다. 단, Fig. 22의 주입율 곡선에서 주입 초반 50초 사이에 현장시험 자료에서 관찰되는 주입율의 증가는 본 연구의 모델링 범위를 벗어난다.

현장시험 자료와 모델링 결과를 개략적으로 비교하면, 균열의 개시 압력, 주입율의 증가 양상, 주입부피 등에서 상당히 일치하는 결과를 얻을 수 있었다. 420초를 전후하여 닫힌 균열이 개방되면서 주입율과 모니터링 압력이 급격히 증가하는 것을 확인할 수 있으며, 420 - 453초 사이에 주입율도 거의 일치하였다. 현장자료에서 이 구간에 단층으로 주입된 물의 부피는 약 7.7×10^-3 m³로서, 수치해석모델에서 계산된 수리간극의 부피증가량 8.0×10^-3 m³과도 근접한 결과를 나타내었다. 이는 단층의 균열이 주입수에 의해 완전히 포화되었다고 가정할 때, 개발된 수치모델이 단층의 수리간극부피 변화를 적절히 모델링하고 있음을 의미한다. 수치모델에서는 음의 주입율, 즉 역류 현상(flow back)이 관찰되었으며, 현장에서도 측정이 이루어지진 않았으나 시험 종료 시 일정량의 역류가 관찰된 것으로 보고되었다.

앵커의 상대변위 해석결과는 현장자료와 비교할 때 수직방향으로 3 - 4배 가량 다소 크게 예측되지만, 방향에 따른 변화 양상과 규모는 대체로 잘 일치하는 것으로 나타났다. 이러한 결과는 Park et al.(2019)이 발표한 연구결과(2단계 해석결과)에 비해 고무적인 결과로, 그들의 연구에서는 420초 이후 단층이 재활성과 함께 과도한 전단변위가 발생하여 모든 방향에서 수백 μm에 달하는 상대변위가 예측되었다. 앵커의 변위벡터를 단층의 전단방향과 법선방향으로 분해하면, 변위 발생의 양상을 좀 더 효과적으로 살펴볼 수 있다(Fig. 24). 앵커의 변위는 단층의 변위와 거의 일치하는 것으로 나타났고, 앵커의 변위는 단층의 전단방향보다 법선방향으로 더 우세하게 발생함을 확인할 수 있다. 420초 이전에 앵커의 변위는 대부분 단층의 법선방향 탄성변형에 기인하였고, 전단파괴 이후에도 단층의 전단방향으로는 10 μm 이내의 작은 변위를 보였다.

Fig. 21.

Variations of monitoring pressures obtained from field experiment and numerical model (best matching case)

또한, 현장의 단층면에서 수리간극이 변화하는 주요메커니즘은 전단파괴에 의한 단층 미끄러짐(hydro-shearing)보다는 단층의 인장방향으로 발생하는 균열이 개방(tensile opening)이 더 밀접한 영향을 미치는 것으로 판단할 수 있다. Yong et al.(2010), Martin and Lanyon(2003)의 연구에서는 현지응력 범위를 σ_V = 6.0–7.0 MPa, σ_H = 4.0–5.0, σ_h = 0.6–3.0로 제시하였으나, 이 범위의 현지응력을 적용하는 경우 모든 해석케이스에서 과도한 전단변위가 발생하였다. Park et al. (2019)의 연구에서는 단층의 방향과 현지응력의 크기를 조합하여, 이론적으로 최소의 초기 전단응력이 작용하도록 파라미터를 가정한 경우에도 과도한 단층 미끄러짐으로 인하여 앵커 변위를 재현하지 못하였다. Table 4의 Best matching case의 경우, 초기 전단응력이 최소화되도록 수직주응력을 σ_V = 5.1 MPa로 하향 적용하였을 때의 해석결과이며, 이때 단층에 작용하는 초기 수직응력과 초기 전단응력은 각각 5.01 MPa, 0.003 MPa로서 초기 전단응력이 굉장히 작은 수준일 때 앵커의 변위를 재현할 수 있었다. 현장시험을 수행한 Guglielimi et al.(2017)은 현지의 응력상태가 시추공 굴착과 응력 재분배에 의해 교란되었을 가능성을 보고한 바 있으며, 현장시험의 해석결과로 볼 때 현지응력 상태가 문헌상 제시된 수치와는 다소 차이가 있을 것으로 추정된다.

Fig. 22.

Variations of injection pressure obtained from field experiment and numerical model (Park et al., 2020)

Fig. 23.

Variations of anchors’ relative displacement obtained from field experiment and numerical model (Park et al., 2020)

Fig. 24.

Comparison between anchor displacement and fault displacement (Park et al., 2020)

5. 결론 및 요약

본 논문에서는 국제공동연구 DECOVALEX-2019 프로젝트 Task B의 연구내용을 소개하고, 국내 참여팀인 한국지질자원연구원의 주요 연구성과를 요약, 제시하였다. 본 연구에서는 TOUGH-FLAC 연동해석 기법을 사용하여 유체 주입(유입)으로 인한 단층재활성 거동을 평가할 수 있는 수치해석기법을 개발하였다. 단층의 수리역학 연계거동 해석 시 일정 두께의 연속체 모델을 사용하는 기존의 접근법과는 달리, 역학적 Interface 요소를 사용하여 단층의 거동을 보다 직접적으로 재현할 수 있는 해석기법을 제안하였으며, 벤치마크 해석과 스위스 Mont Terri 지하연구시설의 단층재활성시험 모델링을 통하여 현장적용성 및 타당성을 확인하였다. 본 연구에서 제시한 TOUGH-FLAC Interface 모델은 현장조건에 최적화된 다양한 수리간극 모델(또는 수리역학 커플링 모델)에 유연하게 적용할 수 있을 것으로 판단되며, 단층과 현지응력 조건에 따른 수리역학적 파괴메커니즘과 거동특성을 예측하는 데에 유용하게 사용할 수 있을 것으로 기대된다.

하지만 개발된 해석기법을 보다 다양한 문제에 적용하기 위해서는 인터페이스 요소와 접촉(contact)의 처리에 있어 세심한 주의와 추가 연구가 필요하다. 현재 DECOVALEX-2019의 공식적인 연구기간은 종료된 상태이나, Step 3에 대한 모델링 연구가 진행중에 있으며, 상호 교차하는 인터페이스(단층) 간 접촉 문제를 처리하기 위한 후속 연구를 수행하고 있다(Fig. 25). 향후 DECOVALEX-2019를 통해 구축한 국외 연구팀들과의 긴밀한 협력관계를 바탕으로, 개발된 수치모델을 지속적으로 개선해 나갈 예정이다.

Fig. 25.

Mesh discretization for Step 3 modeling