Numerical analysis of the hyporheic flow effect on solute transport in surface water

Jun Song Kim

doi:10.3741/JKWRA.2022.55.1.23

Preview

Research Article

Journal of Korea Water Resources Association. 31 January 2022. 23-32
https://doi.org/10.3741/JKWRA.2022.55.1.23

Numerical analysis of the hyporheic flow effect on solute transport in surface water

혼합대 흐름이 지표수 내 용질거동에 미치는 영향 수치모의 분석

Jun Song Kim^a^*

김 준성^a^*

^aAssociate Research Fellow, Environment and Resources Research Division, Korea Research Institute for Human Settlements, Sejong, Korea

^a국토연구원 국토환경·자원연구본부 부연구위원

^{*Corresponding Author}

ABSTRACT

본 연구는 지표수-지하수 연계 2차원 흐름 및 용질거동 수치모의를 통해 혼합대 흐름에 의한 지표수-하상 경계층 용질교환이 지표수 내 용질거동에 미치는 영향을 분석하였다. 먼저, 혼합대 흐름이 부재한 불투수성 하상의 경우 사련 하상파에 의해 발달한 바닥면 근처 재순환 흐름의 포착 효과가 용질 파과곡선 형태를 결정하였다. 반면, 용질이 혼합대로 유입되는 투수성 하상에서는 지표수 대비 유속이 느린 혼합대 흐름에 의해 파과곡선 꼬리 부분 농도가 증가하여 파과곡선 꼬리 연장이 급격히 증가하였다. 또한, 지표수 유속 증가는 바닥면 압력 상승과 함께 혼합대 흐름을 빠르게 하여 파과곡선 꼬리의 연장을 단축시켰다. 이러한 본 연구의 결과는 지표수 내 용질 체류시간 예측에 있어 혼합대 흐름의 중요성을 시사하는 바이다.

Keywords

Hyporheic zone

Solute transport

Ripple bed

Breakthrough curve

Numerical simulation

This paper performs two-dimensional numerical simulation of surface water-groundwater coupled flow and solute transport to investigate the effect of the hyporheic exchange at the sediment-water interface (SWI) on surface solute transport. For the impermeable bed case in the absence of hyporheic flow, the trapping effect of flow recirculation associated with the ripple bed controls the shape of breakthrough curves (BTCs). However, the permeable bed case with hyporheic flow stimulates the extended tailing of the BTCs more significantly due to the elevated concentration of the BTC tailing resulting from slow hyporheic velocity. Also, the increased bottom pressure at the SWI with an increase in surface velocity shortens the BTC tailing because of increasing hyporheic velocity. These results infer that hyporhiec flow is critically important in predicting solute residence times in surface water.

키워드

혼합대

용질거동

사련하상

파과곡선

수치모의

MAIN

1. 서 론
2. 흐름 및 용질거동 모의
2.1 지표수-지하수 연계 흐름해석
2.2 용질거동 해석
3. 수치모의 계산조건 및 검증
3.1 모델 검증
3.2 모델 적용
4. 모의결과 및 고찰
4.1 흐름 모의결과
4.2 용질거동 모의결과
5. 요약 및 결론

1. 서 론

지표수-하상 경계층에서의 흐름교환은 하천, 호소, 연안, 해안 등 다양한 수환경 시스템에서 일반적으로 발견되는 수리·지형학적 특성으로, 흐름교환이 발생하는 경계층 근처 하상층 영역은 혼합대(hyporheic zone)로 정의된다(Findlay, 1995). 혼합대 흐름은 하천의 수질 거동에 있어 큰 영향을 주는 주요 수리 인자로써 용존산소, 질소, 탄소순환 메커니즘을 제어하고, 미생물 서식처를 제공하는 등 수질·생태학 측면에서 하천에 매우 중요한 역할을 제공한다(Boulton et al., 1998). 특히, 혼합대 내에서는 탈질산화(denitrification) 과정을 통해 아산화질소(N₂O)가 생성되기 때문에 혼합대는 하천의 주요한 온실가스 생성 및 배출경로로 알려져 있다(Marzadri et al., 2012). 또한, 혼합대에서 발생하는 지표수와 지하수 간의 흐름교환 특성은 홍수터의 하안식생 분포 변화, 지하수 재충전 및 수위 변동 등에 직접적인 영향을 미치기 때문에 하천복원 및 관리 측면에서 매우 중요하다(Harvey and Gooseff, 2015).

더 나아가, 수질사고로 인해 유해화학물질이 하천에 유입될 경우, Fig. 1과 같이 하천-하상 경계층(sediment-water interface)에서 용질교환이 발생함과 동시에 하류로 이동하기 때문에 지표수 대비 유속이 매우 느린 혼합대의 흐름 특성이 오염물질 체류시간 예측에 지배적인 영향을 미치게 된다(Haggerty et al., 2002). 따라서, 혼합대의 수리·지형학적 특성을 수질모델링에 정확하게 반영하는 것이 매우 중요하다. 하지만, 일반적으로 하천 수질예측에 이용되는 1차원 이송-분산 방정식 기반의 저장대 모델(transient storage model)(Runkel and Chapra, 1993) 등에는 혼합대 효과가 기타 저장대 효과(식생, 사수역, 재순환 흐름 등)와 구분되어 있지 않고, 추적자실험에서 실측된 파과곡선(breakthrough curve)을 이용하여 저장대 면적, 교환계수 등의 매개변수로 현상을 단순화함에 따라 기존 1차원 수질모델은 용질거동에 대한 혼합대 효과를 물리적으로 해석하기에 한계가 있다(Kim et al., 2020).

