1. 서 론

하천과 대수층이 수리적으로 서로 연결되어 있을 때 하천 주변에서 지하수를 양수하면 지하수위 강하 정도에 따라 하천으로의 지하수유출량이 줄어들거나 하천수가 지하수로 유입되어 결국 하천의 기저유량이 감소하게 된다. 이렇게 지하수 양수로 인해 하천의 유량이 감소하는 것을 하천수 고갈 또는 감소(stream depletion)라 일컫는데, 이의 양적 시간적 변화를 모사하기 위해 대수층과 하천의 기하학적, 수리적 특성을 단순화한 조건에 대해 다양한 해석적 모형들이 개발되어 왔다.

균질의 반무한(semi-infinite) 대수층의 바닥까지 하천이 완전 관통한 조건에 대해 지하수 양수 경과 시간에 따른 하천수 감소량을 계산하는 최초의 해석해가 ^Theis(1941)에 의해 개발되었다. 이후 동일한 조건에 대해 형태만 약간 변형된 해석해가 ^{Glover and Balmer(1954)}에 의해 개발되었다. 여기에 ^{Hantush(1965)}는 하천과 대수층 경계면의 지체 영향까지 고려한 해석해로 확장하였고, ^Hunt(1999)는 미소 하폭(line width)과 저투수성의 하상 매질을 가진 하천이 양방향 무한대수층 상단부 일부를 차지한 조건에 대해 하천수 감소량과 지하수위 저하량을 계산하는 해석해를 유도하였다. ^{Baalousha(2012)}는 하천 누수만 있을 때와 양수정만 있을 때의 해를 구하고 중첩원리를 적용하여 ^Hunt(1999) 해석해를 같은 형태로 다시 유도하였고, 또한 유한하폭(finite width)에 대한 해석해도 새롭게 개발하여 하폭의 영향을 직접적으로 고려할 수 있게 개선하였다. ^{Butler et al.(2001)}은 유한 하폭뿐만 아니라 유한 폭의 대수층 조건에 대한 해석적 모형을 개발하였다. 이 중에서 Hunt(1999) 해석해는 미국 USGS에서 전산프로그램화하여 배포하였고 뉴질랜드에서는 지하수 양수에 따른 하천수 감소 영향 평가를 위한 가이드라인(^{Pattle Delamore Partners Ltd and Environment Canterbury, 2000})의 대표적인 방법으로 제시하는 등 국외적으로 활용성이 높으며, 국내에서도 ^Kim(2010)과 ^{Lee et al.(2016c)}이 실제 및 가상 문제에 적용한 사례가 있다.

대수층의 기하학적 조건을 보다 실제적으로 고려하기 위해서 복수의 층으로 구성된 누수대수층(leaky aquifer) 조건에 대한 해석해가 개발되었다. ^Hunt(2003)는 상부에 반대수층(aquitard)이 있는 무한의 반피압대수층(semi-confined aquifer)에서 양수로 인한 하천수 감소량을 산정하는 해석해를 유도하였고, ^Hunt(2008)는 이를 유한의 하폭과 대수층에 대한 해석해로 확장하였다. ^{Butler et al.(2007)}과 ^{Zlotnik and Tartakovsky(2008)}은 수두가 일정한 하부대수층 위에 상부대수층이 반대수층으로 연결되어 있을 때 상부대수층 지하수 양수로 인한 하천수 감소량을 계산하는 해석해를 개발하였다. ^Hunt(2009)는 상부대수층뿐만 아니라 하부대수층의 수두도 미지수로 함께 고려하여 상부층 양수에 따른 하천수 감소량 및 지하수위 저하량을 산정하는 해석해를 유도하였으며, ^{Ward and Lough(2011)}는 동일 조건에서 하부층 양수에 따른 해석해를 개발하였다. 특히 ^Hunt(2009)과 ^{Ward and Lough(2011)} 해석해는 ^{Lee et al.(2016a}; ^2016b)와 ^{Lee and Chung(2020)}에서 다양한 수리조건에 대한 적용성이 평가되었고, 환경부 한강홍수통제소의 하천 인근 지하수 양수 영향을 검토하는 방법(^{HRFCO, 2018}) 중 한 가지로 전산프로그램화 되어 하천수사용관리시스템에 탑재되었다(^{Hong et al., 2019}).

