서 론

최근 수평정시추(horizontal drilling) 및 수압파쇄(hyd-raulic fracturing)등과 같이 생산성을 증대시키는 기술의 향상으로 상업적 생산이 가능해져 셰일가스(shale gas)의 개발이 북미를 중심으로 활발히 진행되고 있다. 셰일가스는 현재 경제적으로 효율적인 개발이 가능한 신석유자원(unconventional resources) 중 가장 부존량이 많은 것으로 평가되고 있으며(Shin et al., 2012), 대규모 연속체 형태로 부존되어 있다. 셰일가스의 개발을 위해서는 생산 가스전의 생산성을 정확히 예측하고 해당 지역에 대한 평가 작업을 우선적으로 수행하여야 한다.

다양한 자원량 평가방법 중 대상 저류층 또는 생산정에 대한 상세 정보가 부족하거나 신속한 분석이 필요한 현장에서는 주로 생산감퇴곡선분석법(Decline Curve Analysis; DCA)을 이용하여 작업이 이루어지며(Lee et al., 2008), 이 방법은 일반적인 유･가스전뿐만 아니라 신가스자원(unconventional gas)에서도 활용되고 있다(Yang and Lim, 2013). 생산감퇴곡선분석법은 과거의 생산 자료를 이용하여 감퇴경향을 파악하는 방법으로 향후 생산 계획 수립을 위해 활용된다. 생산감퇴곡선분석법에서 Arps (1945)의 경험식 중 쌍곡선(hyperbolic)함수를 가장 보편적으로 사용하고 있으며, 감퇴곡선 형태를 결정하는 감퇴지수(decline exponent)를 0에서 1사이의 값으로 가정한다. 그러나 셰일가스 생산 시 발생하는 장기간의 천이유동으로 인해 감퇴지수 값이 1을 초과하는 특징이 나타날 수 있으며, 쌍곡선함수를 적용할 경우 궁극가채량(Estimated Ultimate Recovery; EUR)의 정확한 추정에 어려움이 존재하여 최근에는 셰일가스와 같은 신석유자원에 적용 가능한 다양한 생산감퇴곡선분석법들이 연구되고 있다(Joshi and Lee, 2013). 또한, 셰일가스 지층은 낮은 유체투과도로 인해 수압파쇄를 적용하여 발생된 균열로 생산이 이루어지며, 이로 인해 생산초기의 감퇴율이 큰 반면 후반에는 변화율이 매우 작아 초기 생산거동의 변화가 큰 특징이 나타나기 때문에 단일 예측 값을 산출하는 생산감퇴곡선분석법을 이용하여 생산량을 예측할 경우 생산거동의 불확실성을 고려하기 어렵다(Shin et al., 2014).

이 연구에서는 이러한 불확실성을 최소화하고 신뢰도를 높여주기 위해 대표적인 확률론적 기법인 몬테카를로 시뮬레이션(Monte Carlo simulation)을 적용하여 셰일가스의 궁극가채량을 추정하였다. 생산감퇴곡선분석법 중 일반적으로 사용되고 있는 Arps의 쌍곡선함수와 셰일가스의 생산 시 나타나는 유동형태 변화를 고려하여 생산량 예측이 가능한 Stretched Exponential Production Decline (SEPD)를 몬테카를로 기법에 적용함으로써 두 모델의 궁극가채량에 대한 분포형태를 비교 및 분석하고, 이를 누적확률분포(Cumulative Distribution Function; CDF)로 나타내어 확률에 따른 궁극가채량을 도출하고자 한다.

셰일가스 궁극가채량 추정 방법

셰일가스 지층에서의 생산감퇴곡선분석법 적용

셰일가스는 지층의 두께가 얇고 지질학적으로 광범위하게 부존되어 있는 신석유자원 중 비교적 많은 양을 차지하고 있다. 또한, 셰일가스 지층은 암체가 매우 치밀하고 공극 내 자유가스(free gas) 또는 유기물 및 광물입자에 흡착된 상태로 부존하고 있어 유체 유동에 제한이 있다(Ko, 2007). 근래에는 수평정시추 및 수압파쇄와 같은 기술의 향상으로 셰일가스 지층에 대한 대규모 자극을 발생시킴으로써 인위적인 균열을 생성하여 유체유동을 보다 원활하게 할 수 있다. 이를 통해 생산성이 증대된 셰일가스는 북미를 중심으로 활발한 개발이 이루어지고 있으며, 가스가격을 반영하여 생산이 진행되고 있다. 이러한 생산 가스전의 개발을 위해서는 해당 지역에 대한 자산 가치를 평가하고 잔존하는 자원의 양과 생산수명을 파악하여야 한다. 이를 결정하기 위해서는 해당하는 생산 가스전의 생산감퇴경향과 생산성을 예측하여야 하며, 일반적으로 경험적인 곡선 맞춤(curve fitting)에 의한 표준곡선(type curve)을 도출하거나 혹은 Arps의 경험식과 같은 생산감퇴곡선분석법에 의해 수행된다. 이 기법은 시간에 따른 생산량 변화 관계로부터 회수량을 평가하는 방법으로 생산 자료의 변화에 민감하게 반응하여 자료의 변화가 크지 않을 때 양질의 분석결과를 도출할 수 있다.