혼합대의 흐름 강도 및 구조는 지표수의 유속 및 수심, 하상의 표면 형태 및 투수율(permeability) 등 다양한 하천의 수리·지형학적 인자에 영향을 받는 것으로 알려져 있다(Norman and Cardenas, 2014). 과거 연구들은 하상 표면의 기하학적 구조가 혼합대 흐름 기작에 지배적인 영향을 미침을 실험적으로 밝힌 바 있으며(Packman et al., 2004; Janssen et al., 2012), 해당 연구들은 주로 하천, 해양, 호소 등 자연수계 바닥에서 일반적으로 발견되는 사련하상(ripple bed)을 대상으로 혼합대 흐름 특성을 규명하고자 하였다(Ribberink and Al-Salem, 1994). 사련하상 위의 난류 흐름은 바닥면 근처에 흐름분리(flow separation) 현상을 발생시키는데, 여기서 베르누이(Bernoulli) 효과에 의해 사련 하상파의 고점(crest)를 중심으로 좌측과 우측에 각각 고압력, 저압력 영역이 형성되어 하상 표면을 따라 Fig. 1(a)에 나타낸 바와 같이 국부적인 압력경사가 발생하게 된다(Packman et al., 2004). 이러한 바닥면 압력경사는 하천과 하상 간의 흐름교환을 증진시킴과 동시에 Fig. 1(b)와 같은 혼합대 고유의 2차원적 흐름 구조를 형성하게 되며, 결과적으로 하천-하상 경계층 용질교환 메커니즘을 설명하는 혼합대 플럭스(hyporheic flux)에 큰 영향을 미치게 된다. 위에서 언급한 바와 같이 혼합대는 하천의 온실가스 주요 배출경로로서 Ahn et al. (2019)는 하상의 형태가 N₂O 배출에 큰 영향을 미침을 실험적으로 확인하였는데, 특히 언덕과 비탈로 인해 강한 압력경사와 함께 혼합대 흐름이 증대되는 사구 고점에서 N₂O 배출량이 평탄한 하상 대비 약 3배 가까이 높게 나타났다.

최근까지 지표수-지하수 연계 2차원 수치모델링을 통해 지표수 및 하상과 연관된 수리·지형학적 인자를 변수로 하여 혼합대 모멘텀 플럭스(momentum flux)의 민감도를 분석하고자 하는 연구들이 주를 이루었고, 혼합대 모멘텀 플럭스가 지표수의 유속 변화에 가장 민감한 것으로 나타났다(Cardenas and Wilson, 2007; Ren et al., 2019; Lee et al., 2020; Ren and Zhao, 2020). 그러나, 해당 연구들은 흐름해석 결과를 토대로 혼합대 모멘텀 플럭스만을 계산하였거나, 지표수 영역의 흐름을 고려하지 않고 혼합대 내부에서만 용질거동 모의를 수행함에 따라 하천과 같이 난류흐름을 지닌 지표수에서 용질거동에 미치는 혼합대의 영향 분석이 미흡하다. 따라서, 혼합대 흐름과 지표수 내 용질 체류시간 분포의 인과관계를 물리적으로 규명하기 위해서는 지표수 흐름과 혼합대 흐름을 동시에 고려한 용질거동 해석이 필요하다.

본 연구의 목적은 지표수와 지하수 흐름을 연계한 수치해석 방법을 통해 혼합대가 지표수 용질거동에 미치는 영향을 분석하는 것이다. 먼저 2차원 Reynolds 평균 Navier-Stokes (RANS) 방정식과 Darcy 방정식을 연계하여 지표수와 혼합대(지하수) 내 유동을 모의하고, 이후 연계모의를 통해 얻어진 흐름장과 함께 혼합대 흐름 효과를 고려하여 지표수-하상 경계층에서 2차원 용질교환 및 거동을 해석하고자 한다. 여기서, 하상의 형태는 자연계에서 일반적으로 관측되는 사련하상으로 설정하였고, 과거 연구에서 혼합대 플럭스에 가장 큰 영향을 주는 것으로 나타난 지표수 유속을 변수로 하여 현상을 분석하였다. 더 나아가, 용질이 하상 내로 유입되는 경우와 유입되지 않는 경우를 비교 모의하여 용질 체류시간 분포, 즉 파과곡선 형태에 대한 혼합대 흐름의 기여도를 정량화하는 것을 본 연구의 주요 목적으로 한다.

https://cdn.apub.kr/journalsite/sites/kwra/2022-055-01/N0200550103/images/kwra_55_01_03_F1.jpg

Fig. 1.

Conceptual diagrams of (a) bedform-induced hyporheic exchange and resulting tracer transport and (b) key hydrodynamic processes controlling interfacial mass exchange

2. 흐름 및 용질거동 모의

2.1 지표수-지하수 연계 흐름해석

Cardenas and Wilson (2007)은 2차원 RANS 방정식과 Darcy 방정식을 연계한 혼합대 흐름 해석방법을 제안한 바 있는데, 해당 방법에서는 Fig. 2에 나타낸 바와 같이 RANS 기반의 지표수 흐름 모의를 통해 계산된 하상 압력분포를 경계조건으로 하여 하상 형상에 따른 바닥 압력경사를 고려한 하상 내 혼합대 흐름을 해석하게 된다. 해당 방법은 개수로의 난류흐름 특성과 하상 표면 형태를 동시에 고려하여 혼합대 흐름을 해석하는데 가장 널리 활용되고 있는 수치해석 방법이다(Janssen et al., 2012; Chen et al., 2015; Jin et al., 2020; Lee et al., 2020). 따라서, 본 연구에서는 Cardenas and Wilson (2007)이 제안한 혼합대 흐름해석 방법을 채택하였고, 먼저 지표수의 2차원 난류흐름은 아래의 연속방정식과 RANS 운동량방정식을 이용하여 계산하였다.

https://cdn.apub.kr/journalsite/sites/kwra/2022-055-01/N0200550103/images/kwra_55_01_03_F2.jpg

Fig. 2.

An experimental flume of Janssen et al. (2012) and boundary conditions used for coupled open channel flow-groundwater flow simulation

(1)

\frac{\partial U_{i}}{\partial x_{i}} = 0

(2)

\frac{\partial U_{i}}{\partial t} + U_{j} \frac{\partial U_{i}}{\partial x_{j}} = - \frac{1}{ρ} \frac{\partial p}{\partial x_{i}} + \frac{\partial}{\partial x_{j}} [ν (\frac{\partial U_{i}}{\partial x_{j}} + \frac{\partial U_{j}}{\partial x_{i}})] - \frac{\partial}{\partial x_{j}} u_{i} u_{j}