이상과 같이 지하수 양수에 따른 하천수 감소량을 계산할 수 있는 여러 해석해들이 개발되어왔으며, 비록 하천-대수층-양수정 시스템 조건을 단순화하였지만 투수량계수, 저류계수, 하상수리전도도, 하천으로부터 관정의 이격거리, 양수시간 등 대수층과 하천의 몇 가지 수리특성치만으로 하천수 감소량을 산정할 수 있기에 실무적으로 유용하다. 그러나, ^Hunt(1999), ^Hunt(2009), ^{Ward and Lough(2011)} 등 대부분의 해석해가 하천의 규모가 매우 작다는 가정하에 하폭을 점으로 간주하고 유도되었기 때문에 하폭이 큰 대하천에 적용할 경우 오차가 발생할 수 있다. ^{Lee et al.(2018)}는 미소하폭에 대한 Hunt(1999) 해석해와 유한하폭에 대한 ^{Baalousha (2012)} 해석해를 비교하여 하천-관정 이격거리가 하폭 보다 긴 조건에서 사용되어야 미소하폭 간주에 따른 오차를 줄일 수 있음을 밝혔다. ^{Lee et al.(2018)}의 연구는 단일 대수층 양수로 인한 하천수 감소량을 산정할 때 하폭 크기에 따른 ^Hunt(1999) 해석해의 정확도를 분석한 것으로 이 외에는 하폭의 다양한 크기에 따른 해석해의 적용성을 평가한 사례는 찾아보기 어렵다. 따라서 본 연구에서는 상, 하부 대수층으로 구성된 누수대수층에 대한 ^Hunt(2009) 해석해를 다양한 수리특성 조건에 대해 적용하여 하폭의 크기가 하천수 감소량 계산 결과에 미치는 영향을 분석, 고찰하였다.

2. 재료 및 방법

2.1 Hunt 해석해

^Hunt(2009)는 Fig. 1과 같이 상부, 하부의 두 개 대수층(aquifer)이 투수성이 낮은 반대수층(aquitard)으로 연결되어 있고, 상부 대수층이 하천과 수리적으로 연결되어 있는 조건에 대해서 상부 대수층 지하수 양수로 인한 지하수위 저하량과 하천수 감소량을 산정하는 해석해를 개발하였다.

Fig. 1

A Two-Layered Leaky Aquifer, Stream, and Shallow Well Pumping for Hunt’s Analytical(2009) Solution

Eq. (1)은 지하수 양수량( $Q$)로 나누어져 무차원 형태로 나타내어진 Hunt(2009)의 하천수 감소량( $∆Q$) 산정식으로서 상부와 하부 대수층의 투수량계수 $T_1$, $T_2$ 및 저류계수 $S_1$, $S_2$, 반대수층의 수리전도도 $K_v$와 두께 $B$, 하천과 양수정간의 거리 $L$, 그리고, 하상수리전도도 $K_s$와 하폭 $W$의 곱을 하천바닥층 두께 $M$으로 나눈 값인 하상수리전도성 $\lambda$ 등의 함수로 나타내어져 있다.

여기서, $∆Q$는 하천수 감소량, $Q$는 양수량, $L^{-1}$는 라플라스 역변환, 그리고 기타 변수들은 Eqs. (2,3,4,5,6)와 같다.

^Hunt(2009) 해석해는 2차원 지하수 유동 방정식에 하천과 양수정 위치에 Dirac delta 함수( $\delta$)를 사용한 생성(source)항을 고려하여 해를 구한 것이다. 하폭을 점으로 간주하였지만 하폭의 함수인 하상수리전도성 $\lambda$을 입력변수로 사용하고 있어 불완전하지만 하폭 크기 영향을 간접적으로 반영할 수 있다(^{Lee et al., 2018}). 따라서 하폭이 작은 경우에는 해석해 적용 결과에 큰 문제가 없을 것으로 예상할 수 있으나, 규모가 큰 하천에 적용할 경우에는 하폭의 간접 고려에 따른 오차가 발생할 수 있다. 예를 들어, 하폭이 500 m이고 하상수리전도도가 1 m/d인 경우와 하폭이 100 m이고 하상수리전도도가 5 m/d인 경우는 하폭과 하상수리전도도가 달라 수치해석모형으로 정밀하게 하천수 감소량을 모의하면 서로 다른 결과를 나타내게 되지만, 하상수리전도성 $\lambda$가 같기 때문에 해석해로는 동일한 결과를 얻게 된다.

