Techniques to Estimate Permeability Based on Spectral Induced Polarization Survey

광대역유도분극 탐사에 기초한 유체투과도 예측기법들

Abstract

Permeability-analyzing methods commonly involve small-scale drilling, such as pumping or slug test, but it is difficult to identify overall distribution of permeability of the entire target sites due to high cost and time requirement. Spectral induced polarization (SIP) method is known to be capable of providing distributions of both the porosity and the pore size, the two major parameters determining permeability of the porous medium. The relationship between SIP variables and permeability has been studied to identify the hydrological characteristics of target sites. Kozeny-Carman formula has been improved through many experiments to better predict fluid permeability with electrical properties. In this work, the permeability prediction techniques based on SIP data were presented in accordance with the hydrogeological and electrical characteristics of a porous medium. Following the summary of the techniques, various models and related laboratory experiments were analyzed and examined. In addition, the field applicability of the prediction model was evaluated by field case analysis.

1. 서론

유체투과도(permeability)는 다공성 매질에서 상호 연결된 공극을 통해 유체가 얼마나 잘 흐를 수 있는지에 대한 척도로 저류층 생산, 수자원 관리 등을 계획하는 데 필수적이다. 환경 분야에서는 오염 물질의 유동 범위를 관리하고 정화 전략을 효과적으로 세우기 위해서 유체투과도를 정확히 분석하는 것이 매우 중요하다(Robinson et al., 2018). 현장에서 수리지질학적 특성을 규명하는 것은 가장 어려운 문제 중 하나이다. 매질의 수리지질학적 특성을 파악하기 위해 사용되는 직접적인 방법으로는 양수 시험(pumping test)나 순간수위변화시험(slug test) 같은시추공 시험 자료를 이용하는 것이다. 이러한 측정 방법은 시추공이 위치한 범위 내에서는 비교적 정확한 정보를 제공해줄 수 있으나 실제 매질은 불균질하며 복잡하기 때문에 시추공 주변의 제한적인 자료만으로 현장 스케일에서 수리지질학적 특성에 대하여 파악하기 어렵다. 지구물리탐사법은 공간적으로 비교적 조밀하게 측정할 수 있으며 이에 따라 연속적인 공간적 특성을 구할 수 있다. 매질의 전기적 특성을 이용하여 유체투과도를 예측하는 방법에 대한 가능성에 대한 연구가 꾸준히 이루어져 왔다.

IP 탐사는 전기전자 현상에 기초한 물리탐사 기법 중 하나로 1970~1980년대에 광물 탐사 분야에서 활발이 이용되었으나 2000년대 이후 장비와 해석법의 발전으로 인해 더욱 다양한 분야에 적용되고 있다(부록 A ; Kim et al., 2017). 특히 진동수영역에서 IP 반응을 측정하는 복소전기비저항(complex resistivity; CR) 혹은 광대역 유도분극 탐사(spectral IP; SIP)법은 유체투과도에 큰 영향을 끼치는 수리적 변수와 관계성이 밝혀지면서 큰 주목을 받게 되었다. 공극 내 유체의 특성에 대해서만 주로 분석할 수 있는 전기비저항과 달리 IP 반응은 공극률, 내부 면적과 같은 공극의 기하학적 특성(pore space geometry)과 유체의 염분 등으로 결정되는 공극 내 유체의 물성에 모두 영향을 받기 때문에 실질적인 유체투과도를 결정하는 수리지질학적 변수와의 관계성을 더 잘 설명할 수 있다(Leroy et al., 2017, 2008; Revil et al., 2012; Weller et al., 2015). 이러한 관계성을 기반으로 수리지질학적 변수를 이용한 유체투과도 예측 모델을 허수 전기전도도나 시간 상수 등의 IP 관련 변수로 대체하기 위해 많은 연구가 이루어졌다(Kruschwitz et al., 2016; Weller et al., 2010; Slater et al. 2006; Binley et al. 2005; Börner et al., 1996).

이 논문에서는 먼저 다공성 매질의 수리적 특성과 전기적 특성에 대하여 기술하고 수리지질학적 변수를 이용한 유체투과도 모델 중 가장 대표적인 Kozeny-Carman 모델식에 대하여 정리하였다. 전기적 특성 변수, 특히 IP 변수와 수리적 특성의 관계를 설명하기 위해 실수 전기전도도, 허수 전기전도도, IP 변수와 수리지질학적 변수와의 관계를 기술하였다. 이러한 관계성을 입증하기 위해 이루어진 기존의 실내실험 사례에 대하여 정리하였으며 마지막으로 실제 현장 스케일의 유체투과도 모델로 적용성을 알아보기 위해 현장 사례를 소개하였다.

2. 다공성 매질의 특성

2.1. 다공성 매질 수리적 특성