여기서 U_i와 U_j는 각각 x_i, x_j 방향 시간평균 유속, p는 압력, $ν$ 는 동점성계수, $u_{i} u_{j}$ 는 레이놀즈 응력항으로 다음과 같이 정의된다.

(3)

- u_{i} u_{j} = ν_{t} (\frac{\partial U_{i}}{\partial x_{j}} + \frac{\partial U_{j}}{\partial x_{i}}) + \frac{2}{3} k δ_{i j}

여기서 $v_{t}$ 는 난류점성계수, k는 난류에너지, $δ_{i j}$ 는 크로네커델타(Kronecker delta)를 의미한다. 위에서 언급한 바와 같이 하천의 바닥 표면은 일반적으로 사련 형태로 발달해 있어 바닥 근처 점성류에 의해 역압력경사(adverse pressure gradient)가 발생하여 흐름분리, 즉 재순환 흐름 영역이 바닥면 근처에 형성된다(Motamedi et al., 2014). 따라서, 본 연구에서는 난류모델 중 흐름분리 현상 해석에 강점을 지닌 SST k-omega 모델(Menter, 1993)를 이용하여 Eq. (3)의 $ν_{t}$ 를 계산하였다. SST k-omega에서 $ν_{t}$ 는 다음과 같이 표현된다.

(4)

ν_{t} = \frac{a_{1} k}{\max (a_{1} ω, S F_{2})}

여기서 $ω$ 는 난류운동에너지의 비소산율, a₁은 종결계수, S는 와도(vorticity)의 절대값, F₂는 혼합함수를 의미한다. Eq. (4)를 종결하기 위해서는 다음의 추가적인 수송방정식이 필요하다.

(5)

\frac{\partial k}{\partial t} + U_{j} \frac{\partial k}{\partial x_{j}} = - P_{k} - β^{*} k ω + \frac{\partial}{\partial x_{j}} [(ν + σ_{k} ν_{t}) \frac{\partial k}{\partial x_{j}}]

(6)

\frac{\partial ω}{\partial t} + U_{j} \frac{\partial ω}{\partial x_{j}} = α S^{2} - β ω^{2} + \frac{\partial}{\partial x_{j}} [(ν + σ_{ω} ν_{t}) \frac{\partial ω}{\partial x_{j}}] + 2 (1 - F_{1}) σ_{ω 2} \frac{1}{ω} \frac{\partial k}{\partial x_{i}} \frac{\partial ω}{\partial x_{i}}

여기서 P_k는 생성제어기, F₁는 혼합함수, $α, β, β^{*}, σ_{k}, σ_{ω}, σ_{ω 2}$ 는 종결계수이며, 본 연구에서는 Menter (1993)에서 제시된 값들을 사용하였다.

하상 내 유동은 다공성매질을 통과하는 층류로 가정하여, 다음의 Darcy 법칙기반 흐름방정식을 이용하여 정상상태 2차원 유속을 계산하였다.

(7)

\frac{\partial q_{i}}{\partial x_{i}} = 0

(8)

q_{i} = - \frac{k}{μ} \frac{\partial p}{\partial x_{i}}

여기서 q_i는 Darcy 유속, k는 투수계수, $μ$ 는 물의 점성계수를 의미한다.

2.2 용질거동 해석

지표수-지하수 연계 흐름 모의결과의 흐름장과 2차원 이송-분산 방정식을 결합하여 용질거동을 모의하였으며, 지표수 영역에서 용질의 농도는 다음의 지배방정식을 이용하여 계산된다.

(9)

\frac{\partial C}{\partial t} + U_{i} \frac{\partial C}{\partial x_{i}} = \frac{\partial}{\partial x_{i}} (D_{t} \frac{\partial C}{\partial x_{i}})

여기서 C는 농도, D_t는 $ν_{t} / S C_{t}$ 로 정의되는 난류확산계수이며, SC_t는 난류 Schmidt수를 의미한다. 과거 연구들에 따르면 SC_t를 0.7로 설정하였을 때 추적자실험 측정값과 수치실험 모의값의 오차가 가장 낮게 나타나는 것으로 보고된 바 있어 본 연구에서도 해당 값을 사용하였다. 지하수 영역에서 용질의 농도는 다음의 2차원 이송-분산 방정식을 통해 계산된다.

(10)

\frac{\partial C}{\partial t} + q_{i} \frac{\partial C}{\partial x_{i}} = \frac{\partial}{\partial x_{i}} [(D_{m} + D_{i j}) \frac{\partial C}{\partial x_{i}}]

여기서 D_m는 분자확산계수, D_ij는 역학적 분산계수 텐서(mechanical dispersion coefficient tensor)를 의미한다. 역학적 분산은 다공성매질 흐름구조에 의해 발생하는 수리동역학적 확산과정으로 다음과 같이 정의될 수 있다(Alavian, 1986).

(11)

D_{x x} = D_{L} \frac{q_{x}^{2}}{q^{2}} + D_{V} \frac{q_{y}^{2}}{q^{2}}

(12)

D_{x y} = D_{y x} = (D_{L} - D_{V}) \frac{q_{x} q_{y}}{q^{2}}

(13)

D_{y y} = D_{V} \frac{q_{x}^{2}}{q^{2}} + D_{L} \frac{q_{y}^{2}}{q^{2}}

여기서 q는 $\sqrt{q_{x}^{2} + q_{y}^{2}}$ 로 정의되고, D_L과 D_V는 각각 종방향, 연직방향 분산계수를 의미한다. D_L과 D_V는 다공성매질의 유속, 입자 크기 등에 따라 상이한 값을 지니며, 본 연구에서는 Harleman and Rumer (1963)이 제안한 다음의 경험식을 이용하여 산정하였다.

(14)

D_{L} = 1.6785 {(d U)}^{1.2}

여기서 d는 유사의 평균입경을 의미하며 D_L은 유속을 변수로 하는 멱함수로 나타낼 수 있고, D_V는 D_L의 1/10로 가정하였다(Moltyaner and Killey, 1988; Cardenas et al., 2008). 본 연구에서 주입된 용질을 비점착성 물질로 가정하여 하상 내에서 발생하는 흡착과정은 Eq. (10)에 고려되지 않았다.

