1. 서론

최근 수자원으로서 지하수 사용에 대한 관심이 높아지고 있으며, 이로 인해 대수층 내 지하수위 변화에 대한 많은 연구들이 이루어지고 있다. 지하수위는 강우나 조석, 증발산과 같은 자연적인 현상과 상하수도의 누수, 관개용수, 지하구조물과 같은 인위적인 요인으로 인해 변동이 발생하며 그 중 강우에 의한 지하수 함양이 지하수위 변동의 주요 요인으로 작용한다. 강우로 인한 지하수위 변동은 대수층의 수리전도도, 비저류계수, 대수층의 형태 등 수리지질학적 요인의 영향을 받으며 지표의 투수성은 지표 유출(Surface runoff)과 침투율(Infiltration rate)에 따라 강우의 직접적인 지하수 함양을 결정하는 요인으로 작용한다(Lee and Lee, 2001). 하지만 최근 도시화로 인해 지표의 불투수층 확대, 식생의 감소가 발생하며 토지피복형태가 불투수층으로 변화함에 따라 지표의 투수성이 감소하여 강우의 직접적인 지하수 함양이 저해되고 있다. 이러한 지하수 함양량의 감소는 지하수위의 강하뿐만 아니라 지표유출의 속도를 증가시키며 하천의 한계유출량에 도달하는 시간을 크게 감소시키고 하천의 배출능력을 초과하여 도시홍수와 같은 수자원문제를 발생시킨다(Lee and Koo, 2007). 따라서, 토지피복형태에 따라 강수량 및 지하수위 간 상관관계를 규명하는 연구가 필요한 실정이다.

강수량과 지하수위 간 직접적인 상관관계 분석은 강수가 발생한 시간과 강수가 지하수위 변동에 영향을 주는 시간의 차이인 지연시간을 통해 분석할 수 있다. 지연시간의 산정을 위해 사용되는 시계열 분석법은 대표적으로 자기상관 및 교차상관 분석, 이동평균분석 등이 있으며, 시계열 분석과 함께 지역 특성, 관정의 심도, 한계침투량 등 다양한 인자들을 고려한 연구들이 진행되어 왔다(Lee and Lee, 2002; Yang et al., 2006; Ha et al., 2006; Yang and Kim, 2011). 하지만 시계열 분석법을 이용한 강수량과 지하수위의 상관관계 연구들 중에서도 토지 피복형태를 고려한 연구는 부족한 실정이다. 따라서, 본 연구에서는 교차상관분석 및 이동평균분석을 통해 강수와 지하수위 간 지연시간을 산정하고, Arc-GIS를 이용하여 토지피복형태 및 지형을 분석함으로써 토지피복형태에 따른 지하수위와 강수량의 직접적인 상관관계를 분석하고자 하였다.

2. 연구 방법

1) 연구지역 및 자료수집

연구지역인 연제구는 부산광역시 중심부에 위치하며 2006년부터 15년간 총 면적 대비 공장용지가 62.7% 증가하고 전답지와 임야용지가 각각 40.1%, 6.3% 감소하며 많은 도시화가 이루어진 지역이다(Korea Land and Housing Corporation, 2019). 연평균 강수량은 1567.1 mm로 우리나라 연평균 강수량(1289.8 mm; Korea Meteorological Administration, https://data.kma.go.kr/cmmn/main.do) 보다 높으며, 지하수위는 전체적으로 30-100 m의 해발고도 범위를 가지고 도심을 중심으로 30-50 m 내외의 유사한 지하수 등수위선을 보인다(National GroundwaterInformation Center, http://gims.go.kr).

지연시간 산정을 위해 국가지하수정보센터에서 제공하는 부산광역시 연제구 국가지하수관측망 지하수위 자료와 기상청 강수량 자료를 수집하였다. 지하수위 자료는 연속적인 시계열 분석을 위해 현재 연제구에서 운영 중인 13개의 보조관측망 자료(GW01–GW13) 중 2018년 5월부터 2020년 12월까지 결측 값이 상대적으로 적은 5개의 보조관측망 자료(GW01–GW05)를 선정하였다. 연구지역에 설치된 관측공은 모두 충적대수층이며 강우에 의한 수직적인 지하수 함양과 지하수 유동에 의한 수평적인 지하수 함양이 이루어지는 자유면 대수층이라 가정하고 분석을 진행하였다. 공간분석을 위한 기초자료는 환경부 환경공간정보서비스(Environmental Geographic Information Service, EGIS)에서 제공하는 격자 1 m 공간해상도를 가진 축적 1:5,000의 세분류 토지 피복도(2020) 자료와 국토교통부에서 제공하는 1:5,000의 수치지도(2020) 자료를 사용하였다.

2) Arc-GIS를 이용한 공간분석