본 연구에서는 대수층의 수리특성 조건으로서 1층의 투수량계수 2, 10, 20 m²/d(수평 수리전도도는 0.1, 0.5, 1 m/d), 2층의 투수량계수 2, 10, 20 m²/d(수평 수리전도도 0.02, 0.1, 0.2 m/d)로 각각 3가지 조건을 부여하여 다양한 투수조건에 대해 하천수 감소량을 산정하였다. 연직 수리전도도는 수평 수리전도도의 1/10인 것으로 가정하고 각 층의 대수층 두께를 고려하여 연직 수리전도도의 조화 평균값을 구한 후에 이를 1층과 2층의 누수계수 값으로 활용하였다. 또한 대수층의 저류성을 나타내는 입력변수로서 1층의 저류계수는 0.02, 0.1, 0.2, 2층의 저류계수는 0.001, 0.01, 0.05로 역시 각각 3가지 조건을 부여하였다. 하천의 수리특성으로 하상 수리전도도는 0.001, 0.01, 0.1 m/d, 하천바닥층 두께는 1.0 m, 하폭은 100, 300, 500, 700 m, 하천-관정 이격거리는 100 m 간격으로 100~1,000 m의 값을 부여하였다. 따라서 총 9,720 가지 대수층과 하천의 수리특성 조건에 따라 지하수 양수에 따른 하천수 감소량(이하 하천수 감소비)을 모의하였다. 지하수 양수로 인한 하천수 감소비는 일단위로 지하수 허가기간인 5년 동안 모의하고, 그 평균값을 분석대상으로 하였다.

2.2 MODFLOW 모형의 구축

하천규모가 대수층에 비해 작다는 가정하에 유도된 Hunt 해석해의 적용성을 평가하기 위해서 범용적으로 사용되고 있는 3차원 지하수 유동모형 MODFLOW으로 모의한 수치해와 비교하였다. MODFLOW 모형(^{McDonald and Harbaugh, 1988})은 연속방정식과 Darcy’s law를 결합하여 유도된 Eq. (7)과 같은 지하수 유동방정식을 지배방정식으로 하고 있으며, 함양(recharge), 하천(stream), 주입 또는 양수정(well), 배수관(drain), 수평차수벽(horizontal barrier), 증발산(evt) 등의 경계조건하에 유한차분법으로 미지수인 지하수 수두에 관한 수치해를 구한다.

여기서, $K_{ij}$는 수리전도도 텐서, $h$는 지하수 수두, $W$는 단위체적당 유입 또는 유입량, $S_s$는 다공성 매질의 비저류계수, 그리고 t는 시간을 나타낸다.

Hunt 해석해 유도 조건과 유사하게 상, 하부 두 개의 대수층과 중앙을 가로지르는 하천, 그리고 상부 대수층을 관통한 양수정을 고려하여 Eq. (7)을 육면체로 이산화하였다. 해석해가 수평방향으로 무한 대수층 조건에서 유도된 점을 감안하여 MODFLOW 모델 영역을 가로 10,000 m, 세로 10,000 m로 크게 하여 지하수 양수 모의 결과가 모델 영역 끝단 경계조건의 영향을 받지 않고 내부의 하천 경계조건에 지배적으로 작용하도록 하였다. 격자 한 개의 크기를 100 m×100 m로 하여 총 101개 행과 101개 열로 구성하였다. 대수층은 상부의 1개의 충적층과 그 아래 1개의 암반층으로 총 2개의 레이어를 구성하였고, 충적층은 두께 20 m의 자유면 대수층, 암반층은 두께 100 m를 가지는 피압/자유면 변환 대수층으로 설정하였다. 각 층의 두께는 고정시켰지만 수리전도도 값을 여러가지 값을 부여하여 지하수 유동량에 지배적인 인자인 투수량계수 값을 변화시켰다. 수리전도도, 저류계수, 하상수리전도도, 하폭 등 대수층과 하천의 수리특성값은 앞선 해석해 적용시와 동일한 값을 이용하였다. 경계조건으로서 영역의 외곽 경계는 무흐름으로 간주하여 불투수경계로 설정하였고, 하천 셀에서는 일반수두(general head) 경계조건으로 설정하여 하천수위와 지하수위의 수위차에 따라 하천-지하수 유동량이 계산되도록 하였다.

