1. 서 론
2. 해석개요
2.1 DECOVALEX-2019 Task B Benchmark Model Test (BMT)
2.2 해석모델
2.3 TOUGH-FLAC 연동해석
3. 해석 결과 및 토의
3.1 BMT1-BMT3
3.2 BMT4-BMT7
4. 요약 및 결론
1. 서 론
최근 이산화탄소나 방사성폐기물의 지중처분, 심부지열 개발, 셰일가스 개발 등 심부 지층의 활용기술에 대한 수요가 전 세계적으로 증가하고 있다. 심부 지층은 천부 지층과 비교할 때, 고온, 고응력, 고압의 극한 환경과 제한된 접근성으로 인해 특성 평가와 거동 예측이 매우 어려울 뿐만 아니라, 열-수리-역학-화학적(thermal-hydrological-mechanical-chemical, THMC)으로 연계된 복합적 거동을 대상으로 한다는 점에서 신뢰성 높은 수치모델을 개발하는 일은 매우 중요하다. 상기한 기술 분야에서 환경위해물질의 차폐와 기밀성, 대상 구조물과 주변 지층의 안정성 확보 등을 위해 열-수리-역학-화학적으로 연계된 다양한 복합거동에 대한 고려가 필요하지만, 단층이나 절리와 같은 불연속면의 거동을 고려하는 암반공학 문제에서는 수리-역학적 연계거동 해석 기술의 개발이 우선적으로 해결되어야 할 문제 중 하나이다.
일반적으로 암반을 구성하는 암석의 공극률과 수리전도도는 매우 낮기 때문에 암반 내 유체의 유동은 주로 단층이나 절리와 같은 불연속면을 통해 이루어진다. 따라서 개별적인 불연속면의 수리역학적 구성모델의 개발과 수치모델에의 적용은 가장 기초적이면서 핵심적인 연구라고 할 수 있다. 암반 불연속면의 수치모델링 연구는 주로 이산화탄소 지층처분이나 지열 개발을 고려하는 지층에서 인공적인 유체 주입으로 인한 인접 단층의 재활성화와 관련하여 이루어져 왔다(Gudmundsson, 2004; Vidal-Gilbert et al., 2009; Cuisiat et al., 2010; Cappa and Rutqvist, 2011; Orlic et al., 2011; Morris et al., 2011; Rinaldi et al., 2014; Rutqvist et al., 2015).
수치해석에서 단층을 모사하는 방법은 이를 독립적인 개별 불연속면으로 고려할 것인지, 연속체의 집합으로 고려할 것인지에 따라 크게 두 가지 접근법으로 나눌 수 있다(Park et al., 2018). 전자는 연속체 또는 불연속체 모델 내에 개별 불연속면을 정의함으로써 단층의 개폐(opening and closing)나 미끄러짐(slip)을 직접적으로 재현하는 방법이며, 후자는 단층을 일정 두께를 갖는 연속체의 일부로 가정하는 방법이다. 유한요소법이나 유한차분법과 같은 연속체 모델은 광역 규모에서 연속적인 응력과 변형률을 계산할 수 있다는 장점이 있지만, 소성변형률을 기반으로 파괴면을 결정하기 때문에 단층의 미끄러짐(slip)이나 개폐(opening/closing)를 직접적으로 재현할 수 없고, 해석 메쉬의 형상을 유연하게 적용할 수 없다는 제약이 따른다. 또한, 일반적인 암반공학 문제의 스케일에서 단층 수리간극의 크기는 규모에 비해 매우 작기 때문에 해석 요소의 불량한 종횡비를 개선하기 위해 실제보다 더 큰 간극 두께와 등가 물성을 가정하여야 한다는 단점이 있다. 이에 대한 대안으로서 일련의 절점 위에 인터페이스 요소(interface element)를 정의하여 수직/전단방향 거동을 정의하는 방법이 널리 사용되고 있으나(Vidal-Gilbert et al., 2009; Cappa and Rutqvist, 2011; Orlic et al., 2011; Cuisiat et al., 2010), 연속체 요소와 인터페이스 요소간의 강성차로 인해 해석의 불안정성을 야기할 수 있다는 한계점이 있다. 한편, 개별요소법을 사용하는 불연속체 모델은 보다 직접적인 방식으로 단층의 거동을 재현할 수 있다는 장점이 있다(Gutierrez and Makurate, 1997). Morris 외(2011)는 개별요소 모델을 통해 In Salah 프로젝트의 이산화탄소 주입으로 인한 지층과 인접 단층의 수리역학적 거동을 모사한 바 있다. 한편, Cappa와 Rutqvist(2011)는 TOUGH-FLAC 시뮬레이터를 이용한 지층 내 이산화탄소 주입 모델링 연구에서 FLAC3D를 통해 단층의 역학적 거동을 모사하기 위한 세 가지 방법을 제시하였다: 1) 두께가 없는 불연속 인터페이스 모델(zero-thickness mechanical interface), 2) 유한한 두께의 등가연속체 모델(equivalent continuum model using solid element), 3) 연약면을 포함한 유한한 두께의 등가연속체 모델(combination of solid elements and ubiquitous-joints oriented as weak planes). 그들은 Cappa와 Rutqvist는 주입공과 인접한 단층에서 압력, 유효 수직응력, 전단응력, 전단변위를 비교한 결과, 세 모델이 유사한 해석 결과를 보이며, 단층 모델링 시 작업이 용이한 등가연속체를 사용하여도 무방하다고 보고하였다. 상기한 단층 모델링 기법은 서로 다른 특징과 장단점을 내포하고 있으며, 이에 대한 선택은 관심 영역의 스케일과 가정된 단층 구성모델의 물성 획득 여부, 암반의 상태나 불연속면의 발달 정도 등에 따라 좌우된다(Leijon, 1993; Bohloli et al., 2015; Park et al., 2018).
본 연구는 유체 주입에 따른 단층 활성화와 수리역학적 거동을 모사하기 위한 수치해석 연구로 DECOVALEX(DEvelopment of COupled models and their VALidation against EXperiments)-2019 프로젝트 Task B의 일환으로 수행되었다. DECOVALEX 프로젝트는 방사성 폐기물의 지층처분과 관련하여 공학적 방벽과 지하 암반의 T-H-M-C 해석기법 개발을 위한 국제공동연구로서 한국지질자원연구원은 Task B에 참여하여 연구를 수행하고 있다(Park et al., 2018). Task B의 정식명칭은 ‘Fault Slip Test: Modelling the induced slip of a fault in argillaceous rock’으로 단층 내 유체의 주입으로 인한 수리역학적 거동을 예측할 수 있는 해석모델을 개발하는 데에 최종목표가 있다. 본 연구는 Park et al.(2018)이 보고한 Task B 1단계 연구의 후속연구로서, 그들의 논문에 DECOVALEX-2019 프로젝트의 Task 및 현황에 대해 자세히 소개된 바 있다. 여기에서는 Task B에 새롭게 추가된 Benchmark Model Test(BMT) 시나리오 해석의 결과와 이를 바탕으로 업데이트된 1단계 해석의 결과를 수록하였다. 단층의 수리역학적 거동해석을 위하여 TOUGH-FLAC 연동해석 기법을 이용하였으며, 7가지 BMT 시나리오 해석을 통해 단층의 저류특성(압축률)과 투수특성(투수계수)의 변화, 단층 미끄러짐이 해석 결과에 미치는 영향을 살펴보았다. TOUGH2 코드를 통한 수리해석에서는 여러 상태방정식(equation of state, EOS) 중 물과 공기의 유동을 대상으로 하는 EOS4 모듈을 사용하였고, FLAC3D 해석에서는 불연속 모델인 인터페이스 요소(interface element)를 사용하여 단층의 역학적 거동을 모사하였다. 또한 단층의 역학적 변형(간극의 변화)으로 인한 수리물성 변화(압축률, 투수계수)와 기하학적 변화(해석 메쉬 변형)를 수리해석에 반영할 수 있는 커플링 모듈을 개발해 TOUGH-FLAC 연동해석에 적용하였다.
2. 해석개요
2.1 DECOVALEX-2019 Task B Benchmark Model Test (BMT)
DECOVALEX-2019 프로젝트의 Task B(Fault Slip Test: Modelling the induced slip of a fault in argillaceous rock)는 Swiss Federal Nuclear Safety Inspectorate(ENSI)이 제안하였으며, 총 6개국(스위스, 미국, 독일, 한국, 캐나다, 대만)에서 7개의 연구그룹이 참여하고 있다: Swiss Federal Nuclear Safety Inspectorate(ENSI, 스위스), Lawrence Berkeley National Laboratory (LBNL, 미국), Federal Institute for Geosciences and Natural Resources/Helmholtz Centre for Environmental Research (BGR/UFZ, 독일), Helmholtz Centre Potsdam, German Research Centre for Geosciences(GFZ, 독일), 한국지질자원연구원(KIGAM, 한국), Canadian Nuclear Safety Commission(CNSC, 캐나다), Institute of Nuclear Energy Research(INER, 대만). Task B에 참가하는 연구그룹들은 1년에 2회 개최되는 워크숍과 서면회의를 통해 서로의 연구결과를 공유하고 의견을 교류함으로써 해석코드를 지속적으로 개선하고 있다. Task B의 최종목표는 유체의 주입(유입)으로 인한 단층 활성화와 수리역학적 거동을 모사할 수 있는 수치모델을 개발하는 데에 있으며, 향후 개발된 수치모델을 스위스 Mont Terri 지하연구시설(Underground Research Laboratory, URL)에서 수행된 ‘Fault activation experiment’ 현장시험에 적용하여 실측 결과와의 비교를 통해 그 타당성을 검증할 예정이다. Table 1은 Task B 참가팀들의 수치해석 코드 및 단층 모델링 방법을 정리한 것이다. 한국지질자원연구원은 단층을 불연속면과 등가 연속체로 가정하는 두 가지 접근방식을 통해 수치해석기법을 개발하고 있으며, 본 논문에서는 불연속 인터페이스를 이용한 경우의 해석 결과만을 수록하였다.
Table 1. Numerical models of Task B participating teams
Task B의 추진일정은 2019년까지 크게 세 단계로 구성된다(Park et al., 2018). 1단계 연구(Step 1)는 각 연구그룹이 독자적인 해석 코드를 개발하고 모델링 기법을 수립하는 일종의 준비 단계의 연구로, 1단계 연구에서는 동일한 단층 개념모델에 대한 각 참여그룹의 해석 결과를 상호 비교함으로써 개발된 해석코드를 지속적으로 개선하게 된다. 2단계 연구(Step 2)와 3단계 연구(Step 3)는 각각 1단계에서 개발한 해석코드를 이용해 Mont Terri URL에서 수행된 ‘minor fault’와 ‘major fault’에 대한 주입시험을 모델링한다. 이 단계에서는 현장 조건과 시험 결과가 배포되어 수치해석 결과와의 비교를 통해 각 해석코드/해석기법의 타당성을 검증하게 된다. 당초의 계획에 따르면 2018년 12월 현재 3단계 연구가 진행 중에 있어야 하나, 2017년 4월 영국에서 개최된 제4차 워크숍에서 각 연구팀의 해석결과가 상당히 큰 차이를 보였으며, 근본적인 원인 규명과 해결을 위하여 7가지 BMT 시나리오에 대한 추가 해석의 필요성이 제기되었다(Rutqvist et al., 2017). 따라서 Task B의 연구그룹들은 추가 해석을 통해 해석코드의 완성도를 높이고 서로간의 교차검증을 수행하기로 합의하였다. 본 논문에서는 2018년 10월 제6차 워크숍에서 발표된 BMT 해석 결과에 대해 소개하였다.
2.2 해석모델
상기한 바와 같이, BMT 시나리오 해석은 1단계 연구에서 보인 각 연구그룹의 해석결과의 차이에 대한 근본적인 원인을 규명하고 해석코드의 완성도를 높이기 위해 추가적으로 제안되었다. 따라서 해석모델의 형상, 초기응력, 입력물성 등은 1단계 연구(Park et al., 2018)와 동일하다. 단, 1단계 연구에서는 두 가지 수리간극 모델(FM1, FM2)이 사용되었으나, BMT 해석은 FM2 모델만을 대상으로 한다.
Fig. 1은 BMT 모델의 형상과 좌표계, 모니터링 지점 등의 정보를 보여준다. 모델은 20 × 20 × 20 m3의 정육면체로 중심에 65° 경사의 단층이 관통하며, 단층의 주향은 x축과 동일한 방향이다.
초기 응력은 각 좌표축에 평행한 방향으로 주응력 σx= 3.3 MPa, σy= 6.0 MPa, σz= 7.0 MPa이 작용하며, 초기 지하수압은 0.5 MPa이다. 해석모델의 경계는 각 면의 수직방향으로 변위가 구속되는 이동단(roller) 조건이며, 경계에서의 응력과 압력은 해석 중 일정하게 유지된다. 단층의 주향방향으로 작용하는 주응력 σx는 단층에 수직응력과 전단응력을 유발하지 않으므로 단층에 작용하는 수직응력과 전단응력은 이론적으로 식 (1)과 같이 계산할 수 있다.
| $$\begin{array}{l}\sigma_n=\frac12(\sigma_z+\sigma_y)+\frac12(\sigma_z-\sigma_y)\cos(2\times65^\circ)=6.179\;MPa\\\tau=\frac12(\sigma_z-\sigma_y)\sin(2\times65^\circ)=0.383\;MPa\end{array}$$ | (1) |
단층의 전단강도는 Coulomb 파괴식을 따르며, 이론적으로 단층의 미끄러짐을 유발할 수 있는 이론적인 최소압력(critical injection pressure, Pcrit)은 식 (1)의 수직응력과 전단응력을 이용해 다음 식 (2)와 같이 계산할 수 있다.
| $$P_{crit}=\sigma_n-\frac{(\tau-c)}{\tan\varnothing}=5.231\;MPAa$$ | (2) |
여기서 c는 점착력으로 0이며, ∅는 마찰각으로 22°이다.
BMT 시나리오 해석에서는 암반이 불투수성이며, 주입수의 유동은 단층을 통해서만 이루어지는 것으로 가정한다. 본 연구에서는 물 주입으로 단층과 주변 암반의 역학적 거동과 수리간극 크기 변화를 반영하여 주입수의 유동 특성(유량 및 압력 분포)을 적절히 예측할 수 있는 수치코드를 개발하는 것이 주요 연구내용이라고 할 수 있다. 수리 간극은 다음의 식 (3)과 같은 요소를 통해 구성되며, 수리간극과 역학적 간극은 동일한 것으로 간주된다.
| $$b_h=b_{hi}+\bigtriangleup b_{he}+\bigtriangleup b_{hs}+b_{hc}$$ | (3) |
여기서 bhi는 초기 수리간극, △bhe는 단층의 수직방향 탄성변형에 의한 수직변위, △bhs와 △bhc는 단층의 전단파괴(미끄러짐 또는 활성화) 이후 발생하는 수직변위이다.
△bhe는 단층의 수직강성 kn과 단층에 작용하는 유효 수직응력 변화량 △σnʹ을 이용하여 식 (4)와 같이 계산할 수 있다.
| $$\bigtriangleup b_{he}=\bigtriangleup\sigma_n'/k_n$$ | (4) |
전통적인 암석역학에서는 전단파괴 이후 발생하는 간극의 변화를 단층면의 거칠기에 의한 수직팽창으로 간주하여 왔다. 식 (3)의 △bhs가 이에 해당되며, 전단변위 △us와 팽창각 ψ를 이용하면 다음과 같이 수식화할 수 있다.
| $$\bigtriangleup b_{hs}=\bigtriangleup u_s\times\tan\psi$$ | (5) |
Guglielmi 외(2015)는 프랑스 Tournemire URL에서 수행된 단층 주입시험 및 이에 대한 수치해석 결과를 토대로, 전통적인 팽창각 모델이 현장시험에서 파괴 이후 발생한 간극의 변화량을 과소평가한다고 보고하였다. 그들은 이에 대한 대안으로서 전단 또는 인장파괴 발생 시 일정량의 수직팽창이 발생하는 모델을 제안하였으며, 식 (3)의 △bhc가 이에 해당된다. Task B에서는 △bhc를 ‘creation aperture at rupture’로 지칭하고 △bhs와 구별하여 사용하고 있다.
Task B의 1단계에서는 두 가지 수리간극모델(FM1, FM2)에 대한 수치해석이 제안되었으며, BMT 해석은 FM2 모델만을 대상으로 한다. FM1 수리간극모델에 대해서는 Park et al.(2018)의 논문에 상세하게 기술되어 있다. BMT 해석이 수행된 FM2 모델에서는 단층이 10 μm의 초기 수리간극을 갖는 것으로 가정한다. 수리간극은 탄성변형에 의해 변화하며, 전단파괴 발생 이후에는 팽창각 10°에 의해 전단변위에 따른 수리간극 증가가 발생한다.
Fig. 2는 단층에 적용되는 주입압 조건을 보여주는 것으로 Mont Terri URL의 ‘minor fault’에 대한 현장주입시험 조건과 동일한 것이다. 주입과정은 총 9단계로 구성되는데, 각 주입단계에서는 일정한 주입압을 수십 초간 유지하게 된다. 첫 번째 단계의 주입압은 0.7446 MPa이며, 여덟 번째 단계에서 최대 주입압인 6.302 MPa이 적용된다. 마지막 단계에서는 압력을 낮추어 354초 동안 3.382 MPa의 일정 주입압을 적용한다.
추가된 BMT 해석의 7가지 시나리오는 다음과 같다.
1) BMT1: Uncoupled fluid flow during constant pressure injection
2) BMT2: One-way coupled fluid flow during constant pressure injection
3) BMT3: Two-way coupled fluid flow during constant pressure injection
4) BMT4: Uncoupled fluid flow during step pressure injection
5) BMT5: One-way coupled fluid flow during step pressure injection
6) BMT6: Two-way coupled fluid flow during step pressure injection
7) BMT7: One-way coupled fluid flow during step pressure with activation
BMT1, BMT2, BMT3는 각 팀의 해석모델이 단층의 수리역학적 커플링을 적절히 재현하는지를 검토하기 위한 해석 케이스로, 모두 일정 압력(1 MPa)이 800초간 유지되는 조건이다. BMT1의 경우, 단층에 매우 큰 수직강성(1.0 × 107 GPa/m)을 적용하여 물 주입으로 인해 수리간극의 탄성변형이 거의 발생하지 않는 조건을 가정한다. 즉, 수리간극의 크기는 물 주입에 의해 거의 변화하지 않고 초기값(10 μm)를 유지하게 되며 결과적으로 수리역학적 커플링은 발생하지 않는다(uncoupled analysis). BMT2의 경우 상대적으로 낮은 수직강성(20 GPa/m)을 적용하여 물 주입으로 인한 간극 부피의 변화를 반영한다. 즉, 단층의 높은 압축률(compressibility) 또는 저류계수(storativity)로 인한 영향을 살펴보게 된다. 그러나 역학적 변형(수리간극 변화)으로 인한 투수계수(permeability) 또는 투수량계수(transmissivity)의 변화는 고려되지 않으며, 수리해석으로부터 얻어진 압력 분포를 바탕으로 유효응력 계산을 수행하는 일방향 연동해석(one-way coupling)이 수행된다. BMT3은 BMT2와 입력 물성이 동일하며, 역학적 변형으로 인한 수리간극의 변화를 동시에 고려하는 쌍방향 연동해석(two-way coupling)이다. 따라서 BMT1과 BMT2의 결과, BMT2와 BMT3의 결과를 비교함으로써 각각 단층의 저류계수(storativity)와 투수량계수(transmissivity)의 변화가 해석결과에 미치는 영향을 검토할 수 있다. BMT4, BMT5, BMT6는 각각 BMT1, BMT2, BMT3와 동일한 조건이지만, Fig. 2의 단계적 주입압이 적용된다. BMT1, BMT2, BMT3의 경우, 주입압이 낮아 전단파괴가 발생하지 않을 것으로 예상할 수 있으며(식 (2)), BMT4, BMT5, BMT6에서는 높은 전단강도를 임의로 적용하여 단층 미끄러짐을 허용하지 않는다. BMT7의 경우, BMT5와 동일한 조건에서 단층의 미끄러짐을 유도하고 이로 인한 영향을 함께 검토하게 된다. 참고로, 1단계 연구의 해석조건은 BMT7과 동일한 입력 물성 하에 단층의 미끄러짐을 허용하는 쌍방향 연동해석이다.
Table 2는 1단계 연구의 FM2 수리간극모델과 BMT 시나리오에 대한 수치해석에 사용된 단층과 암반의 역학적 물성 및 주입 유체의 수리적 물성을 나열한 것이다. FM2 수리간극모델에 사용된 단층과 암반의 물성은 스위스 Mont Terri URL에서 수행된 부지조사 결과를 통해 결정된 값으로 현지의 단층과 모암인 ‘minor fault’와 Opalinus Clay의 물성이며, BMT 해석에서는 일부 상이한 단층 물성이 사용된다.
Table 2. Material properties for Step 1 and BMT simulations
상기한 BMT 모델의 수리역학적 연동해석을 위해서는 단층의 수리특성을 결정하는 투수량계수와 저류계수에 대한 이론적 배경을 살펴볼 필요가 있다. 연속체 해석에서 단층 내 유체의 흐름이 Darcy의 법칙을 따르며, 유량과 수리간극의 관계가 평행한 두 평판사이의 유체 흐름 모델인 삼승법칙(Witherspoon, 1980)을 따르는 것으로 가정하면 단층의 수리간극을 통한 유체의 흐름은 대수층 내 유동이론을 적용하여 설명이 가능하다.
대수층을 통과하는 지하수의 2차원 평면 유동을 고려하는 경우 대수층의 유체통과능력은 주로 투수량계수를 통해 다음과 같이 표현된다.
| $$T=Kb=\frac{k\rho_fg}\mu b$$ | (6) |
여기서 T는 투수량계수, K는 수리전도도, b는 대수층 두께, k는 고유투수계수, ρf는 유체의 밀도, g는 중력가속도이며, µ는 점성계수(dynamic viscosity)이다.
단층의 경우, 2차원 평면 유동으로 가정할 수 있으며, 대수층의 두께를 수리간극 bh로 치환하면, Darcy의 법칙과 삼승법칙을 이용하여 투수량계수를 다음과 같이 산정할 수 있다.
| $$T=\frac{\rho_fg}{12\mu}b_h^3$$ | (7) |
따라서 단층의 투수계수는 수리간극 bh와 다음의 관계에 놓인다.
| $$k=\frac{b_h^2}{12}$$ | (8) |
한편, 대수층의 저류계수는 단위수두의 변화에 따라 단면적을 통해 유입되거나 유출되는 물의 부피를 나타내는 지수로, 대수층의 압축성과 유체의 압축성에 의해 수두가 변화함에 따라 물이 저장되거나 방출된다. 이를 수식으로 표현하면 다음과 같다.
| $$S=\rho_fgb(\alpha+n\beta)$$ | (9) |
여기서 α는 대수층의 압축률, β는 물의 압축률, n은 대수층의 공극률을 의미한다.
단층의 경우, 유효응력의 변화에 따른 부피 변화가 수직 방향으로만 발생하며 공극률이 1인 대수층으로 가정할 수 있다. 또한 대수층의 압력 변화량이 유효응력의 변화량과 동일하다면 단층의 압축률은 식 (10)과 같이 수직강성 kn과 수리간극 bh를 통해 표현할 수 있으며, 단층의 저류계수는 식 (11)과 같이 표현할 수 있다.
| $$\alpha=-\frac{dV/V}{d\sigma'}=-\frac{db_h/b_h}{d\sigma'}=\frac1{b_n}\frac1{k_n}$$ | (10) |
여기서 V는 대수층의 부피, σʹ는 유효응력, kn는 단층의 수직강성이다.
| $$S=\rho_f\;g\;b_n\;(\frac1{b_hk_n}+\beta)$$ | (11) |
위에서 살펴본 바와 같이 단층의 수리특성을 결정하는 중요 인자인 저류계수와 투수량계수는 모두 수리간극에 의존적인 변수들이며, TOUGH2 코드에서는 위 인자들의 영향을 압축률과 투수계수 변화를 통해 반영할 수 있다. 본 연구에서는 이를 위하여 TOUGH-FLAC 시뮬레이터에서 수리역학 커플링 관계식을 수치화하였으며, 이와 관련된 자세한 내용은 제2.3절에서 기술하였다.
2.3 TOUGH-FLAC 연동해석
TOUGH-FLAC 연동해석 기법은 Rutqvist 외(2002)가 개발한 해석기법으로 현재까지 방사성폐기물 처분, CO2 지중저장, 심부지열개발, 압축공기에너지저장, 열에너지저장 등 다양한 암반공학 문제에 적용되어 왔다(Rutqvist, 2012; Rutqvist and Tsang, 2012; Rutqvist et al., 2015; Kim et al., 2012; Park et al., 2016; Park et al., 2018). 이는 다상(multi-phase), 다성분(multi-component)의 유체유동/열유동 코드인 TOUGH2와 암반 및 지반의 역학적 거동 해석프로그램인 FLAC3D를 연동하는 방법으로서 구성방정식을 통한 커플링이 아닌 두 프로그램 간 해석 결과의 교환을 통한 연성해석이다. 두 프로그램간의 연동은 TOUGH2의 비선형 반복계산(열-수리 유동해석)과 FLAC3D의 준정적 해석(역학적 거동해석)을 순차적으로 반복하는 방식으로 이루어지며, 두 프로그램의 해석 결과는 각각 TOUGH2 코드에 추가된 별도의 subroutine과 FLAC3D 내의 FISH 함수를 통해 구축된 연계 모듈을 통해 상호 교환된다. 이 과정에서 FLAC3D는 TOUGH2에서 계산된 열-수리적 해석 결과(온도, 압력 등)를 입력받아 역학적 거동을 계산하게 된다.
TOUGH-FLAC 시뮬레이터를 통한 대부분의 선행연구들에서는 단층을 일정 두께를 갖는 연속체의 일부로 가정하여 사용하여 왔다(Rinaldi et al., 2014; Rutqvist et al., 2015). 일반적인 암반공학 문제에서 단층의 수리간극은 연장이나 폭에 비해 매우 작기 때문에, 연속체 모델을 이용하는 경우 해석 메쉬 작성 시 실제보다 더 큰 간극 두께를 가정하게 된다. 따라서 단층 물성(강성, 강도, 수리간극 등)을 요소 두께에 상응하는 등가물성(변형계수, 강도, 투수계수 등)으로 변환하여 입력 자료로 사용하게 되며, 결과 분석 시에도 변형률을 통해 단층의 변위와 수리간극을 역산해야 한다는 단점이 있다. 본 해석 모델의 경우 단층으로 직접적인 물 주입이 이루어지고, 수리간극이 영역별로 다르거나 전단파괴에 의한 갑작스러운 수리간극의 증가가 예상되었으므로, 연속체 모델을 사용하는 경우 수리간극을 등가 투수계수로 변환하는 과정에서 지나친 가정의 중첩이 발생할 수 있을 것으로 판단되었다. 따라서 본 연구에서는 FLAC3D에서 불연속면의 미끄러짐(slip)과 분리(separation)를 재현할 수 있는 인터페이스 요소(interface element)를 통해 보다 직접적인 방식으로 단층의 거동을 모사하고자 하였다.
일반적으로 TOUGH-FLAC 연동해석은 수리해석과 역학적 해석을 두 프로그램에서 별도로 수행하고 순차적으로 해석결과를 교환하는 방식이므로 동일한 해석 메쉬를 필요로 한다. 그러나 본 연구에서는 역학적 해석이 수행되는 FLAC3D에서는 암반과 단층을 모두 모델링하였고, 수리 해석이 수행되는 TOUGH2에서는 단층에 대한 해석 메쉬만을 작성하였다. 이는 단층의 응력과 변위 계산에는 암반과 단층의 거동이 복합적으로 영향을 미치지만, 유체의 유동은 단층을 통해서만 이루어지기 때문이다. 즉, 두 프로그램의 연동해석 시 데이터의 교환은 TOUGH2의 전체 요소와 FLAC3D의 인터페이스 요소 사이에서만 이루어진다.
Fig. 3은 FLAC3D에서 작성된 해석모델의 형상을 보여준다. 단층과 암반의 입력 변수는 Table 2에 나열한 바와 같다. 암반은 탄성 모델이며, 단층은 탄소성 모델인 Coulomb sliding 모델을 따른다. Fig. 4는 TOUGH2 해석에 사용된 FM2 수리간극모델의 초기 해석메쉬를 보여주는 것으로, 주입정 주변의 단층 단면을 확대하여 함께 나타내었다.
TOUGH2의 경우 적분유한차분법(integral finite difference method)을 이용하여 공간이산화(spatial discretization)를 하기 때문에 폭에 비해 두께가 매우 얇은 요소에 대해서도 유연하게 적용할 수 있다. 본 연구에서는 이러한 장점을 활용하여 단층의 두께를 수리간극 크기에 상응하도록 모델링하였다. TOUGH-FLAC 시뮬레이터를 이용한 대부분의 선행연구들에서는 단층을 등가연속체로 모델링하므로 간극 변화에 따른 수리물성 변화에 초점을 맞추고 있으나, 여기에서는 실제 두께의 해석 메쉬를 사용하였으므로, 수리간극의 변화가 발생하면 해석 메쉬를 지속적으로 업데이트할 필요가 있었다. 따라서 본 연구에서는 FLAC3D와 TOUGH2의 연동 시 해석 메쉬의 변형과 입력물성의 변화를 지속적으로 업데이트하기 위한 커플링 모듈을 개발하였다.
Fig. 5는 TOUGH2와 FLAC3D 프로그램 간의 연동해석 과정을 보여주는 것으로(쌍방향 연동해석의 경우), 두 프로그램의 호출(실행)과 데이터 교환은 TOUGH2의 Newton iteration level에서 매 time 스텝마다 순차적으로 반복된다. TOUGH2에서 단층을 구성하는 각 해석 요소의 압력을 계산하여 FLAC3D로 전송하면, FLAC3D에서는 인터페이스 절점의 압력을 업데이트하여 유효응력 기반 준정적 해석을 수행한다. 계산된 인터페이스 절점의 수직변위는 다시 TOUGH2로 전송되어 각 해석 요소의 메쉬 정보와 투수계수, 압축률 등을 업데이트하고, TOUGH2는 다음 time step의 수리해석을 수행한다. 투수계수와 압축률의 경우에는 수직변위(수리간극)을 바탕으로 식 (8)과 식 (10)에 의해 계산된다.
이때 TOUGH2는 요소의 부피 중심에서의 해석 결과를, FLAC3D는 인터페이스의 각 절점에서의 해석 결과를 출력하게 되므로 두 프로그램 간의 데이터 교환을 위한 데이터 변환이 필요하였으며, 본 연구에서는 이를 위해 역거리가중법(inverse-distance weighting)을 이용하여 각 요소 또는 절점에서의 평균값을 대푯값으로 사용하였다. 여기서 가중치는 인접한 요소의 중심 또는 인터페이스 절점의 좌표정보를 이용하여 거리의 제곱에 반비례하는 것으로, 근거리에 있는 데이터에 더 높은 가중치를 부여하게 되므로 불규칙한 형태의 해석요소를 사용하는 경우 산술평균에 비해 더 합리적인 결과를 도출할 수 있다.
Table 3은 각 시나리오 모델의 해석 조건과 수리역학 커플링 시 FLAC3D에서 TOUGH2로 전송되는 업데이트 정보를 보여준다. BMT1과 BMT4 해석에서는 수리역학 커플링을 고려하지 않기 때문에 FLAC3D의 해석 결과는 TOUGH2의 수리해석에 아무런 영향을 미치지 않는다. 실제로 예비해석을 통해 해석 메쉬, 투수계수, 압축률을 모두 업데이트한 경우에도 높은 수직강성으로 인해 수직변위가 거의 발생하지 않아 TOUGH2 해석과 TOUGH-FLAC 연동해석이 동일한 수리적 거동(압력, 유량)을 출력함을 확인하였다. 한편, 일방향연동해석이 수행되는 경우에는 해석 메쉬와 압축률이 지속적으로 업데이트되며, 쌍방향연동해석의 경우 투수계수도 함께 업데이트되는 것으로 설정하였다.
Table 3. Descriptions of BMT simulations and updated information in hydro-mechanical coupling (α: compressibility, k: permeability)
Park 외(2018)의 연구에서는 쌍방향 연동해석 시 해석 메쉬와 투수계수가 역학적 변형에 따라 변화하지만 압축률은 항상 일정한 것으로 가정하였다. 그러나 식 (10)에서 확인할 수 있듯, 압축률은 수리간극에 의존적인 변수이며, 압력(유효응력)의 변화에 따라 단위부피당 유체의 방출량과 저장량을 결정하는 매우 주요한 수리적 인자이다. 본 연구에서는 수리간극의 변화에 따른 압축률의 변화를 추가적으로 고려할 수 있도록 그들의 연구에서 개발된 수치모델을 개선하였으며, 그 영향에 대해서는 제3장의 해석결과에서 자세히 기술하였다.
3. 해석 결과 및 토의
3.1 BMT1-BMT3
Fig. 6 ~ Fig. 9는 BMT1, BMT2, BMT3의 수리해석 결과를 보여준다. Fig. 6은 주입 후 800초 경과 시의 압력분포를 보여주며, Fig. 7은 주입 후 200초와 800초 경과 시 주입지점을 통과하는 주향방향 프로파일 상에서의 압력을 도시한 것이다. Fig. 8과 Fig 9는 각각 시간의 흐름에 따른 물 주입률의 변화와 모니터링 지점 P3에서의 압력 변화를 보여준다.
위 해석결과를 종합적으로 살펴보면, BMT1의 경우 주입 후 수 초 이내에 정상류 흐름(steady state flow)에 도달하는 것을 확인할 수 있다. Fig. 7에서 확인할 수 있듯이, 200초와 800초 경과 시의 압력 분포가 동일하며 이는 200초가 경과하기 이전 이미 정상류 상태에 도달하였음을 의미한다. 주입 초반 주입지점과 주변 단층의 압력 수두차로 인해 일시적으로 높은 주입률을 보이지만 수 초 이내에 일정한 주입률 0.00263 liter/min에 수렴하였고(Fig. 8), 모니터링 지점 P3(경사방향으로 이격거리 1.5 m)에서 압력도 초기상태인 0.5 MPa에서 빠르게 증가한 뒤 약 7초 이후에는 일정한 값을 나타내는 것을 확인할 수 있다(Fig. 9). 해석모델의 경계부까지 전 단층면을 따라 압력 분포가 형성되고 수두 경사는 주입지점으로부터 거리에 따라 지수적으로 감소하는 전형적인 형태의 방사형 흐름을 보인다. 단, 해석영역의 경계부(주입지점에서 10 m 이격)에서도 압력 구배가 0으로 수렴하지 않는 것은 해석영역의 크기가 경계 효과(boundary effect)를 배제할 수 있을 만큼 충분히 크게 설정되지 않았음을 의미한다. 위 모델의 경우, 수리해석 시 일정한 압력(0.5 MPa)이 유지되는 경계 조건이 적용되었는데, 만약 더 큰 해석영역이 설정되었다면 10 m 이격 거리에서 0.5 MPa보다 높은 압력이 발달하였을 것임을 짐작할 수 있다. 유체 주입으로 인한 방사형 유동 해석에서 압력의 분포는 투수층의 투수특성이나 저류특성뿐만 아니라 주입공의 반경에 영향을 받으며(Peng et al., 2002), 단층의 경우 수리 간극의 변형 특성이 압력의 전달 범위에 지대한 영향을 미칠 수 있으므로 차기 연구에서는 해석 모델의 적정 사이즈 결정을 위한 예비해석이 필요할 것으로 판단된다.
단층의 수직강성이 상대적으로 낮은 BMT2의 경우에는 시뮬레이션 종료 시(800초)까지 정상류 상태에 도달하지 않고 물 주입률과 압력이 지속적으로 변화하는 상태를 나타내었다. 해석 종료 시에는 주입공으로부터 약 3 m 이내의 가까운 거리까지만 압력이 전달되었고, 좁은 영역에서 급한 압력 구배가 형성되어 상대적으로 높은 물 주입률을 보였다. BMT1과 BMT2의 해석 종료 시 P3 지점에서의 압력을 비교하면, BMT2에서 매우 느린 속도로 압력의 전달이 이루어지고 있음을 확인할 수 있다. 위와 같은 결과는 단층 압축률의 영향이 해석에 적절히 반영되고 있음을 의미한다. 수리간극의 크기가 일정한 BMT1과 비교할 때, BMT2에서는 낮은 수직강성(높은 압축률)으로 인해 주입공 주변 압력이 증가하면 수리간극의 부피가 팽창하고, 이에 따라 주입수가 원거리로 전달되기보다는 주입공 주변에 저장되는 현상으로 나타났다. 이 경우, 800초 이후로 물 주입을 지속하면 결국 정상류 상태에 도달하여 압력 분포는 BMT1과 동일한 형태를 보이게 될 것이며, 물 주입률은 주입공 주변 수리 간극 부피 증가로 인해 상대적으로 높은 값에 수렴하게 될 것이다.
수리 간극의 변화에 따른 해석 메쉬의 변형, 압축률의 변화, 투수계수의 변화를 모두 고려하는 BMT3의 경우에는 BMT1이나 BMT2와는 다른 압력 분포를 보인다. Fig. 6을 살펴보면, BMT1과 BMT2의 경우 주입공으로부터 거리가 멀어짐에 따라 압력수두 경사가 감소하는 경향을 보이는 반면, BMT3의 경우에는 비교적 일정한 압력수두 경사를 보였다. 투수계수가 일정한 BMT2와 비교하면 압력의 전달 속도가 느리다는 점에서는 유사하나, 동일한 거리에서 더 높은 압력 수준을 유지하고 있음을 확인할 수 있다. 해석 종료 시 주입공으로부터 약 3 m 거리 이내에서는 BMT1보다 오히려 높은 수준의 압력을 보였으며, 이는 P3 지점에서의 압력 변화(Fig. 9)를 통해서도 확인할 수 있다. 물 주입률의 경우에는 주입 초반 높은 수치에서 감소하는 경향을 보이다가 약 10초 경과 시 크게 증가하는 경향을 보였고 해석 종료 시 BMT1이나 BMT2에 비해 10배 이상의 물 주입률을 나타내었다. 이는 주입공 주변에서 단층의 탄성변형으로 인한 수리간극의 부피 증가와 투수계수 증가에 기인하는 것으로 판단할 수 있다.
Fig. 10과 11은 BMT1, BMT2, BMT3의 역학적 해석 결과를 보여준다. Fig. 10은 시간의 흐름에 따른 주입공 주변(P1)의 유효수직응력 변화를 나타낸 것이며, Fig. 11은 주입 후 800초 경과 시 주입지점을 통과하는 주향방향 프로파일 상에서의 수직 변위(탄성변형)를 나타낸 것이다.
식 (1)에서 보인 바와 같이, 초기 응력 조건에 따라 단층에 작용하는 수직응력을 이론적으로 계산하면 6.179 MPa이며, 초기 압력이 0.5 MPa이므로 주입 전 유효수직응력은 5.679 MPa이다. FLAC3D를 통해 구현된 단층의 초기응력 조건도 이론적인 계산과 거의 동일한 수치를 나타내었다.
본 해석의 경우, 물 주입 이외의 외적 변화요인이 작용하지 않으므로 유효 수직응력과 압력의 합인 전응력(total normal stress)이 이론적으로 일정하여야 한다. 즉, 압력의 증가량(△P)은 유효수직응력의 감소량(△σʹ)과 동일할 것으로 기대할 수 있다. 예를 들면 1 MPa의 주입압이 적용되는 경우 단층의 유효수직응력이 5.179 MPa로 감소할 것으로 예상할 수 있으며, BMT1 해석에서는 Fig. 10에서와 같이 주입과 동시에 위 수치와 동일한 유효 수직응력을 나타내었다. 그러나 BMT2와 BMT3에서는 주입공 주변의 유효응력 감소량이 압력 증가량(초기압력 0.5 MPa → 주입압 1 MPa)보다 작은 수치를 보였다(해석 종료 시 BMT2 △σʹ=0.35 MPa, BMT3 △σʹ=0.47 MPa). 이러한 현상은 단층과 주변 암반의 탄성거동과 단층 수직변위의 구속에 따른 전응력의 증가에 그 원인을 찾을 수 있다.
단층의 수직강성과 압력 변화를 고려할 때, 각 해석 케이스에서 주입공 주변에 예상되는 수직변위는 다음과 같이 계산할 수 있다.
| $$\bigtriangleup u_n=\frac{\bigtriangleup P}{k_n}=\frac{0.5\times10^6\;Pa}{10^{16}\;Pa/m}=0.5\times10^{-10}m\;\;\;\;\;\;\;\;\;\mathrm{BMT}1$$ | (12) |
| $$\bigtriangleup u_n=\frac{\bigtriangleup P}{k_n}=\frac{0.5\times10^6\;Pa}{20\times10^9\;Pa/m}=25\times10^{-6}m\;\;\;\;\;\;\;\;\;\mathrm{BMT}2,\;\;\mathrm{BMT}3$$ | (13) |
Fig. 11을 살펴보면, BMT1 해석모델의 주입공(X=0)에서의 수직변위는 식 (12)와 동일하며, BMT2와 BMT3의 경우 이론적인 수치보다 낮은 수직변위를 보인다. BMT1의 경우 높은 수직강성으로 인해 초기수리간극(10 μm)에 비해 수직변위가 거의 발생하지 않았고 주변 암반의 영향을 받지 않는다. 반면, BMT2와 BMT3는 초기수리간극의 약 2배에 달하는 수직변위가 발생하게 되며 주변 암반의 탄성거동으로 인한 반력이 발생하여 수직방향으로의 전응력이 증가하는 것으로 판단된다. 한편, BMT2의 경우 BMT3에 비해 낮은 수직변위가 발생함에도 불구하고 유효수직응력이 더 크게 나타나는 것을 확인할 수 있는데, 이는 BMT2에서 주입지점 주변 압력 구배가 상대적으로 크고 이에 따라 국부적인 변위 구속에 따른 전응력 증가가 발생하기 때문인 것으로 판단할 수 있었다.
3.2 BMT4-BMT7
Fig. 12 ~ Fig. 14는 BMT4 ~ BMT7 모델의 해석결과로서 수리적 거동을 보여준다. Fig. 12는 주입 후 200초와 800초 경과 시 주입지점을 통과하는 주향방향 2차원 프로파일 상에서의 압력 분포를 도시한 것이며, Fig. 13과 Fig. 14는 각각 시간의 흐름에 따른 물 주입률의 변화와 모니터링 지점 P3에서의 압력 변화를 보여준다.
앞서 기술한 바와 같이, BMT4, BMT5, BMT6 해석모델에서는 매우 높은 전단강도를 가정하여 단층 미끄러짐을 허용하지 않으며, 본 연구에서는 이를 위해 마찰각 45°, 점착력 50 MPa을 적용하였다. BMT7의 경우, 마찰각 22°, 점착력이 0인 입력물성을 사용하였으나 단층의 전단파괴는 발생하지 않았다. 따라서 일방향 연동해석을 고려하는 BMT5와 BMT7의 해석결과는 정확히 일치하였으며, Fig. 12 ~ 14에서 ‘BMT5 & BMT7’는 두 해석 케이스의 동일한 해석결과를 의미하는 것이다. 모든 해석 케이스에서 해석 종료 시까지 전단응력의 변화가 거의 발생하지 않았고, 전단변위 역시 무시할만한 수준으로 나타났다.
앞서 BMT1, BMT2, BMT3의 비교를 통해 단층의 수리역학적 거동 특성에 미치는 투수계수(또는 투수량계수)와 압축률(또는 저류계수)의 영향을 검토한 바 있다. BMT4, BMT5, BMT6에서도 위와 유사한 결과를 얻을 수 있었다. BMT4의 경우, 각 주입단계에서 주입 초반에 정상류에 도달하였으며, 전체 단층면에 걸쳐 압력이 빠른 속도로 전달되는 결과를 보였다. 압축률이 큰 BMT5는 매우 느린 유체 유동 속도를 보였으며, 주입공 주변의 일부 영역에서만 압력 구배가 형성되었다. 단층의 압축률이 크고, 수리간극의 증가에 따라 투수계수가 증가하는 BMT6의 경우, 유체의 유동 속도는 느리지만, 좁은 영역에서 높은 압력이 형성되었다. Fig. 12에서 확인할 수 있듯이, 주입공과의 거리에 따라 압력 구배가 급격히 감소하는 BMT4, BMT5와는 달리, BMT6에서는 오히려 위로 볼록한 압력 분포 형태를 보였으며, 주입지점과 1.5 m의 거리에 있는 P3 지점의 압력을 살펴보면, BMT4의 경우 즉각적인 압력전달이 이루어지고 있는 반면, BMT5는 해석이 종료될 때까지 압력 증가가 미미했다. BMT6에서는 P3에서 즉각적인 압력 변화가 관찰되지 않았으나, 주입 후 45초 경과 시부터 압력이 증가하기 시작하다가 60초 이후에는 BMT4의 해석결과보다 훨씬 높은 압력 수준에 도달하였다. BMT6의 물 주입률이 BMT4나 BMT5의 주입률에 비해 1000배 이상 높은 수준임을 확인할 수 있으며(Fig. 13), 이는 주입에 따른 압력 증가로 인해 간극의 부피와 투수성이 동시에 증가하기 때문에 주입공과 가까울수록 많은 양의 물이 저장되어 빠르게 유동함을 의미한다. 한편, 물 주입률 곡선에서 각 주입단계 초반에 보이는 정점(peak)들은 압력변화로 인한 수리간극의 급격한 증가에 따른 것으로서 해석상에서 timestep에 의존적인 것으로 판단된다.
Fig. 15와 Fig. 16은 역학적 거동해석 결과를 보여주는 것으로, Fig. 15는 주향방향 프로파일 상에서의 수직변위를, Fig. 16은 주입지점 P1에서의 응력상태를 보여주는 것이다.
제3.1절에서 관찰한 바와 유사하게, 단층의 수직강성이 커 변형이 거의 발생하지 않는 경우(BMT4), 주입으로 인해 즉각적인 유효수직응력의 감소가 발생하고 압력의 증가량은 유효수직응력의 감소량과 거의 동일하게 나타났다. 단층의 수직강성이 작은 경우(BMT5)에는 주입공에서 수직변위가 크게 발생하고 높은 압력 구배로 인하여 국부적인 변위 구속에 따른 전응력의 증가가 크게 발생한다. 하지만, 수리 간극의 증가에 따른 투수계수의 변화를 함께 고려하는 경우(BMT6) 주입공 주변에서의 압력과 수직변위가 증가하면서 주입공에서의 응력 집중이 완화되는 경향을 보였다. Fig. 16을 살펴보면, BMT5에 비하여 BMT6의 수직응력이 현저히 낮은 값을 보임을 확인할 수 있다.
Fig. 16에는 각 해석모델의 유효수직응력 결과를 바탕으로 마찰각이 22°일 때 이론적으로 계산된 Coulomb sliding 모델의 전단강도 곡선이 함께 도시되어 있다(σʹtan22°). BMT4, BMT5, BMT6의 경우 임의로 높은 전단강도를 할당하여 전단파괴가 발생하지 않았고, BMT7의 경우에는 전응력의 증가로 인해 유효수직응력과 전단저항이 증가하여 파괴상태에 도달하지 않았다. 응력조건을 고려할 때 BMT4와 BMT6에서 만약 22°의 마찰각이 설정되었다면 전단파괴에 도달하였을 것을 예상할 수 있다.
Fig. 17과 Fig. 18은 BMT6과 같은 조건에 22°의 마찰각을 설정한 경우의 해석 결과를 보여준다. 전자는 단층의 전단변위 발생범위(주입 후 453초 경과 시)를 보여주며, 후자는 주변 암반의 변위(주입 후 267초, 420초, 453초, 807초 경과 시)를 보여주는 것으로 주입공 주변에서 전단파괴에 의한 단층 미끄러짐이 뚜렷이 발생하는 것을 확인할 수 있다. 위 조건은 1단계 연구의 FM2 모델에 대한 해석조건과 동일한 것으로 단층 미끄러짐을 포함한 수리역학적 거동에 대해서는 Park 외(2018)의 논문에서 자세히 기술한 바와 같다. 단, Fig. 17과 Fig. 18은 수리간극에 따른 압축률의 변화를 반영한 수정된 해석코드를 통해 도출한 결과이므로 그들의 연구에서 보고한 해석 결과와는 다소 차이가 있다.
4. 요약 및 결론
본 논문에서는 국제공동연구인 DECOVALEX-2019 프로젝트 Task B의 Benchmark Model Test(BMT) 연구내용을 소개하였다. Task B의 주제는 ‘Fault slip modeling’으로 단층 내 주입수로 인해 발생하는 전단파괴 및 수리역학적인 거동을 예측하기 위한 해석코드를 개발하는 데에 그 목적이 있다. BMT 해석은 1단계 연구를 통해 개발된 해석코드가 단층의 수리역학적 커플링을 적절히 반영하고 있는지 검토하기 위한 후속연구로서 주입수의 압력, 단층의 물성, 수리역학적 커플링 고려 여부 등에 따라 7가지의 시나리오 해석으로 구성된다. 본 연구에서는 단층의 거동을 보다 직접적인 방식으로 모사하기 위해 수리 모델(TOUGH2)에서는 실제 수리간극에 상응하는 두께를 갖는 연속체 요소를 이용하였으며, 역학적 모델(FLAC3D)에서는 불연속 인터페이스 요소를 이용하였다. TOUGH의 연속체 요소와 FLAC3D의 불연속 인터페이스 요소의 거동을 연동할 수 있는 모델링기법을 개발하고, 역학적 변형으로 인한 투수계수와 압축률의 변화를 고려하기 위한 커플링 모듈을 추가하였다. 제3장에서 살펴본 바와 같이, 단층 내 주입수의 유동 특성(유량, 압력)은 수리역학적 커플링 반영 정도에 따라 매우 상이한 결과를 보였다. 수리간극의 크기, 투수계수, 압축률의 변화를 모두 고려한 쌍방향 연동해석에서는 비연동해석에 비해 유량이 1000배 이상 증가하였고 주입공 주변으로 높은 압력 분포가 형성되었다. 현재 Task B의 연구결과가 계속 업데이트 중에 있고 본 연구의 해석기법 역시 지속적인 개선과 검증이 필요한 상태이다. 본 연구에서 개발한 해석기법은 Task B에 참여하는 국외 연구팀들과의 의견 교류와 워크숍을 통해 지속적으로 개선하는 한편, 향후 현장시험을 대상으로 하는 2단계와 3단계 연구를 통해 그 타당성을 검증할 예정이다.
