Aims & Scope

정확하고 효율적인 수치 모델링을 수행하기 위하여 기존의 열-수리학적 모델인 TOUGH2 (Pruess et al., 1999)를 기반으로 지열 생산정 내부 및 지층에서 발생하는 복잡한 지하수 유동 및 열 이동의 상호 작용을 적용하기 위한 기법을 개발하였다. TOUGH2는 미국 LBNL에서 개발한 다차원 다상 열-수리학적 수치 모델로 지하수, 환경, 석유공학 등의 분야에서 널리 사용되고 있다. 이 모델은 Integral Finite Difference Method (Narasimhan and Witherspoon, 1976)에 기반하여 모델링 영역을 이산화하기 때문에 불규칙하고 불균일한 격자 요소를 사용할 수 있는 장점이 있다. 또한 이 모델에서는 다양한 유체의 특성을 모듈(module)을 통하여 고려할 수 있다. 즉 온도와 압력의 변화에 따른 지중 유체의 밀도 및 동적점성계수 변화를 고려할 수 있다. 본 연구에서는 지하수의 액체상과 기체상의 유체 특성을 고려할 수 있는 EOS1 모듈을 사용하였다.

지열 생산정은 실제로 가장 외부의 그라우트와 케이싱, 유체 유동 통로, 전기수중펌프(electrical sub-mersible pump, ESP) 및 내부 관로 그리고 내부 유체 통로등으로 구성되고 지층에 비해서 매우 좁고 길며, 심도에 따라서 그 설치 사양이 상이하다. 이러한 복잡한 생산정과 지층의 연결을 위하여 이들 구조를 단순화하고 등가(equivalent)의 격자 요소로 이산화하는 것이 필요하다. 본 연구에서 사용된 TOUGH2는 불규칙하고 불균일한 격자 요소를 실제가 아닌 가상의 연결을 통하여 구성할 수 있기 때문에 이러한 복잡한 격자 요소를 고려할 수 있다.

지열 생산정과 이를 통한 지하수 유동, 열 이동 현상을 수치 모델링 하는데 있어서 다음의 난관이 존재 한다: ① 그라우트, 케이싱, 파이프는 내부적으로 유체 유동은 발생하지 않고 열 이동 중에서도 전도만 일어나야 하고 내/외부의 지층 또는 유체 통로(지하수로 채워진)와도 마찬가지로 열전도만이 일어나야함, ② 지상부의 발전기 작동을 위한 최소한의 압력을 충족하기 위하여 지하 깊은 곳에 설치된 전기수중펌프에서 적절하게 압력의 상승(boost-up)이 일어나야 하는데 중간에 압력이 단속적으로 증가하는 조건을 고려할 수 있는 수치 모델 및 기법이 매우 제한적임. 첫 번째 문제와 관련하여 TOUGH2는 각 격자 요소들의 경계를 통하여 질량 이동(유체 유동) 및 열 이동(열 이송 및 전도)이 반드시 동시에 일어난다. 선행 연구를 통하여 각 격자 요소들의 경계 사이에서 열전도만이 고려되는 모듈을 개발한 바 있으며(Park et al., 2016), 이를 활용하여 해결하였다. 두 번째 문제의 경우 지열 생산정의 내부 유체 유동 통로를 전기수중펌프를 기준으로 두 구간으로 나누고 두 구간의 연계가 수치학적으로 동시에 되도록 TOUGH2의 source code를 수정/개발하여 해결하였다. 먼저 지층과 지열 생산정이 만나는 가장 하부로부터 전기수중펌프까지의 구간(아래 구간)에서는 전기수중펌프로 유입되는 단위 시간당 물의 질량을 고정하고 하부의 스크린구간으로부터 유입된 지열수가 지열 생산정을 따라서 올라오면서 주변 지층과 열교환(열전도)을 하는 조건을 수치 모델링 하였다. 전기수중펌프에서 지상의 지열수 유입구까지의 구간(위 구간)에서는 지상의 지열수 유입구에서 압력(발전기 운영을 위한 최소 압력)을 고정하고 전기수중펌프 위치에서 하부 구간을 통하여 펌프에 유입된 지하수 온도를 유지한 채로 같은 양을 다시 주입하는 경계 조건을 설정하여 수치 모델링 하였다. 이 경우에도 지열수는 상부로 이동하면서 주변의 구조물 및 지층과 열교환(열전도)을 하게 된다. 여기에서 아래 구간에서 계산된 펌프로 유입된 온도가 동시에 위의 구간의 전기수중펌프 위치에 주입되도록 적용되어야 하므로 이를 수치학적으로 선형연립방정식을 구성하는 단계에서 자동으로 반영되도록 source code를 수정/개발하였다. 따라서 이러한 기법을 적용 하면 지상의 지열수 유입구에서 압력과 생산 유량 두 가지가 주요 조건 변수가 되고 이들 변수에 따른 생산정 내에서의 압력과 온도가 연속적으로 계산되게 된다.