Arps의 경험식

생산감퇴곡선분석법의 대표적인 방정식인 Arps의 경험식은 감퇴지수()의 값에 따라 지수(exponential)함수, 쌍곡선함수, 조화(harmonic)함수 3가지로 구분되며(Table 1), 생산정의 유동형태가 경계영향유동(Boundary Dominated Flow)임을 가정하고 있다. 일정한 초기생산량(initial production rate; )에서 시간이 경과함에 따라 감퇴지수가 생산량의 변화에 미치는 양상을 Fig. 1에 나타내었으며 감퇴지수 값이 클수록 장기적인 생산이 이루어짐을 알 수 있다. 지수함수 형태로 감퇴 하는 특성은 통상 균열 저류층 또는 매우 불균질한 저류층에서 나타나며, 궁극가채량을 과소평가하는 경향이 있어 초기에 생산을 집중할 필요가 있는 저류층에 대해 현금흐름을 계획할 경우에 많이 활용된다. 조화함수 형태로 감퇴하는 경우, 수공법(water flooding)을 적용하는 저류층이나 저류층 주변의 대수층에 의해 생산할 때 나타나지만 일반적인 특성은 아니다(KIGAM, 2002). 대부분의 저류층은 감퇴지수가 0과 0.7 사이의 범위를 가지는 쌍곡선함수에 의해 묘사되며(Arps, 1956), 일반적인 저류층에서의 생산감퇴경향을 반영하는 쌍곡선함수는 현재까지도 궁극가채량 추정을 위해 가장 보편적으로 사용되고 있다.

Stretched Exponential Production Decline(SEPD)

수압파쇄를 적용한 저류층에서 유체 유동 형태는 초기의 천이유동이 장기간 나타나며, 경계영향유동이 나타나는 시점이 불확실한 특징을 나타낸다. 이로 인해 감퇴지수가 1을 초과하지 않는 경계영향유동을 가정하는 Arps의 쌍곡선함수를 적용할 경우 궁극가채량이 과대 산정되거나 실제 유량의 거동을 반영하지 못하여 정확한 생산량 예측에 어려움이 있다. 최근에는 셰일가스와 같은 신석유자원에 적용 가능한 다양한 생산감퇴곡선분석법들이 연구되고 있으며, Valko와 Lee(2010)는 셰일가스 생산 시 나타나는 두 가지 구간의 유동특징을 고려하여 생산량을 예측할 수 있는 SEPD를 제시하였고 Can과 Kabir (2012)은 SEPD 모델을 이용하여 확률론적 곡선을 도출하는 연구를 수행한 바 있다. 이 기법은 쌍곡선함수와 초기 생산거동은 비슷하나 감퇴지수가 1을 초과하는 경우를 고려할 수 있으며, 자료의 양이 많을 경우 특별한 보정을 하지 않아도 적합한 모델로 사용될 수 있다. 또한, 생산기간과 생산량에 대한 제한조건 없이 생산량 예측이 가능하고 궁극가채량이 발산하는 경향을 나타내는 쌍곡선함수와 달리 SEPD는 수렴하는 경향을 나타낸다(Fig. 2).

Valko와 Lee(2010)와 Yuan 등(2012)은 SEPD를 이용하여 궁극가채량을 추정 시 36개월의 생산 자료를 이용하는 것을 제시하였으며, 추정기간은 30년으로 설정하였다. SEPD의 생산량 및 궁극가채량을 구하는 방법은 각각 식 (1), (2)와 같으며 SEPD의 주요 변수인 , 는 식 (3)의 연립 방정식을 통해 구할 수 있다. 이는 식 (4)를 통해 계산이 가능하며, 최종적으로 식 (5)를 통해 36개월의 생산 자료에 대한 평균값인 초기생산량 를 도출할 수 있다. 일반적으로 SEPD 이외의 다른 기법에서는 초기생산량 값을 지정하는 경우도 있지만, 셰일가스의 수평정에 대한 생산 자료를 이용하여 SEPD를 적용 할 경우 값을 계산하여 추정함으로써 보다 적합한 예측 값을 도출할 수 있다(Yuan et al., 2012).

                        (1)

                (2)

                                     (3)

    (4)

            (5)

