서 론

지중 유체 주입은 이산화탄소 포집, 활용 및 저장(carbon capture, utilization & storage, CCUS), 인공저류층 지열시스템(enhanced geothermal system, EGS) 구축, 수압 파쇄(hydraulic fracturing) 등 다양한 목적으로 수행된다(Stephens and Jiusto, 2010; Tapia et al., 2018). 최근, 주입정 위치 및 주입량과 지진 발생 시기 및 시점 사이에 대한 통계 분석 연구는 유체 주입 이후 지진 발생률의 급증을 제시하고 있다(Horton, 2012; Ellsworth, 2013; Keranen et al., 2014; Weingarten et al., 2015). 2006년 스위스 바젤(Basel)의 EGS 프로젝트에 의해 규모 2.6의 지진(Deichmann and Giardini, 2009), 2016년 미국 오클라호마(Oklahoma)에서 폐수 주입에 의해 규모 5.1의 지진(Goebel et al., 2017), 2019년 영국 프레스턴(Preston)에서 셰일가스 추출을 위한 수압 파쇄에 의해 규모 2.9의 지진(Cao et al., 2022)이 발생했다. 주입 유발지진은 인적/물적 피해로 이어질 수 있으며(Ellsworth, 2013), CCUS 사업의 경우 재활성화된 단층을 따라 이산화탄소가 누출될 가능성도 존재한다(Rinaldi et al., 2015). 이는, 지역 주민들에게 유체 주입에 대한 불안감을 느끼게 하며, 사회적 수용성을 감소시키는 원인으로 작용한다. 지중 유체 주입에 대한 사회적 반발은 네덜란드의 바렌드레흐트(Barendrecht) 이산화탄소 포집 및 저장(carbon capture and storage, CCS) 프로젝트(Terwel et al., 2012; van Egmond and Hekkert, 2015)나 독일의 바텐팔(Vattenfall) CCS 프로젝트(Dütschke, 2011) 등의 사례처럼, 사업 지연 등의 제한을 가할 수 있다(Brunsting et al., 2011; Kraeusel and Möst, 2012; Lee and Min, 2013; Haug and Stigson, 2016). 원활한 프로젝트 진행을 위해서는 유체 주입과 지구조 응력 사이의 역학적 관계에 대한 상세한 연구를 기반으로 사회적 수용성을 향상할 필요가 있다. 또한, 지중 유체 주입 중/후 발생한 미소지진으로 인해 유체 주입률을 감소시키거나(Cipolla et al., 2012), 주입을 중단하였으나 이후에 더 큰 지진이 발생하는 국외 사례(Yim and Min, 2022)가 있다. 안전한 지중 유체 주입을 위해서는, 유체가 저장되는 층의 수리지질학적 특성(예, 공극률, 투수계수 등)과 지질역학적 특성(예, 지층의 탄성계수, 단층 마찰계수 등)을 고려한 주입-유발지진(단층 미끌림) 발생 시나리오 구축이 필수적이다(Chang et al., 2013).

유발지진 발생의 대표적 메커니즘으로 (1) 주입된 유체가 단층면에 닿음으로써 가해지는 단층면 공극압의 직접적인 증가(예, Nolte et al., 2017; Yeo et al., 2020)와 (2) 유체가 단층에 닿지는 않았지만, 공극 골격의 부피변화에 따른 다공탄성 응력 섭동(예, Yeck et al., 2016; Chang and Yoon, 2020) 등이 제시된다. 다공탄성 이론(Biot, 1941; Segall, 1989; Piombo et al., 2005; Segall and Lu, 2015)은 다공성 매질 내에서의 유체의 흐름과 고체 변형 사이 상호작용에 대한 역학적 모형을 제공하며, 주입-유발지진 연구에 널리 적용되어 왔다(예, Deng et al., 2016; Safari and Ghassemi, 2016; Jin and Zoback, 2018). 본 연구에서는 다공탄성 및 단층 미끌림 결합 수치 모사 방법을 제안하고, 공극압 및 지중응력 변화, 표면 변위에 대한 계산을 수행하고자 한다.