수리전도도(hydraulic conductivity; K [m/s])는 다공성 매질(porous medium) 공극 내에서 지하수의 흐름이 얼마나 용이한지 나타내는 척도이며, 특히 포화 수리전도도(saturated hydraulic conductivity)는 물에 의해 완전히 포화된 매질에서의 수리전도도를 의미한다. 수리전도도에 주로 영향을 주는 요인으로는, 다공성 매질의 특성인 투수율 혹은 유체투과도(permeability; k [mD])와 공극 유체 특성인 유체 밀도(density; \(\rho_{w}\)w [kg/m³]), 유체 점성도(viscosity; \(\mu\) [kg/m·s]), 중력가속도(gravity; g [m/s²)가 있다. 포화 매질의 수리전도도(K)는 이들의 함수로 다음과 같이 나타낼 수 있다.

\(K=k \frac{\rho_{w} g}{\mu}\)       (1)

불포화대에서 수리전도도는 위의 정의에서 지하수의 함수율(water content)를 추가적으로 고려해야 한다. 유체투과도(k)는 암석 내를 유체가 얼마나 쉽게 통과할 수 있는 지를 나타내는 공극 매질 고유의 값이며 공극의 크기, 모양, 분포 상태에 좌우된다. 단위로는 주로 millidarcy를 사용한다. 공극성 매질을 가장 간단하게 나타내는 모델로는 모세관 다발 모델(capillary bundle model)이 있는데 이 모델에서는 모세관을 하나의 유체 경로로 가정하고 한 공극 매질에서 모세관 현상이 일정하게 일어나는 것을 전제로 하여 모세관 망(network of capillaries)을 통한 흐름을 나타낸다(Fig. 1).

Fig. 1. Pore structure models; (a) Porous medium and (b) capillary bundle model.

모세관 다발 모델을 기반으로 한 유체투과도 예측모델 중 가장 많이 사용되는 식은 Kozeny-Carman이 제시한 등가수로방정식(equivalent channel equation)으로(Kozeny, 1927; Carman, 1939) 흔히 Kozeny-Carman 방정식이라 불린다.

\(k=\frac{\phi r^{2}}{a T}\)       (2)

여기서 \(\phi\)는 공극률(porosity), T는 공극 모세관 비틀림계수(pore capillary tortuosity, T = [L_e / L]²), a는 형상계수(shape factor)이며, r [m]은 공극의 수리반지름 혹은 수경(hydraulic radius)으로 단면적 당 공극 부피를 의미한다(이때 유효 통로 길이(L_e)가 샘플 길이(L)보다 값이 커야한다). 수리반지름은 측정 가능한 단위 공극 부피당 표면적(surface area per unit pore volume; S_por [\(\mathrm{m}^{-1}\)])의 역수 관계로 흔히 나타낸다.

\(S_{\text {por }}=\frac{2}{r}=S_{s} \rho_{s} \frac{(1-\phi)}{\phi}\)       (3)

여기서 S_s는 비표면적(specific surface area [m²/kg])이며\(\rho_{s}\)는 광물 밀도이다. 수리반지름(r)과 S_por의 역수 관계를 이용하여 Kozeny-Carman 방정식(식 2)을 다시 정리하면 다음과 같다.

\(k=\frac{\phi}{a S_{\text {por }}^{2} T}\)       (4)

모세관 다발 모델이 모든 매질에서 적용 가능한 것도 아니고 오차 범위도 매우 크다는 한계점이 있지만, 전기적 특성 변수를 이용하여 유체투과도 예측 모델을 만드는데 Kozeny-Carman 식이 기초가 되었기 때문에 매우 중요하다. 이러한 모세관 다발 모델은 동일 직경 r을 가진 수많은 원통형 모세관을 통해 유체의 흐름이 발생한다고 가정한 반면 실제 공극성 매질의 공극의 크기와 연결성은 균일하지 않기 때문에 유체투과도와 지층계수(Formation factor; F)에 기초하여 다음과 같이 계산한 유효 수리반지름(effective hydraulic radius; r_eff [\(\mu\)m])을 이용하여 모델을 구성한다.

\(r_{e f f}=\sqrt{8 F k}\)       (5)

유효 수리반지름은 직접적으로 측정할 수 있는 방법이 없기 때문에 이를 대체할 수 이는 측정 가능한 변수를 상정하는데, 가장 흔히 사용하는 대체 변수는, 수은을 주입하여 모세관압을 측정하는 수은주입 모세관압력(mercury injection capillary pressure; MICP) 방법을 통해 얻은 주공극 크기(dominant pore size; r_dom)이다. 이와 같이 MICP를 이용한 주공극 크기 외에도 모세관뭉치 모델(capillary bundle model; CBM)로부터 다음과 같이 계산하는 공극크기(\(r_{S_{\text {por }}}\))도 있다.

\(r_{S_{\text {por }}}=\frac{2}{S_{\text {por }}}\)       (6)

실제 토양 혹은 암석 시료의 경우 Spor는 BET 공법(Brunauer, Emmett and Teller, 1938)과 같은 가스 흡착 방법으로 얻을 수 있다.

2.2. 다공성 매질의 복소 전기전도도 및 Cole-Cole 모델

지하 매질의 복소 전기전도도(complex conductivity, \(\sigma^{*}\))는 진동수(frequency; f )의 함수로 다음과 같이 나타낼 수 있다.

\(\sigma^{*}(\omega)=\sigma^{\prime}(\omega)+i \sigma^{\prime \prime}(\omega)\)       (7)

\(\sigma^{\prime}\)와 \(\sigma^{\prime \prime}\)는 각각 복소 전기전도도 \(\sigma\)의 실수 성분과 허수 성분이며 는\(\omega\) 각진동수(angular frequency; \(\omega\) =2\(\pi\)f )이다. 낮은 진동수(1000 Hz 이하)에서 다공성 매질의 복소 전기전도도는 1) 상호 연결 되어있는 공극(inter-connected pore) 내 유체를 통한 전기전도도(electrolytic conductivity, \(\sigma_{e l}[\mathrm{~S} / \mathrm{m}]\))와 2) 광물 표면의 복소 표면전도도(surface conductivity, \(\sigma_{\text {surf }}^{*}[\mathrm{~S} / \mathrm{m}]\))로 이루어진 평행회로로 가정하여 나타낼 수 있다(Vinegar and Waxman, 1984).

\(\sigma^{*}(\omega)=\sigma_{e l}+\sigma_{\text {surf }}^{*}(\omega)\)       (8)

실수 전기전도도(\(\sigma^{\prime}\))는 전하 운반체(charge carrier)에 의한 전류와 공극 벽 표면에 형성되어 있는 전기적 이중층(electrical double layer; EDL)의 이온에 의한 전류로, 일반적으로 평행회로로 가정하여 아래와 같이 나타낸다.

\(\sigma^{\prime}(\omega)=\frac{1}{F} \sigma_{w}+\sigma_{\text {surf }}^{\prime}(\omega)\)       (9)