몬테카를로 시뮬레이션(Monte Carlo simulation)

셰일가스 지층은 낮은 유체투과도와 유정완결기술 적용으로 인해 초기에는 균열을 통한 유체유동이 일어나 감퇴율이 크지만, 후반에는 기존 암체에 의한 유동이 발생하여 변화율이 작아져 생산거동의 변화가 큰 특징을 나타낸다. 이러한 셰일가스 지층의 생산 자료를 단일 예측 값만을 도출하는 결정론적 방법의 생산감퇴곡선분석법에 적용할 경우, 생산거동의 변동성으로부터 파생되는 불확실성이 존재하므로 이를 최소화하기 위해 이 연구에서는 대표적인 확률론적 기법인 몬테카를로 시뮬레이션을 적용하였다.

몬테카를로 시뮬레이션이란 자원량 산정에 필요한 저류층 입력변수들의 가능한 모든 범위를 확률 분포로 나타내고, 확률적 의미를 가질 만큼 여러 번 무작위로 각 변수 값을 추출하여 자원량의 예상 분포범위와 이것의 확률을 계산하는 방식이다(Fanchi, 2012). 이러한 몬테카를로 시뮬레이션을 수행하기 위해서는 자료의 특성에 따른 확률분포형태를 결정해야 하며, 자료의 개수가 많은 경우는 정규분포(normal distribution) 또는 로그정규분포(log-normal distribution)를 사용하고 자료의 개수가 적은 경우는 균등분포(uniform distribution) 또는 삼각분포(triangular distribution)를 주로 사용한다(Budzaya, 2013). Fig. 3과 같이 변수에 대한 확률분포형태를 결정한 후 누적확률분포로 나타내어 y축에서 난수를 발생시키고, 이에 해당하는 x축의 값을 추출하여 얻어진 변수 값을 생산감퇴곡선분석법의 식에 대입함으로써 확률론적인 결과 값을 도출할 수 있다. 몬테카를로 시뮬레이션은 모든 종류의 분포에 적용가능하며, 이론적 배경만으로 구하기 어려운 수치들을 쉽게 구할 수 있다. 또한, 입력 값의 확률분포 형태가 분석에 큰 영향을 미치므로 각 모델의 입력 값에 맞는 확률분포를 선택하여 난수를 발생시켜야 하며, 이들 과정을 반복 수행하여 확률론적인 결과 값을 도출할 수 있다.

몬테카를로 시뮬레이션을 결합한 생산감퇴곡선분석법 적용 과정

생산감퇴곡선분석법의 두 모델인 쌍곡선함수와 SEPD에 몬테카를로 시뮬레이션을 적용하여 궁극가채량을 추정한 연구 절차는 Fig. 4와 같다. 우선적으로 셰일가스전의 생산 자료를 이용하여 각 모델의 변수를 추정하였으며, 도출된 각각의 변수에 대하여 분포 형태를 지정하여 난수를 발생 시키고 이에 해당하는 값을 누적확률분포를 통해 구하였다. 이 때, 발생시킨 난수의 수만큼 생산감퇴곡선을 도출할 수 있으며, 생산감퇴곡선 별로 궁극가채량을 도출하여 누적확률분포로 나타내어 10%의 확률에서의 값을 P10, 50%에서의 값을 P50, 90%에서의 값을 P90으로 분류함으로써 확률론적인 결과를 얻을 수 있다.

미국 Barnett shale 지역의 생산 자료 적용

연구대상 지역