수치 모사는 대상 지역의 수리지질학적 특성과 다양한 주입 시나리오(주입 유체의 종류, 주입량 등)에 대한 복합적인 고려가 가능하므로, 시공간적으로 관측이 제한된 지점에서의 결과값(예, 공극압 및 지중응력 변화량 등)에 관한 정량적 추정을 가능하게 해준다(Lee et al., 2008). 기존 다공탄성 연구는 유체 주입에 따른 공극압 및 지중응력의 변화 계산(Altmann et al., 2010; Chang and Segall, 2016)과 지진에 의한 다공탄성 반발 계산(Hughes et al., 2010; McCormack and Hesse, 2018)을 각각 독립적으로 수치 모사했지만, 정확한 결과값 도출을 위해서는 유체 주입과 지질역학적 특성(예, 단층 미끌림)을 결합한 수치 모사가 필요하다(Jha and Juanes, 2014).

본 연구에서는, 유한요소 상용 소프트웨어인 COMSOL Multiphysics^® (이하, COMSOL)를 사용하여 다공탄성(poroelasticity) 방정식을 약형식(weak form)으로 이산화하고, 파이썬(Python) 기반의 유한요소 소프트웨어인 PyLith(Aagaard et al., 2013)를 사용하여 단층 미끌림에 의한 변위를 계산하고자 한다. COMSOL은 수치 모형의 형태, 경계 조건, 물성 등을 자유롭게 조절할 수 있다는 장점이 있다(Lee and So, 2020). 뿐만 아니라, 다양한 솔버(solver)와 요소 차수(element order)를 제공(예, Hoffmann et al., 2009)하기 때문에 유체와 고체의 상호작용에 대한 다공탄성 수치 모사 연구에 적용된 다수의 예가 있다(Pampillón et al., 2018; Wetzler et al., 2019; Hui et al., 2021). 하지만 유한요소법은 연속체 역학을 바탕으로 하므로, 단층 미끌림이나 파괴에 대한 수치 모사 시 어려움이 있다(Zienkiewicz and Taylor, 2005; Agathos et al., 2016; Lee and So, 2019). 이에 반해 PyLith는 영역 내에 존재하는 단층면을 다수의 영역으로 분리하여 내부 마찰 경계로 간주하는 영역분해법(domain decomposition approach)을 이용해서, 불연속면이 미끄러지는 현상인 지진을 엄밀하게 수치 모사할 수 있다(Williams and Wallace, 2015; Kim et al., 2022; Xue et al., 2022). 우리는 공극 유체의 흐름과 다공매질의 연속체적인 변형 사이의 다공탄성 상호작용 계산에는 COMSOL을, 단층의 불연속적인 거동 모사에는 PyLith를 적용했다.

본 연구는 유체 주입/중지 상황 및 단층 미끌림 현상에 대한 수치 모사를 통해 응력 변화량, 변위 등을 계산하였다. 이러한 수치 모사는 단층 미끌림 이후 주입량의 감소 정도를 설정하거나, 주입 중단 이후 지중응력이나 지표 변위 등을 추정하는데 기여할 것으로 판단된다.

연구 방법

지배방정식

다공탄성(Poroelasticity)

Biot(1941)의 다공탄성 이론은 공극 유체 흐름과 암석 변형 사이의 상호작용에 대한 설명으로(Wang, 2000; McCormack et al., 2020), 주입-유발지진 연구(예, Ge and Saar, 2022) 및 지진 후 공탄성 반발 연구(예, Barbot and Fialko, 2010)에 널리 사용되어 왔다. 지층 내 인위적 유체 주입은 공극 내 유체 양을 변화시킴으로써 물질 거동의 변화를 일으킨다.

식 (1)과 (2)는 고체와 유체의 상호작용을 고려하는 등방성 선형 다공탄성 지배방정식의 구성방정식으로, 각각 힘평형방정식(force equilibrium equation)과 유체 질량 보존 방정식(fluid mass conservation equation)이다. p는 공극압, ηf는 유체의 점성도, κ는 투수계수, u는 고체의 변위장(displacement field)을 의미한다. 응력 텐서 σ는 식 (3)에서, 변형률 텐서 ε와 그 대각합이자 부피변형률인 tr(ε)=∇∙u, 배수(drained) 포아송비 ν, 전단계수 G, identity 텐서 I 등의 조합으로 표현된다. 식 (4), (5)를 이용해서 Biot-Willis 계수 α와 비저류계수(specific storage) Sε를 계산할 수 있다.