여기서 \(\sigma_{w}\)는 공극 유체의 전기전도도, \(\sigma_{\text {surf }}^{\prime}(\omega)\)은 복소 표면전도도(\(\sigma_{\text {surf }}^{*}(\omega)\) 의 실수 성분이며 전기적 이중층에서 전하 이동을 의미한다. 지층계수 F는 Archie 식을 이용하여 다음과 같이 구할 수 있다.

\(F=\frac{\sigma_{w}}{\left(\sigma^{\prime}-\sigma^{\prime \prime}\right) / l}\)       (10)

여기서 l은 실험상수로 Börner et al.(1996)은 미고결 시료에 대해서 l = 0.1로 제시하였다.

Fig. 2. Permeability prediction from electrical properties (modified from Slater, 2007).

복소 전기전도도의 허수 성분은 전기적 이중층에서 발생하는 일시적인 전하 축적(분극)을 나타낸다.

\(\sigma^{\prime \prime}(\omega)=\sigma_{\text {surf }}^{\prime \prime}(\omega)\)       (11)

복소 표면 전도도의 실수 성분과 허수 성분 모두 진동수에 따라 변화하며 IP 효과의 세기는 허수 전기전도도에 의해 나타난다.

진동수에 따른 복소전기전도도의 분산 현상을 모사한 모델을 사용하는 데 가장 많이 사용하는 모델은 Cole 모델이다(Cole and Cole 1941; Pelton et al., 1978).

\(\sigma^{*}(\omega)=\sigma_{0}\left(1-\frac{m}{1+(i \omega \tau)^{c}}\right)\)      (12)

여기서 \(\sigma_{0}\)은 DC 전기전도도, m은 충전율(\(m=1- \sigma_{0}/\sigma_{\infty}\), 이때 \(\sigma_{\infty}\)은 높은 진동수에서 전기전도도), \(\tau\)은 진동수 분산에 따른 시간상수, c는 이완상수이며 미고결 혹은 퇴적암에서 0.1~0.6 값을 가진다.

3. IP 변수와 수리학적 변수 관계

3.1. 실수 전기전도도

유체투과도 k와 실수 전기전도도 \(\sigma^{\prime}\)의 관계성은 광물 및 공극 유체의 물리적 및 화학적 성질에 따라 변화하며 이는 전해질 전기전도도(\(\sigma_{e l}\)) 혹은 매질 알갱이 표면과 공극 유체 사이의 경계면 전기전도도(\(\sigma_{\text {surf }}^{*}\))가 \(\sigma^{\prime}\)에 지배적인 영향을 미치는 것임을 보여주는 것이다(Slater, 2007; Attwa and Günther, 2013). 전통적으로 실수 전기전도도(즉, DC 전기비저항)와 k 값의 관계성에 대한 많은 연구가 이루어지면서 k 값을 예측하는 데에 많이 이용되었다. 그 중 하나로 k와 \(\sigma^{\prime}\)의 양의 로그-로그 관계를 나타낸 아래와 같은 모델이 있다(Heigold et al., 1997; Purvance and Andricevic, 2000; Attwa et al., 2009; Attwa, 2012).

\(k_{D C}=a_{D C}\left(\sigma^{\prime}\right)^{b_{D C}}\)       (13)

위 식을 적용하기 위해서는 K 값과 전기적 측정 값을 비교하면서 실험상수인 \(a_{D C}\)와 \(b_{D C}\)값을 유추하는 과정이 필요하다. 하지만 \(\sigma^{\prime}\)값이 전해질 전기전도도(\(\sigma_{\text {surf }}^{*}\))에 의한 공극 부피 통로와 경계면 전기전도도( )에 대한 공극 표면 면적에 모두 민감한 반응을 보이기 때문에 전기적 반응 측정 자료로부터 매질의 특성을 분석하기 쉽지 않아, 실수 전기전도도를 이용해 K 값을 예측하는 데는 명확한 한계가 있다.

3.2. 허수 전기전도도

실수 전기전도도만을 이용할 경우 얻을 수 있는 정보의 한계점은 매질의 전기전도도와 충전효과 특성을 같이 측정하는 IP 탐사법을 통해 해결할 수 있다. 이는 IP 측정값으로 K 값을 예측할 수 있다는 것은 식(11)에서 보여주고 있듯이 허수 전기전도도 \(\sigma^{\prime \prime}\)가 상호 연결된 공극표면적에 대한 특성에 유일하게 관여된 변수라는 사실 때문이다. Börner et al.(1996)은 전기전도도 허수성분과 S_por의 선형관계(\(\sigma^{\prime \prime} \propto S_{\text {por }}\))를 제시하였으며(Fig. 3) 이를 기반으로 IP 변수로부터 수리전도도를 예측하고자 하는 많은 연구가 수행되고 있다(Börner et al., 1996; Slater and Lesmes, 2002; Slater et al., 2006; and Weller et al., 2010b). 이상과 같이 \(\sigma_{w}\)와 \(\sigma_{\text {surf }}^{\prime}\)를 알면 실수 전기전도도(\(\sigma^{\prime}\))로부터 F 값을 얻는 것이 이론적으로는 가능하지만 \(\sigma_{w}\)는 측정을 통해 쉽게 알 수 있는 반면, \(\sigma_{\text {surf }}^{\prime}\)은 직접 알아내기 어렵다.

Fig. 3. Relationship between σ'' and S_por (Weller et al., 2010).

3.3. IP 변수: 시간상수

시간상수 \(\tau\)은 공극 크기 혹은 입자 크기의 특성과 연관이 많으며 이는 전류 송신원을 주었을 때 발생하는 전기적 이중층에 있는 이온의 길이 척도(length scale)를 정의한다(Chelidze and Gueguen 1999; Olhoeft 1985). 다수의 연구에서 광물화된 암석에서 \(\tau\)길이척도가 평균 입자 사이즈 직경에 비례한다는 사실을 입증하였다(Pelton et al. 1978; Wong 1979; Wait 2012). 그러나 최근 연구에서는 \(\tau\)와 임계공극 직경(critical pore diameter) 혹은 공극 간 연결통로 크기에 상관관계가 있다는 것을 보여주고 있다(Fig. 4; Binley et al. 2005; Scott and Barker 2003, 2005).

Fig. 4. Time constant (\(\tau\)) from Cole–Cole relaxation modeling of SIP spectra versus dominant pore size (D₀) for UK sandstone samples from two independent studies (Slater, 2007).

4. 유체투과도 예측모델

유체투과도는 공극의 기하학적 특성에 기초하여 설명할 수 있으며, 이러한 공극의 기하학적 매개변수들은 SIP 변수에 기초하여 추정할 수 있다. 유체투과도를 추정하기 위한 기하학적 모델로는 Kozeny-Carman 모델을 기반으로 기하학적 길이 척도(Geometric length scale)를 이용한 다양한 모델들이 있다. 그 중 대표적인 모델로는 모세관 망을 통한 유체 흐름을 기반으로 하는 PaRiS 모델(Pape et al., 1987)과 Banavar & Johnson(1987)이 공극 크기 분포에 대한 확산 이론(percolation theory)을 기반으로 제시한 Katz와 Thompson(KT)모델(Katz & Thompson, 1987)이 있다. 두 모델에서는 유체 흐름을 조절하는 공극 크기를 나타내는 기하학적 매개변수인 유효 수리반지름 r_eff 을 이용하는데 이를 정의하는 방식에 약간의 차이가 있다. 기하학적 모델의 가장 큰 한계점은 가장 주요한 변수인 r_eff 가 실험실에서 측정해야 알 수 있다는 점이다(Weller and Slater, 2019). 지구물리 길이척도는 지표 탐사 자료나 시추공 자료로부터 얻은 변수를 이용하기 때문에 이러한 문제점을 해결해 줄 수 있다(Robinson et al., 2018). 매질의 전체 부피 특성(total volumetric properties)을 전기적 변수에 기초한 유효 특성(effective properties)으로 변환하는 것은 Kozeny-Carman 방정식(식(2))을 개선하는데에 매우 중요하게 적용되었다(Revil and Cathles, 1999). 즉, 기존 Kozeny-Carman 식에 있는 공극률을 공극 공간(pore space)의 부피 혹은 연결성과 관련된 전기적 특성 변수로 바꾸는 것이다. 그러나 수리반지름(r) 등의 수리학적 변수로 이루어진 Kozeny-Carman 식을 유효 특성 변수로 정확히 바꾸는 것은 간단하지 않다.

Revil과 Cathles(1999)은 전기적 이론으로부터 정의된 수리반지름에 대한 유효 변수일지라도 직접 측정하는 것은 어렵다고 보았으나 최근 연구에서는 이러한 변수에 대한 추가적인 정보를 복소 전기전도도(\(\sigma^{*}(\omega)\)) 자료를 분석하여 얻을 수 있다고 제안하였다. 즉, 단일 진동수를 이용한 IP 측정은 2개의 전기적 변수(해당 진동수에서 \(\sigma^{\prime}\)와 \(\sigma^{\prime \prime}\))를 얻는 반면, SIP 방법에서는 진동수에 따른 \(\sigma^{*}(\omega)\)의 양상을 분석할 수 있으므로 이를 통해 추가적인 전기적 특성 정보를 얻을 수 있다.

4.1. 실험모델(Empirical model)

4.1.1. PaRiS 모델

PaRiS 모델에서 유효 수리반지름(r_eff)은 다음과 같이 정의한다.

\(r_{e f f} \cong \frac{2}{S_{p o r}}\)       (14)

PaRiS 모델에서는 식 (2)의 Kozeny-Carman 방정식에서 공극 모세관 비틀림계수(T)와 공극의 비율을 지층계수(F)로 대체시켜 유체투과도를 다음과 같이 정의할 수 있다.

\(k_{\text {PaRiS }}=\frac{a_{\text {PaRiS }}}{F\left(S_{\text {por }}\right)^{c_{\text {PaRis }}}}\)       (15)

위 식에서 k_PaRiS는 PaRis 모델에 기초하여 예측한 유체투과도(k)임을 나타내며 a_PaRiS와 c_{PaRiS​​​​​​​}는  PaRiS 모델의 실험 상수이다. 결과적으로 이 모델의 기본적인 가정은 위상지연 측정 값이 진동수와 독립적이라는 것이다. 그러나 대부분의 암석은 그렇지 않기 때문에 Pape et al. (1987)은 PaRiS 모델을 이용하여 상수 위상각(constantphase-angle, CPA)을 이용해 kPaRiS 값을 계산했다. 식(14)은 고결된 퇴적 암석에 적합하며 Pape et al.(1987)이 a_PaRiS = 0.00475 값을 제안한 이후, Börner et al. (1996)이 식 (14)가 미고결 암석 시료에 적합하도록 수정하면서 a_PaRiS = 1로 대체하고 지수값인 c_PaRiS는 2.8 < c < 4.6 범위 안에 있다는 것을 보였다. a_PaRiS와 c_PaRiS에 가장 대표적인 상수 값을 적용하면 PaRiS 실험 모델은 다음과 같이 나타낼 수 있다.

\(k_{P a R i S}=\frac{475}{S_{p o r}{ }^{3.1} F}\)       (16)

4.1.2. 공극 부피당표면적(S_por) 예측: 허수 전기전도도 적용 모델

(15)식을 바탕으로 많은 연구에서 유체투과도를 예측하는 모델을 전기적 특성 변수를 이용할 수 있도록 개선시키고자 하였다(Börner et al. 1996; de Lima and Niwas 2000; Slater and Lesmes 2002;). 그 중 대표적인 것이 전기적 비틀림계수(electric tortuosity)가 수리비틀림계수(hydraulic tortuosity)와 같다는 전제 하에 \(T / \phi\)를 지층계수로 교체하는 것이다. Börner et al.(1996)은 넓은 진동수 대역에 걸쳐 위상 지연(phase shift, \(\varphi\))이 나타난 것을 관측했으며 그 중 1 Hz 측정 자료와 수리전도도 k_PaRiS의 관계식을 제시하였다. 허수 전기전도도(\(\sigma^{\prime \prime}\) ;[mS/m])는 공극 표면적에 대한 특성이 유일하게 영향을 미치기 때문에 식 (14)를 다음과 같이 나타낼 수 있다.

\(k_{\sigma^{\prime \prime}}=\frac{a_{\sigma^{\prime \prime}}}{F\left(\sigma^{\prime \prime}\right)^{c_{\sigma^{\prime \prime}}}}\)       (17)

여기서 \(a_{\sigma^{\prime \prime}}\)와 \(c_{\sigma^{\prime \prime}}\)는 허수 전기전도도 기반 모델의 실험 상수이다. Weller et al.(2015)은 이 관계식을 기반으로 다양한 지역의 순수 사암 혹은 셰일질 사암 샘플을 1 Hz 진동수에서 측정하여 다음과 같은 유체투과도 예측실험식을 제시했다.

\(k_{\sigma^{\prime \prime}}=\frac{2.66 \times 10^{-7}}{\sigma^{\prime \prime} 0.66 \mathrm{~F}^{5.35}}\)       (18)

여기서 \(k_{\sigma^{\prime \prime}}\) [m²]는 허수 전기전도도(\(\sigma^{\prime \prime}\))에 기초하여 예측한 유체투과도(k)이다. \(\sigma^{\prime \prime}\)과 S_por의 관계에 대한 선형계수로 고결시료 실험에서는 10 10^–12였으며 미고결시료실험에서는 7.5 10^–12로 나타났다.  Weller et al.(2015)은 식을 허수 전기전도도 대신 Debye decomposition(Nordsiek and Weller, 2008)의 정규화 충전율(m_n)을 이용한 유체투과도 예측식을 다음과 같이 제시한 바 있다.

Fig. 5. Permeability estimation model based on hydraulic parameters and geophysical parameters.​​​​​​​

\(k_{\sigma^{\prime \prime}}=\frac{8.69 \times 10^{-7}}{m_{n}^{0.79} F^{5.35}}\)       (19)

4.2. 역학 모델(Mechanical model)

4.2.1. KT 모델

Johnson et al.(1986)은 제한된 특정 구역의 공극 공간(e.g., 유체유동 결정공극(pore throat))에 기초한 가중 공극 표면적 대비 부피비(a weighted pore volume to pore surface ratio with the weight favoring constricted regions of the pore space(i.e. the pore throats)) 혹은 역동적 연결 공극 반지름(dynamically interconnected pore radius)를 의미하는 전기적 기하학적 변수인 \(\Lambda\)을 도입했다. Revil and Cathles(1999)은 \(\Lambda\)를 유체흐름 공극 표면적당부피라 하여 다음과 같이 정의하였다.

\(\Lambda=1 / S_{p o r}\)       (20)

\(\Lambda\)는 수은을 주입하는 MICP법을 통해 측정하며 S_por은 질소(N₂)와 같은 기체를 주입하여 Brunauer-EmmettTeller(BET)법을 이용해 표면에 흡착된 가스의 양으로 비표면적(S_s)을 측정하여 식(3)을 통해 얻을 수 있다. 

KT 모델에서는 유효 수리반지름(r_eff)이 와 \(\Lambda\)동일하고 본다(Banavar & Johnson, 1987; Katz & Thompson, 1986).

\(r_{e f f} \cong \Lambda\)       (21)

KT모델에서는 유체투과도를 다음과 같이 나타낼 수 있다.

\(k_{K T}=\frac{\Lambda^{2}}{8 F}\)       (22)

위 식에서는 KT 모델에 기초하여 예측한 유체투과도임을 명시하기 위하여 아래 첨자에 모델명을 적시하여 k_KT로 나타내었다.

4.2.2 유체흐름 공극 표면적당부피( ) 예측: 시간상수 적용 모델

Revil et al.(2012)은 KT 모델을 기반으로 Cole-Cole모델 변수 중 시간상수(relaxation time, \(\tau\) [sec])를 이용해 유체투과도를 예측하는 모델을 제시하였다. 이 때 시간상수 \(\tau\)는 입자 크기(R)의 제곱에 비례하며 관계식은 아래와 같다.

\(\tau=\frac{R^{2}}{2 D_{+}}\)       (23)

D_{+​​​​​​​} [m² /s]는 확산계수(Nernst-Einstein relationship)로 \(D_{+}=k_{b} T \beta_{+}^{s} /\left|q_{+}\right|\)에 의해 이온 이동도 \(\beta_{+}^{s}\), 온도 T, 볼츠만 상수(Boltmann’s constant) k_b, 흡착층에서 반대이온의 충전(charge of counterions in the Stern layer) q+와 관계되어 있다. 매질마다 고유의 이온 이동도가 있으며 Revil(2013)과 Revil 등(2012)은 실험을 통해 순수 사암에서 확산계수 D₊:1.3 10^–9 m² /s와 점토질 매질에서는 D₊:3.8 10^–12 m² /s을 제시했다. Revil et al.(2012)은 \(\tau\)이와 비례관계에 있음을 다음과 같은 식으로 제시하였으며

\(\tau=\frac{\Lambda^{2}}{2 D_{+}}\)       (24)

식 (24)를 식 (22)에 대입하면 아래와 같은 유체투과도 모델을 구할 수 있다.

\(k_{\tau}=\frac{\tau D_{+}}{4 F}\)       (25)

식 (25)을 이용해 유체투과도 예측을 위해서는 정확한 D₊을 알아야한다.

5. 실내실험 사례

Börner et al.(1996)은 최초로 r_eff와 지층계수를 결합시켜 지구물리학적 변수를 이용한 유체투과도 예측 모델을 고안했다. 이 연구에서는 \(\sigma^{\prime \prime}\)와 S_por의 직접적인 선형 관계(\(\sigma^{\prime \prime} \propto S_{\text {por }}\))를 제시하여 전기적 특성으로 유체투과도를 유추할 수 있다는 가능성을 보여주었고 이를 기반으로 IP자료를 이용한 유체투과도 예측 모델에 대한 다양한 연구가 진행되었다. 하지만 Börner et al.(1996)에서는 실험 자료가 충분하지 않았다는 한계점이 있는데 이후에 Weller and Slater(2019)는 \(\sigma^{\prime \prime}\)와 S_por가 단순한 선형관계가 아니라 지층계수와도 연관되어 있음을 보여주었다. 대부분의 실내실험에서 단일 진동수를 이용할 때는 주로 1 Hz를 이용하였는데 이는 1 Hz가 SIP 탐사에서 일반적으로 사용하는 진동수 대역의 중간 대역이기 때문에 기준 진동수(reference frequency)로 사용하기에 유리하다(Börner and Schön, 1991; Weller and Slater, 2019). 또한, 일정한 실험 조건이 유지되도록 포화 유체 염도는 1 Hz에서 100mS/m가 되도록 고정시켜 사용하였다. Robinson et al. (2018)과 Weller and Slater(2019)에서는 기존 여러 연구에 활용되었던 고결 및 미고결 시료를 이용해 샘플 개수를 더욱 확장하여 지구물리학적 길이척도를 기반으로 k값을 예측하는 방법이 기하학적 모델 혹은 유효수리반지름을 이용한 방법에 비해 갖는 장점에 대하여 기술하였다(Fig. 6).

Fig. 6. Permeability predictions using geometric length scales (a) PaRiS model per equation (4) (d = 1.023) and (b) KT model with K per equation (16) (d = 0.729) and permeability predicted for the polarization magnitude geophysical model using (c)\(\sigma^{\prime \prime}{ }_{1 \mathrm{~Hz}}\)(d = 0.770) and m_n (d) (d = 0.807) (Robinson et al., 2018).​​​​​​​

Börner et al.(1996)이 단순화한 \(\sigma^{\prime \prime}\)와 S_por의 선형 관계에는 전제조건이 필요한데 바로 1) 유체 구성성분이 일정하고 2) 지층 계수 변동이 매우 작다는 것이다. Weller et al.(2010b)은 9개의 독립적인 연구자료로 쓰였던 사암과 미고결 토양 시료로 구성된 총 114개의 자료를 취합하여 이러한 선형 관계가 합당하다는 것을 확인하였다. 그런데 Niu et al.(2016)은 밀도가 높은, 즉 상대적으로 F와 S_por 값이 큰 점토질 암석에서  \(\sigma^{\prime \prime}\)와 S_por의 비례 관계가 이루어지지 않음을 발견하였다. 이에 따라 \(\sigma^{\prime \prime}\)와 S_por의 단순한 선형 관계를 보강하기 위해 F를 고려한 \(\sigma^{\prime \prime} \propto S_{p o r} / F\)를 제시했으며 이 모델은 점토질 암석에서 잘 적용되었다. Weller and Slater(2019)에서는 이러한 연구 결과를 받아들이고 Kozeny-Carman 방정식에 대입하여 \(\sigma^{\prime \prime} \propto S_{p o r} / F\)를 고려한 새로운 지구물리 변수를 도입한 유체투과도 예측 모델을 제시하였고 기존 모델 보다 더 향상된 결과를 얻을 수 있었다(Fig. 7).