토지피복에 따른 강수량과 지하수위의 상관관계를 알아보기 위해 Arc-GIS v10.5를 사용하여 공간분석을 진행하였다. Arc-GIS의 Buffer와 Clip 기능을 이용해 관측공 반경 250 m 지역을 구분하였으며, 범위 내에서 토지피복지도 대분류항목별 분류기준에 따른 7개의 항목(시가화/건조지역, 농업지역, 산림지역, 초지, 습지, 나지, 수역)으로 분류하였다. Calculate Geometry 기능을 사용하여 각 항목에 대한 토지피복별 면적을 계산하였으며, 7개의 항목 중 도시화의 영향이 크고 불투수층으로 이루어진 시가화/건조지역(Urban)을 도심영향으로, 그 외 도시화의 영향이 적고 투수층으로 이루어진 농업(Agriculture), 산림(Forest), 초지(Grass_land), 습지(Wet_land), 나지(Bare_land)의 경우 식생영향으로 결정하였다. 또한, 토지경사도가 강수량과 지하수위의 상관관계에 미치는 영향을 알아보기 위해 Slope 기능을 이용하여 연구지역의 경사도를 계산한 뒤 Zonal statistics 기능을 이용하여 관측공 반경 250 m의 토지경사도의 최대값, 평균, 표준편차를 계산하였다.

3) 자기상관 및 교차상관분석

자기상관함수(Auto-correlation function)는 주어진 시차에 대한 자기 값의 선형성과 기억효과를 나타낸다(Angelini, 1997; Larocque et al., 1998). 자기상관함수를 시간 t에 따른 시계열 함수를 y_t라고 가정한다면, 식 (1)으로 표현할 수 있다.

\(\rho(h)=\frac{\operatorname{cov}\left(y_{t}, y_{t-k}\right)}{\sqrt{\operatorname{Var}\left(y_{t}\right) \operatorname{Var}\left(y_{t-h}\right)}}, h=0, \pm 1, \pm 2, \ldots\)       (1)

여기서, ρ(h)는 시차 h에 따른 자기상관함수를 뜻하며, cov(y_t, y_t-k)는 y_t, y_t-h간 공분산(Covariance)을 뜻하며, 분모항은 변수 y_t, y_t-k의 표준편차 간 곱을 의미한다. 자기상관함수를 통해 현재 시점의 결과값과 시차 h 이후의 결과값 간 상관관계를 확인할 수 있으며, 시계열 자료가 0에 수렴할 때 까지의 시간이 길수록 강한 자기상관성과 긴 기억효과를 가진다(LeeandLee, 2002).

지하수위와 강수량에 대해 시차 h까지의 상관관계를 파악하기 위해 교차상관함수(Cross-correlationfunction, CCF)를 사용하였다. 시차 h에 따른 교차상관함수 ρ_{x, y}(h)는 식(2)로 표현할 수 있다.

\(\rho_{x, y}(h)=\frac{\lambda_{x, y}(h)}{\sigma_{x} \sigma_{y}}, h=0, \pm 1, \pm 2, \ldots\)       (2)

\(\lambda_{x, y}(h)=\operatorname{cov}\left(x_{p}, x_{t+h}\right)=\frac{1}{n} \sum_{t=1}^{n-k}\left[\left(x_{t}-\mu_{x}\right)-\left(y_{t+h}-\mu_{y}\right)\right]\)       (3)

여기서, μ_x와 σ_x는 각 시계열자료 x_t의 평균과 표준편차, μ_y와 σ_y는 각 시계열자료 y_t의 평균과 표준편차이다. λ_x,y(h)는 시차 h에 대한 시계열자료 x_t와 y_t의 공분산이다. 위 식을 통해 지하수위와 강수량 간 가장 높은 상관관계를 가지는 시차 h를 구하여 지연시간으로 선정하였다.

4) 이동평균분석