νu와 B는 각각 비배수(undrained) 포아송비와 Skempton 계수(Skempton’s coefficient)이다. 우리는 기준 모형(reference model)을 주 대상으로 수치 모사를 진행했으나, 주입층의 수리적/공탄성적 물성 변화에 의한 효과를 정량 분석하기 위해 매개변수 중 주입층의 투수계수와 Biot-Willis 계수를 기준 모형과 다르게 설정한 후 추가 수치 모사를 수행했다. 기준 모형에 대한 자세한 물성 수치는 Table 1에 제시했다. 또한, 본 연구에서는 다공탄성에 관한 다수의 기존 연구(Chang and Segall, 2016; Deng et al., 2016; De Simone et al., 2017; Rathnaweera et al., 2020)와 같이, 초기의 응력 상태를 자기 평형(self-equilibrium) 상태로 가정 후 공극압, 쿨롱 파괴 응력 등의 변화량에 대한 수치 모사를 수행했다. 따라서 중력에 의한 효과를 고려하지 않았으며, 중력항에 의한 영향에 대해서는 이후 부록을 통해 서술하도록 한다.

쿨롱 파괴 응력 변화량(Coulomb Failure Stress Change, CFS)

쿨롱 파괴 응력 변화량은 단층의 미끌림 발생 기준 설정 시 널리 사용된다(King et al., 1994; Rutqvist et al., 2007). 특정 단층의 쿨롱 파괴 응력 변화량 CFS가 음수인 경우, 단층면이 닫히는 상태를 의미하지만, 양수일 때는 단층 미끌림에 가까워진 상태를 의미하며(Shan et al., 2013) 식 (6)으로 표현된다.

τ, σn, μ, p는 각각 전단 응력 변화, 수직 응력 변화, 마찰 계수와 공극압 변화로, σn의 부호가 양수이면 단층이 닫히고 있는 압축상태임을 지시한다. 본 연구에서는 초기 상태를 자기 평형 상태로 가정하고, 중력을 무시한다고 가정함에 따라, 수치 모형의 공극압과 응력 텐서의 초기 조건은 모두 0으로 설정하였다. 따라서 본 연구를 통해 계산된 τ, σn, p 등은 변화량을 의미한다. 우리는 단층 미끌림 발생 여부의 판단 기준인 임계 쿨롱 파괴 응력 변화량 CFS_C를 0.1 MPa로 설정했다. CFS가 1 bar (= 0.1 MPa)인 경우 미끄러짐 임계점에 도달한 것으로 간주하는 기존 지구물리학 연구들(예, Papadimitriou, 2002; Bagge and Hampel, 2016)과 유사한 기준이다. 유체 주입에 의한 공탄성 효과와 공급압의 효과로 인해 단층의 CFS가 CFS_C = 0.1 MPa에 도달하는 시점이 지진 발생 시점이며, 이때부터는 PyLith를 주 소프트웨어로 사용해 지진에 의한 변위를 계산한다.

수치 모형 설정

