암반의 불균질성을 고려한 불포화대 지하수 유동 평가

Numerical Simulation of Groundwater Flow in Feterogenetic Rockmass of Unsaturated Condition

Abstract

본 연구에서는 단열을 포함하는 불포화 암반에서의 지하수 유동 흐름을 예측하기 위한 2차원 수치 모델링을 수행하였다. 특히, 불연속 단열망 모델링(Dscrete Fractured Network, DFN)을 통하여 도출된 투수계수의 분포를 나타내는 k-field를 모델링 입력변수로 적용하였다. 투수계수 차이는 불포화대에서 지하수의 이동속도 차이를 야기시키는 중요한 인자로 적용되었다. 연구지역의 지표 토양층으로부터 지하 대수층까지의 불포화대 깊이를 적용하여 지하수 유동 모델링의 초기조건을 실제와 유사하게 설정하였다. 강우의 지표 침투율은 인공구조물과 자연 토양층의 침투율 차이를 적용하여 실제 불포화 암반을 거쳐서 포화대까지 지하수의 이동흐름을 해석하고 시각화하였다. 특히, 오염물질의 이동 시작점이 될 수 있는 인공구조물의 하부에 모니터링 지점을 설정하여 실제 오염물질원이 지하수에 용해되어 불포화대를 이동할 때의 경로를 예측 하였으며 중력방향인 직하부로 이동하는 것을 확인되었다.

We present the results of two-dimensional numerical simulations predicting the flow of groundwater in a fractured unsaturated zone. We applied the k-field distribution of permeability derived from discrete fracture network (DFN) modeling as the hydraulic properties of a model domain. To model an unsaturated zone, we set the depth from the ground surface to the underground aquifer. The rate of water infiltration into the unsaturated zone was divided into two parts, an artificial structure surface and unsaturated soil zone. The movement of groundwater through the unsaturated zone was simulated with particular emphasis on contaminant transport. It was clearly observed that the contaminants dissolved in groundwater transported vertically from the ground surface to the saturated zone.

서 론

지표와 지하 대수층 사이에 존재하는 불포화 구간 중 지표를 포함한 상부 토양층에 대한 지하수 유동 연구는 활발히 진행되어 왔다(Mualem, 1976; van Genuchten, 1980; Parviz and William, 1984). 반면 토양층 하부로 존재하는 불포화 암반에 대한 연구는 매질 특성 규명에 대한 어려움이 존재하며, 단열·단층 존재 등 불확실성에 대한 문제로 인하여 상대적으로 연구 대상으로 한계가 존재한다. 그러나 지표에서 발생한 오염물질, 강우의 유동 경로, 지표수의 침투 특성 등 암반으로 유입되는 유체에 대한 연구가 계속될 때 불포화 암반의 유체 유동에 대한 불확실성이 저감될 것이다. 불포화 암반에서의 유체 유동에 대한 연구는 인간생활에서 발생하는 다양한 오염물질 등의 지하 침투와 연계된다.

따라서, 우리나라의 지하 암반의 지질 특성 중 단열의 존재를 고려한 연구지역을 선정하여, 불포화대 해석을 위하여 본 연구를 수행하였다. 다만, 기존의 불포화 토양층은 다공성 매질을 대상으로 수행된 물질이동 방정식 적용이 가능하였지만(Bear, 1972), 연구대상 암반은 불포화 암반으로 다공성 매질이 아니기 때문에 토양층 불포화 해석과는 다른 가정이 필요하다. 본 연구에서는 기반암이 존재하며, 암반 내부로 발달한 단열들로 이루어진 불포화 암반의 전체 시스템을 다공성 매질과 동일한 특성을 가진다는 개념으로 모델링을 수행하였다(Liu, Zhang and Ahlers, 2004). 지표에 인공 구조물과 강우 특성을 고려한 침투수량을 산정하였으며, 암반에 존재하는 단열 특성을 반영할 불포화 암반을 대상으로 지하수 침투시 유동 경로와 특성을 파악하기 위하여 TOUGH2 전산코드를 사용하여 본 연구를 수행하였다(Pruess, 2010).

 

연구대상

연구지역은 약 50~60 m 두께로 존재하는 불포화대 암반을 포함하는 임의 지역으로 선정하여, 지하수 유동 관련 불확실성 저감을 위하여 연구를 수행하였다. 연구지역은 산마루 기준, 남동, 북동으로 각각 지하수의 흐름이 형성되어 있다. 산마루 기준으로 남쪽에 위치한 인공 구조물의 지하수 흐름은 동쪽 방향으로 형성되어 있다(KRMC, 2008). 따라서 모델링 영역은 지하수 흐름을 고려하여(Fig. 1)과 같이 북서-남동 방향인 대각선 방향으로 2D 절단면을 선정하여 도메인을 구성하였다.