이동평균분석(Moving average analysis, MA)은 t시점에 따른 평균을 이용해 시계열 자료의 추세를 파악하는 분석법으로 주기 h에 대한 이동평균함수 r(h)는 식 (4)로 표현할 수 있다.

\(r(h)=\frac{x_{t}+x_{t+1}+\ldots+x_{t-(k-1)}}{h}\)       (4)

여기서, x_t+ x_t+1+ … + x_t–(k–1)는 주기 h에 대한 x_t+1을 예측하기 위해 사용된 시계열 자료이다.  본 연구에서는 지하수위와 일평균강수량 간 지연시간을 산정하기 위해 일평균강수량의 주기를 0일부터 1000일까지 설정하여 가장 높은 상관계수를 띄는 일수(argmax)를 지연시간으로 선정하였다. 최적의 지연일수를 계산하기위해 피어슨 상관함수를 목적함수로 설정하여 분석하였으며 시차 h에 대한 피어슨 상관함수는 식 (6)으로 표현할 수 있다.

\(L(h)=\operatorname{argmax} r_{X, r(h),} \quad 0 \leq h \leq 1000\)       (5)

\(r_{X, r(h)}=\frac{\operatorname{cov}(X, r(h))}{\sigma_{X} \sigma_{r(h)}}\)       (6)

여기서, L(h)는 피어슨 상관함수의 목적함수이고, argmax는 r_{X, r(h)}의 최대값을 만들기 위한 주기 h를 구하는 함수이다. X는 지하수위 자료, r_{X, r(h)}는 지하수위 자료와 주기 h에 대한 일평균강수량과 이동평균함수의 피어슨 상관함수를 의미한다. 또한, 강수량 자료의 계절적 주기성을 배제한 최적의 시차 h를 구하기 위해 함수 내 첫 번째 극대점을 최적의 지연일수로 계산하여 분석하였다. 이동평균함수를 차분하여 극대점을 계산하였으며 식(7)으로 표현된다.

\(\operatorname{argmin}\left|r_{X, r(h)}^{\prime}\right|=\left|r_{X, r(h)}-r_{X, r(h+1)}\right|, \quad 0 \leq h \leq 1000\)       (7)

여기서, |r′_{X, r(h)}|는 주기 h에 대한 피어슨 상관함수 차분의 절대값을 의미하며, argmin은 | r′_{X, r(h)}|의 최소값을 만들기 위한 주기 h를 구하는 함수이다. 결과값이 가장 작은 주기(argmin) h를 지연일수로 최종 선정하였다.

3. 결과 및 토의

1) 공간분석 결과

관측공 반경 250 m 내 토지피복별 사용면적을 계산한 결과 5개의 관측공 모두 농업지역의 이용은 거의 나타나지 않았으며 연구지역 특성상 대부분 시가화/건조지역이 관찰되었다(Fig. 1). 5개의 관측공 중 GW01, GW02, GW03, GW05는 도심영향이 각각 96, 88, 71, 78%로 불투수지역이 우세하게 관찰되었으며, GW04의 경우 도심영향이 35%인 반면 식생영향이 65%로 투수지역이 우세하게 관찰되었다(Table 1).

Fig. 1. Land cover map of the study area. (A) Line circles represent the effective radius (250 m) of each well at five sites. (B) The composition of land cover within 250 m radius is shown at each site.

Table 1. Land cover areas (m²) within 250 m radius of each well at five sites

연구지역의 토지경사도는 전체적으로 평탄한 모습을 보이지만 산간지역의 영향으로 서쪽 일부과 남쪽, 동쪽지역에서 높은 경사도를 보였다(Fig. 2). 관측공 반경 250 m의 토지경사도 최대값, 평균, 표준편차를 계산하였으며, GW01은 평균 0.96, 표준편차 0.51의 평탄한 경사도를 보였고, GW02, GW03, GW05는 평균 약 7-10, 표준편차는 8-9로 중간정도의 경사도를 보였다(Table 2). GW04는 평균 17.63, 표준편차 10.80의 상대적으로 높은 경사도와 굴곡진 지형을 보였으며 명확한 지형을 관찰하기 위해 Arc-Scene을 통해 각 관측공의 지형을 3D로 표현하였다(Fig. 2)

Fig. 2. Land slope map and 3D landform within 250 m of each well at five sites.

Table 2. Max, mean and standard deviation of land slope values within 250 m radius of each well at five sites

2) 자기상관 및 교차상관분석 결과