우리는 유체 주입에 의한 공극압 진화 및 단층 미끌림 수치 모사를 위해, Chang and Segall(2016)의 모형 및 매개변수를 바탕으로 수치 모형을 제작했다. 수치 모형을 구성하는 이암, 사암(주입층), 기반암의 두께는 각각 0.6 km, 0.1 km, 2 km이며, 총 25 km (수평방향) × 2.7 km (수직방향)인 2차원 모형이다. 주입된 유체가 지표로 누출되는 현상을 방지하기 위해 사암에 비해 상대적으로 투수성이 낮은 이암을 사암층 위에 위치시켰으며(Shukla et al., 2010), 이는 실제 유체 저장층의 환경과 유사하다(Rutqvist et al., 2010). 기반암 내에 60°의 경사를 갖는 단층은 선 형태의 미끌림면(slip surface)과 이를 둘러싼 면 형태의 단층핵(fault core)으로 표현하였다. 힘평형방정식(식 (1))의 경계 조건으로 상부 경계는 자유면(free surface) 경계 조건, 그 외 경계는 롤러(roller) 경계 조건이 적용되었다. 유체의 질량 보존 방정식(식 (2))에 대한 경계 조건은 하부 경계는 무흐름(no flow) 경계 조건, 그 외 경계는 디리클렛(p = 0) 조건을 부여하였다. 사암층 내, 직경 0.2 m의 주입정을 정의했으며(Fig. 1a의 확대 그림), X = ±0.1 m인 지점의 수직선에서 주입 시나리오에 따라 공극압 값을 조절해 유체의 주입을 모사했다. 디리클렛 조건으로 설정된 주입압력은 Fig. 1c와 같이, 30일 동안 25 MPa의 일정한 압력을 가했다. 본 연구 수행을 위한 수치 모형 및 경계 조건은 Fig. 1a, 수치 모형을 구성하는 격자 구조는 Fig. 1b, 주입 시나리오는 Fig. 1c에 나타내었고, 각각의 층 및 단층핵의 물성 값은 Table 1에 명시하였다. 총 481,932개의 DOF로 이루어진 수치 모형은 삼각형 선형 요소로 구성되었으며, 격자의 최대 크기와 최소 크기는 각각 50 m와 0.3 m이다. 유체 확산 및 단층 미끌림이 발생하는 주입층과 단층핵의 경계는 주변보다 더욱 세밀한 격자 구조를 갖도록 설정함으로써(Fig. 1b) 수치 모사 결과의 정확도를 제고했다. PyLith에서 사용한 수치 모형의 경계 조건과 물성은 COMSOL과 동일하다. 단층 미끌림은 0.5 m, 1 m, 2 m, 3 m 등 다양한 단층 미끌림 길이로 수치 모사하였고, 지표면에서의 수직/수평 변위와 이암층에서의 본 미세스 응력(von Mises stress) 변화를 계산하였다.

Fig. 1.

Model set-up. (a) Numerical model domain with three layers: mudrock, sandstone, and basement. The fault core is 2 km from the injection well at X = 0 km in the sandstone layer. The plane and line represent the fault core and slip surface, respectively. The three points within the slip surface are denoted by Points A, B, and C to track the time evolution of pore pressure and Coulomb failure stress changes. (b) Magnification of the mesh structure near the junction of the sandstone layer and fault. (c) Injection scenario continues with constant injection pressure of 25 MPa for 30 days.

COMSOL과 PyLith 순차 적용 알고리즘

본 연구의 수행 순서는 다음과 같다. (1) COMSOL을 사용하여 일정 기간 동안의 유체 주입과 중지 상황을 수치 모사하여, 시간에 따른 주입 중/후의 공극압, 변위, 응력 등을 계산한다. (2) 단층에 가해지는 공탄성 응력이 사용자가 지정한 임계점에 도달하는 경우, 단층 미끌림에 의한 수직 및 수평 변위를 수치 모형에 삽입한다. 해당 수직 및 수평 변위는 PyLith에 지층 물성(예, 영률, 포아송비 등), 단층 위치, 단층 각도, 미끄러짐 양 등을 입력하고 유한요소 수치 계산하여 구한다. (3) 단층 미끌림 이후의 공극압 및 응력을 계산한다. 각 소프트웨어의 적용 순서 및 변경 시점 등의 알고리즘은 Fig. 2에 제시했다. CFS 분포, 유체 속도장, 공극압 변화, 계산 영역의 변위 등에 관한 시각화 및 후처리는 사용자가 원하는 변량에 관한 시각화 자료를 생성할 수 있는 소프트웨어인 ParaView(Ahrens et al., 2005)를 이용해 수행했다.

Fig. 2.