Fig. 7. Effective hydraulic radius r_eff computed from equation (2) versus geophysical length scales (a) r_{σ ''} and (b) r_{Fσ ''} (Weller and Slater, 2019).

Weller et al.(2015)은 공극율을 이용한 유체투과도 모델과 공극-정규화 표면적(pore-normalized surface area)을 이용한 유체투과도 모델을 각각 고결시료와 미고결 토양시료에 적용하고 비교해보았다(Fig. 8). 이 때 예측 값의 신빙성을 위해 여러 실험실의 자료를 수집하여 총 56개의 사암시료와 35개의 미고결 토양 시료를 이용하였다. 공극률 기반 모델과 표면적 기반 모델의 상대적 중요성은 사암시료 혹은 토양시료인지 여부에 따라 달라졌다. 왜냐하면 고결된 사암 샘플은 지층계수의 영향을 많이 받으며 미고결 토양 샘플은 허수 전기전도도의 영향을 많이 받기 때문이다. 사암에서는 지층계수 혹은 전기전도도 실수 성분이 유체투과도에 주로 영향을 미쳤다. 이는 현장 자료에 적용시 상당히 불리하다. 정확한 지층계수 값을 얻기 위해서는 다양한 염도에서 측정을 하거나 지하수보다 높은 염도에서 측정을 해야하기 때문이다. 반면, 미고결 토양 시료에서는 표면적(여기서는 허수 전기전도도 혹은 정규화 충전율)이 k 값을 예측하는 데 큰 영향을 미친다. 이는 지하수의 화학성분에 대한 보정 계수는 어느정도 필요하겠지만 현장에서 IP 측정을 하고 바로 유체투과도 예측에 이용할 수 있다는 것을 의미한다. 결과적으로 현장에서 전기적 측정자료로부터 유체투과도를 예측할 때 사암보다는 미고결 토양에서 더 신뢰성 높은 결과를 얻을 수 있을 것이다(Fiandaca et al., 2018).