셰일가스의 생산이 활발히 이루어지고 있는 미국의 경우, 2,000년대 이후 약 28,000개의 셰일가스정에서 연간 700 Bcf 이상의 양을 생산하고 있다. 그 중 비교적 오랜 기간 동안 생산을 진행하고 있는 Barnett shale 지역은 미국 텍사스 주 포트워스 분지에 위치하며(Fig. 5), 시추기술 및 유정완결기술의 향상으로 인해 회수가능량이 증가하고 있는 지역이다. Barnett shale의 주요 생산지역은 4,000 km²까지의 면적을 포함하며, 회수가능 자원량 14.7 Tcf에 달하고 약 1,700개 생산정에서의 궁극가채량 중간값(median) 및 최대값(maximum)은 각각 0.52 Bcf, 7 Bcf로 생산성이 매우 높다. 주요 생산 지역에 비해 셰일가스 지층의 두께가 얇고 상, 하부 파쇄 지지대의 일부가 결여된 지역의 면적은 12,500 km²으로 아직까지 탐사·개발이 미진하여 생산이력이 짧으며, 약 134개의 생산정에서 궁극가채량 중간값은 0.28 Bcf, 최대값은 2.4 Bcf이고 회수가능 자원량은 11.6 Tcf로 추정된다(Ko, 2007). 이 연구에서는 Barnett shale지역의 4개 수평정 Well #1～Well #4의 생산 자료를 이용하였으며, 생산기간은 각각 72, 44, 88, 90개월이다(Fig. 6).

생산감퇴곡선 변수 추정

Arps의 경험식 중 하나인 쌍곡선함수와 SEPD를 이용하여 궁극가채량을 추정하기 위해 36개월까지의 생산 자료를 이용하여 변수를 추정하였으며, 그 결과는 Table 2와 Table 3에 나타내었다. 쌍곡선함수의 경우 생산 자료의 감퇴경향을 반영하여 실제 값과 예측 값의 차이인 잔차(residual)를 최소화하는 곡선을 도출하는 비선형회귀분석(nonlinear regression)을 이용함으로써 감퇴지수, 초기감퇴율(), 초기생산량의 값이 도출되었다. 이 때 Table 2에 제시된 바와 같이 쌍곡선함수의 감퇴지수는 최소 1.2517에서 최대 1.7674로 도출되어 이 연구에서 적용한 Barnett shale의 수평정 자료의 경우 감퇴지수가 1을 초과하는 특징을 나타내는 것을 확인할 수 있다. SEPD 경우 식(3)을 이용하여 12개월까지의 누적생산량과 24개월까지의 누적생산량의 비로 구성된 방정식과, 12개월까지의 누적생산량과 36개월까지의 누적생산량의 비로 구성된 방정식을 동시에 만족시키는 , 값을 결정하였다. 도출된 , 값과, 생산 자료에서 취득한 생산기간() 및 생산기간에 따른 생산량() 값을 식(5)에 대입하여 을 도출한 후, 36개월의 생산 자료에 대한 평균값인 을 계산하였다. 생산 자료의 실제 초기생산량 값과 쌍곡선함수와 SPED를 이용하여 추정한 초기생산량 값에 대한 평균 제곱근 오차(Root Mean Square Error; RMSE)는 각각 7.90, 78.20으로, 실제 초기생산량 값과 쌍곡선함수의 초기생산량 값의 차이가 SEPD보다 작게 도출됨을 알 수 있었다.

변수 분포 지정 및 값 생성

몬테카를로 시뮬레이션을 적용하기 위해 우선적으로 앞서 추정한 변수들에 대한 분포형태를 지정해야 한다. 생산 자료가 많을 경우에는 평균과 분산 등의 정보를 필요로 하는 정규분포를 사용하거나 또는 생산 자료 값에 로그를 취한 새로운 값이 정규분포를 따르는 로그정규분포를 사용할 수 있지만(Budzaya, 2013), 이 연구에서는 생산 자료가 충분하지 않아 변수의 정확한 분포추정 어려우므로 자료의 최소값(minimum), 최빈값(mode), 최대값에 대한 비교적 간단한 정보가 필요한 삼각분포를 이용하였다. 이 때 최빈값은 가장 빈도수가 높은 값이며(An, 2005), 앞서 각 모델에서 추정한 변수들의 값을 최빈값으로 지정하였다. 최빈값을 기준으로 한 범위 지정 시 Lin 등(2005)은 5～20%로 범위를 지정하였으며, 이 연구에서는 셰일가스 생산거동의 큰 변화폭으로 인한 불확실성을 최대한 고려하고자 20%로 범위를 적용하였다. 이를 통해 쌍곡선함수와 SEPD의 세 변수들은 Well 별로 다른 형태의 삼각분포가 도출되었고 이를 누적확률분포로 나타내어 난수를 발생시켰으며, 확률론적 결과를 도출하기에 충분히 큰 5,000회를 시행하여 난수에 해당하는 값을 구하였다(Han and Moon, 2012).

생산감퇴곡선 도출

