1. 서 론

위해성이란 인체나 생태계의 구성요소인 수용체가 오염물질에 노출될 경우 발생할 수 있는 악영향의 가능성으로 정의되며 위해성평가는 이러한 가능성에 대한 정량적으로 측정하는 과학적인 과정이다(Ministry of Environment, 2006). 환경 위해성평가는 인체위해성평가와 생태 위해성 평가로 나눌 수 있으며, 그 중 오염부지에 대한 인체위해성평가는 오염원(source), 이동경로(pathway), 수용체(receptor)의 시스템에서 토양에서 대기로 확산되어 나간오염물질의 흡입에 대한 인체 위해성평가, 토양에서 지하수 및 지표수로 이동한 오염물질의 섭취에 대한 인체 위해성평가 및 토양에서 자란 작물을 사람이 음식물로 섭취하였을 때의 위해성평가를 포함한다(Korean Society of Soil and Groundwater Environment, 2008). 미국 환경보호청(USEPA)과 유럽화학물질국(European Chemical Bureau, ECB)은 이미 오염부지 위해성 평가에 대한 지침을 마련하여 적용하고 있으며(Korean Society of Soil and Groundwater Environment, 2008), 국내에서도 이러한 선진국의 위해성평가 절차를 바탕으로 토양환경보전법제 15조의 5규정에 따라 2006년 처음 ‘토양오염 위해성평가지침’을 제정하였으며(Ministry of Environment, 2006), 그 이후 수차례의 개정을 거처 최근 총석유계탄화수소(TPH)가 위해성평가 대상물질로 추가되면서 2018년 ‘토양오염물질 위해성평가 지침’이 일부 개정되어 현재에 이르고 있다.

USEPA는 오염부지를 정화하기 위한 구체적인 지침으로 미국재료시험협회(American Society for Testing and Material, ASTM)의 RBCA(Risk-Based Corrective Action)를 개발하여 사용하고 있는데, 이는 위해성에 근거한 오염부지 정화방법의 체계적인 틀로 보고 있다(ASTM, 2000). RBCA에서는 노출농도 및 이에 상응하는 위해성기반 선별 기준 및 정화 기준을 설정하는데 이송 및 저감 모델(fate and transport modeling)이 사용될 수 있다고 밝혔다

지중환경 내 오염물질의 이동, 희석 및 저감을 모의하기 위해서는 적절한 수학적 모델을 선택하여야 한다. 일반적으로, 해석 모델은 일반적이고 부수적인 수준의 물질이동 분석에 사용되며, 복잡한 수치모델은 보다 구체적인 결과를 원하거나 토양의 불균질성과 같은 조건 때문에 해석해로부터 만족할 만한 결과를 얻을 수 없거나 해석해가 불가능한 경우에 사용한다(Kim and Park, 2007).

불포화대로 표현하는 토양에서 지하수로의 오염물질 이동경로를 파악할 때는 오염원의 농도, 토양의 공극률, 유기탄소함량, 침투율 등이 주요 입력인자들이며 포화대로 표현하는 지하수 대수층에서는 수리전도도, 수리구배 등을 주요 입력인자로 사용한다. 대부분의 위해성평가에서 사용하는 포화대 및 불포화대의 오염농도 예측모델은 보수적인 결정 방식의 하나인 정상상태에서의 해석식을 적용하는 방식이 주를 이룬다.

많은 선행 연구들 중에서는 수치모의 방식이 포화대, 불포화대 중 일부에만 적용되거나 또는 여러 기작이 생략된 단순화된 개념 모델을 사용해 모의하였다. 한국환경정책평가원에서는 2차원 수치모의모델인 2DFATMIC(2-Dimensional subsurface Flow, fAte and Transport of MIcrobes and Chemicals model using a Lagrangian-Eulerian adapted zooming and peak capturing algorithm)을 활용하여 불포화대 및 가변포화대의 토양 구성별 시나리오에 따라 위해노출농도의 해석저감계수를 평가하였다. Ryu(2010)은 BTEX(benzene, toluene, ethylbenzene, xylenes)와 총석유계탄화수소류(total petroleum hydrocarbons, TPH)를 대상으로 불포화대 오염물질 거동을 해석식으로, 포화대 오염물질 거동을 MODFLOW와 MT3D를 이용한 수치모의로 각각 구현하였다. Verginelli와 Baciocchi(2013)는 평가지침에서 사용하는 해석식이 위해 노출농도를 과다하게 평가하고 있다고 판단하여 오염원 고갈, 흡착, 생분해 기작을 고려한 1차원 이송확산 지배방정식에 기반한 해석식으로 오염물질의 거동을 모의하였다. 최근Mazzieri et al.(2016)은 시간에 따른 노출 농도 변화를상정하고 시간과 깊이에 따른 오염물질 이동경로를 예측하기 위해 불포화대를 MODFLOW와 MT3D를 이용하여 수치모의를 수행하였고, 포화대를 유사해석식(semi- analytical model)을 이용하여 해석하였다. 한편 미국 EPA는 부지특이적 경로를 파악할 수 없는 오염 부지 위해성평가 과정에서 활용할 수 있는 미국의 충적 대수층에 대한 희석저감계수(Dilution and Attenuation Factor, DAF) 대표값 20을 불포화-포화대 수치모델인 EPACMPT 모델 및 몬테카를로 모의 결과로 제시하였는데. 이 과정에서 활용한 부지개념모델은 오염물질 이동경로 상의 생분해 및 흡착 등의 오염물질 특성 및 저감 기작에 대한 고려를 최소화하였다(EPA, 1996).