Fig. 8. Measured values of permeability for (a) 56 sandstone and (b) 38 unconsolidated samples, respectively (Weller et al., 2015).​​​​​​​

지구물리학적 변수를 이용할 경우 미고결 토양의 유체투과도를 예측하는 것이 더 효과적이었는데 이는 곧 지층계수(F) 값이 비교적 작고 변동 범위가 크지 않을 때이다(Slater and Lesmes 2002; Maurya et al. 2018). Börner et al.(1996)은 표면적(혹은 공극 크기) 실수 전기전도도가 전해질 전기전도도와 표면 전기전도도 모두에 영향을 받는다는 한계점을 IP와 표면적 특성 관계를 단순화 시키면서 극복할 수 있었다. 하지만 F가 \(\sigma^{\prime \prime}\)에 영향을 미칠만큼 커지게 되면 공극 반지름과 공극망(porous network)의 공극에 모두 IP 측정 자료에 영향을 미치면서 k 예측값이 불명확해지기 때문에 제한된 범위에서만 유체투과도 예측 모델을 적용할 수 있게 되었다.

유효수리반지름(r_eff)은 유체투과도를 예측하는 핵심 변수이기 때문에 r_eff 와 IP 이완 시간(\(\tau\))의 관계성을 알아내고자 하는 연구도 있었으나 유의미한 관계성을 찾지는 못하였다(Weller et al., 2015). 이는 IP 이완시간을 구성하는 확산 계수(D₊)가 일정하다는 전제가 있는데 이 확산계수와 실제 사암 샘플에서는 이 확산 계수가 일정하지 않기 때문에 적용하는 데에 매우 제한적이다.