Fig. 1.Bird eye's view of the study area.

 

수치해석

개념모델

다상 유체 및 열 흐름에 대한 지배 방정식은 유체의 차원 및 구성요소의 수와 관계없이 동일한 수학적 형태를 가진다. 본 연구에서 활용된 TOUGH2 (Transport Of Unsaturated Groundwater and Heat) 전산코드의 지배방정식인 질량, 에너지 보존 방정식에 대한 표현을 동시에 하면 식 (1)과 같다(Pruess, Oldenburg, and Moridis, 1999). TOUGH2는 주 유체의 흐름과 전송모듈이 다른 유체의 속성 모듈과 서로 영향을 줄 수 있는 모듈형 구조로 설정되어 있다. 이러한 특징 때문에 TOUGH2는 다차원, 다상 유동 시스템의 광범위한 처리를 할 수 있다. 밀도, 점도, 엔탈피 등과 같이 혼합유체의 열-물리학적 변수들은 지배방정식에 입력되어 적절한 상태방정식(equation of state, EOS)의 형태로 제공된다.

식 (1)에서 Vn는 유동 시스템에서 임의의 격자, Γn : 표면 경계면, M는 부피당 질량 또는 에너지, k는 질량 구성요소(물, 공기, 수소, 용액 등), F는 질량 또는 열 유속, q는 소스 그리고 n은 표면 요소 dΓn의 법선 벡터를 각각 나타낸다.

모델링에서는 cell-centered 유한체적법을 기반으로 도메인을 구성하는 각 격자의 질량, 에너지 보존 방정식을 해석하였다(Hinds, 2001). 본 모델링에서는 water-air의 물리학적 특성과 불포화대 생성 후 침투되는 지하수의 유동 특성을 고려하기 위하여 EOS7r 모듈을 이용하였고, 모델링 결과의 시각화를 위하여 후처리 프로그램으로 Tecplot을 이용하였다.

단열을 포함한 불포화 암반 매질에서의 다상 유동 특성을 결정짓는 모세관압과 상대투과율 해석 도구로써 van Genuchten 함수를 사용하였다. 함수의 주요 입력변수인 air entry pressures와 λ는 실험값을 사용하였고, 암반 매질의 permeability 값은 문헌 자료를 바탕으로 비례식을 통하여 도출하였다(Leverett, 1940).

경계조건 및 초기조건

영역을 둘러싼 경계 조건으로는 불투수 모암층으로 설정하였고, 불포화대로 침투된 지하수가 이동 후 경계층에서 발생할 수 있는 오류를 최소화하기 위하여 포화대 최상부 도메인의 가장자리 모암 격자의 체적인자(volume factor)에 대하여 일정수두 경계조건(Dirichlet Boundary)을 부여하였다. 인공 구조물 하부부터 지하수대까지 분포하는 약 40~50 m 깊이의 불포화대 초기조건 설정을 위하여 (Fig. 2(a)~(b))와 같이 모델링 초기 포화조건으로부터 깊이에 따른 격자별 수분 포화도 차이를 나타내는 불포화대 초기 조건을 설정 하였다(Fig. 3(a)~(b)).

Fig. 2.(a): Distribution of hydraulic pressure in saturated initial condition, (b): Distribution of water saturation in saturated initial condition.

Fig. 3.(a): Distribution of hydraulic pressure on the unsaturated initial condition, (b): Distribution of water saturation on the unsaturated initial condition.

입력변수

연구지역의 시추공 코어시료 11개 시편을 채취하여 수은 주입 모세관압력 방법(Mercury Injection Capillary Pressure, MICP)을 이용하여 불포화 암반 매질에 대한 공극률, 밀도를 산정하였으며 모델링 입력변수로 평균값을 적용하였다(Table 1).

Table 1.Porosity and density of Matrix on unsaturated zone.

MICP 분석을 통하여 포화도에 따른 모세관압 곡선을 도출하였으며, 이후 van Genuchten 함수 식 (2)에 적용하여 모델링 입력인자로 활용되는 P0, λ 값을 도출하였다. 모델링 입력변수로 Table 2와 같이 평균값을 적용하였다.

Table 2.Saturation-capillary pressure curve and P0, λ.