쌍곡선함수와 SEPD의 변수 분포에 대한 난수발생횟수와 동일한 개수만큼 감퇴곡선을 도출하였으며, 그 결과는 Fig. 7～Fig. 10과 같다. 그래프의 검은색 점은 궁극가채량 추정을 위해 사용한 36개월의 생산 자료이고, 흰색 점은 36개월 이외의 사용되지 않은 생산 자료를 나타낸다. 분석을 실시한 모든 생산정에서 쌍곡선함수가 SEPD 보다 시간 경과에 따라 감퇴율이 점점 작아지는 양상이 나타남을 알 수 있다.

궁극가채량 추정 결과

미국 Barnett shale 지역의 4개의 수평정에 대한 생산 자료를 이용하여 궁극가채량을 추정하고자 쌍곡선함수와 SEPD 모델에 몬테카를로 시뮬레이션을 적용하였으며, 이에 대한 히스토그램(histogram)과 누적확률분포의 결과는 Fig. 11～Fig. 14와 같다. 히스토그램의 범위와 평균을 Table 4에 나타내었고, 자료의 산포정도를 나타내는 표준편차(standard deviation)와 자료가 평균을 중심으로 한쪽으로 치우친 경향을 나타내는 척도인 왜도(skewness)를 Table 5에 표시하였으며 이와 같은 통계치를 통해 히스토그램의 분포형태를 정량적으로 파악할 수 있다. 이 때 왜도 값은 분포모양이 평균을 중심으로 우측으로 긴 꼬리를 가지는 경우 양수(+)로 나타나고, 좌측으로 긴 꼬리를 가지는 경우에는 음수(-)로 나타나게 된다. Table 5를 통해 Well #1, #3, #4의 SEPD 모델은 쌍곡선함수에 비하여 표준편차 값이 더 크게 도출되어 산포가 큰 형태의 확률분포로 나타나는 것을 알 수 있으며, 각 모델의 왜도 값이 양수로 나타나는 것을 통해 히스토그램 분포형태는 모두 평균을 중심으로 우측으로 긴 꼬리를 가진 형태임을 알 수 있다. 특히 SEPD 모델의 왜도 값은 쌍곡선함수에 비해 더 크게 나타나므로 SEPD의 확률분포 또한 보다 평균을 중심으로 우측으로 긴 꼬리를 가지는 형태임을 파악 할 수 있다. 쌍곡선함수와 SEPD 모델 모두 Well #2, #4 보다 Well #1, #3의 히스토그램 분포가 상대적으로 넓게 도출되어 Well #1, #3의 P90과 P10의 값 차이가 더 크게 나타난다. Well #2, #4 보다 Well #1, #3의 초기생산량 추정 값과 사용한 생산 자료의 산포정도가 더 큰 특징을 나타내기 때문에 20%의 범위를 지정한 Well #1, #3의 초기생산량의 삼각분포 범위가 더 넓게 생성되며, 이로 인해 히스토그램의 분포 형태에도 영향을 끼친 것으로 판단된다.

몬테카를로 시뮬레이션을 적용하여 도출된 P10, P50, P90의 값과 결정론적 방법을 이용하여 구한 궁극가채량 값을 Table 6에 정리하였다. SEPD 모델의 P90과 P50 값의 차이와 P50과 P10 값의 차이의 간격이 쌍곡선함수 보다 큰 것을 알 수 있으며, 이를 통해서도 SEPD의 왜도 값이 크게 나타나는 경향과 히스토그램의 분포양상을 예상할 수 있다. 쌍곡선함수와 SEPD 모델의 P50의 값은 Well #1의 경우 각각 2588.60 MMcf, 2002.12 MMcf, Well #2는 1658.86 MMcf, 1140.81 MMcf, Well #3은 5734.92 MMcf, 5650.50 MMcf, Well #4에서는 1776.62 MMcf, 1631.62 MMcf로 도출되었으며, P10, P90 값과 결정론적 방법을 이용하여 구한 궁극가채량 크기 또한 Well #3, #1, #4, #2 순서로 나타났다. 이를 통해 초기생산량 값과 사용한 생산 자료의 산포정도가 클수록 궁극가채량이 크게 추정되는 경향이 나타나는 것으로 판단할 수 있다. Well #2, #4 보다 Well #1, #3의 초기생산량 추정 값이 더 크기 때문에 궁극가채량 또한 크게 도출된 것을 알 수 있다. 또한, 셰일가스 생산 초기 장기간의 천이유동을 고려할 수 있는 SEPD 방법과 달리 경계영향유동을 가정하는 쌍곡선함수 식의 본질적인 특성이 반영되어 도출된 생산감퇴곡선으로 인해 쌍곡선함수를 이용