수치 모델링 대상 지층 및 지열정의 설정은 독일의 그로스 쉐네벡(Groß Schönebeck)의 인공저류층 지열 시스템 실증 현장의 조건을 참조하여 결정하였다(Zimmermann et al., 2010). 해당 실증 현장은 비화산지대에서 수리 자극을 통하여 지열 시스템을 시도한 곳으로 관련 연구 자료들이 많이 공개되어 있어 본 연구에 활용하였다. 현재 지진 안전성 문제로 중단된 국내의 포항 지열 발전 프로젝트와 조건들이 유사한 측면이 있으며, 실제로 연구 개발과 사업 수행에 큰 영향을 준 현장이다. 현재 주입정과 생산정 두 개의 심부 지열정이 완공되어 있으며, 지속적인 수리 자극과 실험 및 연구를 수행 중에 있다(Blöcher et al., 2010, 2015, 2018; Zimmermann and Reinicke, 2010; Jacquey et al., 2016, 2018). 해당 지열정이 시공 되어 있는 깊이까지의 지질학적 조건은 하부로부터 석탄기의 화산암, 하부 페름기의 화산암 및 퇴적암, 역암 등을 포함하는 퇴적암, 상부 페름기의 증발암류, 하부 트라이아스기의 퇴적암, 중부 트라이아스기의 증발암을 포함하는 퇴적암, 상부 트라이아스기의 석회암을 포함하는 퇴적암, 쥐라기의 미고결 퇴적암, 백악기의 미고결 퇴적암, 제3기의 미고결 퇴적암 및 현생의 퇴적물로 구성된다. 전반적인 지층들은 층상으로 분포하며, 암염돔 구조 및 구조 운동에 의한 완만한 습곡 등이 분포한다(Moeck et al, 2005; Blöcher et al., 2010). 그 중에서 하부 페름기의 화산암 및 퇴적암이 분포하는 -4,354 m ~ -3,878 m가 주요 지열 저류층이다(Reinicke et al., 2005; Zimmermann et al., 2010).

수치 모델링 영역은 대상 지층에 대하여 반경 500 m (x = 0 ~ 500 m), 높이 5,000 m (z = 0 ~ -5,000 m)의 이차원 방사상 축대칭 격자로 설정하였다(그림 1) 격자 요소는 수평 방향으로 20개, 수직 방향으로는 451개가 고려되었다. 수평 방향으로는 지열정에 가까울수록 격자 간격이 줄어드는 방식으로 구성하였다. 수직 방향으로는 지열정이 설치된 영역에서 10 m 간격으로 세밀하게, 하위의 지층에서는 조밀하게 구성하였다. 대상 지층의 물성은 기존의 문헌(Freeze and Cherry, 1979; Domenico and Schwartz, 1990; Fetter, 1994) 및 해당 지역의 연구 결과(Reinicke et al., 2005; Blöcher et al., 2010)를 참조하여 결정하였으며, 특히 열전도도의 경우에는 측정된 지온경사로부터 구간별로 계산된 결과(Ollinger et al., 2010)를 이용하였고 표 1에 정리되어 있다. 한편 지열 생산정은 하부에서 경사를 가지고 기울어져 설치되어 있으나 본 연구에서는 수치 모델링의 안정성과 효율성을 확보하기 위하여 지층의 경우 축 대칭 격자로 구성하고 수직인 지열정이 각 심도에서 만나는 방식으로 구성하였다. 수치 모델링 영역의 최상부와 최하부는 각각 지하수와 열 모두 통과하지 않는 경계조건을 부여하였고 측면 방향의 끝에는 (x = 500 m) 지하수 압력과 지온이 모두 초기 조건에 고정된 경계를 부여하였다. 수치 모델링은 총 10일의 기간에 대하여 수행하였다.

Fig. 1. 
Schematic diagram of the geologic formation (modeling domain), geothermal production well, and grid elements used in the numerical simulations. The horizontal coordinate axis x is exaggerated 10 times.