지하수위의 자기상관성을 알아보기 위해 5개의 관측공을 대상으로 자기상관분석 결과를 Fig. 3에 제시하였다. 분석 결과 GW01과 GW05는 각각 131일과 146일로 비교적 긴 시차를 가지며 0으로 수렴하였으며, GW02와 GW04는 각각 87일과 58일로 짧은 시차를 보이며 0으로 수렴하였다. GW03의 경우 300일 이내에 0에 수렴하지 않았으며, 자기상관분석 결과가 0에 수렴하는 시차가 길수록 강한 자기상관성을 나타내기 때문에 GW01과 GW03, GW05의 지하수위는 GW02와 GW04에 비해 더 안정적인 지하수위의 변동을 가지는 것을 알 수 있다. 이는 외부적인 요인인 강수의 함양으로 인한 지하수위의 변동이 적다는 것을 의미한다.

Fig. 3. Auto-correlation coefficient (ρ) and lag (day) for groundwater levels at five sites (A: GW01, B: GW02, C: GW03, D: GW04, and E: GW05).

교차상관분석 결과 지하수위와 강수량 간 가장 높은 상관관계를 가지는 최적의 지연시간을 Fig. 4에 제시하였다. GW01과 GW03이 71일과 67일로 비교적 긴 지연시간을 보였고 GW02와 GW05는 35일과 42일로 중간정도의 지연시간을 보였다. GW04는 4개의 관측공에 비해 4일로 가장 짧은 지연시간을 보였는데 이는 GW04의 지하수위의 변화가 강수량과 단기간에 반응한다는 것을 의미하며, 강수의 수직적인 함양이 지하수위에 직접적인 영향을 미친 것을 알 수 있다.

Fig. 4. Cross-correlation coefficient (ρ) and lag-time (day) between precipitation and groundwater levels at five sites (A: GW01, B: GW02, C: GW03, D: GW04, and E: GW05).

3) 이동평균분석 결과

추가적으로 지하수위와 강수량 간 가장 높은 상관관계를 가지는 최적의 지연시간을 산정하기 위해 이동평균분석을 실시하였다. 일평균 강수량을 이용하여 생성한 0일부터 1000일까지 이동평균 강수량을 각 관측공의 지하수위 변동과 비교했을 때 가장 높은 상관계수를 가지는 이동평균 시차를 최적의 지연시간으로 선정하였다(Fig. 5). GW01, GW03, GW05는 각각 153일, 161일, 148일로 긴 지연시간을 보였고, GW02는 94일로 중간정도의 지연시간을 보였으며, GW04는 67일로 가장 짧은 지연시간이 나타났다. 각 관측공에서 최적의 시차를 가지는 이동평균 강수량과 지하수위 변동 그래프를 비교하였을 때 모두 유사한 패턴을 보이며, 이를 통해 지하수위 변동과 이동평균 강수량 간 상관성을 확인할 수 있다(Fig. 6).

Fig. 5. Correlation coefficient (ρ) and lag-time (day) from moving average analysis for groundwater levels at five sites (A: GW01, B: GW02, C: GW03, D: GW04, and E: GW05).

Fig. 6. Trends of groundwater levels and precipitation simulated from moving average analysis from May 2018 to December 2020 (A: GW01, B: GW02, C: GW03, D: GW04, and E: GW05).

4) 토지피복형태 및 지형과 지연시간의 상관관계 분석

각 관측공에서 교차상관분석과 이동평균분석을 통해 산정한 강수량과 지하수위 간 최적의 지연시간과 관측공 반경 250 m 의 토지경사도와 토지피복형태를 정리하였다(Table 3). 분석 결과 GW01, GW03, GW05는 대부분 불투수층의 토지피복형태(Urban: 71%–96%)로 인해 긴 지연시간(CCF: 42–71일, MA: 148–161일)이 나타났으며, 이는 도심영향의 토지피복형태가 강수의 수직적인 지하수 함양에 부정적인 영향을 준 것으로 해석된다. GW04는 관측공 중 유일하게 불투수층보다 투수층의 토지피복형태가 많이 존재하는 관측공(Urban: 35%)으로 다른 관측공에 비해 짧은 지연시간(CCF: 4일, MA: 67일)을 보이며 식생영향의 토지피복형태가 강수의 지하수 함양에 긍정적인 영향을 주는 것을 확인할 수 있다. 이러한 연구결과는 토지피복형태로 인한 지표의 투수성이 강수의 지하수 함양에 영향을 미치기 때문에 강수량과 지하수위 간 지연시간을 결정하는 중요한 요소로 작용하는 것으로 판단된다. GW02의 경우 도심영향의 토지피복형태(Urban: 88%)가 나타났지만 다른 도심영향의 관측공에 비해 상대적으로 짧은 지연시간(CCF: 35일, MA: 94일)이 나타났는데 이는 GW02 관측공 주변에 위치한 도심공원에 의해 긍정적인 영향을 받은 것으로 판단된다.