3. 수치모의 계산조건 및 검증

3.1 모델 검증

위에서 설명한 혼합대 모델링 방법을 적용하기에 앞서 Janssen et al. (2012)의 실험자료를 이용하여 흐름모의 결과를 검증하였다. 해당 실험에서는 Fig. 2에 해당하는 길이 1.5 m 수로에서 높이 0.02 m, 길이 0.2 m (앞단 0.15 m, 뒷단 0.05 m)의 사련 구조가 연속적으로 발달한 두께 0.09 m 하상 위의 수심 0.1 m, 평균유속(U₀) 0.07 m/s 개수로 흐름을 측정하였으며, 흐름분리 영역이 사련하상 표면에 연속적으로 발생함을 바닥면 근처에서 흐름방향 기준으로 역방향 유속결과(-)를 통해 확인하였다. 해당 실험에서 투수계수는 k = 1.5 × 10^-11 m²으로 설정하였는데 이는 사련파가 관측되는 고운 모래(fine sand) 하상에 해당하는 수치이며(Freeze and Cherry, 1979; Raudkivi, 1997), 해당 투수계수는 유사 평균입경 0.174 mm, 표준편차 0.15 mm에 상응한다(Janssen et al., 2012).

Fig. 2에 나타낸 바와 같이 지표수 흐름해석의 RANS 모의에서 유입부에 실험조건과 동일한 평균 유속 0.07 m/s, 유출부에는 압력 경계조건을 입력하였다. 수표면은 실험조건의 Froude 수가 0.1 이하로 매우 낮음에 따라 자유수면(free surface)이 발생하지 않음을 가정하여 대칭(symmetry) 경계조건, 바닥면은 점착(no slip) 경계조건을 설정하였다. Darcy 방정식을 이용한 지하수 흐름해석에서는 Cardenas and Wilson (2007)이 제시한 방법과 같이 RANS 모의에서 얻어진 바닥면 압력분포를 지표수-지하수 경계의 경계조건으로 입력하였고, 입출부 경계에는 무경사(zero gradient) 경계조건, 바닥에는 점착 경계조건을 설정하였다. 그리고, 지표수와 지하수 영역의 격자수 모두 20,000개(흐름방향 400개 × 연직방향 50개)의 사각망으로 동일하게 구성하였고, RANS 모의에서 바닥면 근처에서 발생하는 흐름분리 현상을 정밀하게 모의하기 위해 무차원 벽면 거리 $y^{+}$ (= $y_{p} u * / ν$ )를 1 이하로 설정하였다(Chen et al., 2015; Lee et al., 2020). 여기서, $y_{p}$ 는 바닥에서 첫 번째 격자까지의 거리, $u *$ 는 전단유속(shear velocity)을 의미하며, 계산시간 간격( $∆ t$ )은 Courant 수 0.4 이하를 만족하도록 설정하였다(Lee et al., 2020; Li et al., 2020).

Janssen et al. (2012) 실험에서는 하상 내부의 유속은 측정하지 않았고, 바닥면 근처의 흐름분리 현상을 확인할 수 있는 지표수 흐름의 종방향 유속분포와 혼합대 흐름 구조 및 강도를 결정하는 바닥면 압력분포를 측정하였다. 본 연구의 RANS 모의결과와 해당 실험의 측정값들을 비교하였으며, 이를 Figs. 3 and 4에 도시하였다.

실험값과 비교결과, RANS에서 모의된 종방향 유속이 실험값에 매우 근사하게 나타나는 것을 확인할 수 있으며, 특히 바닥면 근처에서의 음수값의 종방향 유속을 정확하게 모의함에 따라 사련하상표면 형상에 의한 재순환 흐름(흐름분리) 발생을 적절히 재현한 것으로 사료된다(Fig. 3). 더 나아가, 혼합대 흐름해석의 경계조건으로 이용되는 바닥면 압력 역시 사련 하상파 고점 사이에서 발생한 역압력경사를 포함하여 압력분포의 전체적인 경향을 정확하게 모의하였다(Fig. 4). RANS에서 취득된 바닥면 압력분포를 Darcy 방정식과 결합하여 모의한 혼합대의 유속장은 Fig. 1(b)에 나타낸 바와 같으며, 바닥면 압력경사에 의해 혼합대 내에 2차원적 흐름 구조가 발생함을 확인할 수 있다.

https://cdn.apub.kr/journalsite/sites/kwra/2022-055-01/N0200550103/images/kwra_55_01_03_F3.jpg

Fig. 3.

Simulated longitudinal velocity profiles against observation of Janssen et al. (2012)

https://cdn.apub.kr/journalsite/sites/kwra/2022-055-01/N0200550103/images/kwra_55_01_03_F4.jpg

Fig. 4.

Comparison of measured (Janssen et al., 2012) and simulated bottom pressure

3.2 모델 적용

Janssen et al. (2012)의 실험값과 검증된 흐름 모의결과를 바탕으로 용질거동 모의를 수행하였다. 검증사례의 경우, 유출입부 경계조건으로 인해 계산영역 전체를 완전 발달된 흐름(fully developed flow)으로 가정하기 어렵다. 따라서, 계산영역을 Fig. 5와 같이 단축시킴과 동시에 해당 영역에서 완전 발달 흐름을 모의할 수 있는 주기(periodic) 경계조건을 지표수와 지하수 영역 유출입부 경계에 적용하여 지표수-지하수 연계 흐름해석을 하였다. 주기 경계조건은 유출부를 통해 빠져나간 유동이 다시 유입부로 유입되는 것으로 가정함에 따라 계산영역을 무한한 연장을 지닌 수로로 해석할 수 있다. 본 연구에서는 지표수 영역의 평균유속을 0.07 m/s과 0.12 m/s로 설정하여, 두 가지 다른 지표수 흐름조건에 대해 용질거동을 모의하였다. 여기서, 수심, 사련의 길이 및 높이는 Janssen et al. (2012)의 실험수로와 동일하게 설정하였으나, 하상 두께는 0.7 m 약 8배 증가시켰다. Ren and Zhao (2020)는 Janssen et al. (2012)와 동일한 흐름 및 지형 조건에서 지하수(하상) 영역 바닥면의 점착 경계조건에 의해 혼합대의 유동 구조가 변화함을 확인하였고, 하상 두께를 0.7 m 이상으로 증가시켰을 때 혼합대 흐름에 대한 바닥 경계조건 영향이 나타나지 않음을 보고하였다. 또한, 용질거동 모의를 통한 체류시간 분포 해석의 계산시간을 단축하기 위해 하상재료는 얇은 자갈로 가정하여 투수계수를 k = 5 × 10^-6 m²으로 증가시켰다. 위의 모의조건을 기반으로 한 지표수 평균유속(지표수 영역 유입부 경계조건) 0.07 m/s, 0.12 m/s에 대한 혼합대 흐름 모의결과는 Fig. 5와 같다.