본 연구는 오염부지 위해성 평가시 활용할 수 있는 포화대 및 불포화대 오염물질 이동 경로 평가를 동일한 수치모의 방식을 적용해 수행하고자 하였다. 이를 위해 본 연구에서는 불포화대에 위치해있는 오염원에서 포화대를 거쳐 노출되는 용질 수송 과정의 모의를 MODFLOW (McDonald and Harbaugh, 1988)와 MT3D(Zheng, 1999)를 이용하여 유한차분 모델을 구축하고 대표적인 오염물질인 벤젠을 대상으로 노출 경로 시나리오를 적용하여 지하수에서의 노출농도를 모의하였다.

이 과정에서 오염부지 위해성평가시 오염물질 노출이동경로 평가를 위한 현장 개념 모델 구축 방식 및 매개 변수 적용 방식을 다양화한 여러 시나리오 적용을 통해 수치 모의의 현장 적용성을 알아보고 각 시나리오별로 다르게 적용한 불포화대 및 포화대 모델 연계 방식이 위해성농도 도출에 미치는 영향을 알아보고자 하였다 마지막으로 기존 위해성 평가에서 사용하는 해석식과 기존 선행연구에서 채택했던 단순화된 개념모델이 수치모델로 구현되는 경우 주의할 점을 알아보고자 하였다.

2. 연구방법

2.1. 개념 모델

Fig. 1은 오염 부지 내 오염물질 이동 및 노출경로를 평가하는데 사용하는 부지개념모델의 수평 수직 방향의 단면이다. 본 모델의 불포화대 모델은 토양 매질 내 흡착/용해/휘발의 평형 상태이며, 정상상태의 수분 침투율과 함께 오염물질이 정속도로 지하수면을 향해 수직방향으로 이동한다고 가정한다. 오염원 자체의 두께는 고려하지 않으며 오염물질 오염경로에서 오염물질의 총질량 감소 및 이에 의한 최고 농도값(peak concentration)의 변화, 그리고 부산물의 생성을 야기할 수 있는 생분해기작은 고려하지 않는다. 관정이 오염원 근처에 가까이 위치하여 오염부지를 벗어난 이후에 관정까지의 거리를 0으로 설정한다. 또한 오염부지에 위치한 오염원은 고갈되지 않아 시간이 지나도 농도가 일정한 것으로 가정한다.

Fig. 1. Conceptual description of the contaminant migration to groundwater

2.2. 해석해(Analytical Solution)

일반적이고 보수적인 위해성평가에 사용하는 노출농도 산정식은 토양에서 용출되는 오염물질이 지하수로의 이동과정에서 전혀 희석 및 저감이 일어나지 않고 그대로 지하수 포화대에 전달됨을 가정한다(EPA, 1996). 이 설정에서는 비포화대 오염물질의 농도는 깊이에 따라 변하지 않는 값이며 포화대 오염물질 농도는 지하수의 혼합대에서 고르게 분산된다는 가정에 따라서 희석되는 시간을 고려하지 않는 정상상태의 해석식으로 계산한다. 이 해석해는 일반적으로 ASTM 모델로 불린다(ASTM, 2000).

포화대 해석식의 도출은 혼합대에서의 물수지에 기초한다. 지하수면에 유입되는 오염물질 농도는 포화대의 담수흐름과 섞이면서 농도가 낮아지는 현상을 아래의 식으로 나타낼 수 있다.

\(\begin{equation} Q_{L} C_{L}=\left(Q_{L}+Q_{w}\right) C_{w} \end{equation}\)       (1)

여기서 Q_L는 지하수 포화대로 오염물질과 함께 유입되는 수직방향 침투 유량 [L³/T], Q_w는 지하수 포화대로 유입되는 수평방향 담수 유량 [L³/T], C_L는 토양유출수 농도[M/L³], C_w는 지하수노출농도 [M/L³]이다.

이 때 지하수 유입 유량은 침투율과 면적의 곱이며, 포화대의 담수 흐름은 단위면적당 유량과 포화대 단면적의 곱으로 나타내진다.

\(\begin{equation} (I L W) C_{L}=(I L W+K i d W) C_{w} \end{equation}\)       (2)

여기서 K는 수리전도도 [L/T], i는 지하수 수리구배 [L/L], I는 침투율 [L/T], d는 오염토양유래 지하수 오염원깊이 [L], L은 지하수흐름방향과 평행한 오염원 길이 [L], W는 지하수 흐름과 직각인 오염부지 길이 [L]이다.

이 때 포화대의 지층 깊이 중에 오염물질 혼합대가 차지하는 깊이는 아래 식 (3)과 같은 경험식으로 계산할 있으며(EPA, 1996) 국내 토양선별지침에서는 d의 값에 대표값 2m를 제시하고 있다(Ministry of Environment, 2006).