Flow chart for simulating pore pressure evolution, poroelastic stressing, and fault slip in the fluid injection model. COMSOL (blue) and PyLith (green) compute the effects of poroelasticity and fault slip, respectively. To build a detailed poroelastic model, we select an option to customize the partial differential equations in the first step. In the discretized domain with the numerical conditions (i.e., boundary and initial conditions) and properties (e.g., permeability and fluid viscosity), the pore pressure and poroelastic stress fields are redistributed by injection for 30 days. When the Coulomb failure stress change of the fault overcomes the pre-defined critical value CFS_C, the main numerical software is switched from COMSOL to PyLith to simulate the fault slip and the resulting coseismic stress distribution.

결과 및 토의

COMSOL을 이용한 유체 주입에 의한 공극압 및 응력 상태의 변화와 지표에서의 변위 계산

유체 주입 및 중지에 따른 이암층, 사암층(주입층), 기반암층에서의 공극압 변화 p(Fig. 3a~3d)와 수직 변위 v(Fig. 3e~3h)의 시계열 진화를 수치 모사했다. 초기 30일 간 주입층에 25 MPa의 주입압력을 주고(Fig. 3a and 3b 및 Fig. 3e and 3f), 200일(Fig. 3c and 3d 및 Fig. 3g and 3h)까지의 p와 v의 진화를 관찰하였다. 유체 주입이 활발한 시기에는 사암층의 p가 주입정 근처에서 가장 큰 값을 갖으며, 투수율이 높게 설정된 단층핵에 p의 상승이 국지화되고 음의 Y 방향으로 확산된다. 주입이 중지된 후에는, 전체적으로 p가 시간에 따라 완화하지만, Fig. 3b와 Fig. 3d의 단층 뿌리 부분의 p를 비교해보면, 단층의 뿌리 부분은 주입이 중지되어도 p가 상승하는 구간이 존재함을 파악할 수 있다. 주입 시나리오(Fig. 1c 참조)에 따라 주입정에서 최대 압력(25 MPa)에 도달한 시기에는 주입정 바로 위인 X = 0 km 지점에서 ~3 cm로 가장 높은 표면 수직 변위를 보여준다(Fig. 3f). 주입정에 가해지는 압력이 주입 시나리오(Fig. 1c)에 따라 주입이 중단되면, 표면 수직 변위의 크기는 점점 감소한다(Fig. 3g and 3h).

Fig. 3.

Time evolution (10 to 200 days) of (a-d) pore pressure and (e-h) vertical displacement owing to fluid injection. The residual pore pressure changes after injection stops on the 30th day. Vertical displacements along the surface (dotted lines) are exaggerated by 5,000 times to emphasize surface deformation.

단층 안정성을 정량화하는 쿨롱 파괴 응력 변화 CFS와 지중 유체 주입 사이의 관계를 이해하기 위해, 주입 중/후 시점에 따른 모형의 전체적 CFS분포 및 단층 내 CFS분포를 Fig. 4에 도시했다. 대상 단층이 양의 CFS구역에 존재할 경우 단층의 파괴 가능성을 지시하므로, 단층과 그 주변부 CFS분포 양상에 관한 정량적 추정은 주입-유발지진 위험성을 최소화하기 위한 지중 유체 주입 시 필수적이다. 본 연구에서는 0.1 MPa (예, Papadimitriou, 2002; Bagge and Hampel, 2016)을 단층 미끄러짐을 위한 임계 쿨롱 파괴 응력 변화 CFS_C로 설정하였다. 유체 주입이 지속되는 30일 동안 CFS는 사암층을 통해 주입정에서부터 수평방향으로 증가하며, 이암과 기반암에 비해 높은 투수율을 갖는 단층핵에서도 CFS가 증가했다(Fig. 4a and 4b). 유체 주입 중단 이후(즉, 30일 이후), 공극압이 확산함에 따라(Fig. 3c and 3d 참조) 사암층 주변부(X = ±6 km)와 단층핵의 뿌리 부분(Y = -2 km)에서는 100일 및 200일 시점에도 양의 CFS가 잔류한다(Fig. 4c and 4d). 우리는 단층핵 내 3개의 지점(Points A, B, and C)에서 CFS의 시간 진화를 비교했다(Fig. 4e). 3개 지점의 CFS 최대값은 A의 경우 ~120일에 ~3.5 MPa이며, B와 C는 300일 이후에도 꾸준히 증가하는 것으로 계산되었다. A, B, C 지점의 CFS는 각각 ~13일, ~21일, ~32일 시점에 임계점 CFS_C에 도달하며(Fig. 4e), 단층핵 내의 CFS 분포 양상(Fig. 4f)도 사암층과 근거리인 점에서 더 빠르게 CFS_C에 도달함을 보여준다.