앞에서 언급한 바와 같이 수치 모델인 TOUGH2의 특성을 활용하여 지열정의 복잡한 형태 및 구조를 등가의 격자 조합으로 구성하였다(그림 2). 외부의 그라우트(grout)와 케이싱(outer pipe), 케이싱 내의 유체 유동 통로(inside of outer pipe), 내부 관로(inner pipe), 내부 유체 유동 통로(inside of inner pipe)를 구성 요소로 구분하고 각각의 심도에서 외경과 내경을 이용하여 등가 부피 및 접촉 면적과 연결 관계를 구성하였다(표 2). 이들은 수직적으로는 일차원의 수직 격자들의 연결로 구성되며, 수평적으로 각 요소들끼리 연결된다(그림 2). 각 요소들과 지층의 연결 관계에서 모두 수직, 수평적으로는 열전도가 고려되었고 유체가 실제로 유동하는 케이싱 내의 유체 유동 통로, 내부 유체 유동 통로는 수직적으로 지하수 유동과 열 이동의 전도 및 대류가 모두 고려되었다. 최하부의 스크린 구간에서는 지층과 케이싱 내의 유체 유동 통로에서는 지하수 유동과 열 이동이 모두 일어나며, 전기수중펌프 구간에서는 앞에서 기술한 바와 같이 펌프를 중심으로 상하부를 분리하고 펌프 하부까지 도달한 지하수가 온도를 유지한 채로 펌프 상부에 다시 공급되는 형태의 수치 모델링 조건을 적용하였다. 이때 지상의 지열수 유입구에서 압력(운영 압력)은 1.52 MPa (15 atm)으로 고정하고 펌프 주변에서의 유출 유입량, 즉 발전 생산 유량은 20 kg/sec로 설정하였다.

Fig. 2. 
Schematic diagram of the geothermal production well and its discretization and assemblage for IFDM (Integral Finite Difference Method) grids.

지열 생산정의 구성 요소인 그라우트, 케이싱, 내부 관로는 심도에 따라서 다른 직경을 가지며 이들의 열적 물성은 표 2에 정리되어 있다(Lindeburg, 2013). 또한 케이싱 내의 유체 유동 통로와 내부 유체 유동 통로의 경우 수치학적인 불안정성을 최소화하기 위하여 유체의 유동이 층류로 발생하여 Darcy의 법칙을 따른다는 가정을 하여 Darcy-Weisbach 방정식에 의하여 산정된 유효고유투수계수(effective intrinsic permeability)를 적용하였다(Olson and Wright, 1990; Sutera and Skalak, 1993; Munson et al., 1998).

발전을 위한 지열수의 생산이 시작됨에 따라서 스크린 구간과 인접한 지층 내에서는 압력의 하락이 발생하고 압력이 감소하는 영역이 시간이 지남에 따라서 점차적으로 지열정에서 멀어지는 방향으로 확산된다(그림 3). 생산 시작 후 1분 이내에는 지열정 중심의 20 m 내외의 영역에서 압력의 감소가 발생하고(그림 3b) 1시간 후에는 해당 영역이 약 200 m까지 확장되며(그림 3c), 10일 이후에는 영역의 경계에 해당하는 500 m 부근까지 도달한다(그림 3d). 지층 내 온도 분포의 초기 상태에서 등온도선의 수직 방향 간격이 심도에 따라서 구간마다 차이를 보여준다(그림 4a). 이는 앞에서 설명한 바와 같이 그로스 쉐네벡의 현장에서 측정된 지온 경사로부터 계산된 각 지층의 열전도도를 물성값으로 설정하였기 때문이다. 지층 내의 압력이 주로 스크린 구간을 중심으로 지열 저류층에서 변화가 발생한 것과는 달리 지층 내 온도의 분포는 주입정을 따라서 전체적으로 변화한다(그림3과 그림 4 참조). 생산 시작 후 6시간이 경과하면 생산정 주변의 지층의 온도가 소폭 증가하기 시작하며(그림 4b), 1일 이후에는 온도의 상승이 더 커지고 그 영역이 지열정에서 멀어지는 방향으로 확장된다(그림 4c). 10일 이후에는 지열정 바로 주변에서 온도가 현저하게 상승한 영역이 나타나고 이 영역이 지열정으로부터 수평으로 약 5 m 정도 멀어진 영역까지 도달한다(그림 4d). 또한 지열정 주변의 온도 상승이 점차적으로 둔화되고 10일 이후에는 변화가 매우 느리게 나타나게 된다. 이것은 스크린구간을 통하여 지열정으로 유입된 지하수가 상부로 상승하면서 주변의 지층과 열전도에 의한 열교환이 발생하기 때문이다. 또한 시간이 지남에 따라서 주변 지층의 온도가 증가하고 해당 영역이 수평으로 확장하면서 수평 방향의 온도 경사가 점차적으로 완화되어 전도에 의한 열교환 속도가 급격하게 감소하기 때문이다. 전반적으로 지층 내의 온도는 전도에 의한 열교환으로 수직적으로는 지열정과 만나는 모든 영역에서 변화한다. 반면 그 변화가 압력에 비해 상대적으로 더 느리고 수평 방향으로 훨씬 더 좁은 지역까지만 확장된다.

