1. 서 론
2. 포화/불포화 흐름 및 염분 수송 합성 모델링
2.1 합성 모델의 개념
2.2 지배 방정식
2.3 수치모델
3. 지하수 및 토양 염류화 평가 적용 사례
3.1 가상 해안 농경지
3.2 대수층 특성 및 수치모델의 구성
3.3 지하수 및 토양 염류화 결과
4. 결 론
1. 서 론
많은 해안 지역의 농경지는 지하수와 토양수의 염분농도가 높고 지하수면이 높아서 작물의 경작에 어려움을 겪고 있다(Kim, 2014). 지하수를 관개용수로 사용하는 경우 지하수위 염분 농도는 작물의 생산성을 감소시킨다. 지하수를 관개용수로 사용하지 않아도 염 지하수가 모세관 현상으로 상승하여 토양의 염분도를 높이고 그에 따라 작물의 생육에 악영향을 미친다. 염분 농도가 높지 않은 경우에도 지하수위가 높으면 작물의 생육에 역시 부정적 영향을 미친다. 우리나라에는 남서해안을 중심으로 약 1,624여개의 간척지구가 있다(Korea Rural Community Corporation, 1999). 그런데 간척지는 지반 표고가 낮기 때문에 전술한 3개의 문제점들이 더욱 심각하게 대두되고 있는 실정이다. 기후변화에 동반하여 나타나는 해수면 상승은 더욱 문제를 심화시킬 것으로 우려되고 있다.
지난 100년(1906~2005년)간 전 지구의 평균온도는 약 0.74°C 상승하였으며 우리나라의 경우에도 1970년대에 비해 평균기온이 0.7°C 상승한 것으로 알려져 있다. IPCC (Intergovernmental Panel on Climate Change)에 따르면 향후 100년간 지구의 평균온도는 2.8°C 증가할 것으로 예측하고 있으며, 우리나라의 경우에도 약 4°C가 상승할 것으로 예측 하고 있다. 이러한 기후변화는 극단적인 기후 현상의 빈도와 강도는 증가시키고, 단기적인 환경 변화뿐만 아니라 장기간에 걸친 CO2의 변화, 온도 상승, 강수량의 변화, 해수면 상승 등을 유발할 것으로 현재의 과학적 관측 자료를 통해 설명 되어 지고 있다(IPCC, 2007; Oh at al., 2012). 또한 IPCC (2007)는 해안시스템을 기후변화에 가장 취약한 시스템으로 평가하고 있으며, 최빈국 기후변화 적응 지원 분야도 해수면 상승에 의한 해안 지역 적응이 최우선 순위를 차지고 있다. 그러나 해수면 상승에 의한 사회 경제 시스템의 취약성은 개도국은 물론 선진국도 매우 높은 것으로 나타나고 있다. 해수면 상승은 해안침식, 홍수 및 호우피해의 증가, 저지대의 침수, 염수 침입(Saltwater intrusion)에 의한 지하수 오염 등 많은 부작용을 수반한다(Kang et al., 2005; Ferguson and Gleeson, 2012). 우리나라 해수면의 상승률은 1.9 mm/yr (1964~2006)정도로 세계 평균값(약 1.8 mm/yr)을 약간 상회하고 있으며, 특히 제주지역의 경우 매년 약 5.1 mm씩 지난 40년간 약 22 cm가 상승하여 그 증가속도가 매우 빠른 것으로 나타났다(http://www.nori.go.kr). 우리나라 해안별로 상승률을 살펴보면 남해안 3.4 mm/yr, 동해안 1.4 mm/yr, 서해안 1.0 mm/yr, 제주연안 5.1 mm/yr의 상승률을 보이고 있으며 최근 10년 사이 상승속도가 더욱 빨라지고 있다(Park, 2009).
해수면 상승에 대한 지하수 해수침투 연구는 일찍이 시작된 바 있다. Essink et al. (2010)은 MOCDENS3D 모델을 이용하여 네덜란드의 뫼즈쉘트 강의 델타지역에서 지하대수층의 염분화 및 간척지 지표수의 염분화를 이송 확산에 대해 모의하였고, Narayan et al. (2007)은 SUTRA모델을 이용하여 호주 버드킨의 델타지역에서 염수침투를 모의하였다. Werner and Si mmons (2009)은 경계면모델을 이용하여 가상의 해안대수층에서 해수면상승에 따른 염수침투를 모의하였다.
관개로 인하여 진행되는 농경지 불포화 토양의 염류화는 잘 알려진 현상이며, 많은 연구가 수행되었다(Dehaan and Taylor, 2002; Kitamura et al., 2006 등). 그러나 저자들은 해수면 상승으로 발생한 지하수 해수침투로 인한 토양 염류화에 대한 발표된 연구결과를 발견하지 못하였다. 관개로 인한 토양 염류화는 통상 연직 1차원 토양 수분의 흐름과 염분 농도에 대한 해석으로 평가가 가능하지만 해수면 상승으로 인한 지하수 해수침투 및 토양 염류화 평가에는 3차원 포화/불포화 흐름해석이 필요한데 방정식의 비선형성이 커서 고해상도의 시공간 격자가 요구되고, 그에 따른 계산량이 방대하여 현재의 컴퓨터 성능으로도 많은 시간이 걸린다.
농업분야의 기후변화 적응대책 정책 수립을 위해서는 해수면 상승으로 인한 해안 경작지에서 예상되는 피해 범위를 정량적으로 산정하는 것이 필요하다. 본 연구에서는 다수의 해안 경작지에 대한 해수면 상승의 영향을, 과도하지 않은 계산량으로, 평가하기 위한 합성 수치모델링 방법을 적용하였다. 합성 수치모델은 포화대의 지하수 해수침투 평가에는 3차원 모델을, 지하수면 위의 불포화 토양층에 대해서는 연직 1차원 모델을 적용한다. 합성 수치모델을 이용하면 다수의 해안 농경지에 대한 해수면 상승이 지하수와 불포화 토양의 염분도 그리고 불포화 토양의 두께에 미치는 영향을 정량적으로 평가하여 영농 정책 수립의 근거 자료가 될 것으로 기대한다. 개발된 합성모델의 적용성을 제시하기 위하여 가상의 간척 농경지에 적용하였다.
2. 포화/불포화 흐름 및 염분 수송 합성 모델링
2.1 합성 모델의 개념
포화 흐름은 불포화 흐름의 특수한 경우이다. 따라서 3차원 불포화 흐름 모델을 이용하여 포화대와 불포화대의 흐름과 용질의 수송 현상을 모의할 수 있다. 그런데 포화대의 흐름과는 달리 불포화대의 흐름은 비선형성이 커서 포화대의 흐름보다 모의하기가 더욱 까다롭다. 게다가 염분 농도가 지하수 밀도에 큰 영향을 미치는 지하수 해수침투의 경우에는 포화대의 흐름을 모의하는 것도 밀도가 일정한 경우의 포화대 흐름을 모의하는 것보다 훨씬 어려워진다. 따라서 해수침투가 포함된 포화/불포화 흐름 모델링은 지하수 모델링에서 가장 까다로운 분야라고 생각할 수 있다.
본 연구에서는 포화흐름과 불포화흐름을 구분하여 모의하는 합성 모델링 방법으로 지하수 해수침투 현상과 토양 염류화 현상을 모의하였다. Fig. 1(a)에는 가상 해안 지역의 연직단면을 나타냈다. 해안선의 내륙으로부터 담수 지하수가 유출되고 바다쪽에서는 해수가 침투하며 그 사이에는 염분농도가 변하는 천이대가 위치한다. 지하수면 아래의 흐름은 포화흐름으로 통상적인 해수침투 모델링으로 다루어진다. 지하수면과 지표면 사이에 존재하는 불포화대의 흐름은 지표면에서 발생하는 강우의 침투와 지하수면의 조건에 의하여 결정된다. 합성 모델링의 기본 개념은 불포화대의 흐름이 거의 연직방향이라는 점에 기초하며 절차는 다음과 같다. (i) 포화대의 흐름과 염분 수송 현상을 3차원 모델로 모의하고 (ii) 지하수면의 염분농도가 기준치를 초과하는 경우 해당 지점에서 불포화 흐름을 연직 1차원 모델로 모의한다(Fig. 1(b)). (iii) 농경지의 토양 염류화는 토양층의 뿌리심도에서 염분농도가 기준치를 초과하는지 여부로 판단한다.
2.2 지배 방정식
지하수 밀도가 염분의 농도에 따라 변하는 포화/불포화대의 지하수 흐름과 염분 수송 현상에 대한 지배방정식은 흐름방정식과 용질수송방정식으로 구성된다. 흐름방정식은 다음과 같다(Voss, 1984):
(1)
여기서 
은 공극률, 
는 염분 농도(염분질량/지하수질량), 
는 압력, 
는 유체의 점성계수, 
는 유체의 밀도, 
는 중력가속도 벡터 그리고 
는 유체의 Source항이다. 
는 비압력저류계수(specific pressure storage), 
는 포화도(흙속의 간극 부피당 물의 부피 비율)로 지하수염류화 모델에서는 1의 값을 가지며, 토양 염류화 모델에서는 
의 범위를 가진다. 
는 고유투수계수, 
은 상대고유투수계수(relative intrinsic permeability)로 지하수염류화 모델에서는 
이 되며, 토양 염류화 모델에서는 다음의 Van Genuchten (1980)식에 의해 산정된다.
(2)
여기서 
는 토양 물성치이고 
는 유효 포화도이다. 유효 포화도 
는 포화도 
와 다음과 같은 관계식으로 나타난다.
(3)
여기서 
는 잔류포화도, 
는 모세관 압력, 
은 매질의 간극비이다. 유효 포화도 
는 압력의 함수로 다음의 Van Genuchten (1980) 관계로 나타난다.
(4)
여기서, 
는 air entry pressure, 
와 
는 토양 물성치이다.
용질 수송 방정식
염분 확산에 대한 항이 포함된 일반적인 염분 수송 방정식은 다음과 같다.
(5)
여기서 
은 분자확산계수, 
유입/유출되는 지하수의 염분 농도, 
는 단위텐서이며 확산텐서 
는 다음의 일반 좌표계에서 다음 행열로 나타낼 수 있다.
(6)
여기서 
는 대칭행렬이며, 주대각선 요소는 다음과 같다.
(7a)
(7b)
(7c)
여기서 
는 
방향 Darcy 흐름률 벡터의 성분, 
과 
는 각 각 종방향과 횡방향 확산지수이다. 그리고 주대각선 밖의 요소는 다음과 같다.
(7d)
가변 밀도 불포화 흐름에 대한 Darcy 흐름률 벡터는 다음과 같이 정의된다.
(8)
1차원 연직 모델의 경우 지배방정식은 위 식에서 중력 방향 성분식들로 구성된다.
2.3 수치모델
미국 지질조사소(United States Geological Survey, 이하 USGS)에서 개발한 SUTRA (Voss, 1984)는 포화/불포화 흐름 및 염분 수송 방정식을 유한요소법으로 이산화한 수치모델이다. 그러나 포화/불포화 흐름은 염분 농도의 영향을 받는 지하수 밀도와 불포화 흐름 해석에 필요한 각종 관계식에 큰 비선형성이 포함되어 있다. 따라서 방정식의 수치해를 구하기 위해서는 높은 해상도의 시공간 격자가 필요하므로 계산량이 방대하여 현장 규모의 포화/불포화 흐름 모델링을 수행하는 것이 현실적으로 어렵다.
본 연구에서는 포화 지하수 해수침투 해석을 위한 3차원 모델링에 SUTRA를, 토양층의 연직 1차원 모델링을 위하여 미국 환경청(EPA)에서 개발한 VADOFT 모델을 선정하였다(Huyakorn and Buckley, 1987). 합성 모델의 절차는 먼저 3차원 지하수 해수침투 문제에 대한 SUTRA 모델링을 수행하고 그 결과를 분석하여 토양 염류화 평가가 필요한 최상부 절점들에 대한 연직 1차원 불포화 모델링을 VADOFT를 수행하는 순서로 진행된다. 두 모델을 수행하고, 지하수 해수침투, 작물 뿌리 심도에서의 토양 염류화, 불포화층의 두께 등을 평가하는 프로그램을 작성하여 적용하였다.
3. 지하수 및 토양 염류화 평가 적용 사례
3.1 가상 해안 농경지
합성 모델의 시범 적용을 위하여 육지부 497 ha와 담수호 12 ha로 총 면적은 509 ha로 구성된 가상의 간척 농경지를 선정하였다. 간척지의 지표 표고는 0 m (해안선)에서 10 m의 분포를 보이며 평균 표고는 2.28 m이고, 해안선의 길이는 4,250 m이다(Fig. 2).
3.2 대수층 특성 및 수치모델의 구성
본 연구의 가상의 해안 농경지에는 간척지에서 나타나는 일반적인 수리지질 특성을 적용하였다. 지하수 함양율은 227 mm/year (연 평균 강수량 1,356 mm의 16.7%), 담수의 밀도는 1,000 kg/m3, 해수의 밀도는 1,024 kg/m3를 적용하였다. 대수층의 두께는 80 m로, 투수계수는 10.7 m/d, 유효공극률은 0.3, 종분산지수는 10 m, 횡분산지수는 1 m로 가정하였다. 담수호의 관리 수위는 -0.5 mMSL이며, 수질은 담수이다.
KMA (Korea Meteorological Administration, 2012)는 기후변화 시나리오 RCP 4.5와 8.5의 경우에 우리나라 해역에서 예상되는 해수면 상승자료를 Fig. 3과 같이 제시한 바 있다. 1996년 평균 해수위를 기준으로 2100년에는 RCP 4.5 시나리오에서는 72.5 cm, 8.5 시나리오에서는 99.6 cm가 상승하는 것으로 예측되었다(Fig. 3).
대상 지역의 지하수 해수침투 모의를 위한 SUTRA 모델은 절점 개수는 40,436개, 비정형육면체 유한요소는 34,890개, 모델의 대수층은 10개 층으로 구분하였다. 유한요소망은 총 11개의 곡면으로 구성되고 포화 지하수 흐름 격자망의 최상부 곡면에는 3,676개의 절점이 존재한다. 따라서 토양 염류화 평가 대상 후보 절점의 수는 3,676개이다. 절점간의 수평 간격은 40~60 m의 범위이며 층의 두께는 최소 5 m에서 최대 10 m로 불균일하다. 담수 지하수가 유출되는 대수층 상부에 최소 두께의 층을 배치하고 하부 방향으로 층의 두께를 증가시켰다.
해안선에는 해수위와 바다의 염분농도(0.0357 kg/kg)을 경계조건으로 지정하였으며, 내륙의 경계에는 불투수 조건을 지정하였다. 해수면 상승 기간에 대한 부정류 모델링의 초기조건으로 해수위 0 m의 경계조건에 대한 정상상태 결과를 사용하였다. 해수면 상승 시나리오에서는 해수위의 시간 변화를 time-dependent 경계조건으로 적용하였다. RCP 4.5와 8.5 시나리오에서 1560개의 time step을 사용하여 2100년까지 모의하였고, 각각의 time step에서 수렴하기 위한 연산방법은 ORTHOMIN method (Vinsome, 1976)를 사용하였다. 2050년까지는 해수면 상승 경향에 큰 차이가 없어 2050년까지는 시나리오 구분없이 동일한 부정류 해수위 경계조건을 적용하였다. 지하수 함양율을 담수호에 속한 절점을 제외한 모든 절점에 227 mm/yr로 연중 균일하게 적용하였다. 담수호에 속한 절점에는 담수호 수위를 일정 수두로 지정하였다. 지하수 함양량과 담수호 수위는 모의기간 동안 일정하게 유지되는 것으로 가정하였다.
3차원 SUTRA 유한요소망에서 최상부면의 절점 중 염분농도가 기준값을 초과하는 경우 해당 절점과 지표면 사이의 토양층에 대하여 연직 1차원 불포화 흐름 모델링을 수행하였다. 대상 농경지의 불포화 토양층의 두께가 작아 불포화 흐름은 정상상태로 가정하였다. 토양층은 총 40개 요소, 41개 절점으로 이산화하였다. SUTRA에서 결과된 지하수압과 염분농도를 토양층의 바닥 절점에, 지하수 함양율(담수)을 토양층의 상부 절점에 경계조건으로 지정하였다. 토양 염류화 판단을 위하여 작물의 뿌리심도에서 보통 작물의 적정 생육한계 염농도를 사용하였다. Korea Rural Community Corporation (2015a)은 보통작물의 생육에 알맞은 염농도는 뿌리심도(지표하 1 m)에서 2,500 
S/cm (0.0016 kg/kg) 이하이며, 6250 
S/ cm (0.004 kg/kg) 이상에서는 심각한 장해를 받는다고 지적하였다.
불포화대 흐름 모델링에는 표토층의 토양 특성이 고려되어야한다. 대상 농경지의 토양은 Fig. 2(b)의 특성 분포를 가지는 것으로 취급하였으며, 토양통별 특성 변수의 값은 Table 1에 제시하였다. 불포화대의 공극율은 0.45를, 종확산지수는 0.1 m를 동일하게 적용하였다.
3.3 지하수 및 토양 염류화 결과
지하수 해수침투는 밀도의 영향으로 대수층의 하부에서 가장 많이 일어나고, 대수층 상부에서 가장 적게 일어난다. 해수침투 평가에서는 통상 최대 침투 지점이 중요하게 다루어지지만, 본 연구에서는 토양 염류화 평가가 포함되므로 대수층 상부면의 해수침투가 중요하다. 대수층 상부면 지하수의 시기별, 시나리오별 염분농도 분포를 Fig. 4에 도시하였다. 대상 지역의 해수침투에 가장 큰 영향을 미치는 인자는 담수호인 것으로 나타난다. 담수호의 관리 수위가 해수면보다 낮기 때문에 지하수가 담수호로 유출되고 있으며 그에 따라 해수침투가 촉진되고 있다. 초기상태 1996년 대비 해수면 상승에 따라 해수침투가 심화되는 것으로 나타난다. 1996년 지하수염류화 결과를 기준으로 농어촌공사에서 제공하는 해수침투보고서(Korea Rural Community Corporation, 2015b)에 수록되어있는 간척지 중 유사한 지역에서 지하수위 및 해수침투를 비교를 통해, 적용한 매개변수의 타당성을 검토하였다.
3차원 SUTRA 유한요소망에서 최상부면에 속한 절점은 3,606개이다. 초기 상태(1996년)에서는 960개의 절점에서, 2050년에는 1,103개, RCP 4.5 2100년에는 1,228개, RCP 8.5 2100년에는 1,247개의 절점에서 염분농도가 담수보다 높게 나타났으며, 해당 절점에 대하여 1차원 불포화 모델링을 수행하고, 지표면 아래 1 m 지점의 염분 농도를 계산하였다. 불포화 토양의 두께가 1 m 미만인 경우에는 지하수면의 염분농도로 염류화를 판단하였다. Fig. 5에는 뿌리 심도에서 보통작물의 적정 생육한계 0.0016 kg/kg를 초과하는 지역을 도시하였다. 역시 담수호 인접 지역에서 토양 염분농도가 높게 나타나고 있다.
Fig. 6에는 RCP 8.5 시나리오에서 2100년 시점의 불포화 모델링 대상 절점 1,247개 중 한 개 절점(위치는 Fig. 5에 원으로 표시)에 대한 불포화 모델링 결과를 도시하였다. 해당 지점의 지표표고는 3.2 m이며 1996년 상태에서는 MSL 상부 토양층에서 압력수두가 모두 음의 값을 가지며, 포화도도 1 미만으로 불포화 상태이다. 뿌리 심도에서 염분 농도는 0.00110 kg/kg으로 적정 생육이 가능한 염도이다. 해수면이 약 1 m 상승하는 RCP 8.5 시나리오에서 2100년에는 불포화 토양의 두께가 약 0.9 m 감소하였으며 뿌리심도의 염분농도가 0.0032 kg/kg으로 상승하여 적정 생육 기준을 초과하였다. 따라서 고려된 지점은 해수면 상승으로 인하여 토양의 염분농도가 높아져서 영농 피해가 우려된다.
기후변화 시나리오 하에서 해수면 상승이 시범 적용 농경지에 미치는 영향을 Table 2에 제시하였다. 1996년 대비 2100년(RCP 8.5)에 지하수 해수침투 면적은 39.7 ha, 염류화 토양 면적은 29.9 ha, 그리고 불포화층의 두께가 1 m (뿌리심도) 이하인 지역의 면적은 7.8 ha 증가할 것으로 예상된다. 위 결과를 도출하는 데 소요된 컴퓨터 시간은 3차원 지하수 해수침투 모델링에 약 2시간, 연직 1차원 불포화대 모델링에는 수 분 정도가 소요되었다.
4. 결 론
기후변화로 인한 해수면 상승은 해안 농경지의 지하수 해수침투, 토양의 염류화 및 지하수위 상승 피해 우려를 가중시키고 있다. 피해 규모를, 과도하지 않은 수준의 계산량으로, 정량적으로 평가하기 위해서 3차원 포화 지하수 해수침투 모델링과 연직 1차원 불포화 지하수 흐름 모델링을 연계하여 적용하는 합성모델링 방법을 적용하였다. 해수침투 3차원 모델링에는 USGS의 SUTRA 모델을, 연직 1차원 불포화 모델링에는 미 환경청이 개발한 VADOFT 모델을 선정하고, 두 개의 수치모델을 연동시키는 합성 모델을 개발하였다. 개발된 합성모델은 지하수 해수침투 면적, 토양의 염류화 면적 그리고 불포화 토양 두께 감소를 정량적으로 평가할 수 있다. 모델의 적용성을 보이기 위하여 가상의 간척 농경지에 대하여 RCP 4.5와 8.5 기후변화 시나리오에서 2100년까지 예상되는 해수면 상승 자료를 적용하고 농경지의 예상 변화량을 도출하였다. 기준 시점 대비 지하수 해수침투 면적은 39.7 ha, 염류화 토양 면적은 29.9 ha 그리고 불포화층의 두께가 1 m (뿌리심도) 이하인 지역의 면적은 7.8 ha 증가할 것으로 예상되나 해당 농경지의 경우는 높은 지하수위가 영농에 가장 큰 영향을 미치는 것으로 판단된다.
본 연구에서는 2100년까지 기후변화 시나리오 별로 예상되는 해수면 상승 예측 자료를 경계조건으로 취급하여 농경지에 미치는 영향을 평가하였다. 그러나 상승하는 해수위 경계조건에 지하수계가 반응하는 데 상당한 시간이 걸리므로 2100년 시점에 대한 예측 결과에 2100년까지의 해수면 상승량이 완전하게 반영된 것이 아니라는 점은 언급될 필요가 있다. 2100년 시점의 해수위를 일정수위조건으로 적용하고 2200년까지 추가 모델링을 수행하는 경우에도 해수침투 면적 등이 지속적으로 증가하는 것으로 나타나 해수면 상승의 실제 영향은 훨씬 더 오랫동안 지속되고 가중된다는 점을 고려하면 해수면 상승에 대한 적절한 적응 대책이 시급함을 알 수 있다. 본 연구에서 개발된 합성 모델은 우리나라 서남해안에 산재한 많은 농경지에 대한 기후변화 대책 수립에 필요한 근거자료를 산출하는 데 기여할 수 있을 것으로 판단된다.