Fig. 4.

Coulomb failure stress change (CFS) induced by fluid injection. (a-d) Evolution of CFS with time from (a) 10 to (d) 200 days. (e) CFS changes during 300 days measured at Points A (pink), B (green), and C (orange). The gray dotted line indicates the threshold, 0.1 MPa. (f) Variations of CFS changes with time within the fault core.

CFS는 공극압 관련 항 μp와 다공탄성 응력(poroelastic stressing) 관련 항 τ-μσn으로 구성되어 있다. 두 항 모두 더 큰 양수일수록 CFS값도 더 크다. Fig. 5a~5d는 μp 분포를 각각 10일, 30일, 100일, 200일 시점에서 보여준다. 주입이 진행되는 동안 주입정을 중심으로 사암층을 따라 μp가 확산되며, 단층을 통해서도 μp가 확산되는 현상이 나타난다(Fig. 5a and 5b). 주입이 종료된 이후에도 μp는 사암층과 단층 내에 여전히 존재한다(Fig. 5c and 5d). 양의 τ-μσn 값은 단층에 작용하는 수직 응력보다 전단 응력의 크기가 더 큼을 의미하고, 이는 단층이 미끄러지기에 용이한 상태임을 지시한다(Fig. 5e~5h). 유체를 주입하는 30일 동안 주입정 주변의 사암층에서는 음의 τ-μσn 값이 나타나 사암층이 압축 중임을 지시하고, 반대로 단층에서는 양의 τ-μσn값이 나타나 단층의 미끄러짐 가능성을 암시한다(Fig. 5e and 5f). 이후, 주입이 종료된 후, 100일과 200일에서 사암층과 단층에서는 주입 중 시점과 유사하게 각각 음의 τ-μσn값과 양의 τ-μσn값이(Fig. 5g and 5h) 나타난다.

Fig. 5i와 5j에 각각 Points A, B, and C지점에서 시간에 따른 μp와 τ-μσn 값 변화 그래프를 도시했다. 주입 기간(~30일) 중 초반부에는 μp값이 음수이며, 중반부터는 증가하다가 결국 양의 값을 보인다(Fig. 5i). 이는, 주입 유체의 도달 시점과 μp가 관련이 있음을 지시한다. 주입층인 사암층과 상대적으로 가장 가까운 Point A에서는 주입 시작 후 ~15일부터 μp가 상승하며 ~130일 시점에서 최대 μp값 ~2 MPa가 나타났다. 반면, 사암층에서 가장 멀리 떨어진 Point C에서는 주입 후 ~30일부터 상승하는 추세를 보였으며, 300일 이후에도 μp값은 지속적으로 상승했다. τ-μσn는 주입 시작과 동시에 상승한다. ~13일까지는 사암층에서 멀리 있을수록 더 큰 값을 보이지만, 이후 상대적으로 사암층과 가까이에 위치한 Point A가 가장 큰 값을 보이고 ~115일 시점에서는 최대 τ-μσn값이 ~1.6 MPa에 이른다. ~15일 시점에서 Point B는 Point C보다 더 큰 값을 갖기 시작하며, ~280일 시점에 1.2 MPa의 최대 τ-μσn값을 갖는다. 사암층으로부터 가장 멀리 떨어진 Point C는 300일 이후에도 τ-μσn값이 상승하는 경향을 보였다.

Fig. 5.

Evolution of (a-d) μp and (e-h) τ-μσn, which are terms related to pore pressure and poroelastic stressing, respectively. (i, j) Daily variations of each value measured at Points A (pink), B (green), and C (orange) in the slip surface.