Fig. 3. 
Spatial distribution of (a) pressure at the initial state, and difference in pressure from the initial value after (b) 1 minute, (c) 1 hour, and (d) 10 days in the geologic formation.

Fig. 4. 
Spatial distribution of (a) temperature at the initial state, and difference in temperature from the initial value after (b) 6 hours, (c) 1 day, and (d) after 10 days in the geologic formation.

지열 생산정 내부에서 압력과 온도의 수직적 분포 변화는 그림 5에 도시되어 있다. 압력의 경우 지열수 생산 시작 후 비교적 빠른 시간 내에 전 영역에서 급격한 변화를 보이며, 약 1분 후까지 최대로 감소하였다가 점차적으로 소폭 회복되어 10일 이후에는 거의 정상 상태에 도달하게 된다(그림 5a). 압력의 회복은 지열수의 생산이 진행되면서 주변 지층과의 열교환이 지속적으로 이루어지고 생산정 내부가 비교적 높은 온도의 지하수로 채워지게 되면서 감소된 밀도에 의해 발생하는 것으로 판단된다. 압력의 수직적 분포에서 구간에 따라서 그 기울기가 달라지는 현상이 나타나는데 이것은 지열수가 통과하는 영역의 직경 변화로 단면적과 유효고유투수계수가 동시에 변화하기 때문이다. 지상의 발전 시설에 필요한 압력을 제공하기 위하여 1,100 m 부근에서 전기수중펌프의 작동에 의한 압력의 급격한 증가가 발생하며, 해당 압력의 증가가 앞에서 설명한 수치 모델링 기법의 적용으로 정상적으로 계산됨을 확인 할 수 있다.

Fig. 5. 
Vertical profiles of (a) pressure and (b) temperature of the geothermal fluid in the production well.

온도의 경우 스크린 구간의 상대적으로 뜨거운 지하수(스크린 구간 내에서 약 141℃ ~ 149℃)가 시간이 지남에 따라서 상부로 이동하면서 전체적으로 증가한다(그림 5b). 초기에는 주변 지층과의 온도차가 커서 전도에 의한 열손실로 상부에서 빠르게 온도가 상승하지 못하나 주변의 지층과의 열교환이 진행되면서 전도에 의한 열손실은 점차적으로 감소함에 따라서 지열정 내의 온도가 점차적으로 상승한다. 10분 후의 수직적 온도 분포가 직선이 되지 않고 중간에 변곡점이 나타나게 되는데 이것은 주변 지층의 열적 특성(열전도도)의 차이에 따른 열손실량의 차이 때문이다(그림 5b). 생산정 내 온도의 상승은 점차적으로 느려지며, 1일 이후에는 10일 이후와 거의 다르지 않은 수직적 분포를 보여준다. 10일 이후에는 지층 내에서 수평 방향의 온도 경사가 거의 완화되어 주변 지층으로의 열손실 속도가 매우 작아지게 되며, 특징적으로 스크린구간 최상부의 온도가 140.9 ℃인 반면 생산정 내에서의 온도는 상부로 갈수록 소폭 상승하여 지상부에서 최종적으로 145.7℃에 도달한다. 주변 지층과의 열교환에 의한 열손실이 비록 작은 값이지만 지속적으로 이루어지고 있음에도 온도가 상부로 갈수록 소폭 증가하는 현상은 물의 열역학적 특성으로 설명될 수 있다.

