In this study, a three-dimensional hydrogeological model was developed to numerically deliniate contaminant transport and distribution in the area surrounding the Dalseong mine. The model’s accuracy was assessed by comparing its output with the field data. The spatial distribution of contaminants predicted by numerical simulation showed a good similarity to the field data obtained by electrical resistivity surveys, long-term groundwater observation, and the findings from previous hydrogeochemical studies. These results suggest that a 3D hydrogeological model can be used as an effective tool for predicting contaminant transport in groundwater, even in situations where field data is limited.

Keywords: 3D hydrogeological model, Groundwater flow, Groundwater contamination, Numerical modeling

한국광해광업공단(Mine Damage Safety)에 따르면, 2021년 광해 실태조사 결과 전국 5,475개 휴·폐광산 중 3,300개 광산에서 광해가 발생하는 것으로 파악되고 있다. 이 중, 달성광산은 1938년부터 텅스텐, 구리, 금, 은 등의 광산으로 개발되어 운영되다가, 1975년 이후 휴·폐광된 광산으로, 환경 오염 및 지역 사회의 경제적 피해를 야기하고 있는 대표적인 사례이다.
Jung(1996)은 달성광산 인근의 하상 퇴적물 분석을 통해 산성 광산 배수에 의한 오염 정도를 조사하여 오염 원인을 규명하고자 하였고, Lee et al.(2003)과 Choo et al.(2007)은 수리지구화학적 연구를 수행하여 해당 지역에서 폐석 침출수 및 갱내 유출수로 인한 오염이 정화시설에서 거의 정화되지 않고, 슈베르트마나이트(schwertmannite)의 침전이 오염 농도의 감소에 영향을 미칠 수 있음을 제시하였다. 최근 Kim et al.(2020)은 황산염의 동위원소 및 아연과 카드뮴의 비율(Zn/Cd)을 이용해 달성광산 주변에서 토양 및 지하수 오염의 주요 원인이 광미(mine tailings)임을 확인하였고, 후속 연구(Kim et al., 2023; Kim et al., 2024)에서는 구리, 납, 아연, 황 및 산소 동위원소 분석을 통해 오염원과 오염 메커니즘을 규명하고자 하였다. 오염 정화 또는 오염 저감을 위한 노력으로 1996년 산성 광산 배수 처리를 위한 반응조 시스템이 실험적으로 운영·평가되었고(Cheong et al., 1998), 1998년에는 5개의 소택지로 구성된 자연정화 처리시설이 설치되어 현재까지 운영되고 있다(Lee et al., 2003). 이와 더불어, 2012년부터 2015년까지는 슬래그 반응조를 적용한 정화 실험이 진행된 바 있다(MIRECO, 2015; Kim et al., 2018; Kim et al., 2021).
이와 같이 달성광산 주변 지역의 오염원인 규명 및 정화를 위한 연구는 지구화학적 혹은 수리지구화학적 측면에서 지속적으로 이루어져 왔으나, 오염 거동의 주요 기작인 지하수의 유동과 이를 통한 오염 경로 규명에 대한 연구는 찾아보기 어렵다. 이는 대수층 매질의 불균질성과 오염물의 생화학적 반응 등 다양한 요인으로 인하여 오염물의 전파 경로를 파악하고 추적하는 것이 매우 어렵기 때문이다. 따라서, 보다 상세하고 정밀한 오염 경로 규명을 위해서는 3차원 수치 모델링을 통해 수리지질학적 이해를 뒷받침할 필요가 있으며, 그 결과를 바탕으로 지하수와 오염물 거동의 관계를 해석하여야 한다.
대수층의 수리학적 특성, 수문학적 현황 및 지질학적 매질 분포 등 부지 특성에 대한 이해를 바탕으로 지하수 유동, 용질 이동, 열 전달 등 다양한 현상을 평가하거나 예측하기 위해서, 또는 오염지역에서의 오염 확산 경로와 범위 등을 예측하거나, 역방향 모델링을 통한 오염원 및 오염 경로 추정을 위해서, 3차원 수리지질 모델을 활용할 수 있다. Artimo et al.(2004)은 핀란드 Hitura 니켈 광산 주변 대수층 오염의 원인을 밝히기 위해서 3차원 모델링을 통해 오염 유출의 경로를 파악한 바 있고, Tomiyama et al.(2020)은 일본 야마가타현 Yatani 폐광에 대해 3차원 모델링을 통해 산성 광산 배수의 유출 경로를 모의한 바 있다. 또한 Stollberg(2012)는 독일 Bitterfeld-Wolfen 지역의 폐석탄광에 대해 역방향 모델링을 통해 오염의 근원지를 판별하였다. 국내에서는 광산 개발 및 운영 등이 지하수계에 미치는 영향을 예측하기 위해서(Lee et al., 2023; Kim et al., 2018; Yoon, 2011), 또는 오염 정화 및 저감 방법 적용의 효용을 예측하기 위해서(Park and Kihm, 2021) 3차원 수치모델링 연구가 수행된 바 있다. 이와 같이 국내에서는 정화의 관점에서 오염 범위 및 정화 효율성 평가 등에 3차원 수치 모델링이 활용되고 있을 뿐, 오염 확산 경로 규명에 대한 활용은 미진하다. 따라서 본 연구에서는 국내의 대표적인 폐광으로 인한 오염 사례인 달성광산에 대해 기존의 문헌자료 및 현장 조사를 바탕으로 3차원 수리지질 모델을 구축하고, 수치모델링을 통해 지하수 유동과 이에 따른 오염 경로를 규명하고, 현장 자료와의 비교를 통해 본 기법의 활용성을 제시하고자 한다.