Fig. 6a and 6b는 단층핵과 사암층(주입층)의 물성 조합(투수계수 및 Biot-Willis 계수)에 따라 단층의 최하단 Point C지점(Fig. 1 참고)의 CFS가 임계점 CFS_C에 도달하는 시간(critical time)을 도시했다. Fig. 6a의 가로축은 단층 투수계수의 로그값, 세로축은 사암층 투수계수의 로그값이다. 이 때 사용된 사암층의 Biot-Willis 계수는 기준 모형(reference model)의 물성(Table 1 참조)을 이용해 계산하고 고정하였다. Fig. 6a는 단층핵의 투수계수가 ~1 × 10^-14.5 m²보다 클 경우, 단층핵과 사암층의 투수계수가 클수록 더 빠르게 CFS_C에 도달함을 보여준다. 예를 들어, 단층의 투수계수가 1 × 10^-13 m²이고, 사암층의 투수계수가 6 × 10^-13 m²인 경우, Point C에서의 CFS는 10일 이내에 CFS_C에 도달하지만, 단층과 사암층의 투수계수가 작을 경우(예, 각각 1 × 10^-14 m², 6 × 10^-15 m²인 경우)에는 100일 이후에 도달한다. Fig. 6b의 가로축은 단층 투수계수의 로그값, 세로축은 사암층 Biot-Willis 계수이다. 이때 사용된 사암층의 투수계수는 기준 모형 값인 6.4 × 10^-14 m²으로 고정하였다. Fig. 6b의 도달시간에 대한 등고선의 형태는 CFS의 CFS_C 도달 시간에는 주입층의 Biot-Willis 계수보다는, 단층의 투수계수가 더 큰 영향을 미침을 암시한다. 예를 들어, 단층의 투수계수를 1 × 10^-14 m²으로 고정하고, 주입층의 Biot-Willis 계수가 0.1부터 1.0까지의 범위 내에 존재하는 경우, CFS_C에 도달하는데 약 26일부터 약 46일의 시간이 소요되지만, 주입층의 Biot-Willis 계수를 0.6으로 고정 후, 단층의 투수계수를 1 × 10^-15 m²부터 1 × 10^-13 m²로 설정할 경우, CFS_C까지 도달하는데 걸리는 시간은 약 11일부터 약 326일 사이다.

Fig. 6.

Maps of critical time at which the CFS value reaches CFS_C (i.e., 0.1 MPa) at Point C (see Fig. 5a). (a) The critical time map varied with the permeability of sandstone (horizontal axis) and the fault (vertical axis). (b) Variation of critical time with the Biot-Willis coefficient of sandstone and the permeability of the fault when sandstone’s permeability is 6.4 × 10^-14 m². The dotted lines connect points with the same critical times. Other material properties are identical to the reference model (see Table 1).

PyLith를 이용한 단층 미끌림 계산

단층 최하부(Fig. 5a의 Point C지점)에 가해지는 CFS값이 임계값인 CFS_C에 도달하면, PyLith를 이용해 0.5, 1, 2, 3 m의 역단층 미끄러짐을 수치 모사했다. 단층의 가장 아래 부분을 기준으로 삼은 이유는 단층의 asperity가 단층 최하부에 존재한다고 가정하여, 가장 늦은 미끄러짐 시점을 추정하기 위해서다. 만약 기준을 Point A(Fig. 5a)로 설정하면 더욱 이른 미끄러짐 시점을 추정할 것이다. Fig. 7a~7d는 단층 미끌림 크기에 따른 본 미세스(von Mises) 응력 분포를 보여준다. 본 미세스 응력은 편차응력텐서를 대표하는 양으로 재료의 파괴기준 설정에 주로 사용되며, 높은 본 미세스 응력은 높은 파괴 가능성을 지시한다. 직관적인 응력 분포 및 변위 관찰을 위해 수직 방향으로 300배 확대하여 표현하였다. 단층의 미끄러진 길이가 증가할수록 더 큰 변위와 본 미세스 응력이 도출되었음을 파악할 수 있다. 정량적인 판단을 위해서 Fig. 7e에 단층 미끌림 크기에 따른 지표 수평 변위(실선)와 지표 수직 변위(점선)를 도시했다. 단층이 0.5 m 미끄러진 경우, 최대 수평 변위는 -0.03 m, 수직 변위는 0.11 m이다. 단층이 3 m 미끄러진 경우, 최대 수평 및 수직 변위는 각각 -0.21 m와 0.67 m로 단층 미끄러짐 양이 클수록, 수평 및 수직 변위가 증가했다. Fig. 7f는 단층 미끌림 크기에 따른 사암층 내 수평 측선(Y = 0.4 km)에서 추출한 본 미세스 응력 프로파일이다. 미끌림에 의한 최대 본 미세스 응력은 단층 주변에서 나타나며, 단층의 미끌린 양이 증가할수록 더 큰 본 미세스 응력이 발생한다.