Fig. 5의 흐름장을 Eqs. (9) and (10)의 2차원 이송-분산 방정식과 결합하여 용질거동을 모의하였다. 수치모의 기반의 추적자실험을 위해 Fig. 5의 흐름 구조가 하류에서 균일하게 반복된다고 가정하여 유로연장 2 m 직선수로에 해당 값들을 투영시킨 후, 최상류에서 0.2 m 떨어진 두 번째 사련 최저점(trough)에서 용질을 점으로 순간 주입하였다(Fig. 5(a)). 지표수의 전단흐름(shear flow)에 의한 분산 효과를 최소화하고, 혼합대 흐름이 체류시간 분포에 미치는 효과를 보다 명확하게 하기 위해 용질을 지표수-하상 경계에서 주입하였다. Eq. (10)의 D_m은 추적자로 일반적으로 사용되는 소듐아이오다이드(NaI)의 분자확산계수인 1.1 × 10^-5 cm²/s를 입력하였으며, D_L은 d를 8 mm로 가정하고 Fig. 5의 지하수 영역의 평균유속을 Eq. (14)에 대입하여 계산하였다. 위에서 설명한 모의조건을 바탕으로 혼합대로 용질이 유입되는 경우와 유입되지 않고 지표수 영역에서만 용질이 거동하는 경우를 비교 모의하여, 혼합대가 용질거동에 미치는 영향을 정량적으로 분석하였다.

https://cdn.apub.kr/journalsite/sites/kwra/2022-055-01/N0200550103/images/kwra_55_01_03_F5.jpg

Fig. 5.

Computational domains with boundary conditions for model application cases, and simulated open channel pressure and groundwater (hyporheic) velocity magnitude fields for (a) U₀ = 0.07 m/s and (b) U₀ = 0.12 m/s. Solid lines represent hyporheic flow patterns

4. 모의결과 및 고찰

4.1 흐름 모의결과

용질거동 수치모의에 앞서 지표수 유속 변화에 따른 흐름해석 결과를 분석하였다. Fig. 6에 도시한 바와 같이 지표수 평균 유속이 약 1.7배 높아짐에 따라 하상 경계에서 최대 압력이 10배 이상 상승하였고, 결과적으로 혼합대 내 평균유속 강도가 4.6 × 10^-4 m/s에서 13.8 × 10^-4 m/s로 약 3배 가까이 증가함을 확인할 수 있다. Fig. 6의 혼합대 내 연직방향 유속분포를 이용하여 지표수-하상 경계층에서 혼합대의 모멘텀 플럭스를 계산한 결과, 지표수 유속이 증가함에 따라 모멘텀 플럭스가 2배 이상 높아지는 것으로 나타났다. 여기서 양의 방향과 음의 방향의 플럭스는 질량보존 법칙에 따라 동일한 값을 지니게 되는데, 본 연구에서는 Cardenas and Wilson (2007)이 제안한 방법에 따라 혼합대 내 연직방향 유속을 지표수-하상 경계층에 따라 적분하여 절반값을 취한 후 수로 길이로 나누어 혼합대 플럭스를 계산하였다.

혼합대 내 유속 강도와 달리 지표수 유속 변화에 대한 혼합대 내 유선(streamline) 형태의 민감도는 크게 확인되지 않았다(Fig. 5). 이는 과거 연구에서 보고된 바와 같이 혼합대 흐름 구조가 지표수의 흐름조건 즉, 유속 경계조건이 아닌 하상의 표면 형태에 지배적인 영향을 받기 때문으로 설명된다(Lee et al., 2020; Ren and Zhao, 2020). 또한, 본 연구의 흐름 모의결과, Fig. 5에서 확인된 바와 같이 하상의 깊은 영역에서 유속이 0에 가까운 정체점(stagnation point)이 형성되었는데, 이러한 흐름 현상은 해당 지점에서 종방향 및 연직방향 유속의 합이 0에 접근하여 흐름이 상쇄됨에 의한 것이며, 정체구간의 발생 깊이는 과거 연구에서 보고된 바와 같이 혼합대의 깊이와 유사하게 나타났다(Anderson, 2003; Marzadri et al., 2016). 본 연구에서 정체구간 깊이로 대변되는 혼합대의 깊이는 지표수 유속 증가에 의해 크게 변화하지 않았는데, 이는 위에서 언급한 바와 같이 혼합대의 유선 형태가 지표수의 유속 변화에 민감하지 않기 때문이다.

https://cdn.apub.kr/journalsite/sites/kwra/2022-055-01/N0200550103/images/kwra_55_01_03_F6.jpg

Fig. 6.

Simulated (a) bottom pressure distribution at the sediment-water interface (SWI) and (b) hyporheic fluxes across the SWI for model application cases

4.2 용질거동 모의결과

지표수-지하수 연계 흐름 모의결과의 유속 및 난류점성계수 분포를 2차원 이송-분산 방정식에 결합하여 용질거동을 모의하였다. 지표수 용질거동에 미치는 혼합대의 영향을 정량화하기 위해 먼저 지표수와 하상 간의 용질교환이 발생하지 않는 불투수성 하상을 대상으로 모의를 수행하였다. Fig. 7은 용질거동 모의를 바탕으로 파과곡선을 도시한 것이며, 하천 추적자실험과 유사하게 각 농도의 도달시간은 지표수 영역 유출부(x = 2 m)에서만 측정되었다.