물의 압력-온도 조건에 따른 등엔탈피(iso-enthalpy)곡선(International Formulation Committee, 1967)과 생산정 내에서의 압력-온도의 변화 양상이 그림 6에 도시되어 있다. 생산정 내에서 압력-온도의 변화가 등엔탈피 곡선에 평행하게 진행한다면 열교환에 의한 열손실/이득이 없는 상태의 의미한다. 압력-온도 변화가 등엔탈피값이 더 낮은 곡선으로 이동하면 열손실, 더 높은 곡선으로 이동하면 열이득 현상이 일어난다. 지열수 생산 시작 후 10분 경과한 경우나 30분 경과한 경우 압력-온도의 변화 곡선이 등엔탈피 곡선을 가로질러서 급격하게 낮은 쪽으로 향하게 되는데 이 상태는 주변 지층으로의 열교환에 의한 열손실이 크게 발생함을 의미한다. 반면에 최종적으로 10일이 지나면 압력-온도의 변화 양상은 등엔탈피 곡선과 비교하여 거의 평행하지만 소폭 더 낮은 방향으로 이동하게 되는데 이것은 주변 지층으로의 열교환에 의한 열손실이 있으나 그 양이 매우 작음을 의미한다. 이 상태에서는 압력이 감소함에 따라서 온도는 소폭 상승하는 양상을 보여주고 있으며, 앞에서 언급한 온도 상승의 원인이 이에 기인한다.

Fig. 6. 
Pressure and temperature transitions in the production well with the iso-enthalpy curves of water.

인공저류층 지열 시스템에서는 지열 저류층의 생산성을 확보하기 위하여 수리 자극 등의 방법으로 지층의 고유투수계수의 증가를 시도하게 된다. 이때 효율적인 계획 수립 및 설계를 위해서는 충분한 지열수 생산이 가능한 고유투수계수를 사전에 설계하는 과정이 필요하다. 기존의 방법에서는 이러한 설계 과정에서 일반적으로 석유 생산 등의 분야에서 활용되는 생산성지수와 주입정 내에서의 압력 강하(well loss)를 모두 활용하여 생산 지하수량에 대한 압력의 강하량을 산정하게 된다(Zimmermann et al., 2010). 이 과정에서 물의 밀도 등의 물성이 하나의 상수로 가정하게 되는데 앞 절에서 설명된 바와 같이 실제 지열 생산정 내에서의 압력과 온도는 큰 범위 내에서 변화하게 되고 따라서 기존의 방법을 사용 시에 부정확한 결과를 초래할 수 있다.

본 연구에서 제시된 수치 모델링 기법을 활용하게 되면 지상의 지열수 유입구에서 압력과 생산 유량에 대응하는 생산정 내에서의 지하수 압력과 온도의 분포를 예측할 수 있다. 따라서 주어진 조건하에서 지열 저류층의 고유투수계수를 변화시키면서 반복적인 수치 모델링을 수행하면 해당 지열수 생산 조건을 만족할 수 있는 최소한의 고유투수계수를 정밀하게 산정할 수 있다. 또한 일반적인 지하수 부존량 등을 계산할 때와는 다르게 전기수중펌프에서 압력의 경계값을 대기압이 아니라 해당 온도에서의 포화수증기압(saturated vapor pressure)으로 설정하여야 한다. 이는 지하수의 압력이 해당 온도에서의 포화수증기압 이하가 되면 지하수가 기체상으로 변화하여 전기수중펌프가 정상적으로 작동할 수 없기 때문이다. 스크린 구간과 접하는 지열 저류층의 평균 고유투수계수를 변화시키면서 시행착오(trial-and-error) 수치 모델링을 수행하였고 그 결과가 그림 7에 도시되어 있다. 현재 주어진 지열수 생산 조건에서 전기수중펌프에 유입되는 지하수 온도는 평균 145.4 ℃이며, 이때의 포화수증기압은 4.172 × 10⁵ Pa이다(International Formulation Committee, 1967). 해당 압력보다 높은 압력을 유지하는 최소한의 조건은 지열 저류층의 평균 고유투수계수가 2.26 × 10^-15 m²일 때로 분석되었다. 다만 본 연구에서 시도된 수치 모델링 기법과 방법을 적용하여 인공저류층 지열 시스템의 목표 고유투수계수를 산정할 때에는 그림 7에 나타난 바와 같이 해당 구간에서 평균 고유투수계수 변화에 압력의 변화가 매우 민감하게 반응하므로 주의를 기울여야 한다.

Fig. 7. 
Average intrinsic permeability of the target depth of the geologic formation with respect to the pressure at the inlet of the electrical submersible pump.