1. 서 론
유체 주입은 셰일 가스 개발, 인공저류층 지열시스템(Enhanced Geothermal System, EGS), 이산화탄소 지중저장과 같은 사업들에 수반되는 주요 공정이다(Ellsworth, 2013, Kim et al., 2019). 유체 주입중 발생하는 미소지진은 유체의 유동 범위를 파악할 수 있게 하지만 주변에 피해를 야기하는 큰 규모의 지진(Large Magnitude Earthquake, LME)이 세계 여러 지역에서 다수 보고되어 이에 관한 메커니즘 규명은 과학기술계의 큰 해결과제이다. (Davis and Frohlich, 1993, Ellsworth, 2013) 유체 주입 중에 발생하는 LME가 문제가 된다고 판단되면 유체 주입을 중단시키지만(Majer et al., 2012, Kim et al., 2018) 그 이후에 이전보다 더 큰 규모의 LME가 발생하는 사례가 종종 보고되었다(Charléty et al., 2007, Baisch et al., 2009, Mukuhira et al., 2017, Park et al., 2020). 예를 들어 스위스 바젤(Basel)에서 수행된 EGS 사업의 경우 유체 주입 중단 중에 발생한 LME에 의해 사업이 중단되었으며, 우리나라 포항 EGS현장에서도 규모 5.4의 지진이 유체 주입 중단 이후에 발생한 바 있다(KGC, 2019).
EGS현장에서 유체 주입 중단 후에 발생하는 큰 규모의 유발 지진은 몇 가지 특징을 가지고 있다. 이러한 큰 규모의 지진들의 진원 메커니즘(focal mechanism)을 분석하면 인장 메커니즘보다는 전단 메커니즘을 따르며 전단면은 지진 군집(seismic cloud)이 이루는 면과 평행하다(Deichmann et al., 2009, Charléty et al., 2007). 주입 중단 후에 지진의 발생이 점점 천부로 이동하는 반면, 큰 규모의 지진들은 주입점보다 깊거나 얕은 심도에서 모두 일어났다(Schoenball et al., 2013, Mukuhira et al., 2017). 스위스 바젤과 프랑스 슐츠에서 유체 주입에 따른 유발지진의 진원 메커니즘과 위치가 자세히 조사되었는데, 이에 따르면 LME는 주입점이 포함된 전단면 상에 존재하지 않았고 여러 면 구조를 통해 연결된 것으로 파악되었다. 또한 주입 중단 후에 LME는 주입 중에 발생한 지진 군집의 가장자리에 주로 발생하는 것으로 관측되었다(Dorbath, 2009, Kraft et al., 2014, Mukuhira et al., 2017).
주입을 중단한 상태에서 주입을 하는 상태보다 LME의 규모가 더 커지는 문제에 대하여 많은 연구에 의해 다양한 설명들이 제시되어왔다. 유체 주입을 중단한 후에 공극압의 구배는 완만해지지만 공극압은 계속 전파하여 전단파괴가 발생할 수 있는 영역이 커지는 것으로 설명한 사례가 있으나(Baisch et al., 2009, Mukuhira et al., 2017) 이러한 설명은 개념적인 설명으로 지진의 규모 증가를 설명하는 데는 한계가 있다(Fig. 1 참조). 또한 전단이 동시에 발생하는 영역이 어느 정도 이상이 될 때 단층면 전체로 전단이 불안정 전파(Unstable rupture propagation)하게 되는 것을 큰 지진이 발생하는 메커니즘으로 보는 파괴역학적 설명도 있다(Galis et al., 2017). 통계적으로 많은 지진이 일어나면 규모가 큰 지진이 일어날 확률이 올라간다는 사실을 이용하여 주입이 중단되고도 지진은 계속 발생할 수 있기 때문에 큰 규모의 지진이 발생 가능성이 사라지는 것은 아니라는 설명도 있다(Segall and Lu, 2015, Langenbruch et al., 2020). 그러나 이러한 설명들은 유체 주입 중단과 지진 규모의 물리적 관계를 직접 설명하지 않는다. 또한 수치해석을 이용하여 현장을 모사한 연구들은 주입 중단 이후에도 지진이 발생할 수 있음을 보였지만 규모가 더 커지는 것은 모사하지 않았다(McClure, 2015, De Simone et al., 2017, Kim et al., 2022).
본 연구의 목적은 유체 주입을 중단하는 것이 어떻게 유발 지진 규모에 영향을 미치는지 수리역학적으로 분석하는 것이다. 수리적 해석해와 수리역학 복합 수치해석을 이용하여 주입을 중단하지 않은 경우와 중단한 경우를 비교 분석함으로써 개념적으로 제시된 공극압 구배의 효과를 물리적으로 확인하고 유발 지진의 규모에 미치는 영향을 확인하고자 한다.
Fig. 1.
A schematic diagram of pore pressure distribution and shear slip area increase after the suspension of fluid injection (Mukuhira et al., 2017)
2. 수리적 해석해를 이용한 분석
2.1 모델 설명
유체 주입에 의한 유발 지진은 역학적인 측면에 대한 고려없이 수리적인 측면에서 공극압의 확산으로 설명되기도 한다. 지진기반 저류층 특성화(Seismicity-based reservoir characterization, SBRC)는 유체 주입에 따라 유발된 미소지진의 시공간적 분포를 공극압의 확산으로 해석하여 수리적 확산도(hydraulic diffusivity)를 구하는 방식으로써 공극압에 대한 확산 메커니즘이 현상 관측결과를 잘 설명하는 모델임을 알 수 있다(Shapiro et al., 2015). 본 연구에서는 아래와 같은 쿨롱파괴조건식(식 (1))을 만족할 경우 전단이 발생하면서 전단 팽창에 의해 수리적 간극이 커지는 것을 반영한 수리모델을 이용하여 유체 주입 중단의 효과를 분석하고자 한다.
여기서 p는 공극압, pcri는 임계 공극압, σn은 수직응력, τ는 전단응력, f는 마찰계수이다.
Fig. 2는 불투수성 암반 내에 단일 절리면 상에서 축대칭 확산에 의해, 전단이 발생한 영역을 수리적으로 자극된 영역(stimulated zone)으로, 아직 전단이 발생하지 않은 영역을 수리적으로 자극되지 않은 영역(unstimulated zone)으로 구분한 수리적 해석해의 개념적 모델을 보여준다. 절리상에서 어떤 한 영역의 공극압이 임계 공극압에 도달하여 전단이 발생하면, 정적마찰력보다 운동마찰력이 작기 때문에 준정적 상태에서는 운동 마찰력과 평형을 이룰 때까지 전단이 진행된다. 정적 마찰력과 운동 마찰력을 각각 임계 공극압의 상한과 하한으로 대응시켰고, 수리 자극된 영역의 경계 부분(r = R)에서 임계 공극압 상한에 도달할 때 수리 자극된 영역의 경계를 임계 공극압 하한까지 확장시킴으로써 시간에 따른 전단의 발생과 수리 자극 영역의 팽창을 모사하였다. 수리 자극된 영역이 확장된 직후에는 늘어난 경계 부분에서 아직 공극압이 임계 공극압 상한에 도달하지 못하기 때문에 일정 시간이 지난 뒤에서야 다시 임계 공극압 상한에 도달하면서 전단이 발생하고 수리 자극된 영역이 확장된다(Fig. A1). Table 1은 본 해석에 쓰인 입력값을 나타내고 Fig. 3은 유체를 일정 유량으로 주입하고 118시간이 지났을 때의 해석해에 의한 공극압의 구배를 보여준다(Appendix; Yim, 2018). 파선으로 된 상자로 나타낸 바와 같이 수리 자극된 영역과 수리 자극되지 않은 영역의 경계가 위치함에 따라 공극압의 공간적 기울기에 급격한 변화가 발생함을 알 수 있다. LME는 주로 지진 군집의 경계부에 위치했다는 관측결과에 기반하여 주입 중에는 수리 자극된 영역에 높은 공극압이 유지가 되지만 경계부에서 새로 생기는 수리 자극 영역의 증가분만을 전단 영역(shearing area)으로 가정하고 시간에 따른 전단 영역의 변화를 관찰하였다. 결과에서는 일정 유량으로 주입하다가 5일째에 중단하는 경우(shut-in)와 중단 없이 계속 주입하는 경우(continuous injection)를 비교하였다.
Table 1.
Input parameters for hydraulic analysis
2.2 결 과
수리적인 해석에 의하면 수리 자극된 영역의 경계로 나타내어지는 압력 선두(pressure front)는 유체 주입을 중단하더라도 수 시간 이후까지 임계 공극압 이상에 도달하면서 계속 진행하였으나 그 진행 속도가 감소하였다(Fig. 4). 장시간 동안 경계부에서 공극압이 증가할 수 있는 이유는 수리 자극되지 않은 영역의 작은 수리적 간극으로 인해 저류율과 투수율이 수리 자극된 영역에 비해 현저히 작으므로 유체 주입을 중단하더라도 수리 자극되지 않은 영역으로 흘러나가는 유량이 상대적으로 적어 수리 자극된 영역의 공극압이 유지될 수 있기 때문이다. 유체 주입이 중단되면 공극압을 바깥방향으로 확산시키는 동력이 사라지는 것이므로 압력 선두의 진행 속도 감소, 즉 전단 영역 발생이 보다 늦게 일어나게 된다. 이는 유체 주입을 중단하더라도 지진이 계속 발생하지만 그 발생 빈도는 감소한다는 관측 사실의 설명이 될 수 있다. 또한 압력 선두에서 발생하는 전단 영역의 크기는 유체를 주입함에 따라 계속해서 증가하는데 이는 지진 규모와 유체의 주입부피가 양의 관계를 가진다는 사실(McGarr, 2014)과 연관될 수 있다. Fig. 4는 유체 주입을 중단하는 것이 시간에 따른 전단 영역의 증가를 가속화하는 것을 보여준다. 전단 영역의 크기의 증가는 본 해석에서 수리 자극 영역의 경계에서 공극압 구배의 경사가 완만해진다는 의미이며, 유체 주입 중단이 공극압 구배를 더 빨리 완만하게 만들기 때문에 전단 영역의 증가가 가속화되는 것이고, 이는 Mukuhira et al.(2017) 등이 제시한 Fig. 1의 개념적 설명을 수리적으로 보여주는 것이다(Fig. A1). 한편 주입을 유지한 경우 약 7일이 경과했을 때 발생할 전단 영역의 크기가 5일째에 주입을 중단함으로써 약 5.6일에 발생한 것을 보여주며 이를 발생가능한 지진의 규모와 연관지어 생각할 수 있다.
수리적인 해석만으로도 관측 사실을 정성적으로 모사할 수 있으나 본 해석에는 몇 가지 한계점이 있다. 먼저 수리 자극된 영역을 확장시킬 때 갑작스런 수리적 간극의 증가에 의한 공극압의 감소와 같은 동적 거동은 무시가 되었다. 또한 절리 외에 암반을 통한 확산은 무시하고 균질한 단일 절리만을 모사하였다. 그리고 역학적 거동은 고려되지 않았기 때문에 전단에 의한 응력 재분배 등은 고려되지 않았다.
3. 수리역학 복합 수치해석을 이용한 분석
3.1 모델 설명
다공탄성학적 해석을 통해 공극압의 직접적인 변화뿐만 아니라 공극압이 주변 암반의 응력 상태를 변화시키는 것이 유체 주입에 따른 유발 지진의 특성들을 모사하는 데 큰 영향을 미친다는 것이 여러 수치해석 연구에서 보여진 바 있다(Segall and Lu, 2015, De Simone et al., 2017, Chang et al., 2020, Kim et al., 2022). De Simone et al.(2017)은 전단에 의한 주변 응력의 변화가 크다는 것을 보였고 Kim et al.(2022)은 전단에 의한 영향을 통해 포항에서 발생한 유발지진의 발생 특성을 모사할 수 있었다. 따라서 수리역학 복합 해석은 유체 주입에 의한 유발 지진을 분석할 때 필요하며, 앞 절에서 수행된 수리적 해석의 한계를 극복하여 균열의 수리역학적 특성, 방향, 응력집중 등 역학적 영향을 복합적으로 고려하기 위해, 본 절에서는 수리역학 복합 해석을 수행하였다.
본 연구에서는 불투수성 암반 내 절리를 따라 유체 유동이 발생한다는 가정하에 2차원 개별요소법을 이용한 소프트웨어인 UDEC을 사용하였다(Itasca, 2011). 바젤에서 조사된 진원 메커니즘과 위치 정보를 참고하여 규모 3.4의 지진이 발생한 단층면(Case 1 단층)을 모사하였고 비교를 위해 주입 지점에서 멀어지며 발달했지만 같은 임계 공극압을 가지는 단층면(Case 2 단층)을 가정하여 Fig. 5와 같은 단순화된 두 모델을 구성하였다. Table 2는 수치해석에 입력된 변수로 Basel 현장의 사례를 참고하였으며 필요 시 일반적인 값을 가정하였다. 주입공정의 경우 초기 4시간 동안은 급격한 공극압의 증가를 막기 위해 작은 유량으로 주입한 뒤, 그 이후에는 유량을 증가하여 주입하였다. 비교를 위해 주입 후 16시간 후에 유체주입을 중단한 경우(Shut-in)와 유지한 경우(continuous injection)를 비교하였다. Case 2 단층은 초기 응력에 의해 Case 1 단층과 같은 1.5MPa의 임계 공극압을 가지지만 방향은 다른 단층으로 여러 단층의 기하학적 관계가 어떻게 유발 지진에 영향을 미칠 수 있는지 확인하기 위한 가상의 단층이다. 수치해석상에서 절점에서 최소 전단 속도(minimum shear velocity)이상 일 때 이를 본 해석의 지진으로 간주하였고, 여러 절점에서 지진이 발생한 경우 서로 이웃한 노드끼리는 하나의 지진으로 보았다(Yim, 2018). 각 지진에 대응하는 환산 규모는 식 (2)에 따라 계산하였다(Hanks and Kanamori, 1979).
여기서 G는 단층의 전단 강성, d는 전단 변위, A는 전단이 발생한 면적으로 절점에 해당하는 길이와 현재 2차원 모델의 두께(L)의 곱으로 나타내었다.
Table 2.
Input parameters for UDEC numerical simulation
| Parameters | Value | Reference | |
| Flow rate, q |
0.139E-5 m3/s/m (0-4h) 3.333E-5 m3/s/m (After 4h) | ||
| Residual hydraulic aperture, amin | 1.00 μm | ||
| Maximum hydraulic aperture, amax | 120 μm | ||
| Aperture at zero normal stress, a0 | 120 μm | ||
| Viscosity, μ | 0.1421cp | ||
| Dilation angle, φ | 5° | ||
| Static friction angle, θs | 30° | De Simone et al., 2017 | |
| Dynamic friction angle, θd | 28° | ||
| Cohesion, c | 0 MPa | ||
|
Stiffness characteristics (dkn/dσ'n) | 50000 | ||
| Minimum shear velocity, vmin | 10 μm/s | ||
| Open hole length, L | 400 m | Häring et al., 2008 | |
| Boundary stress | σxx | 120.0 MPa | Häring et al., 2008 |
| σyy | 75.0 MPa | Häring et al., 2008 | |
| σxy | 11.0 MPa | Häring et al., 2008 | |
| Initial pore pressure | p0 | 46.0 MPa | Häring et al., 2008 |
3.2 결과 및 논의
3.2.1 결과
Case 1 해석 결과에 따르면 유체 주입 중단 후에 발생한 지진의 규모는 1.44로 주입 중에 발생한 지진의 규모인 0.86보다 0.58만큼 더 컸다(Fig. 6(a)). 최대 규모 지진이 발생한 위치는 교차점으로부터 약 900 미터 떨어진 지점으로 Fig. 6(d) 빨간 상자로 표기되어 있다. 이 경우에 상대적으로 큰 규모의 지진이 발생할 수 있었던 이유는 단층에 임계 공극압의 기울기가 공극압의 기울기와 비슷한 상태에서 공극압이 임계 공극압에 도달하여 여러 노드에서 한 번에 전단이 발생할 수 있었기 때문이다. 반면 유체 주입을 유지한 경우 큰 지진 대신 작은 지진이 여러 번 발생하였다. Fig. 6(c)와 Fig. 6(e)에서 임계 공극압이 주변보다 2~3 MPa정도 급격히 감소한 것을 볼 수 있는데, 이는 한 노드에서 발생한 지진에 의해 그 옆 노드에 응력이 집중되어 나타난 현상이다. 이로 인해 응력이 집중된 영역에서만 그 다음 지진이 발생하게 되어 지진의 규모가 제한되었다. 응력 집중으로 인해 규모가 제한되는 효과를 극복하고 Case 1의 유체 주입을 중단한 경우에서 큰 규모의 지진이 발생할 수 있었던 것은 주입 지점을 지나는 단층에서 발생한 전단에 의해 단층의 응력이 재분배되었기 때문이다(Fig. B1). Case 1의 경우 응력 재분배로 인해 Fig. 6(b)와 같이 교차점에 가까운 영역에서 초기 임계 공극압 47.5MPa 이상으로 임계 공극압이 증가했다. 이로 인해 Fig. 6(c)와 같이 단층 내에 공극압이 초기 임계 공극압 이상이 될 수 있었다. 유체 주입 중단 이후에도 Fig. 6(d)의 0~500m 구간과 같이 초기 임계 공극압 이상의 공극압이 계속 유지될 수 있었던 것은 500m 이후에는 전단으로 인해 수리적 간극이 증가하지 않아서 저류율과 투수율이 작게 유지되고 그로 인해 전단이 일어난 영역에서 유출되는 유량이 적기 때문이다. 그러나 교차점에서 먼 구간(Fig. 6(d)의 1000m이후)에서는 주입 지점을 지나는 단층의 전단에 의한 응력 재분배의 효과가 작기 때문에 여전히 초기 임계 공극압을 유지하였다. 결과적으로 응력 재분배는 넓은 영역에서 한 번에 전단이 일어날 수 있도록 Fig. 6(c)에 임계 공극압의 분포와 공극압의 분포를 조성하였다.
교차한 단층이 주입지점으로부터 멀어지며 발달한 Case 2는 유체 주입 중단 여부에 무관하게 큰 규모의 지진이 관측되지 않았다(Fig. 7(a)). 이는 주입 지점을 지나는 단층의 전단에 의한 응력 재분배의 효과가 Case 1과 다르기 때문이다. Case 2에서 단층은 주입 지점을 지나는 단층과 비슷한 방향으로 놓여있기 때문에 응력 재분배의 효과는 단층의 임계 공극압을 감소시켜 불안정화시킨다. 그로 인해 공극압의 증가가 아직 주입 지점을 지나는 단층에서 주로 일어나고 있음에도 불구하고 Fig. 7(b)의 뾰족한 임계 공극압 분포에서 알 수 있듯이 전단이 이미 발생하였음을 알 수 있다. 전단에 의한 응력 집중으로 인해 뾰족하게 감소한 임계 공극압 이상으로 공극압이 증가하기 어렵게 되어 Case 2에서는 단층내의 공극압이 대부분의 구간에서 초기 임계 공극압 이하로 유지되고, Fig. 6(d)에서와 같이 큰 영역에서의 전단이 발생할 수 없다. 결국 응력 집중이 연속적으로 일어나는 효과만이 주요하게 영향을 미쳐서 지진의 규모가 제한되었다. 이를 통해 유체 주입 중단 후에 발생하는 지진 규모 평가에서 공극압의 압력 구배가 완만해지는 효과와 더불어 응력 재분배에 의한 효과도 고려되어야만 함을 알 수 있다.
3.2.2 수리적 해석과 수리역학 복합 해석의 비교
복합거동을 고려하지 않고 수리적 해석만을 수행할 경우, Case 1과 Case 2에서 각 단층면이 초기 상태에서 같은 임계 공극압을 가지며 수리적 특성은 큰 차이가 없기 때문에 두 경우에서 지진 양상은 유사해야 한다. 그러나 수리역학 복합 해석은 두 경우에서 지진 규모와 발생 양상이 크게 달랐으므로 수리역학 복합 해석의 필요성을 알 수 있다. 이는 수리적 해석과 달리 수리역학 복합 해석에서는 전단 팽창, 응력 재분배, 절리의 방향 등을 고려할 수 있었기 때문으로 생각된다. 본 연구에서 수행한 수리역학 복합 해석은 유체 주입을 급격히 중단함으로써 바젤의 사례와 같이 오히려 큰 규모의 지진이 발생할 경우가 있음을 보였지만 더불어 다양한 단층과 물리적 기작이 관여되는 복잡한 암반 거동을 고려해야 하는 필요성을 보여주기도 한다. 따라서 유체주입 시 발생한 유발지진에 대한 대응을 위해 급격한 중단, 점진적 중단, 점진적 배수, 전면배수 등 다양한 공정이 유발지진에 미치는 영향 등에 대한 연구가 필요하며, 특히 3차원 균열의 형상, 단층 물성의 불균질성, 암반의 투수성 등을 추가로 반영하여 해석할 필요가 있다.
4. 결 론
본 연구는 유체 주입 중단과 지진 규모의 관계를 분석하기 위해 불투수성 암반에 대하여 수리적 해석과 수리역학적 복합 수치해석을 실시하였다. 수리적 해석을 위해 수리 자극된 영역과 수리 자극이 되지 않은 영역을 나누어 확산 방정식을 적용하여 해석해를 구하였다. 이를 통해 유체 주입 중단이 공극압 구배를 완만하게 함과 동시에 수리 자극된 영역의 경계 부근에서 수 시간 동안 상당한 임계 공극압에 도달하면서 추가적인 전단을 일으킬 수 있음을 확인하였다. 수리역학 복합 해석은 2차원 개별요소법을 이용하여 불투수성 암반 내 절리 유동을 모사하는 수치해석을 수행하였다. 단순화된 균열모델을 대상으로 유체 주입 중단 후에 최대 규모의 지진이 발생하는 것을 모사하였다. 유체 주입을 유지하는 경우에는 큰 규모의 지진 대신 여러 번의 작은 규모의 지진으로 나뉘어서 발생하였다. 이를 통해 급작스러운 유체 주입 중단에 따른 지진의 발생양상 및 규모에 미치는 영향은 보다 발전된 모델링으로 검토할 필요가 있음을 알 수 있다. 또한 공극압 구배와 더불어 단층 간의 기하학적 관계와 응력 재분배 등에 의해 임계 공극압의 구배가 변하는 것 또한 유발 지진의 규모 평가에 중요하게 고려되어야 한다.
Appendix A
수리적 해석해를 구하기 위한 미분방정식과 경계조건은 식 (1)과 같다.
a는 절리 간극, Kn은 절리의 수직 강성, μ는 유체 점성, p cri는 임계 공극압, Q는 주입 유량, r w는 시추공 반경을 나타낸다. 0~R 영역을 1번 영역으로 수리 자극된 구역(stimulated zone), 나머지 영역을 2번 영역으로 수리 자극되지 않은 구역(unstimulated zone)으로 나누었다. 시간에 대한 라플라스 변환을 적용하면 압력분포는 식 (2)와 같이 얻어진다. I, K는 수정된 베셀 함수(modified Bessel function)을 나타낸다.
유체 주입에 따른 해는 식 (2)에 라플라스 역변환을 적용하면 구할 수 있고, 유체 주입 중단은 t shut을 유체 주입 중단 시점이라 할 때 식 (3)와 같이 해를 중첩하여 구할 수 있다. 본 연구에서는 라플라스 역변환을 위해 수치적인 방법으로 Talbot’s method (Abate and Whitt, 2006, Verruijt, 2013)을 사용하였다.
수리 자극된 구역이 시간에 따라 확장되는 것을 모사하기 위해 주입이 계속됨에 따라 r = R 지점에서 공극압이 임계 공극압의 상한 값에 도달했을 때 (Find t i in Eq. (4)), 임계 공극압의 하한 값을 가지는 지점을 새로운 R (Find R i in Eq. (4))로 정하였다(Fig. A1). 이를 위한 알고리즘은 아래와 같고 확장되는 순간 R의 증가분을 그 시점에서 전단된 영역으로 가정하였다.