불투수성 하상에 대한 모의결과, Fig. 7에 실선으로 나타낸 바와 같이 지표수 유속이 낮은 경우에는 상대적으로 낮은 지표수 영역의 최대유속으로 인해 최대농도 도달시간이 지표수 유속이 높은 경우 대비 약 2배 빨라졌으나, 최저농도(C/C_max = 10^-5) 도달시간은 지표수 유속 증가에 따라 약 1.4배 지체되는 것으로 나타났다. 이는 지표수 유속 증가와 함께 흐름분리 영역 내 재순환 흐름에 의한 관성효과가 커져 흐름분리 영역에 포착(trapping)된 용질이 난류확산에 의해 유속이 빠른 주 흐름 영역으로 쉽게 이동하지 못하고, 흐름분리 영역 내에서 상대적으로 긴 시간 체류하였기 때문이다(Kim et al., 2020). 본 연구에서는 혼합대, 재순환 흐름 등 저장대 효과가 반영된 non-Fickian 거동을 설명하기 위해 파과곡선 꼬리의 연장을 최저농도 도달시간과 최대농도 도달시간의 비( $t_{T_{c}}$ = t_max/t_c)로 정의되는 절단시간(truncated time, t_T)을 통해 정량화하였다(Drummond et al., 2012). non-Fickian 혹은 비정상(anomalous) 거동은 파과곡선 꼬리 부분에서 비정상적으로 높은 용질 농도의 검출로 정의된다(Kim and Kang, 2020). 불투수성 하상의 파과곡선 절단시간을 계산한 결과, 지표수 유속이 증가함에 따라 절단시간이 약 3배 증가하였다. 이는 위에서 언급한 바와 같이 벽면 근처 재순환 흐름의 용질 포착 효과 강화에 따른 최저농도 도달시간 지체에 의한 것으로 설명된다.

다음으로 혼합대 흐름, 즉 지표수와 하상 간의 흐름교환이 반영된 투수성 하상에 대해 용질거동 모의를 수행하였다. 불투수성 하상의 파과곡선과 비교하였을 때, Fig. 7에 점선으로 나타낸 바와 같이 투수성 하상의 경우 두 가지 다른 지표수 유속 경계조건에서 모두 파과곡선 꼬리가 급격하게 길어짐을 확인할 수 있다. 이러한 경향은 난류확산과 사련하상 표면의 압력경사(Fig. 6(a))가 동반하는 이송펌핑(advective pumping) 효과로 불리는 혼합대 고유의 2차원적 흐름 구조(Fig. 5)에 의해 지표수 흐름 대비 유속이 매우 느린 혼합대로 용질이 유입되고(Packman et al., 2004), 유입된 용질이 혼합대 내에 장시간 체류한 후 지표수 영역으로 다시 돌아오는 과정이 하류를 따라 지속적으로 반복된 결과이다. 투수성 하상에서는 불투성 하상과는 달리 지표수 유속이 증가함에 따라 파과곡선 꼬리가 짧아지는 현상이 나타났고, 이는 4.1에서 설명한 바와 같이 지표수 유속 증가와 함께 혼합대 흐름의 평균 유속이 증가하여 혼합대 내에 정체된 저농도의 용질이 빠르게 지표수 방향으로 이송되었기 때문이다. 투수성 하상의 파과곡선 절단시간 계산 결과, 지표수 평균유속 0.07 m/s, 0.12 m/s에 대해 절단시간이 불투수성 하상의 결과값 대비 각각 약 5배와 14배 증가하였다. 이는 Fig. 7에 도시한 바와 같이 지표수 유속이 증가함에 따라 바닥면 근처 재순환 흐름과 혼합대 흐름에 의한 최저농도 도달시간 편차가 감소한 것에 기인한다. 더 나아가, 투수성 하상에서는 불투수성 하상과 달리 지표수 유속 변화에도 불구하고 혼합대 흐름에 의해 결정되는 최저농도 도달시간과 최대농도 도달시간의 증가비가 유사하여 유사한 값의 절단시간이 계산되었다. 위의 결과를 통해 혼합대 내 흐름 특성이 체류시간 분포 형태, 특히 파과곡선 꼬리 연장을 결정함에 따라 혼합대 흐름이 지표수 용질거동 메커니즘에 지배적인 영향인자로 작용함을 확인할 수 있다.

https://cdn.apub.kr/journalsite/sites/kwra/2022-055-01/N0200550103/images/kwra_55_01_03_F7.jpg

Fig. 7.

Simulated breakthrough curves (BTCs) measured in the surface flow region at the channel outlet for model application cases. Black boxes indicate an increase in BTC tailing induced by flow recirculation and hyporheic flow

5. 요약 및 결론

본 연구에서는 2차원 RANS 방정식과 Darcy 방정식을 연계한 지표수-지하수 연계흐름 모의를 수행하였고, 모의된 흐름장과 2차원 이송-분산 방정식을 이용하여 지표수-하상 경계층에서 혼합대 흐름교환이 지표수 용질거동에 미치는 영향을 수치해석적으로 분석하였다. 먼저, 지표수의 흐름 특성과 하상 표면 형태에 큰 영향을 받는 바닥면 압력분포와 밀접한 인과관계를 지니는 혼합대 흐름 특성을 정확하게 구현하기 위해 사련하상을 대상으로 지표수 영역에서 RANS 모의를 수행하였고, 유속분포와 바닥면 압력분포 수치해석 결과와 기존 연구의 실험값을 비교하였다. 검증결과, SST k-omega 모델이 결합된 RANS 수치모의가 사련 하상파에 의해 발생하는 바닥면 근처 역압력경사와 재순환 흐름을 재현함과 동시에 바닥면 압력분포를 정확하게 모의하였다. RANS에서 얻어진 바닥면 압력분포를 Darcy 방정식의 경계조건으로 하여 하상 내 혼합대 흐름을 해석하였고, 지표수 유속이 빨라짐에 따라 바닥면 압력이 증가하여 혼합대 흐름이 가속화되었다.