2.1. 연구 지역
본 연구 지역은 대구광역시 달성군 가창면 상원리에 위치한 달성광산 인근으로, 대구지역 지하수기초조사(MOCT, 2004)에서는 신천상류 유역으로 분류되어 조사가 수행되었던 지역이다. Fig. 1은 본 연구 지역을 나타나고 있으며, 분수령에 해당하는 유역간 경계를 붉은색 실선으로 표시하였다. 그림에서와 같이 연구 지역의 북쪽 및 동쪽에는 각각 신천하류 유역 및 남천 유역이 인접하고 있다. 주요 오염원으로는 달성광산(Dalseong mine), 광미적치장, 달성광산에서 유래한 배출수를 처리하기 위한 자연정화 처리시설(passive treatment system)이 있으며 Fig. 1에서는 노란색으로 표시하였다. 2개의 하천이 달성광산의 남쪽과 북쪽에서 연구 지역의 서쪽으로 흘러가며, 남쪽 하천의 상류에는 저수지(Sangwon reservoir)가 있다. 하천 및 저수지는 Fig. 1에서 파란색으로 표시하였다.
2.2. 현장 조사 및 모니터링
2.2.1. 전기비저항 탐사
전기비저항 탐사기법 중 하나인 쌍극자 기법을 3개 측선에 대해 5 m 및 10 m의 전극 간격으로 수행하여, 충적층 및 암반층에서의 개략적인 지하수 유동 경로 및 심도를 파악하고자 하였다. 연구 지역 내 전기비저항탐사 측선은 Fig. 2에 노란색 점선으로 표시하였으며, 측선별 탐측 결과는 Fig. 3에 도시하였다.
Line-A의 전극 간격 5 m 탐사 결과에서는 지표하 7~15 m 인근에서 저비저항대를 확인할 수 있었고, 10 m 간격 탐사에서는 지표하 20~30 m 인근에서 저비저항대를 확인할 수 있었다. Line-B에서는 5~25 m, 그리고 Line-C에서는 5~30 m에서 전극 간격 5 m 및 10 m의 모든 탐사 결과에 대해서 저비저항대가 관측되었다. 전기비저항 탐사 결과를 바탕으로 현장 상황을 고려하여 지하수 장기 관측을 위한 NX 규격의 25 m 깊이 관정(BH)을 설치하였으며, 그 위치는 Fig. 2에 빨간색 점으로 표시하였다. BH 지점은 Line-C의 좌측 저비저항대를 관통하는 위치이다.
2.2.2. 지하수 관측
총 25 m 깊이 중 유공관의 범위는 지표하 5~25 m까지이며, 대수층의 발달 정도와 관측정 설치 후 초기 전기전도도(EC) 수직분포를 고려하여 지하수 유동 심도로 판단되는 8 m, 13 m, 22 m 각각에 Schlumburger사의 CTD-diver를 설치하고, 관측정 보호시설 내부의 지표면에 Schlumburger사의 BARO-diver를 설치하였다(Fig. 4 참조). 2022년 1월 24일부터 동년 9월 21일까지 10분 간격으로 지하수의 수위, 수온 및 EC를 자동 관측하였으며, 그 결과는 Fig. 5와 같다.
Fig. 5에서 파란색 막대는 일 강수량, 검은색 실선은 수위, 노란색 실선, 주황색 실선, 회색 실선은 각각 8 m, 13 m, 22 m에서의 수온 및 EC의 경시 자료를 나타낸다. Fig. 5(a)의 지하수위는 지표로부터 지하수면까지의 심도를 나타낸다. 그림에서 보는 바와 같이 2022년 8월 중순 이전에는 수위와 EC가 거의 일정하게 유지되었으나, 2022년 8월 중순 이후부터 지하수위가 평시 대비 약 60 cm 상승하고 있으며, EC와 수온에서도 변화가 나타나고 있다. 수위의 급격한 상승에 대한 원인은 파악되지 않았으나, 지속적인 강수로 인해 배후 산지의 갱도 내에 물이 차올라 상류의 수두가 상승하여 관측정에서의 수위가 상승했을 가능성이 있다. 일반적으로 강우로 인해 함양이 증가되면 지하수의 EC가 감소할 것으로 기대되지만, 이와는 반대로 13 m 및 22 m 심도에서의 EC는 증가하는 경향을 보인다(Fig. 5(b) 참조). 특히, 13 m 심도에서의 EC 값은 8 m 심도의 EC와 마찬가지로 초반에는 약 800 μS/cm까지 감소하였다가, 그 이후에는 최대 약 1,400 μS/cm까지 상승하고 있다. 또한 수온의 경우 Fig/ 5(c)에서와 같이, 13 m 심도의 수온은 수위의 갑작스러운 상승과 동일한 시기에 증가하여 22 m 심도의 수온에 도달하고 있다. 이러한 관측 결과를 바탕으로, 8 m 심도에서는 강수 함양으로 인한 오염 물질의 희석으로 인해 EC가 감소하는 반면, 13 m 심도에서는 지하수 유동에 의한 오염 물질의 유입으로 EC가 증가하였다고 해석할 수 있다.
지하수위의 갑작스러운 상승과 EC의 변동이 강수에 의한 것인지 파악하기 위해, 지하수위 및 EC를 일강수량에 대해 교차 상관분석을 수행하였으며, 그 결과는 Fig. 6과 같다. 교차 상관분석에 사용한 자료는 강수의 경우 일강수량 자료를 사용하였고, 지하수의 경우 10분 단위의 측정 자료 중 각 일자의 24시 정각에 해당하는 자료만을 추출하여 사용하였다. 지하수위와 강수의 상관관계를 분석한 결과, Fig. 6(a)와 같이 강수 이후 2일째에 교차 상관계수가 -0.258로 그 절대값이 가장 크며 신뢰구간을 벗어나므로, 둘 사이에 상관관계가 있음을 알 수 있다. 이는 지하수위의 상승이 강수로 인한 것일 수 있고 그 경우, 강수 후 수위 상승까지 약 이틀간의 시차가 있음을 의미한다. Fig. 6(a)에서 상관계수가 음수로 표현된 것은 지하수위를 지표로부터의 심도로 표현했기 때문이다.
강수와 심도별 EC에 대해 교차 상관분석을 수행한 결과, 8 m 심도의 경우(Fig. 6(b) 참조) 시차가 5~6일일 때 가장 높은 음의 상관성을 보인다. 이는 지표 인근 심도임을 고려할 때 강수 유입으로 인한 희석으로 EC가 낮아짐을 의미한다. 13 m 및 22 m 심도의 경우(Fig. 6(c) & 6(d) 참조) 17일에 가장 높은 양의 상관성을 보이고 있는데, 이는 강수로 인해 오염 물질의 유입이 증가하였음을 나타낸다. 또한 시차의 차이를 보이는 이유는 8 m 심도의 경우 관측정 부근 상류에서의 함양에 의한 것이며, 13 m 이하에서는 배후지 상류에서의 함양에 의한 효과를 지시한다.
2.2.3. 지하수위 측정
연구 지역 내에 위치한 지하수 이용 시설 중 조사가 가능한 11개소 및 관측정 1개소에 대해 2022년 4월 21일 하루 동안 지하수위를 측정하고, 지하수 이용 여부를 확인하였다. 지하수위 측정 위치 및 수위 측정 결과는 Fig. 7 및 Table 1과 같고, 연구 지역 내 양수 현황이 없음을 확인하였다. 측정된 지하수위는 모델 보정에 활용하였다.
2.3. 모델링
연구 지역의 지하수 유동 및 오염 물질 이동을 모델링하기 위해 DHI사의 FEFLOW(Diersch, 2014)를 활용하여 수치모델링을 수행하였다. FEFLOW는 3차원 유한요소 모델링 소프트웨어로, 포화대 및 불포화대에서의 지하수 유동, 오염 물질 이동 및 열 전도 등을 모의할 수 있는 도구이다.
2.3.1. 모델 설계
달성광산을 중심으로 반경 약 1 km의 구역을 포함하도록 유역 경계에 해당하는 능선을 따라 Fig. 8(a)와 같이 모델 영역을 설정하였다. 격자는 삼각망을 활용하여 작성하였으며, 하천, 광산, 광미적치장 및 관측 지점과 같은 관심 대상에 대해서는 더욱 세밀하게 격자를 구성하였다. 모델 영역의 측면 경계는 배수 분지의 능선을 따라 설정하였으므로 무흐름 경계로 간주하였다. 영역 내부의 남쪽 하천은 상원지로부터 물이 공급되어 수위가 일정하게 유지되고, 북쪽의 하천 또한 현장 조사시 수위가 일정하게 유지되는 것을 확인하였으므로, 제1형 고정 수두 경계(Dirichlet boundary)를 설정하였다(Fig. 8(b) 참조). 모델 보정에 활용한 지하수위 관측 자료는 갈수기에 해당하며, 하천 수심 또한 깊지 않으므로 하천 수위는 하천 지형 고도와 동일하게 설정하였다. 수직적으로는 총 14개의 층으로 모델을 구성하였는데, 1~4층은 각각 5 m 두께로, 5~13층은 갱도 간의 고도 차이(MIRECO, 2016)를 고려하여 20 m 두께로 설정하였으며, 모델의 바닥 경계는 EL.-60 m의 평면을 가정하였다.
2.3.2. 대수층 수리특성
대구지역 지하수 기초조사 보고서(MOCT, 2004)에 따르면, 신천상류 유역 내 연구 지역의 대수층은 크게 암반층과 충적층으로 구분되며, 해당 유역의 암종별 수리 전도도는 Table 2와 같이 정리할 수 있다. 각 층별 수리 전도도의 평균 값을 모델의 초기 수리 전도도 값으로 설정하였다.
대구지역 지하수 기초조사 보고서(MOCT, 2004)에 수록된 신천상류 유역의 함양률은 14.3%이며, 달성 방재기상관측소의 20년 평균 연강수량은 1,164.18 mm/year이다. 따라서, 166.5 mm/year의 함양량을 모델의 입력값으로 적용하였다.
2.3.3. 모델 보정
모델의 보정을 위해 Table 1의 실측 수위와 수치 모의를 통해 얻은 계산 수위가 근사해질 때까지 대수층의 수리 전도도를 반복적으로 변경하는 시행 착오법을 수행하였다. 최종 결과에 대한 RMSE가 4.53으로 산정되어 보정이 적절히 이루어졌음을 확인하였다. 보정 완료된 최종 수리 수두 분포 및 관측 수두와 계산 수두의 산점도는 Fig. 9(a)와 같고, 보정 완료 후 수리 전도도 분포는 Fig. 9(b)와 같다. 1~4층까지는 0.007342~4.037 m/d, 5~14층은 0.007342 m/d의 단일 값으로 보정되었다.
2.3.4. 입자 역추적 모의
Jung(1996), MIRECO(2016), Kim et al.(2020) 등의 여러 문헌을 살펴보면, 연구 지역 내 많은 위치에서 오염이 보고되고 있다. 따라서, 오염이 보고된 지점 중 하류 2곳을 설정하여 오염이 유래한 경로를 알아보고자 하였다. 첫 번째 지점(P1)은 두 하천 사이의 합류부 부근이며, 두 번째 지점(P2)은 자연정화 처리시설 인근에 위치한다(Fig. 10 참조). 역추적 결과 P1 지점에서 관측된 오염은 광미적치장에서 유래하고, P2 지점에서의 오염은 달성광산에서 유래했을 가능성을 보인다.
2.3.5. 오염 거동 모의
광산 배수 및 광미 침출수가 대수층으로 유입되는 것을 모의하기 위해 흡착, 지연, 분해 등 자연 저감의 영향을 받지 않는 임의의 오염 물질을 가정하여 오염원으로 설정하였다. 광산에 대해서는 광산 면적의 전체 층에 걸쳐서 오염원을 설정하였으며, 광미적치장에 대해서는 지표에 오염원이 위치하는 것으로 설정하였다. 각 해당 오염원의 위치는 Fig. 11의 파란색 원과 같다. MIRECO(2016)에 따르면 EC와 오염물 농도 사이의 선형적인 상관관계가 존재한다. 하지만, 절대적인 오염원 농도에 대한 정보가 부재하여 본 연구에서는 MIRECO(2016)에서 보고된 달성광산 지역 내 최대 EC 측정값인 2,400 μS/cm를 2,400 mg/L로 가정하여 적용하였다. 분산지수는 Zech et al.(2015)와 Schulze-Makuch(2005)를 참조하여 최대 종분산지수 10 m, 최소 종분산지수 1 m, 최대 횡분산지수 1 m, 그리고 최소 횡분산지수 0.01 m로 설정하였다. 1974년 달성광산이 휴광한 이후 오염 유출이 시작되었다고 가정하여, 총 50년이 경과한 시점에서의 오염 분포를 모의하였다.
Fig. 12는 오염 거동 모의 결과를 지표 및 지표하 20 m 심도에 대해 평면적으로 나타낸 그림이다. 지표 인근(Fig. 12(a))에서는 광미적치장에서 유래한 유선은 대다수가 남쪽 하천의 하천 합류부 인근으로 향하고, 광미적치장 인근에서 상대적으로 오염 물질이 고농도로 모의되었다. 지하 20 m(Fig. 12(b))에서의 일부 유선은 광미적치장을 지나는 경우 지표 인근의 유선과 유사하게 연장되어 남쪽 하천의 합류부로 향하고 있고, 그 외의 대다수 유선은 자연정화 처리시설을 경유하여 남쪽 하천으로 향하고 있다. 오염 물질의 분포 또한 광산 주변, 광미적치장 하류, 자연정화 처리시설 주변에서 상대적으로 주변보다 높은 농도를 보인다.
Fig. 13은 관측정 인근에서의 용질 이송·확산 모의 결과를 상세히 표현한 그림이다. Fig. 13(a)는 지표로부터 20 m 아래(5층) 관측정 인근에서의 용질 이송·확산 모의 결과를 나타낸 것이다. 오염은 동쪽에서 서남쪽으로 자연정화 처리시설과 그 남쪽의 관측정(BH) 사이를 가로질러 분포한다. Fig. 13(b)는 BH와 MLSW 두 관정을 지나는 Fig. 13(a)의 선분 A-B를 따른 단면도로, 단면 좌하단의 지표하 20 m 부근에서 오염 물질의 농도가 더욱 높게 분포하는 것으로 모의되었다. 이러한 모의 결과를 고려할 때, BH와 자연정화 처리시설 인근 지표하 20 m 부근에서 오염 물질의 주요 이동 경로일 가능성이 높다고 할 수 있다.