일반적으로 수치해의 정확도를 평가하기 위해 해석해와 비교, 검증을 수행하지만, 본 연구에서는 역으로 해석해의 정확도를 평가하기 위해서 정교하게 모의된 수치해와의 비교를 수행하였다. 이를 위해서 수치해의 수렴 조건을 보수적으로 부여하고, 수치모의시 입력과 출력 인자간의 질량보존 결과를 엄밀하게 평가하였다. 따라서 MODFLOW 모의시 수렴 조건으로 허용오차를 0.001 m로 하였고, 전체 모의 조건 중에서 매 연산시간 종료마다 출력되는 모의영역 내외로 유입·출입되는 총지하수량간의 불일치율(percent discrepancy)이 1 %를 초과한 결과는 비교 분석 대상에서 제외하였다.

3. 결과 및 고찰

Figs. 2 and 3은 각각 하폭이 100 m와 500 m일 때 주어진 다양한 수리특성조건에 대해 Hunt 해석해와 MODFLOW 모형으로 계산한 천부 지하수 양수 기간 5년 평균 하천수 감소비를 1:1로 비교하여 나타낸 것이다. Fig. 2는 하폭이 100 m이고 하천과 관정간의 이격거리가 100 m와 300 m일 때 해석해와 수치해를 도시한 것으로 두 값간의 결정계수가 거의 1.0로서 서로 일치된 결과를 나타낸 것을 알 수 있다. 정교하게 모의된 수치해를 참값으로 간주하고 해석해와의 상대제곱근오차(Root Mean Squared Error, RMSE)와 최대오차(Maximum Error, ME)를 계산한 결과 하천-관정 이격거리가 100 m일 때 각각 0.006, 0.018로 오차가 무시할 정도로 매우 작으며, 이격거리가 300 m 일 때 역시 RMSE와 ME가 각각 0.003, 0.008로 오차가 매우 작다.

Fig. 2

Comparison of Analytical and Numerical Solutions for Streamflow Depletion (W=100m)

Fig. 3은 하폭이 500 m일 때의 해석해와 수치해를 비교한 것으로 결정계수가 0.9를 초과하여 역시 두 해가 서로 잘 일치하고 있으나, Fig. 2에 비해 하폭 증가의 영향으로 두 해간의 산포도가 증가하고 결정계수가 줄어들었으며, 하천과 관정간의 이격거리가 짧은 100 m인 경우 강한 하천경계조건 영향으로 더욱 그러한 양상이 나타나고 있다. 이러한 수치해와 해석해간의 차이가 발생하는 이유는 수치해는 하천 횡단면에 대해 셀별로 다르게 계산된 지하수위를 이용하여 하천수와 지하수의 상호 교환량을 산정하지만 해석해는 하나의 등가 지하수위를 이용하기 때문이다. 따라서 하폭이 증가하여 하천셀 갯수가 많아질수록 수치해와 해석해의 차이가 커질 가능성이 높다. RMSE와 ME가 하천-관정 이격거리가 100 m일 때는 각각 0.040, 0.105, 이격거리가 300 m일 때는 각각 0.029, 0.059의 오차가 발생하였다. 따라서 미소 하폭에 대해 유도한 해석해로 하천수 감소량을 산정할 때 하폭 증가에 따라 오차가 커지므로 사용에 유의해야 한다.

Fig. 3

Comparison of Analytical and Numerical Solutions for Streamflow Depletion (W=500m)