혼합대 흐름교환 과정이 고려된 투수성 하상에서 용질거동을 모의한 결과, 바닥면 압력경사가 동반하는 이송펌핑 효과에 의해 용질이 지표수 대비 유속이 느린 혼합대 공극흐름을 경험함에 따라 체류시간이 증가하여 파과곡선 꼬리가 길게 연장됨을 확인할 수 있었다. 여기서, 지표수 유속 증가는 바닥면 압력 상승과 함께 혼합대 흐름의 평균유속을 증가시켰고, 이로 인해 혼합대에서 지표수 방향으로의 용질 이송이 가속화되어 최저농도 도달시간이 단축되었다. 반면, 하상 내부로 용질 유입이 없는 불투수성 하상에서는 하상 표면에 발생하는 흐름분리 영역 내 재순환 흐름에 의한 용질 포착 효과가 파과곡선 꼬리의 형태를 결정하였고, 용질이 지표수 대비 흐름이 매우 느린 혼합대로 유입되지 않아 투수성 하상과 비교했을 때 파과곡선 꼬리의 연장이 과소 산정되는 것으로 나타났다. 결론적으로 본 연구의 결과를 통해 하천과 같은 지표수 내 용질거동 예측에 미치는 혼합대 흐름의 중요도를 정량적으로 확인할 수 있었다. 특히, 혼합대 흐름은 하천의 온실가스 배출과 직접적으로 연관된 질소순환 메커니즘에 지대한 영향을 미치고, 유해화학물질들이 하천으로 유입될 경우 혼합대 내 흐름이 위험농도 잔존시간에 대한 제어인자로 작용하기 때문에 향후 하천 수질관리 및 수질사고대응시스템 구축에 있어 혼합대 흐름의 물리적 특성을 효과적으로 고려 가능한 수질모델 개발 및 국내하천 적용이 필요할 것으로 판단된다.

Acknowledgements

본 연구는 국토연구원의 연구비 지원에 의하여 수행되었습니다.

References

Ahn, M., Kim, Y., Ji, U., Gu, J., Ko, J., Bae, I., and Kang, H. (2019). "Measurement and analysis of nitrous oxide emissions over time around a dune in the experimental flume." Journal of Korean Society of Environmental Engineers, Vol. 41, No. 4, pp. 228-234. 10.4491/KSEE.2019.41.4.228

Alavian, V. (1986). "Dispersion tensor in rotating flows." Journal of Hydraulic Engineering, Vol. 112, No. 8, pp. 771-777. 10.1061/(ASCE)0733-9429(1986)112:8(771)

Anderson, E.I. (2003). "An analytical solution representing groundwater - surface water interaction." Water Resources Research, Vol. 39, No. 3. 10.1029/2002WR001536

Boulton, A.J., Findlay, S., Marmonier, P., Stanley, E.H., and Valett, H.M. (1998). "The functional significance of the hyporheic zone in streams and rivers." Annual Review of Ecology and Systematics, Vol. 29, No. 1, pp. 59-81. 10.1146/annurev.ecolsys.29.1.59

Cardenas, M.B., and Wilson, J.L. (2007). "Dunes, turbulent eddies, and interfacial exchange with permeable sediments." Water Resources Research, Vol. 43, No. 8. 10.1029/2006WR005787

Cardenas, M.B., Wilson, J.L., and Haggerty, R. (2008). "Residence time of bedform-driven hyporheic exchange." Advances in Water Resources, Vol. 31, No. 10, pp. 1382-1386. 10.1016/j.advwatres.2008.07.006

Chen, X., Cardenas, M.B., and Chen, L. (2015). "Three dimensional versus two‐dimensional bed form induced hyporheic exchange." Water Resources Research, Vol. 51, No. 4, pp. 2923-2936. 10.1002/2014WR016848

Drummond, J.D., Covino, T.P., Aubeneau, A.F., Leong, D., Patil, S., Schumer, R., and Packman, A.I. (2012). "Effects of solute breakthrough curve tail truncation on residence time estimates: A synthesis of solute tracer injection studies." Journal of Geophysical Research: Biogeosciences, Vol. 117, No. G3. 10.1029/2012JG002019

Findlay, S. (1995). "Importance of surface-subsurface exchange in stream ecosystems: The hyporheic zone." Limnology and Oceanography, Vol. 40, No. 1, pp. 159-164. 10.4319/lo.1995.40.1.0159

Freeze, R.A., and Cherry, J.A. (1979). Groundwater. PrenticeHall, Englewood Cliffs, NJ, U.S.

Haggerty, R., Wondzell, S.M., and Johnson, M.A. (2002). "Power law residence time distribution in the hyporheic zone of a 2nd order mountain stream." Geophysical Research Letters, Vol. 29, No. 13, pp. 18-1-18-4. 10.1029/2002GL014743

Harleman, D.R.F., and Rumer, R.R. (1963). "Longitudinal and lateral dispersion in an isotropic porous medium." Journal of Fluid Mechanics, Vol. 16, No. 3, pp. 385-394. 10.1017/S0022112063000847

Harvey, J., and Gooseff, M. (2015). "River corridor science: Hydrologic exchange and ecological consequences from bedforms to basins." Water Resources Research, Vol. 51, No. 9, pp. 6893-6922. 10.1002/2015WR017617

Janssen, F., Cardenas, M.B., Sawyer, A.H., Dammrich, T., Krietsch, J., and de Beer, D. (2012). "A comparative experimental and multiphysics computational fluid dynamics study of coupled surface - subsurface flow in bed forms." Water Resources Research, Vol. 48, No. 8. 10.1029/2012WR011982

Jin, G., Zhang, Z., Li, R., Chen, C., Tang, H., Li, L., and Barry, D.A. (2020). "Transport of zinc ions in the hyporheic zone: Experiments and simulations." Advances in Water Resources, Vol. 146, 103775. 10.1016/j.advwatres.2020.103775