6. 현장 적용 사례

5장에서 언급했듯이 낮은 진동수를 이용하는 물리탐사(전기비저항 및 IP) 자료로부터 매질의 수리지질학적 특성을 파악하고자 할 때, 상대적으로 F 값이 작고 그 변동범위가 작은 미고결 토양 지대에서 예측이 유리하기 때문에 현장 적용 사례는 아직까지 대부분 토양으로 이루어진천부 지역에서 이루어졌다. IP 변수를 기반으로 예측한 유체투과도 분포 모델을 검증하기 위해 주로 입자 크기 분포나 양수시험(pumping test) 등을 통한 유체투과도 계산 값과 비교하며 일반적으로 그 차이가 상용로그 스케일에서 1(10배) 이내이면 수용 가능한 오차 범위로 본다(Hördt et al., 2007; Fiandaca et al., 2018).

Hördt et al.(2007)은 지표 전기비저항과 IP 자료를 역산하고 식(15)와 Börner et al.(1996)의 모델 변수를 활용하여 Krauthausen 지역에서 수리전도도 분포를 예측하였다. 이 때 진동수 지수 c에 대한 정보가 없었기 때문에 Börner et al.(1996)가 제안했던 미고결 토양에서 변동 범위(2.8 < c < 4.6)로 설정하였다. 이 지역에서는 S_por 값이 비교적 작았기 때문에(\(\leq 10 \mu \mathrm{m}^{-1}\)) c 값의 변화에 큰 영향을 받지 않고 k 예측 값도 오차범위 내에서만 변동했다. 또한, 복소 전기전도도 모델을 역산할 때 일반적인 거칠기 제한조건 대신 집중 제한조건(focusing regularization)을 주어 지하수층 경계가 좀 더 명확하게 역산되도록 하여 수리전도도 예측의 정확도를 향상시킬 수 있었다.