P0는 매질의 기체투과 초기 압력(gas entry pressure), λ는 van Genuchten 함수의 경험적 매개 변수 값(λ = m = 1 − 1 / n), Slr는 액상의 잔류 포화도(residual liquid saturation), Sls는 액상 포화도(full liquid saturation), Sgr는 기상의 잔류 포화도(residual gas saturation)를 각각 나타낸다.

모델링 입력변수는 (Table 3)과 같이 요약하였다.

Table 3.Input parameter.

또한, 불포화대에 존재하는 단열의 특성을 반영하기 위하여 단열망 모델링 결과를 통하여 얻어진 투수계수 분포도(k-field)를 적용하였다. k-field로 구성된 도메인 전체 단열분포에 대한 통계처리를 통하여 투수계수 범위를 10구간으로 설정하였으며(Table 4), 각 구간별 평균값을 해당 그리드에 적용하였다. 이와 같은 방법으로 구간별 투수계수 차이에 따른 도메인을 구성하는 단열 10구간 및 퇴적암 암반 매질 등을 포함하는 총 9,057개의 요소로 모델링 전영역을 구성하였으며(Fig. 4), 600년의 기간에 대해 해석하였다.

Table 4.Variation of permeability in modeling domain.

Fig. 4.Grids in modeling domain and monitoring points for sensitivity analysis.

지표 침투수량

불포화대로 침투하는 지표수량은(Fig. 4)와 같이 크게 두 영역(I, II)으로 구분하여 설정하였다. 도메인 영역 중심에 존재하는 인공구조물의 길이에 해당하는 약 310 m의 (I)영역과 대기에 직접 노출된 (II)영역으로 구분하여 지표 침투수량 조건을 설정하였다.

(II)영역의 침투수량은 장기(1985.1.1~2014.12.31, 30년) 수리·수문 분석을 통하여 산출하였다. 강우량, 기온, 최대 일조량 등의 분석을 통하여 증발산량, 직접유출량, 기저유출량, 연간 침투수량을 산정하였다(Table 5). 특히, 연간 침투 수량은 식 (3)과 같이 물수지 분석법을 적용하여 도출하였으며, (II)영역에 적용된 침투수량으로 469.76 mm/yr 값을 적용하였다.

Table 5.Hydraulic characteristics during 30years (1985.1.1~2014.12.31).

인공구조물의 하부인 (I)영역의 침투수량은 (II)영역의 침투수량을 기반으로 인공구조물 침투 후 최종적으로 불포화대로 직접 침투하는 양으로, 불포화대 지하수 유동 모델링(I)영역의 침투수량으로 32 mm/yr의 값을 가정하여 적용하였다.

 

결과 및 토의

모델링결과

불포화대가 형성된 초기조건에서 인공구조물 하부로 일정한 비율로 강우의 침투가 이루어지며, 포화대 도달하기까지의 불포화 암반 매질에서 포화도 변화가 시작된다. 모델링 결과 검토 시 불포화대의 포화도 변화에 초점을 맞추어 결과 분석을 하였다. 600년간의 모델링에 대한 포화도 변화 결과는(Fig. 5(a)~(l))와 같다. (Fig. 5(a))에서 나타낸바와 같이 해석기간 1년 후 포화도 변화는 불포화대-포화대 경계지점에서 포화도가 다소 증가한 것을 관찰할 수 있다. 약 30년 후까지 경계면에서의 포화도가 증가하며 50년 이후 지표 부근의 포화도가 증가하기 시작한다. (Fig. 5(e))의 해석 기간 100년 후 포화도 변화는 불포화대에 불균질하게 분포하고 있는 단열의 방향, 위치와 일치하는 포화도 증가 양상을 더욱 뚜렷하게 확인할 수 있다. 이로 인하여 불포화대로의 침투수가 암반매질이 아닌 보다 높은 투수계수를 가지는 단열 방향으로 침투와 이동이 우세하게 일어나고 있는 것으로 판단된다.

Fig. 5.(a): Saturation changes in the unsaturated zone after 1 year, (b): Saturation changes in the unsaturated zone after 20 years, (c): Saturation changes in the unsaturated zone after 30 years, (d): Saturation changes in the unsaturated zone after 50 years, (e): Saturation changes in the unsaturated zone after 100 years, (f): Saturation changes in the unsaturated zone after 150 years, (g): Saturation changes in the unsaturated zone after 200 years, (h): Saturation changes in the unsaturated zone after 250 years, (i): Saturation changes in the unsaturated zone after 300 years, (j): Saturation changes in the unsaturated zone after 400 years, (k): Saturation changes in the unsaturated zone after 500 years, (l): Saturation changes in the unsaturated zone after 600 years.