Fig. 4는 주어진 수리특성 조건에 대한 해석해와 수치해의 차이를 하폭 및 하천-관정 거리에 따라 종합하여 나타낸 것이다. 하나의 하천-관정 거리별로 4가지 하폭에 대한 결과를 함께 나타내었기에 약간의 간격을 두고 도시하였다. 하폭이 100 m일 때에는 하천-관정 이격거리와 상관없이 해석해와 수치해간의 차이가 거의 없으며, 하폭이 점차 증가함에 따라 두 해간의 차이의 변동폭과 최대치가 커지고 그 양상은 하천-관정 이격거리가 짧을수록 분명해지고 있다. 하폭이 작으면 하폭 방향을 따라 지하수두 분포의 차이가 적지만, 하폭이 증가하면서 하상퇴적층내 지하수두의 비균질적 분포로 인해 하천수와 지하수간의 상호 교환량 크기가 달라지기 때문이다. 최대 차이는 주어진 조건에서 하폭이 가장 클 때인 700 m이고 하천-관정 이격거리가 가장 짧을 때인 100 m일 때 0.13 만큼 발생하였다. 전반적으로 Hunt 해석해를 이용하여 하천수 감소비를 산정할 때 발생하는 오차는 크지 않으나, 오차가 0.05 이하로 보다 정확한 결과를 얻으려면 대략 하폭 만큼 관정이 하천에서 이격되어야 하는 것으로 분석되었다.

Fig. 4

Difference of Stream Depletions Computed Using Analytical and Numerical Solutions According to Stream-Well Distance

우리나라 지하수법에서는 하천구역 경계로부터 300 m 이내를 ‘하천 인근’으로 정의하고, 이 지역에 대해서는 하천의 수량에 미치는 지하수 양수의 영향을 사전에 검토한 후 지하수개발·이용 허가를 해주고 있다. 따라서 Fig. 4의 전체 모의 결과 중에서 하천과의 이격거리 300 m 이내의 결과만을 추출하여 Fig 5와 같이 하폭 크기에 따라 박스플럿(box-plot) 형태로 나타내었다. 하폭이 증가함에 따라 변동폭이 점차 증가하고 있는 것을 확인할 수 있다. 주어진 조건에 대해 해석해가 수치해에 비해 최대 0.13 만큼 과다,–0.07 만큼 과소 산정되었으나, 중앙값이 거의 0에 근접하고 1분위와 3분위간 폭이 0.05 이내로서 전반적으로 해석해의 오차가 크지 않다. 따라서 하천 인근 지하수 양수에 따른 하천수 감소량을 산정할 때 해석적 모형을 사용해도 큰 무리는 없을 것으로 판단된다. 다만, 하폭이 500 m 보다 클 때에는 최대 0.1을 넘는 오차가 발생할 수 있음을 감안하고 실무에 활용해야 할 것이다.

Fig. 5

Boxplots of Difference of Stream Depletions Computed Using Analytical and Numerical Solutions for below 300m of Stream-Well Distance

하천으로부터 관정의 이격거리를 제곱하고 이를 대수층의 수리확산계수( $D=T/S$)로 나눈 것을 하천고갈인자(Stream Depletion Factor, SDF)이라 하는데, 대수층과 하천의 다양한 수리특성 조건에 대한 해석해와 수치해의 차이를 Fig. 6과 같이 1층의 SDF값에 따라 나타내었다. 또한 하상수리전도성 $\lambda$, 2층의 수리확산계수 $D_2$, 1층과 2층의 수리확산계수의 비 $D_1/D_2$, 그리고 1층과 2층간의 누수계수 $K_v$의 범위에 따라 Figs. 6(a)~6(d)와 같이 회색과 검정색 표식으로 구분하여 나타내었다. SDF 값이 작을수록 해석해와 수치해의 차이 및 변동폭이 크게 나타나고 있다. 이는 SDF 값이 작을수록 즉 하천-관정 이격거리가 작거나 수리확산계수가 클수록 하천 경계조건 영향이 양수정에 도달하는 시간이 빨라진데서 기인한다. 물론 SDF 값이 클수록 하천수감소비가 줄어드는 양상(^{Lee et al., 2016b})을 일반적으로 보이기 때문에 오차 또한 작게 산정된 원인도 포함되어 있다. SDF 값이 약 400일 이하에서는 해석해와 수치해의 차이가 0.1을 넘는 것으로 나타났고, SDF 값이 약 3,000 일 이상에서는 두 해의 차이가 0.05로 매우 작아 Hunt 해석해로도 정밀하게 모의한 수치해에 준하는 결과를 얻을 수 있는 것으로 분석되었다. Fig. 6(a)와 같이 하상수리전도성이 $\lambda <1m/d$일 때 해석해와 수치해간 차이가 대부분 양의 값을 보여 해석해가 과다 산정되고 있다. 즉, 하상의 수리전도도가 작거나 하폭이 작을 때 해석해가 과다 산정되는 경향을 주로 보이고 있고, 대수층의 수리확산계수가 작고 하천-관정 이격거리가 작아 SDF가 작을수록 더욱 이러한 특성을 보이고 있다. Fig. 6(b)는 해석해와 수치해의 차이를 1층의 SDF와 2층의 수리확산계수 $D_2$에 따라 나타낸 것으로 $D_2<100m^2/d$로 작을 때 해석해가 과다 산정된 경우가 대부분이며, $D_2\ge 100m^2/d$ 일 때는 다른 수리특성치와의 복합적 작용으로 과다 혹흔 과소 산정되고 있다. Fig. 6(c)는 1층과 2층의 수리확산계수의 상대비에 따라 오차를 나타낸 것으로 $D_1/D_2\ge 5$ 일 때 해석해가 과다 산정된 경우가 많은 것으로 분석되었다. Fig. 6(d)는 1층과 2층간의 누수계수 $K_v<0.0004m/d$로 작을 때 양의 오차가 발생하는 것으로 분석되었다. 하부 2층의 수리확산계수가 상대적으로 작거나 층간 누수계수가 매우 작아 2층의 수리특성 영향이 상부 1층으로 지체되어 전달될 때 즉, 1층의 수리특성치가 지배적인 조건이 된 경우에는 해석해가 과다 산정되는 것을 알 수 있다.