Hördt et al.(2009)은 총 4개 지역(Hasloh, Lüdingworth, Kappelen, Krauthausen)에서 얻은 2차원 IP 탐사 자료를 기반으로 유체투과도 분포를 예측했으며 이를 종합하여 IP 자료 기반 유체투과도 예측에 대한 가능성을 재고하고 고려해야 할 사항을 정리하였다. Hasloh와 Lüdingworth 지역에서는 크기 10배 오차 범위 내에서 유체투과도 분포 예측이 잘 이루어졌다. 반면, Kappelen에서는 IP 자료 기반으로 예측한 유체투과도 값이 다소 과소평가 되게 나타났다. 이는 해당 지역에서 얻은 자료의 위상값이 1 mrad 이하로 매우 작았는데 IP 자료 측정에 사용한 장비(SIP256)의 잡음 수준이 1 mrad였기 때문에 정확도가 떨어진 것으로 판단되었다. 입자 크기가 큰 경우 위상 지연이 매우 작아지기 때문에 정확도에 영향을 미칠 수 있다. Krauthausen 지역에서 얻은 IP 자료로 예측한 유체투과도 값은 입자 크기 기반 유체투과도 값보다 매우 크게 나왔는데 이는 이 지역에 입자 크기가 불균질했기 때문에 IP기반 유체투과도 예측이 잘 이루어지지 않는 것으로 추측하였다.