Fig. 5(g)에서 나타난 바와 같이 해석기간 200년 후에는 모델링 전 영역에 대한 포화도 증가 양상 중에서 특히 460~550 m 사면에서 포화도 증가가 두드러지게 나타나며, 침투수량 조건과 지하수 유동 방향으로 인하여 사면 부근에서 침투수량의 증가가 집중되는 것을 확인할 수 있다. 하지만, 600년 후에도 포화도는 0.35 정도를 유지하고 있으며, 이 결과는 이 시점에도 완전 포화되지 않음을 나타내는 것으로 판단된다.

Fig. 5를 통해 단열 발달 구조에 대하여 확인할 수 있고, 발달된 단열을 통한 중력방향으로의 지하수 유동이 이루어짐을 확인할 수 있다. 포화도 변화 양상을 확인해보면, 지점별로 투수계수의 차이로 인한 최대 도달 포화도의 차이는 있지만 약 150~200년 사이에 모든 모니터링 지점에서 포화도가 최대에 도달하는 것으로 확인된다. 그리고 최대치 도달 후 최대 포화도 상태를 유지하는 경향을 보이고 있다. 이러한 현상은 일정한 강우조건에서 불포화 상태 매질의 포화도가 점점 증가하다가 하부로 이동하는 지하수량과 지표로부터 유입되는 침투수량이 평형에 도달하였기 때문에 포화도는 일정한 수준으로 유지한다고 판단된다.

A, C, E 경사면 지점에서 수평방향으로 10 m, 20 m 정도 암반 내부로 이동한 지점에서의 시간에 따른 포화도 변화를 분석하였으며, 모니터링 위치는 (Fig. 6)에 나타내었다.

Fig. 6.Monitoring points for a comparison with water saturation.

경사면 A 지점으로 부터 내부로 10 m, 20 m씩 이동한 A', A"의 포화도 최대치는 0.32, 0.34로 경사면 최대 포화도인 0.37보다 낮은 값을 가지는 것으로 나타났다(Fig. 7). 또한 C로부터 내부로 10 m, 20 m씩 이동한 C', C"의 포화도 최대치는 모두 0.28로 경사면 노출 지점 최대 포화도인 0.35보다 낮은 값을 가지는 것으로 나타났다(Fig. 8).

Fig. 7.Comparison with water saturation at the monitoring point A, A' and A''.

Fig. 8.Comparison with water saturation at the monitoring point C, C' and C''.

E로부터 내부로 10 m, 20 m씩 이동한 E', E"의 포화도 최대치는 0.28, 0.29로 경사면 노출지점인 E에서의 최대 포화도 0.35 보다 낮은 값을 가지는 것으로 나타났다(Fig. 9). 이와 같이 경사면으로부터 내륙으로 이동한 지점의 포화도 분석결과 같은 해수면 기준 높이에서 경사면에서 내부로 들어갈수록 포화도가 감소하는 경향이 나타났다.

Fig. 9.Comparison with water saturation at the monitoring point E, E' and E''.

특히, 인공구조물에서 유입되는 지하수의 유동 특성을 파악하기 위하여 (Fig. 10(a))와 같이 인공구조물 하부를 40 m 간격으로 8지점을 설정하여 지하수의 이동경로를 분석하였다. 인공구조물 하부로부터 유입된 지하수는 불포화대를 지나면서 투과 정도에 따른 이동 속도 차이가 나타났다. 하지만 8지점 모두 중력방향으로 이동을 하고 시간이 지남에 따라 점차 하부로 이동하여, 포화대에 도달한 후 서에서 동쪽방향으로 이동하는 것을 확인할 수 있다(Fig. 10(b)~(e)).

Fig. 10.(a): Groundwater flow in the unsaturated zone after 10 years, (b): Groundwater flow in the unsaturated zone after 30 years, (c): Groundwater flow in the unsaturated zone after 100 years, (d): Groundwater flow in the unsaturated zone after 300 years, (e): Groundwater flow in the unsaturated zone after 600 years.

민감도 분석

적용된 입력변수 중 불포화 환경에서 지하수 이동 특성을 결정짓는 주요인자로 예상되는 van Genuchten 함수 입력변수 λ와 Po에 대하여 민감도 분석을 수행하였다. 불포화영역 내에서 지하수 이동 모델링 시 지정한 320m 지점의 깊이 20m 지점을 A로, 440m 지점의 깊이 40m 지점을 B로 설정하였고, 민감도 모니터링 지점은 Fig. 4에 나타내었다.

단열매질의 van Genuchten 입력변수 λ