Kim, J.S., and Kang, P.K. (2020). "Anomalous transport through free-flow-porous media interface: Pore-scale simulation and predictive modeling." Advances in Water Resources, Vol. 135, 103467. 10.1016/j.advwatres.2019.103467

Kim, J.S., Seo, I.W., Baek, D., and Kang, P.K. (2020). "Recirculating flow-induced anomalous transport in meandering open-channel flows." Advances in Water Resources, Vol. 141, 103603. 10.1016/j.advwatres.2020.103603

Lee, A., Aubeneau, A.F., and Cardenas, M.B. (2020). "The sensitivity of hyporheic exchange to fractal properties of riverbeds." Water Resources Research, Vol. 56, No. 5, e2019WR026560. 10.1029/2019WR026560

Li, B., Liu, X., Kaufman, M.H., Turetcaia, A., Chen, X., and Cardenas, M.B. (2020). "Flexible and modular simultaneous modeling of flow and reactive transport in rivers and hyporheic zones." Water Resources Research, Vol. 56, No. 2, e2019WR026528. 10.1029/2019WR026528

Marzadri, A., Tonina, D., and Bellin, A. (2012). "Morphodynamic controls on redox conditions and on nitrogen dynamics within the hyporheic zone: Application to gravel bed rivers with alternate-bar morphology." Journal of Geophysical Research: Biogeosciences, Vol. 117, No. G3. 10.1029/2012JG001966

Marzadri, A., Tonina, D., Bellin, A., and Valli, A. (2016). "Mixing interfaces, fluxes, residence times and redox conditions of the hyporheic zones induced by dune-like bedforms and ambient groundwater flow." Advances in Water Resources, Vol. 88, pp. 139-151. 10.1016/j.advwatres.2015.12.014

Menter, F.R. (1993). "Zonal two-equation k-ω turbulence model for aerodynamic flows." Proceedings of AIAA 24th Fluid Dynamic Conference, AIAA Paper 93-2906, Orlando, FL, U.S. 10.2514/6.1993-2906PMC275754

Moltyaner, G.L., and Killey, R.W.D. (1988). "Twin lake tracer tests: Transverse dispersion." Water Resources Research, Vol. 24, No. 10, pp. 1628-1637. 10.1029/WR024i010p01628

Motamedi, A., Afzalimehr, H., Singh, V.P., and Dufresne, L. (2014). "Experimental study on the influence of dune dimensions on flow separation." Journal of Hydrologic Engineering, Vol. 19, No. 1, pp. 78-86. 10.1061/(ASCE)HE.1943-5584.0000754

Norman, F.A., and Cardenas, M.B. (2014). "Heat transport in hyporheic zones due to bedforms: An experimental study." Water Resources Research, Vol. 50, No. 4, pp. 3568-3582. 10.1002/2013WR014673

Packman, A.I., Salehin, M., and Zaramella, M. (2004). "Hyporheic exchange with gravel beds: Basic hydrodynamic interactions and bedform-induced advective flows." Journal of Hydraulic Engineering, Vol. 130, No. 7, pp. 647-656. 10.1061/(ASCE)0733-9429(2004)130:7(647)

Raudkivi, A.J. (1997). "Ripples on stream bed." Journal of Hydraulic Engineering, Vol. 123, No. 1, pp. 58-64. 10.1061/(ASCE)0733-9429(1997)123:1(58)

Ren, J., and Zhao, B. (2020). "Model-based analysis of the effects of rippled bed morphologies on hyporheic exchange." Journal of Hydrologic Engineering, Vol. 25, No. 6, 04020023. 10.1061/(ASCE)HE.1943-5584.0001931

Ren, J., Wang, X., Zhou, Y., Chen, B., and Men, L. (2019). "An analysis of the factors affecting hyporheic exchange based on numerical modeling." Water, Vol. 11, No. 4, 665. 10.3390/w11040665

Ribberink, J.S., and Al‐Salem, A.A. (1994). "Sediment transport in oscillatory boundary layers in cases of rippled beds and sheet flow." Journal of Geophysical Research: Oceans, Vol. 99, No. C6, pp. 12707-12727. 10.1029/94JC00380

Runkel, R.L., and Chapra, S.C. (1993). "An efficient numerical solution of the transient storage equations for solute transport in small streams." Water Resources Research, Vol. 29, No. 1, pp. 211-215. 10.1029/92WR02217

Journal of Korea Water Resources Association ISSN:2799-8746(Print) 2799-8754(Online) 한국수자원학회 논문집

Preview

Numerical analysis of the hyporheic flow effect on solute transport in surface water

ABSTRACT

MAIN

Fig. 1.

Conceptual diagrams of (a) bedform-induced hyporheic exchange and resulting tracer transport and (b) key hydrodynamic processes controlling interfacial mass exchange

Fig. 2.

An experimental flume of Janssen et al. (2012) and boundary conditions used for coupled open channel flow-groundwater flow simulation

(1)

(2)

(3)

(4)

(5)

(6)

(7)

(8)

(9)

(10)

(11)

(12)

(13)

(14)

Fig. 3.

Simulated longitudinal velocity profiles against observation of Janssen et al. (2012)

Fig. 4.

Comparison of measured (Janssen et al., 2012) and simulated bottom pressure

Fig. 5.

Computational domains with boundary conditions for model application cases, and simulated open channel pressure and groundwater (hyporheic) velocity magnitude fields for (a) U0 = 0.07 m/s and (b) U0 = 0.12 m/s. Solid lines represent hyporheic flow patterns

Fig. 6.

Simulated (a) bottom pressure distribution at the sediment-water interface (SWI) and (b) hyporheic fluxes across the SWI for model application cases

Fig. 7.

Simulated breakthrough curves (BTCs) measured in the surface flow region at the channel outlet for model application cases. Black boxes indicate an increase in BTC tailing induced by flow recirculation and hyporheic flow

Acknowledgements

References

Computational domains with boundary conditions for model application cases, and simulated open channel pressure and groundwater (hyporheic) velocity magnitude fields for (a) U₀ = 0.07 m/s and (b) U₀ = 0.12 m/s. Solid lines represent hyporheic flow patterns