Fiandaca et al.(2018)은 덴마크 Grindsted 쓰레기매립지에서 시추중검층(logging while drilling; LWD) IP 탐사 자료를 이용해 수리전도도 분포를 유추하였다. 이 매립지 침출수에는 각종 유기물 및 무기물 오염물질이 포함되어 있다. 오염지역에서 전기적 특성을 이용해서 수리전도도 분포를 예측할 때는 지하수를 제외한 오염물질이 전기적 특성을 가지면 안된다는 전제조건을 내재하고 있다. 즉, 전기비저항이 낮거나 IP 효과가 있으면 측정 자료에 강한 영향을 주기 때문에 수리전도도로 인한 IP영향을 보기 힘들게 된다. 시간영역 IP자료를 측정하고 변수를 개선한 Cole-Cole 모델이 사용되었다. 유체투과도는 보정하지 않은 실내실험 실험식을 이용해 계산되었으며 IP 자료를 통해 유추된 k 값은 slug 테스트 방법과 입자 크기 분석으로 추정한 k 값과 비교하였다(Fig. 9). IP 자료를 통한 k 값 추정은 slug 테스트 측정값과 매우 유사하게 나타났으며 입자크기 분석으로 추정한 k 값과도 10배(1 order of magnitude) 이내로 차이가 발생하였다. 이 수치는 실내실험 수준과 비슷한 예측 정확도이다. 또한, Fiandaca et al.(2018)은 현장 스케일에서 IP 자료로부터 유체투과도를 추정할 때 이완시간 기반으로 하는 것보다 허수 전기전도도로 할 것을 추천하였다. 왜냐하면 복소 전기전도도 자료에서 이완시간을 구하기 위해 역산을 하면서 불확실성이 증가하기도 하며 시추공 자료가 아닌 지표 자료의 경우 시간상수 분포도를 구하는 것이 더 어려워지기 때문이다.

Fig. 9. Cross plots between measured k (k_meas) (from grain size analyses and slug tests, Bh1-3, and Sc1) and IP-derived k (k_IP) for El-logs E1, E2, and E3. Dashed lines show deviation of 1 order of magnitude (Fiandaca et al., 2018).​​​​​​​

같은 현장에서 Maurya et al.(2018)은 총 16개 측선으로 지표 탐사를 수행하고 얻은 시간영역 IP 자료를 통해 매질의 수리전도도 분포를 영상화하였다. IP 자료로부터 계산한 유체투과도 값은 slug 테스트와 입자 크기 분석으로부터 얻은 결과와 매우 좋은 상관관계를 보였다. IP 변수로 얻은 k 값 중 몇 개 자료는 다소 큰 오차를 보였는데 이는 낮은 유체투과도를 가진 층 사이에 낀 얇은 모래층 자료로 확인되었으며 역산할 때 주는 평활화 제한 조건(smoothness constraint)으로 인한 결과로 보인다. 또한, 지하수 전기전도도와 역산 전기전도도 값이 양의 상관관계를 보였으며 이 두 변수로부터 계산한 지층계수의 평균 값은 약 F = 5.1로 나왔다. 유체투과도와 유체 전기전도도를 영상화하였을 때 점토 혹은 갈탄(lignite)의 지질학적 분포와 물의 전기전도도 분포가 훨씬 잘 구분되었다. 최종적으로 이러한 정보를 기반으로 3차원 지질모델을 구성하였으며 예측 유체투과도 모델과 유체 전기전도도 모델도 도시하였다(Fig. 10).

Fig. 10. (a) 3-D Geological model, (b) 3-D permeability model and (c) 3-D water electrical conductivity model (Maurya et al., 2018)​​​​​​​

7. 결론

유체투과도를 예측하기 위해 기존에 사용하던 직접 시추 방법을 물리탐사가 대체할 경우 비교적 짧은 시간 안에 넓은 범위의 지하 매질 특성을 파악할 수 있다는 장점이 있다. 특히, 유체투과도를 결정하는 중요한 요인 중 하나인 수리학적 반지름(r) 및 공극 부피당 표면적(S_por)이 전기적 이중층(EDL)에 영향을 미치게 되는데 이를 기반으로 IP 탐사를 이용한 유체투과도 예측법이 관심을 받게되어 많은 연구가 이루어지게 되었다. 지구물리학적 변수를 이용한 유체투과도 예측모델은 IP 변수와 수리학적 변수의 관계를 이용해 유체투과도 모델을 재구성하여 만들어졌다. 수리학적 변수를 이용한 유체투과도 예측모델로 는 실험적 모델인 PaRiS 모델과 역학적 모델인 KT 모델이 있다. PaRiS 모델은 다공성 매질의 공극 부피당표면적 S_por을 이용하며 S_por와 허수 전기전도도 \(\sigma^{\prime \prime}\)는 비례관계를 이용해서 허수 전기전도도를 통한 유체투과도 예측 모델을 구할 수 있다. KT 모델은 S_por의 역수인 를 이용하였으며 시간상수 \(\tau\)로 대체하여 유체투과도 예측 모델을 구할 수 있다.

기존에 허수 전기전도도와 S_por가 단순 비례한다는 전제하에 실내실험을 진행하여 허수 전기전도도 기반의 유체투과도 예측모델을 구성하였다. 하지만 기존 연구들의 시료 개수가 부족하다는 문제점을 해결하기 위해 연구진들이 실험 자료를 취합하여 더 많은 시료로 관계를 분석한 결과 지층계수(F)도 반비례 관계가 있다는 것이 밝혀지면서 모델을 \(\sigma^{\prime \prime} \propto S_{p o r} / F\)로 개선하면서 더 정확한 모델을 구성할 수 있게 되었다. 현장 탐사에서 IP 변수를 이용해 유체투과도를 예측할 때 알갱이 입자 크기 분포가 비교적 균질하고 상대적으로 F 값이 작고 그 변동 범위가 작은 미고결 토양 지대에서 예측이 유리하기 때문에 대부분 토양으로 이루어진 천부 지역에서 이루어졌다. 현장 스케일에서 IP 변수로 유체투과도를 예측하는 기법은 아직 모든 현장에서 적용될 수 있는 것은 아니지만 모델을 좀 더 개선시키는 방안 등을 통해 실내 실험과 비교적 유사한 오차범위 내로 유체투과도 예측을 할 수 있었다. 전기비저항 및 IP 자료를 이용해 유체투과도를 현장 스케일로 예측하게 되면 복잡한 실내실험 과정을 거치지 않고 비교적 빠른 시간 내에 넓은 지역의 수리지질학적 분포 특성을 파악할 수 있기 때문에 다양한 분야에서 유용하게 적용될 수 있을 것으로 판단된다.