\(\begin{equation} d=\left(0.0112 L^{2}\right)^{0.5}+d_{a}\left(1-\exp \left[\frac{-L I}{K i d_{a}}\right]\right) \end{equation}\)       (3)

여기서 d_a는 지하수 포화대의 수직방향 두께 [L]이다.

식 (2)를 정리한 아래 식 (4)은 EPA(1996)의 토양 선별지침에서 대수층의 혼합 및 희석 효과를 평가할 수 있도록 제시하였으며 국내에서도 지하수 노출 농도를 결정하기 위해 현행 토양위해성평가지침에서 목표정화수준의 산정에 이용한다.

\(\begin{equation} C_{w}=C_{L} /\left(1+\frac{K i d}{I L}\right) \end{equation}\)        (4)

이 식은 유도되는 과정에서 지하수 흐름과 직각인 오염부지 너비가 소거되어 Fig. 1과 같이 수평 수직 방향의 2차원 개념 모델로 표현된다.

2.3. 수치모델(Numerical Model)

2.3.1. 지배방정식

불포화와 포화대의 오염물질 이동경로를 균질한 매질의 지하수 포화 유동 방정식 및 용질 이송확산방정식을 바탕으로 계산한다.

3차원의 균질한 다공매질을 통과하는 포화대 대수층 유동의 지배방정식은 다음과 같이 편미분 방정식으로 표현할 수 있으며 2차원 또는 1차원 포화대 유동은 식의 좌변에서 해당 좌표계의 항에 대해서만 표현한다.

\(\begin{equation} \frac{\partial}{\partial x}\left(K_{x x} \frac{\partial h}{\partial x}\right)+\frac{\partial}{\partial y}\left(K_{y y} \frac{\partial h}{\partial y}\right)+\frac{\partial}{\partial x}\left(K_{z z} \frac{\partial h}{\partial z}\right)+W=S_{s} \frac{\partial h}{\partial t} \end{equation}\)       (5)

여기서, h [L]는 수두, K [LT⁻¹]는 수리전도도, S_s [M^{− 1}LT²]는 비저류계수, t [T]는 시간, W는 지하수계로 유입되거나 유출되는 단위체적당 유량, source/sink [T⁻¹], x, y, z는데카르트 좌표로서 x, y는 수평방향 z는 수직방향이다.

용질이동과 거동에 관련된 모델들은 이류, 분산, 확산흡착 등을 예측하는데 사용하며 다음 식과 같이3차원 편미분 방정식으로 표현한다. 주로 이송-확산 방정식에 기반하며 본 연구에서 생분해는 고려하지 않는다.

\(\begin{equation} \begin{aligned} &R \frac{\partial C}{\partial t}=\\ &-V_{x} \frac{\partial C}{\partial x}-V_{y} \frac{\partial C}{\partial y}-V_{z} \frac{\partial C}{\partial z}+\frac{\partial}{\partial x}\left(D_{x} \frac{\partial C}{\partial x}\right)+\frac{\partial}{\partial y}\left(D_{y} \frac{\partial C}{\partial y}\right)+\frac{\partial}{\partial z}\left(D_{z} \frac{\partial C}{\partial z}\right) \end{aligned} \end{equation}\)       (6)

여기서, \(\begin{equation} V_{x}, V_{y}, V_{z} \end{equation} \)[L/T]는 지하수의 평균 선형 유속, C [ML⁻³]는 용질 농도, \(\begin{equation} D_{x}, D_{y}, D_{z} \end{equation}\)[L²/T]는 수리 분산지수를 나타낸다. 이 때 각 격자에서의 평균 선형 유속을 계산하기 위해 단위면적당 유량(specific discharge)을 공극률로 나눈다.

식 (6)에서 사용되는 지연계수는 다음과 같이 식 (7)으로 계산하며 지연계수의 값은 오염물질의 이동 속도가 지연되는 정도에 영향을 끼친다. 본 연구에서는 선형 등온흡착(linear isotherm)을 가정한다.

\(\begin{equation} R=1+\frac{\rho_{b} k_{o c} f_{o c}}{\theta} \end{equation}\)       (7)

여기서, ρ_b[ML⁻³]는 토양용적밀도, k_oc[L³M⁻¹]는 유기물 분배계수, f_oc[dimensionless]는 토양유기탄소비율, θ는공극율이다.

2.3.2. 불포화대 유한차분 모델

Fig. 2은 지하수 불포화대를 모의 영역으로 한 경계 조건을 나타낸 것으로 오염원에서 용존되는 토양유출수가불포화대를 거쳐 지하수로 유입되는 과정을 모의하였다. 오염원이 지하수와 수평방향을 이루는 길이는 50 m이며 오염원에서 지하수면에 이르는 불포화대 깊이는 5m이고포화대 깊이는 10 m이다. 본 연구에서는 토양유출수 C_L가불포화대 경로에 유입되는 시점의 토양유출수 농도를 1mg/L로 설정하고 상대적으로 저감 및 희석되는 농도 분포를 관측한다. 본 연구에서 적용한 불포화대 모델은 1차원 이송분산 모형이기 때문에 유한차분 격자를 짜지 않더라도 이송방정식만을 풀어내는 다양한 방법을 사용하여 지하수 유입농도를 계산해낼 수 있으나 본 연구는 유한차분 수치모의 활용예를 제시하기 위하여 1차원 격자를 구성하여 모델의 결과값을 제시하고자 하였다.