본 연구에서는 3차원 수리지질 모델을 활용하여 달성광산 지역에서의 오염 경로를 도출하고, 이를 현장 관측 결과와 비교하였다. 지하수 유동 및 오염 이송․확산에 따른 오염 물질 분포를 모의한 결과, 상류로부터 오염 물질의 주요 이동 경로가 관측 결과와 유사한 것을 확인할 수 있었다. 즉, 모델링 결과 지표하 20 m 부근에서 높은 농도로 모의되어 해당 심도가 주요 오염 이동 경로일 가능성이 높은데, 이는 저비저항대 심도가 20 m로 관측된 전기비저항탐사 결과와 유사하며, 뿐만 아니라 가장 높은 EC 값을 나타내는 구간이 심도 22 m인 지하수 장기 관측 결과와도 유사하다. 또한 연구 부지 하류 지역에서의 역방향 입자 추적 모의 결과 다양한 연구에서 오염원으로 판명된 광미적치장 및 달성광산이 수리지질학적 측면에서 해당 부지의 상류 함양지에 해당함을 확인할 수 있었다.
본 연구에서 발생한 다양한 불확실성 중 다음의 2가지를 검토하고자 한다. 첫째, 현장 조사 결과 확인된 오염 경로의 심도와 모델링을 통해 추정한 오염 경로의 심도가 정확하게 일치하지 않는데, 이는 모델링의 측면에서 보자면 부지에 대한 3차원 수리지질 모델을 구축할 때 지질의 수직적인 분포를 단순화하였기 때문일 수 있다. 반면, 현장 조사의 측면에서 보자면 전기비저항탐사의 가탐 심도 제한과 관측공 설치 심도의 제한 등으로 인해 모델링 결과에 대응하여 비교하기에는 충분한 심도까지 확인하지 못하였기 때문일 수도 있다. 둘째, 용질 이동 모델링에서 오염 물질 농도 분포와 이동 시간은 대수층의 수리 전도도, 유효공극률 및 분산지수의 영향을 크게 받는데, 기존 문헌에 수록된 값들은 현장 시험 위치가 본 연구의 모의 영역과 일치하지 않는다. 이러한 자료의 한계를 반영하여 단순화한 모델을 활용하였으므로 모의 결과를 특정한 현장 시료의 분석 결과에 대응하여 비교하는 것은 적절하지 못하다. 따라서 본 연구에서는 절대적인 비교보다는 각 조사 결과 자체에서 확인되는 상대적인 농도 차이의 분포를 비교하는 것으로 결과의 적정성을 판단하였다.
Moore and Doherty(2005)에 따르면 불확실성은 크게 2가지로 구분할 수 있는데, 첫째는 자료가 많을 경우 자료 자체에 내재된 오류의 누적에 따른 불확실성이며, 둘째는 개념화 과정 중 지나친 간략화로 인해서 보정 과정에서 현장 조사자료를 바탕으로 반영되었어야 할 불균질성이 무시되어 발생하는 불확실성이다. 전자는 자료상 요인의 불확실성이라고 하며, 후자는 개념적 요인의 불확실성이라고 한다. 이러한 측면에서 볼 때 상술한 오염 심도의 불확실성과 농도분포 및 이동시간의 불확실성은 모두 기존 문헌자료, 즉 관측 지점의 공간적 분포도가 성긴 자료만을 바탕으로 모델링을 수행하였기에 발생한 개념적 요인의 불확실성에 해당한다. 결론적으로 본 연구에서의 3차원 수리지질 모델링 결과의 불확실성을 낮추기 위해서는 모델 규모에 적합한 수준으로 현장에서의 수리지질 조사가 선행되어야 한다. 모델의 규모에 맞게 적절한 공간적 분포를 얻을 수 있도록 모델 영역 내에서 여러 현장조사 지점을 선정하고 시험을 수행하여 그 결과를 모델에 적용한다면, 모델의 복잡성을 적절한 수준으로 높이고, 이를 통해 결과의 불확실성을 낮출 수 있을 것이다.
또한, 기존 문헌 정보만을 바탕으로 구성한 단순한 모델임에도 현장 조사 결과 및 기존 수리지구화학적 연구들의 결과와 상당히 일치하는 결과가 도출되는 것으로 볼 때, 3차원 수리지질 모델이 현장 조사 자료가 부족한 상황에서도 지하수 유동에 대한 유용한 단서를 제공할 수 있음을 시사한다. 알려진 자료를 활용한 개략적인 모델링을 통해 부지의 상황을 이해하고, 지하수 유동 경로나 수리지질학적으로 불확실성이 높은 구간을 파악한다면, 오염 조사 및 관리의 측면에서 초기 비용과 시간을 감소시키는 데에도 기여할 수 있을 것이다.

본 연구는 2024년도 정부(산업통상자원부)의 재원으로 사용후핵연료관리핵심기술개발사업단 및 한국에너지기술평가원의 지원을 받아 수행된 연구입니다(과제번호: 2021040101003B).