Fig. 7.

von Mises stresses and vertical displacements with difference the amount of fault slip (a-d). The developments of von Mises stress and vertical displacements with increasing fault slip displacements of 0.5, 1, 2, and 3 m. The surfaces of the model are deformed proportionally to the vertical displacements (which are exaggerated here by a factor 300). (e) Horizontal (solid lines) and vertical (dashed lines) displacements along the surface with different slip amounts. The red, orange, green, and blue indicate the slip amount of 0.5, 1.0, 2.0, and 3.0 m respectively. The gray lines represent the location of the fault contacted by the sandstone layer (X = 2 km). (f) von Mises stress within the mudstone layer (Y = 0.4 km).

토의 및 결론

본 연구에서는 공극 유체의 질량 보존 방정식과 공극 골격의 힘평형방정식을 결합한 다공탄성 모사에 유용한 COMSOL과, 영역분해법 바탕의 단층 미끄러짐 모사에 특화된 PyLith를 순차적으로 적용하여 유체를 지중 주입하는 수치 모형을 제작하였으며, 주입에 의한 공극압, 지중응력 변화, 단층 미끌림에 의한 변위 및 응력장 변화를 계산하였다. 유체 주입압력을 시간에 대해 일정하게 유지하면, 주입정 주변의 공극압 및 쿨롱 파괴 응력 변화, 수직 변위량이 증가하는 양상을 보이는데, 이는 유체 주입이 장기간 지속될 경우 지중 암반에 가해지는 응력이 증가함으로써, 불안정한 암반 상태로 이어질 수 있음을 의미한다. 주입이 종료된 후에도 주입 중 발생한 공극압과 응력이 주입정의 근거리에서부터 원거리까지 확산 및 전달되며, 단층의 경우에는 유체 주입 중단 이후의 쿨롱 파괴 응력 변화 값이 주입 중에 비해 오히려 더 큰 값을 보이는 것을 확인하였다. 이는, 유체 주입의 중단이 지중 암반의 안정적인 상태를 담보하는 것이 아니며, 주입 중지 후에도 지속적인 모니터링이 필요함을 제안한다.

사암층(주입층)과 단층 물성의 특성(예, 투수계수와 Biot-Willis 계수 등)은 유체 주입에 의해 단층의 쿨롱 파괴 응력 변화량이 임계점에 도달하는 시간을 결정하는 요인이다. 사암층과 단층의 투수계수가 모두 높은 경우, 단층의 최하부인 C지점에 가해지는 쿨롱 파괴 응력 변화량이 임계점에 도달하는데 걸리는 시간은 10일 이내지만, 사암층과 단층의 투수계수가 낮은 경우에는 200일 이후에 도달한다. 사암층의 Biot-Willis 계수도 쿨롱 파괴 응력 변화량 임계점 도달 시간에 영향을 미치지만, 단층의 투수계수에 비하면 그 효과가 작음을 확인했다. 이는, 주입정 위치 탐색에 있어 주입층과 주변 지층의 물성, 특히 투수계수를 파악하는 것의 중요성을 지시한다. 유체 주입에 따른 응력 변화와 단층 미끌림에 대한 다양한 시나리오 설정 및 다공탄성 수치 모사는 주입층 주변의 지층 파괴 가능성을 파악하고, 유체 누출 위험성을 저감하는 데 중요할 것으로 판단된다.