Fig. 2. Computation domain and boundary conditions used in the numerical model describing unsaturated zone.

모델에서 일정한 토양투수도를 가정하고 있으며 토양 내 수분함량에 따른 지하수 수리전도도 변화를 고려하지 않고 있다. 또한 불포화대에서 지속적 오염물질 누출로 인한 지하수 유동의 정상상태를 상정하였다. 따라서 일정한 유속을 가진 이송방정식이 포화수리전도도를 활용하여 모의되는 것을 전제로 하며 이에 따라 포화대 모델을 기반으로 활용한다. US EPA(1996)에서도 정속 토양 투수도에 의한 토양 내 이송만을 고려한 바 있으며 Mazzieri(1996)역시 MODFLOW를 활용하여 불포화대를 모의하였다. 본불포화대 모델 구성은 유체가 등속으로 흘러가는 컬럼 실험과 유사한 모델 구성이다. Mazzieri(1996) 연구와 마찬가지로 본 연구에서 채택한 불포화대의 지하수 등속 유동, 포화수리전도도를 가정하는 방식은 불포화대 영역 내 오염물질 거동을 실제보다 단순화하여 모의하는 방식이다. 모델이 부지특이성을 모사하기 위해서는, 최소한의 입력변수와 여러 적절한 가정(assumptions)들을 내포한 단순한 모델에서부터, 실제 상황을 모사하기 위한 많은 입력변수와 늘려가는 방식으로 연구를 진행할 수 있으며, 본 연구의 불포화대 모델은 단순한 형태의 범용 모델은 불포화대모델의 시작점으로서의 가치가 있다고 판단하였다. 향후 고려할 수 있는 부지특이적 입력변수들은 시변성 강수량, 토양 수분의 변화에 따른 수리전도도 변화, 토양의 pH에 대한 고려, 다양한 생물학적 저감 기작, 불포화대에서 지하수면까지의 계절적 변동 등, 불포화대 모델에서 고려할 수리지질학적 변수 및 화학적 생물학적 변수 등이 포함될 수 있다.

모델의 모의영역을 행방향길이 50 m, 열방향길이 5m로 설정하였다. 층은 단위 길이를 표현하기 위해 1m로 설정하였다. 이 연구의 불포화대 모델에서의 z 방향은 격자 구성 과정에서 중력방향을 의미하기 보다 모델에서 침투한 유량의 흐름 방향으로 정의된다. 5 m 길이의 모델열방향 길이는 0.25 m 크기의 21개의 격자로 구성하였다, 모델 한쪽 경계는 오염원에서 유출되는 토양유출수가 등속도의 토양투수도에 의해 토양 내로 유입되는 현상을 모사하기 위해 침투 유량을 경계조건에 설정하는 Neumann경계 조건을, 나머지 다른 경계는 수두값이 지정되는 Dirichlet 경계조건을 설정하였고 양 옆의 경계는 흐름이 없는 경계로 설정하였다. 이 연구에서는 수두를 지중의 임의의 기준면에서부터 지하수면까지의 높이로 정의하며 지하수위로 표현하기도 한다.

2.3.3. 포화대 유한차분 모델

Fig. 3은 지하수 포화대를 모의 영역으로 한 경계 조건을 나타낸 것으로 지하수로 유입된 오염물질이 포화대에서 이동하면서 지하수 관정에서 인체로 노출되는 과정을 모의하였다.

Fig. 3. Computation domain and boundary conditions used in the numerical model describing saturated zone.

오염원이 지하수와 수평방향을 이루는 길이는 50 m이며 대수층의 깊이는 10 m이다. 모델 상부 지하수 불포화대를 통과한 0.1 m/yr의 단위유량이 오염물질과 함께 모델 안으로 수직방향으로 유입되고 모델 좌측에서는 수리구배 0.0186 m/m에 의해 지하수가 유입된다. 두 방향으로 유입된 지하수는 모델 영역을 통화해 오른쪽 경계로 배출되며 것으로 가정하며 오른쪽 경계에서 시간에 따른 오염물질 농도 변화를 관찰하였다. 모의영역을 행방향 길이 50 m, 층방향길이 10 m로 설정하였다. 이 연구의 포화대모델에서는 z 방향을 중력방향으로 정의한다. 포화대 모델의 오염물질 이동을 1차원으로 모의하는 경우 Fig. 3의 개념모델을 격자크기 ∆x, ∆y, ∆z를 5×1×10m로 분할하여 11개 행과 1개의 열, 1개의 층으로 구축하였다. 포화대 오염물질 이동을 2차원으로 모의하는 경우 5× 1× 1m 크기의 격자크기로 11행과 1개의 열, 10개의 층로구성하였다. 포화대 오염물질 이동을 3차원으로 모의하는 경우 모델을 5×1×1m 크기의 격자크기로 11행, 8개의 열, 10개의 층으로 구성하였다. 좌측과 우측의 경계 조건은 모두 수두값으로 고정하며 각각 10.98 m, 10 m이다.