Table 3. Lag-time (day) calculated using CCF and MA for each groundwater monitoring well at five sites. Land cover composition and mean (± standard deviation) of slope values within 250 m radius at five sites.

관측공 주변의 지형을 알아보기 위해 관측공 반경 250 m 의 토지경사도를 분석하였다. 토지경사도가 급할수록 강수의 지표 유출이 강하게 일어나 강수의 수직적인 지하수 함양이 줄어들고 지연시간이 길어질 것으로 추측되었다. 그러나 도심영향의 GW01, GW03, GW05의 관측공에서 토지경사도가 다양하게 관측되었음에도 토지경사도 차이에 따른 지연시간의 차이(CCF: 29일, MA: 13일)는 크게 나타나지 않았다. 반대로 평탄한 지형을 가진 GW01 (Slope: 0.96 ± 0.51)이 더 높은 경사도를 가진 GW05 (Slope: 6.98 ± 7.80)보다 더 긴 지연시간을 보였으며, 식생영향의 GW04에서 가장 높은 경사도 (Slope: 17.63 ± 10.80)가 나타났지만 토지경사도가 낮은 다른 관측공보다 짧은 지연시간을 보였다. 이는 강수와 지하수위간 지연시간을 결정하는 요인 중 토지피복형태로 인한 영향이 강하게 작용하여 토지경사도의 영향이 미비한 것으로 추정된다. 또한 각 관측공별 토지경사도의 차이는 크지만 경사도의 최대값이 17.63으로 경사가 급하지 않기 때문에 토지경사도에 따른 영향이 크게 나타나지 않는 것으로 판단된다.

4. 결론

본 연구에서는 부산광역시 연제구 지하수관측공 5개를 대상으로 시계열분석을 통해 강수량과 지하수위 간 최적의 지연시간을 산정하고, Arc-GIS를 이용해 토지피복형태를 분석함으로써 토지피복형태가 강수와 지하수위 간 지연시간에 미치는 영향을 평가하였다. 연구결과 토지피복형태가 도심영향이 많은 불투수층으로 이루어진 GW01, GW03, GW05의 경우 지연시간이 길게 나타났으며, 식생영향이 많은 투수층으로 이루어진 GW04는 짧은 지연시간을 보였다. 토지경사도에 따른 지연시간의 차이는 미비하였으며, 이는 토지피복형태가 강수와 지하수위 간 지연시간에 영향을 강하게 미치는 요소이기 때문에 토지경사도에 따른 영향이 작용하지 못한 것으로 판단된다.

본 연구 결과를 통해 부산광역시 연제구의 토지피복형태에 따른 지하수위와 강수량의 지연시간을 규명함으로써 도시 내 강우의 자연적인 지하수 함양을 위한 도시 녹지공간 조성과 수자원 문제 예방을 위한 관리계획의 자료로 활용될 수 있을 것으로 기대된다. 하지만, 본 연구에서 지하수위 관측자료의 결측 값으로 인해 분석된 관측공의 수는 5개로 비교군과 대조군이 적기 때문에 추후 결측값을 보강한 자료를 통해 추가적인 연구가 필요하다. 또한 부산광역시 연제구 외 다른 연구지역을 대상으로 분석을 진행하여 각 지역특성에 알맞은 도시 별 녹지공간 조성과 수자원 문제 예방을 위한 대책 수립에 도움이 될 것으로 사료된다.

사사

이 성과는 정부(과학기술정보통신부)의 재원으로 한 국연구재단의 지원을 받아 수행된 연구임(No. 2021R1 C1C1009177).

이번 연구는 2021년도 교육부의 재원으로 한국기초 과학지원연구원 국가연구시설장비진흥센터의 지원을 받아 수행된 연구임(2021R1A6C101A415).