Fig. 6

Difference of Stream Depletions Computed Using Analytical and Numerical Solutions According to Hydraulic Properties

4. 결 론

본 연구에서는 하천-누수대수층-지하수 양수정 경계치 문제에 대해 ^Hunt(1999) 해석해와 MODFLOW 수치모형을 이용하여 대수층과 하천의 다양한 수리특성 조건별로 지하수 양수량 대비 하천수 감소량(하천수 감소비)을 산정하고 상호 비교를 통해서 해석해의 적용성을 평가하였다. 상부와 하부대수층의 투수량계수 및 저류계수, 층간 누수계수, 하천으로부터 관정의 이격거리, 하상퇴적층의 수리전도도, 하폭 등의 값을 변경하면서 천부대수층 양수에 따른 5년 평균 하천수 감소비를 해석해와 수치해로 각각 산정하고 그 차이를 분석하여 다음과 같은 결론을 얻을 수 있었다.

(1)하폭이 100 m로 작을 때에는 해석해와 수치해간의 결정계수가 1.0에 가깝고 최대 차이가 0.018로 두 해가 거의 일치하였으나, 하폭이 증가함에 따라 두 해간의 차이가 커져 하폭이 700 m, 하천-관정 이격거리가 100 m일 때 최대 0.13의 차이를 나타내었다. 정교하게 모의된 수치해를 참값으로 보고 미소하폭 조건에 대해 유도된 Hunt 해석해의 오차를 평가한 결과 하천-관정 이격거리가 하폭 이상되어야 그 오차가 0.05 이하가 되는 것으로 분석되었다.

(2)상부 1층의 하천고갈인자(SDF) 값에 따라 해석해 오차를 평가한 결과 하천-관정 이격거리가 작거나 수리확산계수가 커서 SDF 값이 작아질수록 오차가 증가하는 경향을 나타내었다. SDF < 400 d인 조건에서 해석해의 최대오차가 0.1을 초과하였으며, SDF > 3,000 d에서는 오차가 0.05 이하로 작게 나타났다.

(3)하상수리전도성이 $\lambda <1m/d$, 2층의 수리확산계수가 $D_2<$ $100m^2/d$, 1층과 2층의 수리확산계수비 $D_1/D_2\ge 5$, 1층과 2층간의 누수계수가 $K_v<0.0004m/d$ 일 때 해석해가 전반적으로 다소 과다하게 산정되는 경향을 나타내었다. 즉, 하상수리전도도와 층간 누수계수 값이 작아 양수에 대한 반응이 지체되거나, 하부 2층에 비해 상부 1층의 수리확산도가 큰 경우 해석해가 과다 산정될 가능성이 높다.