불포화대와 포화대를 구성하는 토질은 loam이고 침투율는 0.1 m/yr이며 불포화대와 포화대의 수리전도도는 모두 27.12 m/yr로 설정하였다. 수리지질 매개변수값 및 오염물질 특성별 할당값은 Mazzieri et al.(2016) 값을 참고하였으며 그 중에서도 수리전도도 및 침투율값은 실제 부지에서 보여주는 값의 범위 중에서 다소 낮은 편에 속한다. Table 1에서 불포화대 및 포화대 모델에 사용되는 모든 입력변수들의 기본값을 볼 수 있다

Table 1. Summary of numerical model parameters

2.3. 시나리오 구성

본 연구는 7가지 시나리오에서 불포화대와 포화대의 오염물질 이동 및 거동 모델 구성방식을 검토하였는데 가장 단순한 포화대 및 불포화대 모델의 조합부터 점차 모델의 차원 및 매개변수의 개수를 늘려가며 여러 개의 시나리오를 구성하였다.

불포화대에서의 가장 단순한 가정은 오염물질이 희석및 저감, 확산이 일어나지 않은 채 토양 유속과 같은 속도로 수직방향으로 이동하는 것으로 본 연구에서는 이와 같은 설정을 지연침투(delayed infiltration)로 명명하였다. 지연침투는 수치모의 없이 불포화대 침투 유속으로 지연시간을 계산하였다. 또한 불포화대 오염물질을 1차원 유한차분 방식(Finite Different Method, FDM)으로 모의하였으며 설정에 따라 흡착의 유무를 구분하여 모의하였다. 시나리오 상에서 포화대 개념모델은 모두 유한차분 방식으로 모의하였으며 시나리오 설정에 따라 각기 1차원, 2차원 3차원 오염물질 거동으로 모의하고 흡착 기작의 유무를 구분하였다.

Table 2. Description of modeling scenarios

시나리오 1에서 불포화대 모델을 활용한 수치 모의를 수행하지 않으며 지연침투로 포화대로 유입된 오염물질은 1차원 유한차분 방식으로 모의한다. 시나리오 2는 불포화대에서의 지연침투 이후 포화대에 유입된 오염물질을 2차원 유한차분 방식으로 이송 및 확산을 모의한다. ASTM 모델 입력 변수 또는 국내 토양오염지침에서 논의하는 오염물질 및 오염 부지 입력데이터만을 가지고 구축할 경우에는 포화대가 2차원 형태인 것에서 착안한 시나리오이다. 포화대에서 2차원 모델로 이송 확산을 모의한다. 시나리오 3는 불포화대에서 오염물질의 수직방향 이송을 1차원 유한차분 모델로 모의하며 후 포화대에서 2차원 모델로 이송 확산을 모의한다. 시나리오 4는 불포화대에서 오염 물질의 수직방향 이송 및 흡착을 1차원 유한차분 모델로 모의 후 포화대에서 2차원 모델로 이송 확산 및 흡착을 모의한다. 시나리오 5는 불포화대에서의 오염물질 희석 및 저감, 확산이 일어나지 않은 채 토양 유속과 같은 속도로 수직방향으로 이동 후 포화대에서 3차원 모델로 이송 확산을 모의한다. 시나리오 6은 불포화대에서 오염물질의 수직방향 이송 및 흡착을 1차원 유한차분 모델로 모의하며후 포화대에서 3차원 모델로 이송 확산을 모의한다. 이 시나리오는 EPA(1996)가 EPACMPT 모델을 활용해로 미국 전역의 부지특이적 희석저감계수의 대표값 20을 도출할 당시 적용한 설정과 가장 유사하다. ASTM 모델 해석식 또는 국내 토양오염지침에서 논의하는 해석식에 입력하는 오염물질 및 오염 부지 입력하는 자료 외에, 지하수흐름 방향과 교차하는 y 방향의 오염부지 길이를 임의로 설정하게 된다. 시나리오 7는 불포화대에서 오염물질의 수직방향 이송 및 흡착을 1차원 유한차분 모델로 모의하며후 포화대에서 3차원 모델로 이송 확산 및 흡착을 모의한다. 이 시나리오는 EPACMPT모델에서 모의의 편의를 위해 생략했던 흡착 및 생분해 기작 중 흡착 과정을 포함하기 때문에 부지특이적인 오염경로를 평가하는데 근접한 설정이다.

3. 연구결과

3.1. 토양에서 지하수로의 이동(soil to groundwater
pathway)

본 연구에서 불포화대 경우는 5m 깊이의 불포화대를 토양오염물질이 통과하여 지하수위로 도달하는 경로를 모의하였다. Fig. 4는 불포화대 모델을 거쳐 모델경계를 통과하는 오염물질의 농도를 시간에 따른 파과곡선(breakthrough curve, BTC)으로 나타낸 그래프이다. 그래프의 x 축은 오염원에서의 유출 이후의 시간이고 y축은 토양에서 유출된 이후의 농도로서 단위는 mg/L이다. 그래프에서 보여지는 세 가지의 곡선은 i) 토양유출에 의한 불포화대에서의 오염물질의 농도변화가 일어나지 않는 지연 침투(시나리오 1, 2, 5), ii) 수치모의에서 1차원 이송, 분산을 모의한 경우(시나리오 3, 6), iii) 1차원 이송, 분산 및 흡착도 고려하였을 경우(시나리오 4, 7)로 구분하였을 경우의 모의 결과이다.