불포화대 모델링에서 단열의 입력변수는 단열망 모델링을 통해 도출된 투수계수를 제외하고는 암반 매질의 입력변수와 동일하게 적용하였다. 그러나, 불포화 환경에서 공극 크기 분포와 관련된 λ값의 민감도 분석을 위하여 λ값을 기준 모델링 0.65에서 ±0.1 적용, 0.55, 0.75로 설정하였다. A, B 지점에서의 각 λ 민감도에 따른 모델링 결과는 (Fig. 11, 12)와 같다.

Fig. 11.Sensitivity analysis of λ at the monitoring point A.

Fig. 12.Sensitivity analysis of λ at the monitoring point B.

λ의 민감도 분석 결과, 0.75일 때 시간변화에 따른 가장 많은 양의 지하수가 이동하였음을 알 수 있다. A, B 두 지점 모두에서 λ값이 커질수록 모니터링 지점에서의 최대 유동량이 증가하였다. λ값은 공극분포의 균질성과 관련된 인자로써, 암반 매질 균질성이 증가함에 따라 이동성이 증가하는 것으로 해석된다. 즉, 암반매질 내의 균질한 공극구조로 인하여 지하수의 유동 속도가 빨라지므로, 모니터링 지점에서 확인되는 시점도 빨라지는 것으로 판단된다.

λ값 변화에 따라 지하수 유동 시간과 최대 유동량에 차이는 나타나지만 최초 모니터링 지점에서 확인되는 시간은 0.55, 0.65, 0.75 일 때 각각 26년, 18년, 10년이고, 최대 포화량은 0.9, 0.91, 0.93으로 차이는 크지 않은 것으로 확인된다.

단열매질의 van Genuchten 입력변수 Po

불포화대 모델링에서 단열의 입력변수는 단열망 모델링을 통해 도출된 투수계수를 제외하고는 암반 매질의 입력변수와 동일하게 적용하였다. 그러나 불포화 환경에서 매질로의 기체 투과 초기압력과 관련된 Po값의 민감도 분석을 위하여 Po값을 기준 모델링 7.91E-5에서 10배, 5배, 1/5배, 1/10배를 적용한 7.91E-4, 3.95E-4, 1.58E-5, 7.91E-6으로 설정하였다. A, B 지점에서의 각 Po 민감도에 따른 모델링 결과는 (Fig. 13, 14)와 같다.

Fig. 13.Sensitivity analysis of Po at the monitoring point A.

Fig. 14.Sensitivity analysis of Po at the monitoring point B.

Po 값의 민감도 결과, 실험값 입력 데이터 기준 압력의 1/10의 값을 가질 때 가장 높은 이동성을 나타내었으며, 최초 변화가 시작된 시점이 20년으로 가장 빠른 변화가 확인되었다. Po 값이 작아짐에 따라 이동성과, 변화 확인 시점이 늦어지는 경향을 나타낸다. 특히 Po 민감도 분석에서 고려해야할 점은 실험값 기준으로 압력이 작아지는 경우에는 암반매질 사이의 공극으로 작은 압력으로도 침투가 일어나기 때문에 A×1/10, A×1/5, A 세 경우에서는 지하수 유동 또한 비교적 활발하게 이루어진다. 하지만 Po 값이 5배, 10배로 증가된 경우에는 이동성이 급격히 떨어지는 현상을 확인할 수 있다. 암반매질 사이의 공극으로 침투하기 위하여 큰 압력이 필요하며, 이러한 이유로 지하수의 이동성이 급격히 떨어지는 것으로 판단된다.

 

결 론

암반을 포함하는 불포화대에서의 침투 지하수 유동 예측을 위하여 본 연구를 수행하였다. 모델링 수행시 지표면과 포화대 사이의 실제 불포화대가 존재하는 두께 및 포화도를 반영한 불포화대를 생성하였으며, 불포화대에 존재하는 단열 특성을 고려하여 지하수 유동 모델링을 수행하였다. 모델링 결과, 불포화대에 존재하는 단열 발달 경로를 따라 지하수가 이동하는 현상을 확인하였으며, 인공구조물 하부로 부터 유입되는 지하수의 이동 경로를 분석함으로써 불포화 대에서는 중력방향으로 이동, 지하수대에 도달하는 현상을 확인하였다. 본 연구는 암반을 포함하는 불포화대내에서 지하수 유동예측을 위한 하나의 방법 제시일 뿐 불포화 암반을 해석하는 적합한 방법인지에 대한 검증은 다양한 해석 방법과 비교 검증을 통하여 이루어져야 한다고 판단된다.