1. 서 론
2. 저류함수법의 이론적 배경
2.1 저류함수법의 분류
2.2 운동파법의 변환에 의한 2항 저류함수법 유도
3. 운동파법과 2항 저류함수법의 비교
3.1 유역의 운동파법과 2항 저류함수법 비교
3.2 하도의 운동파법과 2항 저류함수법 비교
4. 2항 저류함수법의 적용
4.1 적용 대상 지역
4.2 매개변수 추정
4.3 홍수사상 모의
5. 결 론
1. 서 론
저류함수법은 비선형성을 가지는 실제 강우-유출 관계를 해석하기 위한 방법으로 Kimura(木村俊晃; Kimura, 1961)에 제안되었다. 1972~1973년 준비기간을 거쳐 1974년 우리나라 한강홍수통제소에 도입되었다. 이후 홍수예보 대상 유역의 확장에 따라 현재는 4개 홍수통제소 12개 유역의 홍수예측모형으로 이용되고 있다. 한강유역에서 저류함수 매개변수의 정확한 초기치 산정을 위한 최초의 시도는 1973~1977년에 걸쳐 진행되었으며 일본의 경험식을 한강유역에 적합하도록 계수를 조정하는 것이었다(Japan International Cooperation Agency, 1977). 이후에도 한강홍수통제소는 정확한 저류함수 매개변수 추정을 위한 연구를 지속해 왔으며(HRFCO, 1978; 1980; 1982; 1983), 현재 사용되고 있는 저류함수법의 매개변수는 제한된 수문사상 분석에 의해 작성된 2012년도 경험식에 의해 추정된 것이다(HRFCO, 2012).
정확한 강우-유출 해석을 위해서는 정확한 저류함수 매개변수 산정이 매우 중요하다. 저류함수의 매개변수를 산정하는 가장 정확한 방법은 유역의 경우 실측 강우량과 유출량, 하도의 경우 실측 유입량과 유출량을 분석하여 매개변수를 구하는 것이다(HRFCO, 1991; Chung and Park, 2013). 그러나 분할유역 및 분할하도에서 실측자료가 가용한 경우는 그리 많지 않다. 이러한 상황을 극복하기 위하여 실측자료가 가용한 유역 및 하도에 대해서 저류함수 매개변수를 구하고, 구해진 매개변수와 유역 및 하도의 특성인자(경사, 길이 등)를 회귀 분석하여 경험식을 도출하고 이를 이용하여 타 분할유역 및 분할하도의 매개변수 초기치를 추정하게 된다. 경험식의 결정계수가 높지 않은 경우에는 추정된 매개변수의 정확도가 높지 않을 수 있다는 단점이 있다.
경험식을 이용하는 대신에 유역에 대한 1항 저류함수를 운동파의 지배방정식에 근사시켜 유역의 물리적 특성으로부터 매개변수를 추정하는 방법이 지속적으로 시도되고 있다. 이에 대한 최초의 시도는 Kuribayashi and Sadamichi (1969)에 의한 것으로 추정되며 이후 일본에서는 Hashimoto (1978), Nagai (1988), Sugiyama and Kadoya (1988), Sugiyama et al. (1988), Sugiyama et al. (1997), 국내에서는 Nam and Kim (1996), Choi et al. (2010a; 2010b; 2010c) 등의 연구가 있었다. 아쉽게도 국내의 연구는 실무에 적용되지는 못하였다. 하도에 대한 1항 저류함수를 정류 운동파의 지배방정식에 근사시켜 하도의 물리적 특성으로부터 매개변수를 추정하는 방법은 도네가와(利根川, Tone River) 공식(Japan Institute of Country-ology and Engineering, Small and Medium River Planning Study Group, 1999)으로 알려져 있다. 2항 유역 저류함수의 운동파 모형 근사화는 Hoshi and Yamaoka (1982), Katayama et al. (2005b), Gotoh and Saga (2010), Foundation of Hokkaido River Disaster Prevention Research Center (2006) 등의 연구가 있으며, 2항 하도 저류함수의 운동파 방법 근사화는 Miyahara et al. (1999), Foundation of Hokkaido River Disaster Prevention Research Center (2006) 등의 연구가 있다.
일본에서는 2항 저류함수법이 부분적으로 Hokkaido Development Bureau(北海道開発局) 등의 홍수예측 실무에 적용하고 있으나(Katayama et al., 2005a), 우리나라에서는 아직까지 연구가 진행된 바가 없다. 운동파법은 유역 및 하도의 물리변수에 기초한 방법으로 강우-유출 해석에 많이 이용되고 있다. 운동파법의 변환에 의한 2항 저류함수법을 채택하면 물리변수를 기반으로 매개변수를 추정할 수 있어 매개변수의 추정이 용이하며, 유출계산의 정확도가 향상될 수 있을 것으로 기대된다. 본 연구에서는 2항 저류함수법의 운동파 방법 재현성을 가상 유역 및 하도에 대하여 검토한다. 또한 실제 유역 및 하도에 대한 2항 저류함수법의 매개변수를 물리적 인자로부터 도출한 후, 실제 홍수사상에 적용하여 2항 저류함수법의 현업화 가능성에 대하여 검토한다.
2. 저류함수법의 이론적 배경
2.1 저류함수법의 분류
비선형 저류함수법은 Kimura (1961)에 의해서 제안되었다. 저류함수법의 운동량방정식은 유출량과 저류량 관계를 나타내는 저류방정식이며, 이를 연속방정식과 연립하여 유역에서는 강우-유출, 하도에서는 유입량-유출량을 해석하는 방법이다. 저류함수법은 크게 1항과 2항 저류함수법으로 나눌 수 있다. 각각 선형과 비선형 형태가 있으나, 실제 자연현상에서는 유출량과 저류량의 관계가 비선형이기 때문에 본 연구에서는 비선형 방정식만을 다루었다.
2.1.1 1항 저류함수법
1항 저류함수법은 저류방정식에 미분항이 포함되지 않은 특징이 있다. 1항 저류함수법의 지배방정식은 Eqs. (1)~(4)와 같이 표시된다.
유역에 대한 연속방정식 :
| $$ds/dt=r(t-T_l)-q$$ | (1) |
유역에 대한 운동량방정식 :
| $$s=k_1q^{p_1}$$ | (2) |
여기서 s는 저류고[mm], q는 직접 유출고[mm/h], r은 유효우량의 강도[mm/h], k1 및 p1은 각각 저류상수 및 저류지수, Tl은 지체시간[h]이다.
하도에 대한 연속방정식 :
| $$ds/dt=q_0-q_s$$ | (3) |
하도에 대한 운동량방정식 :
| $$s=k_3q_s^{p_3}-T_lq_s$$ | (4) |
여기서 s는 저류량[m³], q0는 상류단 유량[m³/s], qs는 하류단 유량[m³/s], k3 및 p3은 각각 저류상수 및 저류지수, Tl은 지체시간[h]이다.
한강 유역의 경우, 유역에 대한 1항 저류함수법 매개변수의 추정식으로 , p1=0.6, (여기서 n은 유역의 조도계수, i는 경사면 평균 경사[-], A는 유역면적[km²], Au는 ‘시가화/건조지역’의 면적[km²], 는 유효우량의 평균 강우강도[mm/hr])이 제안된 바 있다(HRFCO, 2012). 또한 하도에 대한 1항 저류함수법의 매개변수의 추정식으로 (R2 = 0.560), (R2 = 0.463), (R2 = 0.344) (여기서 L은 하도 길이[km], i는 하도 경사[-])이 제안된 바 있다. 이러한 매개변수 추정식은 타 유역 및 하도의 매개변수를 유역 및 하도 특성의 함수로 추정한 것이기 때문에 정확하다고 보기는 어려울 뿐만 아니라, 추정식의 결정계수도 높지 않다.
2.1.2 2항 저류함수법
2항 저류함수법은 미분항이 포함된 방정식으로 Eqs. (5)~(8)과 같으며, Prasad (1967)의 2항 저류방정식인 를 Hoshi and Yamaoka (1982)가 선형 미분항을 비선형으로 수정한 것이다. 비선형 미분항은 시간에 따른 실제 비선형 저류량을 통제하여 첨두량 및 첨두 발생시간을 실제 현상에 맞추어 주는 장점이 있다.
유역에 대한 연속방정식 :
| $$ds/dt=r-q$$ | (5) |
유역에 대한 운동량방정식 :
| $$s=k_{11}q^{p_1}+k_{12}dq^{p_2}/dt$$ | (6) |
여기서 s는 저류고[mm], q는 직접 유출고[mm/h], r은 유효우량의 강도[mm/h], k11, k12 및 p1, p2은 각각 저류상수 및 저류지수이다.
하도에 대한 연속방정식 :
| $$ds/dt=q_0-q_s$$ | (7) |
하도에 대한 운동량방정식 :
| $$s=k_3q_s^{p_3}+k_4dq_s^{p_4}/dt$$ | (8) |
여기서 s는 저류량[m³], q0는 상류단 유량[m³/s], qs는 하류단 유량[m³/s], k3, k4 및 p3, p4은 각각 저류상수 및 저류지수이다.
2.2 운동파법의 변환에 의한 2항 저류함수법 유도
운동파법은 유역 및 하도의 물리변수에 기초한 방법으로 강우-유출 해석에 많이 이용되고 있다. 운동파법의 변환에 의한 2항 저류함수법을 채택하면 물리변수를 기반으로 매개변수를 추정할 수 있는 장점이 있다. 운동파법에서 2항 저류함수법으로 변환하는 절차는 Foundation of Hokkaido River Disaster Prevention Research Center (2006)의 연구에 나타나 있으며, 본 연구에는 여기에 제시되어 있는 절차를 따랐다.
2.2.1 유역의 2항 저류함수법 유도
Fig. 1에 보인 바와 같이 경사가 일정한 경사면에서 흐름의 운동파 지배방정식은 Eqs. (9)~(11)과 같으며, 적용할 초기조건 및 경계조건은 Eq. (12)와 같다.
연속방정식 :
| $$\frac{\partial h}{\partial t}+\frac{\partial q}{\partial x}=r$$ | (9) |
운동량방정식 :
| $$q=\alpha h^m$$ | (10) |
저류량방정식 :
| $$s=\int_0^Lh(x,t)dx$$ | (11) |
초기조건 및 경계조건 :
| $$h(x,\;0)=0,\;\;h(0,t)=0$$ | (12) |
여기서 h는 수심[m], q는 단위폭당 유출고[m/s], t는 시간[s], x는 흐름 길이에 대한 좌표, r은 단위폭당 유효강우강도[m/s], s는 저류량[m³], L은 경사면의 길이[m]이며, 와 m은 운동파 공식별로 Eq. (13)과 같은 계수를 갖는다.
여기서 i는 유역경사[-], n은 Manning의 등가조고, c는 Chezy 계수, ks는 표층토의 투수계수, 는 표층토의 유효공극이다. 본 연구에서는 Manning 공식을 적용하였다. Manning 공식에서 마찰계수는 과 같이 정의된다.
Eqs. (9)~(12)에 Buckingham II-정리에 의한 무차원 변수 (여기서 이며, 은 평균 유효강우강도[m/s]이다)를 대입하면 Eqs. (14)~(17)과 같은 무차원 운동파방정식이 구해진다.
연속방정식 :
| $$\frac{\partial H}{\partial T}+\frac{\partial Q}{\partial X}=R$$ | (14) |
운동량방정식 :
| $$Q=H^m$$ | (15) |
저류량방정식 :
| $$S=\int_0^1H(X,\;T)dX$$ | (16) |
초기조건 및 경계조건 :
| $$H(X,\;0)=0,\;\;H(0,\;T)=0$$ | (17) |
Eqs. (14)~(17)에 나타나 있는 무차원 운동파방정식의 해는 Eqs. (18) and (19)와 같이 무차원 유역 저류함수법으로 변환될 수 있다.
| $$S=K_{11}Q^{p_1}+K_{12}\frac d{dT}(Q^{p_2})=K_{11}Q^{p_1}+K_{12}p_2Q^{p_2-1}\frac{dQ}{dT}$$ | (18) |
| $$\frac{dS}{dT}=R-Q$$ | (19) |
여기서 K11, K12는 무차원 저류계수이다.
경사면 흐름에 대한 실제 영역에서의 저류함수법은 Eqs. (20) and (21)과 같다.
연속방정식 :
| $$ds_h/dt_h=r_h-q_h$$ | (20) |
여기서 sh는 유역저류량[mm], th는 시간[h], rh는 강우강도[mm/h], qh 단위폭당 유출량[mm/h]이다.
운동량방정식 :
| $$s_h=k_{11}q_h^{p_1}+k_{12}\frac d{dt_h}\left(q_h^{p_2}\right)$$ | (21) |
여기서 k11, k12는 저류상수, p1, p2는 저류지수이다.
경사면에서 Manning의 운동파 공식을 적용하면 저류상수와 저류계수는 Eqs. (22)~(25)와 같이 구해진다.
| $$k_{11}=\frac m{m+1}\left(\frac{10^{3m-6}}{3.6}\right)^{1/m}\left(\frac L\alpha\right)^{1/m}$$ | (22) |
| $$k_{12}=\frac{m^{1/5}}{10}\left(\frac{m+1}m\right)^2k_{11}^2\left(\overline{r_h}\right)^{(2/m-1-p_2)}$$ | (23) |
| $$p_1=1/m$$ | (24) |
| $$p_2={(1/m)}^{3/2}$$ | (25) |
여기서 는 평균 강우강도[mm/h]이다.
또한 Fig. 2에 보인 직사각형 유역 길이와 면적의 관계는 Eq. (26)과 같이 Hack (1957) 공식으로 표현될 수 있다.
| $$L=\varepsilon W,\;\;A=WL_s,\;\;L_s=\eta A^\xi$$ | (26) |
여기서 L은 평균 경사면 폭[km], W는 경사면 폭[km], A는 경사면 면적[km²], Ls은 주흐름의 경사면 폭[km], , 및 는 경사면 형상계수이다. 경사면 길이를[km] 단위계에서[m] 단위계로 변환하면 이 구해진다.
Eq. (26)에 Hack 공식의 일반 계수 = 0.6, = 1.35, = 0.6을 적용하면 Eqs. (27)~(31)의 관계가 구해진다.
| $$k_{11}=2.8235f_cA^{0.24}$$ | (27) |
| $$k_{12}=0.2835k_{11}^2\left(\overline{r_h}\right)^{-0.2648}$$ | (28) |
| $$p_1=0.6$$ | (29) |
| $$p_2=0.4648$$ | (30) |
| $$f_c={(n/\sqrt i)}^{0.6}$$ | (31) |
유역의 무차원 운동파방정식에 의해 구해진 해로부터 구한 무차원 2항 저류함수법의 매개변수는 Fig. 3에 보인 바와 같다. k12는 = 3.0 mm/h를 적용한 것이며 홍수사상에 따라 달라진다.
2.2.2 하도의 2항 저류함수법 유도
Fig. 4와 같은 경사가 일정한 하도에서 흐름의 지배방정식은 Eqs. (32)~(34)와 같으며, 적용할 초기조건 및 경계조건은 Eq. (35)와 같다.
연속방정식 :
| $$\frac{\partial a_s}{\partial t}+\frac{\partial q_s}{\partial x}=0$$ | (32) |
운동량방정식 :
| $$a_s=\alpha q_s^m$$ | (33) |
저류량방정식 :
| $$s_s=\int_0^La_s(x,\;t)dx$$ | (34) |
초기조건 및 경계조건 :
| $$q_s(x,\;0)=0,\;\;q_s(0,\;t)=q_0(t)$$ | (35) |
여기서 as는 흐름 단면적[m²], qs는 유량[m³/s], t는 시간[s], x는 흐름 길이에 대한 좌표[m], ss는 저류량[m³], L은 하도 길이[m]이며, q0는 상류단 유량 [m³/s], 와 m은 하도 상수로 Eq. (36)과 같은 관계를 갖는다.
| $$\alpha={(n/\sqrt i)}^mb^{1-m}$$ | (36) |
여기서 n은 Manning 조도계수, i는 하도 경사[-], m은 상수[-], b는 하폭[m]이다. 장방형 하도 단면에서 m은 3/5의 값을 갖는다.
Eqs. (32)~(35)에 Buckingham II-정리에 의한 무차원 변수 (여기서 이며, 은 상류단 평균유량[m³/s]이다)를 대입하면 Eqs. (37)~(40)과 같은 무차원 운동파방정식이 구해진다.
연속방정식 :
| $$\frac{\partial A_s}{\partial T}+\frac{\partial Q_s}{\partial X}=0$$ | (37) |
운동량방정식 :
| $$A_s=\alpha Q_s^m$$ | (38) |
저류량방정식 :
| $$S_s=\int_0\;S_s(X,\;T)dX$$ | (39) |
초기조건 및 경계조건 :
| $$Q_s(X,\;0)=0,\;\;Q_s(0,\;T)=Q_0(T)$$ | (40) |
실제 영역의 하도 저류함수법은 Eqs. (41) and (42)와 같이 표시된다.
연속방정식 :
| $$\frac{ds_s}{dt}=q_0-q_s$$ | (41) |
운동량방정식 :
| $$s_s=k_3q_s^{p_3}+k_4\frac d{dt}\left(q_s^{p_4}\right)$$ | (42) |
상류단 경계조건은 Fig. 5와 같이 삼각형 유량수문곡선으로 가정하였으며, 기저시간은 tr[h]이며 첨두 발생시간은 ta[h]이다.
Eqs. (37)~(40)의 무차원 운동파방정식의 해는 Eqs. (43) and (44)와 같이 무차원 하도 저류함수법으로 변환될 수 있다.
| $$S_s=K_3Q_s^{p_3}+K_4\frac d{dT}\left(Q_s^{p_4}\right)$$ | (43) |
| $$\frac{dS_s}{dT}=Q_0-Q_s$$ | (44) |
여기서 K3, K4는 무차원 저류계수이다.
하도의 무차원 운동파방정식 의해 구해진 해로부터 구한 무차원 2항 저류함수법의 최적 매개변수는 Fig. 6과 같다.
저류상수 k3와 k4는 Eqs. (45) and (46)과 같이 주어진다.
| $$k_3=K_3\alpha L\left({\overline q}_0\right)^{m-p_3}$$ | (45) |
| $$k_4=K_4\alpha^2L^2\left({\overline q}_0\right)^{2m-p_4-1}$$ | (46) |
3. 운동파법과 2항 저류함수법의 비교
3.1 유역의 운동파법과 2항 저류함수법 비교
Fig. 7에 보인 바와 같이 무차원 시간으로 2시간 동안 1의 무차원 강우가 발생한 경우에 대하여 운동파법과 2항 저류함수법으로 모의하였다. 두 방법에 의한 유량수문곡선 및 유량-저류량 관계곡선은 매우 유사함을 알 수 있으며, 유량-저류량 관계곡선이 고리형태를 보이고 있어 강우-유출의 비선성을 확인할 수 있다.
3.2 하도의 운동파법과 2항 저류함수법 비교
Fig. 8에 보인 바와 같이 상류단 유입수문곡선이 단일 첨두인 경우와 이중 첨두인 경우에 대하여 운동파법과 2항 저류함수법으로 하류단 유입수문곡선을 모의하였다. 하도 길이 L = 19.81 km, 하폭 b = 100 m, 하도 경사 i = 1/200, m = 0.6, Manning 조도계수 n = 0.035을 적용하였다. 이는 가상 조건이며, 무차원화 하였을 때 단순한 숫자로 무차원 유량, 첨두 발생 및 기저시간이 도출되도록 하기 위한 것이다. 두 방법에 의해 모의된 하류단 유량수문곡선이 매우 유사함을 알 수 있다.
4. 2항 저류함수법의 적용
4.1 적용 대상 지역
환경부 한강홍수통제소는 홍수특보 지점을 확대하기 위한 노력을 기울이고 있다. 본 연구의 적용 대상 지역은 이들 지점 중의 하나인 한강 수계 홍천강 상류유역으로 Fig. 9에 보인 바와 같다.
적용 대상 유역은 Fig. 10에 보인 바와 같이 표준유역 101401- 101406으로 구성되며, 분할하도는 4개이며, 유출량 계산 지점은 홍천군(굴운교) 지점으로 수위관측소 코드는 1014060이다. 이는 현업 홍수예측모형에서 사용되는 유출모식도와 동일하다.
4.2 매개변수 추정
유역 저류함수법의 매개변수를 추정하기 위하여 1초(약 30 m) DEM과 환경부 중분류 토지피복지도를 이용하였다. 토지이용 및 토지피복별 등가조도(Hashimoto, 1978)는 Table 1에 보인 바와 같다. 이를 이용하여 격자별로 을 구하여 유역평균조도를 산정하였으며 2.22~2.45의 범위에 있다.
Table 1. Equivalent roughness by land use and land cover
하도의 매개변수는 m = 0.6, = 0.5를 적용하여 무차원 매개변수를 구한 후, 하도의 조도계수 n, 하도의 경사 i와 하폭 b를 적용하여 실제 영역의 매개변수로 변환하였다. 하도 관련 변수는 해당 하천의 하천기본계획 또는 구글어스(Google Earth)를 이용하여 구하였다.
4.3 홍수사상 모의
분석 대상 홍수사상의 기간은 ‘2017/07/02~07/04’와 ‘2018/08/28~08/30’이다. 두 기간에 대한 모의 결과는 Fig. 11에 보인 바와 같다. 1항 저류함법은 홍수통제소 현업 홍수예측모형에 탑재된 것이다. 모의 결과를 살펴보면 두 홍수사상 모두에서 RMSE와 NSE (Nash-Sutcliffe Efficiency)가 개선되었음을 알 수 있다.
5. 결 론
정확한 유출해석을 위해서는 정확한 저류함수의 매개변수의 산정이 중요하나, 현재 사용되고 있는 저류함수법의 매개변수는 제한된 수문사상 분석에 의해 작성된 2012년도 경험식에 의해 추정된 것으로 정확도가 낮은 실정이다. 운동파법은 유역 및 하도의 물리변수에 기초한 방법으로 강우-유출 해석에 많이 이용되고 있다. 본 연구에서는 운동파법을 변환하여 구한 2항 저류함수법의 매개변수의 물리적 추정 방법을 제시하였다. 운동파법을 직접 적용하여 구한 모의결과와 2항 저류함수법을 적용하여 구한 모의결과는 서로 잘 일치함을 알 수 있다. 운동파법을 홍수예측에 적용하기 위해서는 유역 및 하도에 대한 분포형(distributed) 격자망을 작성하여야 하고 계산시간도 오래 소요되는 불편함이 따른다. 그러나 2항 저류함수법을 적용하면 기존의 집체형(lumped) 유출모식도를 그대로 이용하고도 운동파법의 정확도를 유지할 수 있는 장점이 있다. 특히 유역유출에서 중요한 물리변수인 등가조도는 토지이용 및 토지피복을 이용하여 쉽게 구할 수 있으며, 하도유출의 물리변수는 하천기본계획 보고서 및 지형도 등에서 쉽게 구할 수 있어 매개변수의 추정이 매우 용이한 장점이 있다. 현업 홍수예측모형인 1항 저류함수법(경험식에 의한 매개변수 추정)과 비교한 결과에서도 2항 저류함수법이 더 정확한 해를 제시하는 것으로 나타났다. 미계측 유역 및 하도에서도 물리변수를 이용하여 매개변수를 용이하게 추정할 수 있기 때문에 저류함수법의 활용도가 증가될 수 있을 것으로 기대된다.