Fig. 4. Model simulated results for temporal change of contaminant concentration of unsaturated zone for i) delayed infiltration, ii) flow & transport only, iii) flow & transport with sorption.

Fig. 5. Model simulated results for temporal change of contaminant concentration of saturated zone based on scenario 1.

세 가지로 구분한 경우 중에 우선 지연 침투에 의한 오염물질 유입시간의 단순 지체를 가정한 경우는 오염물질이 유속의 흐름과 같은 속도로 불포화대를 통과하는 것을 가정한다. 본 연구에서 채택한 모델에서의 유속은 0.33 m/yr로서 침투율 0.1 m/yr을 공극률로 나누었을 때의값이며 불포화대 5m 깊이를 통과하여 지하수면에 도달시간은 15년이 걸린다. 1차원 이송확산 방정식으로 불포화대 오염물질 경로를 모의했을 경우 상대농도의 중심값인 0.5가 모델 경계에 도착하는 시간은 약 15년으로 앞서모의한 단순 지체와 같고 최고 농도가 경계에 도착하는 시간은 약 30년이 걸린다. 용질농도의 중간값 도달시간은 비교적 분산에 영향을 받지 않는 값으로서 흡착이 없는 경우에는 유속과 같은 속도로 토양오염물질이 이동을 비교하는 지표로 활용할 수 있다. 마지막으로, 오염물질인벤젠의 흡착계수를 입력하여 1차원 이송확산 및 흡착을 모의한 경우 오염물질이 지하수면에 유입까지의 시간이 지연되어 중간농도값이 경계에 도달하기까지의 시간은 약56년 정도이다.

3.2. 지하수 포화대에서의 이동(groundwater pathway)

실측값을 모델 입력값으로 사용하는 경우 지하수위 측정에 의해 수두구배를 결정하는 것이 일반적이다. 모델의양 경계에서 계산한 수리구배 차가 0.0186인 것을 고려하여 모델 좌우 경계를 고정수두로 설정한 모델은 오염물질의 침투율 값이 유의미하지 않을 정도로 작은 값을 띄어수리구배 분포에 큰 영향을 미치지 않는 경우로 한정할 수 있다. 불포화대에서 이동한 토양유출수가 15년 경과 후 포화대로 유입되자마자 거의 바로 포화대 경계에서 오염물질이 관측되기 시작하며 약 50년 이후에는 0.65 mg/L 정도의 농도로 유지된다.

Fig. 6는 시나리오 2의 설정에 의해 오염물질이 불포화대를 통과해 지하수면에 도달한 이후에 포화대로 유입되는 경로를 2차원 유한차분 모델로 모의한 결과로서 시간에 따른 오염운의 농도 분포를 색깔로 표현하였다. 포화대 모델에서의 격자 내 오염물질 농도가 0인 경우는 푸른색으로, 1인 경우는 붉은색으로 나타내었고 2차원 단면을 다수의 수치모의용 지층(numerical layer)으로 구분한 것을 검은색 격자로 표현하고 있다. 토양유출수가 불포화대에 유입되는 시점을 0년으로 하였으며 그림 4는 0년, 20년, 40년, 60년 이후의 포화대에서의 농도분포를 나타내었다. 토양오염물질이 불포화대를 지나는 시간이 15년 이기 때문에, 오염물질 유입은 15년 이후에 일어난다. 20년, 40년, 60년의 포화대 단면에서는 오염운이 지하수 흐름의 방향대로 오른쪽으로 이동하면서 동시에 수직방향으로도 확산되는 것을 볼 수 있다. 60년 이후에는 거의 오염운의 형태가 변화하지 않고 일정한 농도구배를 보이고 있다. 오염운이 유선형의 곡면의 형태로서 오른쪽 상향 경계에서는 유입 농도와 거의 유사한 농도 수준을 보이는 반면수리구배에 의해 지하수가 유입되는 유입되는 좌측과 하단에서는 거의 오염물질이 존재하지 않는다.

Fig. 6. Numerical prediction for the contaminant plume transport based on scenario 2.

Fig. 7은 포화대에서의 수치모의에서 모델 횡단면의 우측 경계에서의 파과 곡선을 검출 깊이별로 구분한 그래프이다. 앞서 보였던 Fig. 6에서 붉은 색 오염운이 있던 지하수면 아래 깊이 1m, 초록색 오염운이 있던 5m, 푸른색 담수 영역으로 보였던 9m 지점에서 관측되는 노출농도의 파과곡선을 각각 그래프로 나타내었을 때 깊이가 깊어질수록 점차 농도가 낮아져 9m에서는 상대농도가 약 0.15로서 현저히 낮은 농도를 유지하는 것을 나타낸다. 

Fig. 7. Model simulated results for temporal change of con- taminant concentration according to a) scenario 2, b) scenario 3, c) scenario 4.

Fig 7(a)에서 1m 지점의 중간 농도값 0.5 mg/L가 지나는 시간은 약 25년 정도이며 토양유출수가 불포화대를 통과해 15년 후에 지하수면에 유입된 후 포화대를 이동해 10년 후에 경계조건으로 배출되는 경로를 따르는 결과를 보여주고 있다.

Fig. 7(b)는 시나리오 3을 바탕으로 2차원 포화대에서의 검출 농도를 검출 깊이별로 구분하여 그래프로 나타낸 것이다. 1차원 불포화대 수치모의 결과로 유입된 시변성 지하수 유입농도로 인해 시나리오 2의 모의 결과 보다 기울기가 완만해졌지만 그 차이는 크지 않으며 중간농도의 도달 시점 및 최고농도 수치가 거의 유사함을 알 수 있다.

Fig. 7(c)은 흡착을 고려하였을 경우의 이송확산 시나리오4에 기반한 모의 결과에서 포화대에서의 검출 농도를 검출 깊이별로 그래프로 나타낸 것이다. 지하수면 아래 깊이 1m에서의 오염물질 파과곡선은 유입되는 농도 분포곡선의 형태와 유사하지만, 모델에 유입되기 시작하는 시간 및 경로를 통과하는 시간이 상대적으로 길어졌다. 파과곡선의 형태 또한 상당히 넓게 분포하고 있는 것을 알수 있다. 지하수면 아래 깊이 5m의 농도는 1에 가까운반면 지하수면 아래 깊이 9m에서 모의 시간(250년)동안 거의 농도가 검출되지 않는다. 흡착계수로 인해 수치모델에서는 경계조건까지 도달하는 속도가 느려지는 효과뿐만 아니라 분산효과도 작아졌기 때문에 오염운 분포의 최고농도에서 최저농도까지의 농도 구배가 급해지는 결과가 도출되었다.

현재까지 시나리오 2,3,4 모의 결과 이송, 확산, 흡착에 의한 각각의 효과는 파과 곡선의 지연이나 기울기의 변화로는 나타내어지나 최고 농도 도달 수준은 토양유출수 농도와 같은 1이었다. 모의 결과는 수직방향의 2차원 이송-확산-흡착 모델에서는 깊이별 농도차이가 존재하나 지하수면 근처에서는 시간이 지나도 희석, 저감되지 않는 농도가 존재하게 된다는 것을 알려준다.

Fig. 8는 시나리오 5의 모의 결과로서 60년 후에 오염운이 평형 상태를 이루고 있을 때의 평면과 단면의 격자별 농도 분포이다. 시나리오 5에서부터 7까지는 오염운이 이동하는 경로에서 모든 방향의 분산을 충분히 구현하기 위하여 3차원 포화대 모델을 모의하였다. 이차원 모델인 시나리오 2,3,4에서는 y 방향 모델 두께를 1m 인 단위넓이를 사용하였었는데 3차원 모델로 변형하면서 y 방향의 두께를 8m로 확장하고 열의 개수를 8개로 증가시켰다. 오염원이 위치해 있는 영역을 4-5열의 함양으로 할당하고 나머지 1-3열, 6-8열에서 오염원이 종방향으로 충분히 이동할 수 있도록 하였다. 모델 모의 시간으로 60년이 경과했을 무렵의 오염운 분포는 그림 9에서와 같이 나타나는 데 Fig. 8(a)에서 오염원이 오른쪽 경계로 이동할수록 오염운이 종방향으로 점차 확장하는 것을 색구분으로 보여주고 있으며, Fig. 8(b)의 단면을 살펴보면 2차원 모의 결과와 비교하였을 때 깊이별 농도 수준이 낮아진 것을 확인할 수 있다.

Fig. 8. Numerical prediction for the contaminant plume after 60 year simulation (a) plan view (b) crosscut of the model.

분산이 오염물질 경로에 미치는 영향을 알아보기 위하여 횡분산지수의 값을 0.01, 0.1, 1 m의 범위 내에서 모의하였다. 분산지수는 오염물질이 지하수와 혼합되는 과정에서 오염운이 주변으로 퍼지면서 농도가 감소하는데 영향을 미치는 정도를 나타내는 값으로서 값이 클수록 높은 혼합을 의미한다. 분산지수는 횡분산지수와 종분산지수로 구분할 수 있으며 본 연구에서는 유속 방향의 분산정도를 고려한 횡분산지수 값을 달리하면서 분산이 오염물질 경로에 미치는 영향을 알아보았다. Fig. 9은 분산의 영향에 따른 포화대 오염물질 농도변화를 파과곡선으로 보여준다. 최고농도 수준이 걸린 시간은 유사하나 분산 정도가 커지면서 최고 농도값이 각각 다르게 나타나면서, 지하수 흐름 방향 이외에 종분산으로 인한 희석 효과가 큰 영향을 나타내는 것을 볼 수 있다.

Fig. 9. Model simulated results for temporal change of contaminant concentration at 1 m depth below of three-dimensional saturated model for i) lower mixing condition with longitudinal dispersity αL=0.01 m⁻¹, ii) default mixing condition with longitudinal dispersity αL=0.1 m⁻¹, iii) higher mixing condition with longitudinal dispersity αL=1 m⁻¹.

Fig. 10에서는 최고위해농도 수준 및 노출지점 도달 시간에 대해 시나리오별 결과를 비교했다. 이 때 노출지점이란 포화대 모델인 MODEL 2의 우측 경계에 해당하는 지점이다. 시간지표인 x축은 시나리오별로 최고 농도의 0.5 mg/L 수준이 도달하는 시간이고 농도지표인 y축은 모의 결과의 최고농도를 나타낸다. 본 연구에서 주어진 매개변수값을 활용하여 지침에서 사용하는 해석식(equation 4)에서 계산한 위해노출농도값은 0.5 mg/L이며 모든 수치모의 시나리오에서는 최고 농도 수준이 0.5 보다 높은 수치를 보여주고 있다.

Fig. 10. Maximum contaminant concentration and time to reach at half of the maximum concentration for all scenarios.

그래프에서는 시간지표 및 농도지표 결과가 시나리오 2와 3, 시나리오 5와 6에서 유사하게 보여진다. 이와 같은 시나리오 사이의 유사한 결과는 동일한 포화대 모의 방식을 유지하는 경우 불포화대의 오염물질 거동 방식의 구분이 최종위해노출농도에 미치는 영향은 크지 않다는 것을 보여준다. 그러나 불포화대에서 흡착을 고려한 이송확산을 모의하는 경우인 시나리오 4 및 시나리오 7은 도달시간이 유의미하게 지연되기 때문에 흡착이 중요한 저감기작이 될 경우에는 불포화대에서 흡착을 고려한 이송확산이 전체 노출이동경로에서 주요한 영향을 주는 기작이 된다. 반대의 경우 흡착 및 자연저감이 불포화대에서 유의미하게 일어나지 않는 환경에서는 불포화대 수치모의를 생략할 수도 있다고 보인다.

최고농도 수준을 비교했을 때 포화대의 모의방식에 따라 결과가 구분되는 것을 알 수 있다. 최고 농도 수준이 오염물질 1차원 이송확산모델(시나리오 1), 3차원 이송확산 모델(시나리오 5,6,7), 2차원 이송확산모델(시나리오 2, 3, 4) 순으로 높아진다. 3차원 모델이 부지 및 오염특성을 반영한 현장특이적인 위해노출농도 계산을 가능하게 하는 모델로서 현장상황을 제일 유사하게 모의한다고 보았을 때, 2차원 이송확산 모델은 부지내 오염물질 종방향 분산을 제한하여 오염물질의 위해성을 과다하게 평가할 수 있는 오류가 있고, 1차원 이송확산 모델은 부지 내 수직방향 오염 구배를 무시한 균질한 오염물질 혼합을 가정함으로써 최고 농도를 과소평가할 수 있는 오류를 내포하고 있음을 그림으로부터 알 수 있다.

이송확산 모델에서 흡착을 고려한 경우(시나리오 4와 7)은 오염물질 도달 시간이 지체되는 정도가 유의미한 것을 볼 수 있다. 위해성평가가 지향하는 현장특이적 평가방식은 오염물질 및 토성의 특성을 모델에 반영하고자 하는 방식으로서 본 연구의 흡착 모의 결과를 살펴보면 시간이 중요한 지표가 된다는 것을 볼 수 있다.

해석식으로 도출한 노출농도 수준은 0.5 m/L이며 본 연구의 시나리오에서 도출한 최고노출농도는 모두 해석식값보다 높았다. 이는 토양오염지침에서 계산하는 해석식방식이 보수적이라고 알려져있지만, 일부 부지 상황에 따라 실제 농도보다 위해성을 과소 평가할 위험이 존재하는 경우가 발생할 수 있다는 것을 보여주는 결과이다.

4. 결 론

본 연구는 오염부지 위해성평가 시 오염물질 노출이동 경로 평가를 위한 수치모델 적용에 관한 연구를 진행하였다, 포화대의 모델을 1차원에서 3차원까지 다양한 차원의 수치모의 결과 해석 및 기존 위해성평가에서 사용하는 산식과의 비교를 통해 본 연구의 결론은 다음과 같다.

1. EPA에서 제시하고 국내 위해성 평가에서 사용하는 해석식 및 경험식은 일반적이고 보수적인 평가 방식으로 알려져 있다. 그러나 유입된 오염물질과 대수층 지하수의 혼합이 잘 일어나지 않는 부지 상황을 모사한 포화대 모델의 수치모의 결과에서는 보수적인 해석식에서 산출된 수치보다도 높은 위해노출농도가 예측됨을 본 연구결과로부터 알 수 있었다.

2. 부지특이적 오염물질 이동경로 예측을 위한 모델 운영 시 주로 수리전도도, 오염물질의 특성 등을 다양화하는 데 주력하지만, 매개변수 설정에 앞서 현실보다 과장되거나 왜곡된 농도 분포 예측을 방지하기 위해서는 오염운의 이동경로를 제한하지 않는 적절한 영역 구성 및 경계 조건 구성이 선행되어야 한다.

향후 본 연구결과는 부지특이적 위해성 평가 기술 개발을 위한 다양한 포화대 모의예제 발굴에 활용될 수 있을 것으로 판단한다. 또한 수많은 가정에 기반하여 간략히 평가하는 불포화대에서의 오염물질 이동 경로의 부지특이적인 평가를 위해 불포화대 모델 구성 및 매개변수 선정에 대한 정밀한 후속 연구가 필요하다.

사사

본 연구는 환경부의 재원으로 한국환경산업기술원의 지 중환경오염위해관리기술개발사업의 지원을 받아 연구되었 습니다(과제번호 2018002450004)