유사 이송

편집 역사 토론
유체역학 전산유체역학 시뮬레이션 마지막 수정: 2026-09-01 04:13:27

1. 개요[편집]

유사 이송
Sediment transport
양식 구분소류사(bed load) · 부유사(suspended load) · 세류사(wash load)
구동 무차원수실즈 수 τ* = τb / ((ρs−ρ)gd)
부유 판별루스 수 P = ws / (κ u*)
하상변동엑스너 방정식 (1925)
대표 소류사식마이어페터-뮐러 · 아인슈타인 · 반 레인
대표 코드HEC-RAS Sediment · Delft3D · TELEMAC-GAIA
핵심 난점흐름은 초 단위, 하상은 연 단위 — 시간규모가 4자릿수 다르다

물만 푸는 사람은 바닥이 고정된 세상에 산다. 모래가 끼는 순간, 경계조건이 미지수가 된다.

유사 이송은 흐르는 물이 하상과 하안의 토사를 움직여 운반하고, 그 결과 하상 자체가 변형되는 현상 및 그 정량화 이론이다. 얕은 물 방정식이 물만 푸는 것이라면, 유사 이송은 여기에 “바닥이 움직인다”는 조건을 얹은 것이고, 그 순간 문제는 고정 경계 유동에서 이동 경계 자유경계 문제로 승격한다.

핵심 감각 하나만 먼저 잡고 가자. 하상이 변형되는 속도는 물이 흐르는 속도보다 네 자릿수쯤 느리다. 홍수 첨두는 몇 시간이면 지나가는데 여울 하나가 이동하는 데는 몇 년이 걸린다. 이 시간규모 분리가 이 분야의 수치적 축복이자 저주다 — 축복인 이유는 흐름과 하상을 떼어 풀 수 있어서고, 저주인 이유는 프루드 수가 1에 가까워지는 순간 그 분리가 통째로 무너지기 때문이다.

2. 소류사와 부유사[편집]

같은 모래알이라도 어떻게 움직이느냐에 따라 완전히 다른 수식을 쓴다.

  • 소류사(bed load) — 굴러가고, 미끄러지고, 짧게 튀어(saltation) 하상 근처를 벗어나지 않는다. 입자 무게가 하상에 직접 전달되므로 접촉하중(contact load)이라고도 한다. 유량은 단위 폭당 부피 유량 qb [m2/s]q_b\,[\mathrm{m^2/s}] 로 쓰고, 국소 하상 전단응력의 대수적 함수로 모형화한다. 관성이 짧아 흐름 변화에 거의 즉시 반응한다 — 이 “즉시”라는 성질이 뒤에서 결정적으로 중요해진다.
  • 부유사(suspended load) — 난류 와가 침강을 이겨 물기둥 안에 떠다닌다. 무게가 유체에 전달되므로 하상 전단응력에 기여하지 않는다. 이송 자체는 대류-확산 방정식 형태의 편미분방정식으로 풀어야 하며, 기억(memory)이 있어 흐름 변화에 지연되어 반응한다.
  • 세류사(wash load) — 너무 잘아서 애초에 하상에 없는 성분(실트·점토). 흐름이 얼마든 다 떠서 나가므로 이송량이 하상이 아니라 상류 공급량으로 결정된다. 유사 이송 공식으로 예측하려 들면 안 되고, 실제로 예측하지 않는다. 하상변동에도 거의 기여하지 않으나 저수지 퇴사와 수질에서는 이쪽이 주인공이다.

경계를 가르는 것은 루스 수다. 침강속도 wsw_s 와 마찰속도 u∗=τb/ρu_*=\sqrt{\tau_b/\rho} 의 비로,

P=wsκ u∗,κ≈0.41P=\frac{w_s}{\kappa\,u_*},\qquad \kappa\approx0.41

대략 P≳2.5P\gtrsim2.5 면 소류사, 0.8≲P≲2.50.8\lesssim P\lesssim2.5 면 부유사, P≲0.8P\lesssim0.8 이면 세류사로 본다.1 경계값은 문헌마다 조금씩 다르고, 하나의 하천 단면에서 입경별로 세 양식이 동시에 존재한다는 것이 실무의 현실이다. 침강속도는 잔 입자면 스토크스 유동 해 ws=(s−1)gd2/(18ν)w_s=(s-1)gd^2/(18\nu) 로 충분하지만(s=ρs/ρ≈2.65s=\rho_s/\rho\approx2.65), 모래 크기부터는 항력계수가 레이놀즈수에 의존해 경험식(루비, 딧리치, 반 레인)을 쓴다.

3. 한계 소류력과 실즈 곡선[편집]

입자가 언제 움직이기 시작하는가. 이 질문에 대한 표준 답이 실즈 수(Shields number)다. 하상 전단응력을 입자의 수중 중량으로 무차원화한 것으로,

τ∗=τb(ρs−ρ) g d\tau_*=\frac{\tau_b}{(\rho_s-\rho)\,g\,d}

이며, 분자는 입자를 굴리려는 힘, 분모는 붙잡는 힘이라고 읽으면 된다. 실즈(A. Shields, 1936)는 이 값의 임계치 τ∗c\tau_{*c} 가 입자 레이놀즈수 Re∗=u∗d/ν\mathrm{Re}_*=u_*d/\nu 의 함수임을 실험으로 정리했다. 이것이 실즈 곡선이다.

곡선의 모양에는 경계층 물리가 그대로 들어 있다. Re∗≲1\mathrm{Re}_*\lesssim1 에서는 입자가 점성저층에 완전히 잠겨 있어 τ∗c∝Re∗−1\tau_{*c}\propto\mathrm{Re}_*^{-1} 로 떨어지고, Re∗≈10\mathrm{Re}_*\approx10 부근에서 τ∗c≈0.03\tau_{*c}\approx0.03 의 최소값을 지나며, Re∗≳400\mathrm{Re}_*\gtrsim400 의 완전거친 영역에서는 τ∗c≈0.045 ⁣− ⁣0.06\tau_{*c}\approx0.045\!-\!0.06 로 평평해진다. 벽함수에서 보던 그 세 영역과 같은 이야기다.

여기 실무자를 골탕 먹이는 함정이 하나 있다. 실즈 곡선은 음함수다. 양 축에 모두 u∗u_* 가 들어 있어서, 입경 dd 를 알고 임계 전단응력을 구하려면 반복해야 한다. 그래서 무차원 입경

d∗=d[(s−1)gν2]1/3d_*=d\left[\frac{(s-1)g}{\nu^2}\right]^{1/3}

을 가로축으로 다시 그린 판(얄린 도표)이나 소울스비-화이트하우스류의 명시적 맞춤식을 쓴다. 코드 안에서는 후자가 사실상 표준이다.

그리고 이 “한계”라는 말은 생각보다 무르다. 실측 τ∗c\tau_{*c} 의 산포는 쉽게 2배가 나고, 하상 경사·입도 분포·노출도(hiding-exposure)·생물막에 따라 다 달라진다. 혼합입경 하상에서는 잔 입자가 굵은 입자 뒤에 숨어 늦게 움직이고 굵은 입자는 더 노출되어 일찍 움직이는 은폐 효과 때문에, 입경별 τ∗c\tau_{*c} 가 단일입경 값에서 크게 벗어난다. 이것이 뒤에 나올 armoring의 뿌리다.

4. 소류사 공식[편집]

τ∗>τ∗c\tau_*>\tau_{*c} 일 때 얼마나 움직이는가. 무차원 소류사량(아인슈타인 수)

qb∗=qb(s−1)g d3q_{b*}=\frac{q_b}{\sqrt{(s-1)g\,d^3}}

을 τ∗\tau_* 의 함수로 쓰는 것이 표준 형식이다.

마이어페터-뮐러(Meyer-Peter & Müller, 1948). 취리히 ETH 수로 실험에서 나온 국룰.

qb∗=8 (τ∗−τ∗c)3/2,τ∗c=0.047q_{b*}=8\,(\tau_*-\tau_{*c})^{3/2},\qquad \tau_{*c}=0.047

여기서 반드시 짚어야 할 것 — 원논문의 τ∗\tau_* 는 총 전단응력이 아니라 하상형상 저항을 뺀 입자저항분이며, 그 보정(리플 팩터 μ=(n′/n)3/2\mu=(n'/n)^{3/2})을 빼먹고 계수 8을 그대로 쓰는 코드가 흔하다. 웡과 파커(Wong & Parker, 2006)가 MPM 원자료를 평탄하상 조건으로 재분석해 얻은 계수는 8이 아니라 약 4(qb∗≈3.97(τ∗−0.0495)3/2q_{b*}\approx3.97(\tau_*-0.0495)^{3/2})다. 즉 보정 없이 8을 쓰면 소류사량을 계통적으로 2배쯤 과대평가한다.2 브리프든 교과서든 “MPM = 8”만 외우고 있다면 이 대목을 확인하고 넘어가야 한다.

아인슈타인(H. A. Einstein, 1942·1950). 결정론적 문턱값 자체를 거부한 계보다. 난류 양력이 요동하므로 입자의 운동은 확률이고, 단위 시간당 이탈 확률과 평균 이동거리(입경의 약 100배)를 곱해 이송량을 낸다. 그래서 아인슈타인 공식에는 임계 전단응력이 명시적으로 안 나오고, 아주 작은 τ∗\tau_* 에서도 이송량이 0이 아니다. 브라운이 단순화한 아인슈타인-브라운 식은 강한 이송 구간에서 qb∗∝τ∗3q_{b*}\propto\tau_*^{3} 의 3제곱 의존을 준다. 문턱 이론과 확률 이론이 오늘날까지 공존하는 것은 어느 쪽도 결정적으로 이기지 못했기 때문이다.

반 레인(van Rijn, 1984). 이송단계 T=(τ∗′−τ∗c)/τ∗cT=(\tau_*'-\tau_{*c})/\tau_{*c} 와 무차원 입경 d∗d_* 로 쓴 반경험식.

qb∗=0.053 T2.1d∗0.3q_{b*}=0.053\,\frac{T^{2.1}}{d_*^{0.3}}

부유사 기준농도와 한 세트로 제시되어 있어 소류사와 부유사를 일관되게 처리할 수 있다는 것이 실무 채택 이유다. Delft3D 계열의 기본값이 여기서 나온다.

이 밖에 총유사량을 한 방에 주는 식들(엥글룬-한센, 애커스-화이트, 양(Yang), 토팔레티)이 있고, HEC-RAS는 이들을 메뉴로 늘어놓고 사용자에게 고르게 한다. 같은 단면에 같은 유량을 넣고 공식만 바꾸면 이송량이 한 자릿수씩 벌어지는 일이 드물지 않다. 이 바닥에서 “공식 선택 = 최대 불확실성”이라는 말이 나오는 이유이고, 결국 현장 실측으로 보정하지 않은 유사 모델은 정성적 도구로만 취급하는 것이 정직하다.

5. 부유사 — 루스 분포와 비평형[편집]

정상 균일류에서 난류 확산과 침강이 균형을 이루면

wsc+εsdcdz=0w_s c+\varepsilon_s\frac{\mathrm{d}c}{\mathrm{d}z}=0

이고, 여기에 대수법칙에서 나오는 포물선형 와점성 εs=κu∗z(1−z/h)\varepsilon_s=\kappa u_*z(1-z/h) 를 넣어 적분하면 루스 분포가 나온다.

c(z)ca=[h−zz⋅ah−a]P\frac{c(z)}{c_a}=\left[\frac{h-z}{z}\cdot\frac{a}{h-a}\right]^{P}

aa 는 기준 높이(보통 2d2d 또는 0.05h0.05h), cac_a 는 그 높이의 기준농도다. PP 가 작으면 농도가 수심 전체에 고르게 퍼지고, 크면 하상 근처에만 몰린다 — 루스 수의 물리적 의미가 이 지수 하나에 다 들어 있다. 문제는 cac_a 를 이론으로 못 준다는 것이다. 스미스-매클린, 반 레인, 자이서만-프레드쇠 같은 경험식이 난립하고, 이 기준농도 하나의 선택이 총 부유사량을 좌우한다.

실제 하천은 정상 균일류가 아니다. 홍수가 오르면 흐름이 먼저 세지고 부유사 농도가 뒤늦게 따라 올라가며, 내려갈 때도 마찬가지로 늦게 떨어진다. 그래서 유량-유사량 관계가 히스테리시스 루프를 그린다.3 이를 담으려면 평형 가정을 버리고 농도의 이송방정식을 직접 푼다.

∂(hcˉ)∂t+∇ ⁣⋅ ⁣(hucˉ)=∇ ⁣⋅ ⁣(hεs∇cˉ)+E−D\frac{\partial(h\bar c)}{\partial t}+\nabla\!\cdot\!\left(h\mathbf{u}\bar c\right)=\nabla\!\cdot\!\left(h\varepsilon_s\nabla\bar c\right)+E-D

EE 는 하상 픽업(연행), D=wscbD=w_s c_b 는 침강 플럭스다. 두 항의 차이가 하상변동에 그대로 들어간다. 이 정식화의 수치적 성격은 강성(stiff)이다. 적응길이 La∼α uh/wsL_a\sim \alpha\,uh/w_s 가 격자 크기보다 훨씬 짧아지면 E−DE-D 가 서로 상쇄되는 큰 두 수의 차가 되어, 양해법으로 풀면 시간스텝이 h/wsh/w_s 수준으로 잘린다. 그래서 픽업/침강 항은 관례적으로 음해로 처리한다.

6. 엑스너 방정식[편집]

하상변동의 지배식은 결국 토사에 대한 질량보존 한 줄이다. 공극률 pp 인 하상 퇴적층에 대해

(1−p)∂zb∂t+∂qb∂x=D−E(1-p)\frac{\partial z_b}{\partial t}+\frac{\partial q_b}{\partial x}=D-E

이것이 엑스너 방정식(Exner, 1925)이다. 소류사만 있으면 우변은 0이고, qbq_b 가 흐름을 통해 zbz_b 에 의존하므로 이 식은 비선형 이류방정식이다. 즉 하상 프로파일에 충격파(하상 계단)와 팽창파가 생길 수 있고, 리만 문제의 언어가 그대로 적용된다. 모래의 pp 는 대개 0.4 정도이며, 이 계수 하나가 하상변동 속도를 곧바로 스케일링한다.

여기서 이 분야 전체를 설명하는 스케일 하나를 뽑을 수 있다. 유량이 준정상이라 보고 qb=Aumq_b=Au^m 으로 두면, 엑스너를 선형화해서 얻는 하상파 전파속도는

cb=11−p⋅uh⋅dqbdu⋅11−Fr2=m qb(1−p) h (1−Fr2)c_b=\frac{1}{1-p}\cdot\frac{u}{h}\cdot\frac{\mathrm{d}q_b}{\mathrm{d}u}\cdot\frac{1}{1-\mathrm{Fr}^2} =\frac{m\,q_b}{(1-p)\,h\,(1-\mathrm{Fr}^2)}

이다. 여기서 세 가지가 한꺼번에 읽힌다.

  • 평범한 하천에서 하상은 물보다 4자릿수 느리다. qb=10−4 m2/sq_b=10^{-4}\,\mathrm{m^2/s}, h=2 mh=2\,\mathrm{m}, m=3m=3, Fr=0.5\mathrm{Fr}=0.5 를 넣으면 cb≈3×10−4 m/sc_b\approx3\times10^{-4}\,\mathrm{m/s}, 즉 하루에 약 29 m 다. 물의 파속 gh≈4.4 m/s\sqrt{gh}\approx4.4\,\mathrm{m/s} 의 만분의 일이 안 된다.
  • Fr→1\mathrm{Fr}\to1 에서 발산한다. 분모의 1−Fr21-\mathrm{Fr}^2 가 그 자리에 있다. 상류(常流)-사류 전이 근처에서는 하상파와 수면파의 속도가 같은 자릿수가 되고, 뒤에 볼 약결합이 물리적으로 무의미해진다.
  • 사류에서는 부호가 뒤집힌다. Fr>1\mathrm{Fr}>1 이면 cb<0c_b<0 — 하상 기복이 상류로 거슬러 올라간다. 급류에서 관찰되는 역사구(antidune)가 정확히 이것이며, 사구(dune)와 역사구를 가르는 것이 프루드 수라는 관측 사실이 이 한 줄에서 나온다. 프루드 수가 개수로 유사 이송에서도 왕 노릇을 하는 이유다.

7. 흐름-하상 결합의 수치[편집]

약결합(decoupled)은 한 시간스텝 동안 하상을 고정한 채 얕은 물 방정식을 풀고, 그 유동장으로 qbq_b 를 계산해 엑스너를 한 스텝 진행한 뒤 다음으로 넘어간다. 위 스케일 논증이 이 전략의 정당화다. 구현이 쉽고 기존 수리 코드에 유사 모듈을 얹기만 하면 되며, 하천 실무의 압도적 다수가 이 방식이다.

강결합(coupled)은 (h,hu,zb)(h,hu,z_b) 를 하나의 쌍곡선계로 놓고 3×3 야코비안의 고유구조로 리만 문제를 푼다. 비용은 비싸지만 다음 상황에서 약결합이 실제로 깨진다.

  • Fr≈1\mathrm{Fr}\approx1 근처 — 위에서 본 공진. 하상파 고유값이 수면파 고유값에 근접해 분리 오차가 O(1)O(1) 이 된다. 분할 기법이 왜 실패하는지, 쌍곡성이 어떻게 도와주는지는 코르디에·르·모랄레스(2011)의 논문 제목 자체가 요약하고 있다.
  • 이송량이 클 때 — 댐 붕괴 유사류, 토석류. qbq_b 가 커지면 하상 고유값이 더 이상 작지 않다.
  • 급경사 하상 — 하상 경사가 안식각(모래 약 32°)을 넘으면 물리적으로 사면 붕괴가 일어나야 한다. 엑스너 자체는 이를 모르므로 별도의 사면활동(avalanching) 연산자를 매 스텝 돌려 경사를 안식각 이하로 눌러 준다. 이걸 안 넣으면 협착부나 세굴공 주변에서 수직 벽이 서고, 그다음 스텝에 이송량이 폭주한다.

이 밖에 유사 이송 코드를 실제로 깨뜨리는 것들.

  • 웰밸런스드성. 정지수 위에서 하상이 울퉁불퉁하면 수치 플럭스와 하상 경사 원천항이 정확히 상쇄되어야 하는데, 그렇지 않으면 가만히 있어야 할 물이 흐르고 그 가짜 유속이 곧바로 가짜 이송량을 만든다. 유사가 얹히면 웰밸런스드 도식은 선택이 아니라 필수가 된다.
  • 수치확산. 부유사 이송은 페클레 수가 큰 이류 지배 문제라 1차 풍상차분을 쓰면 수치분산으로 농도 전선이 뭉개진다. 그렇다고 고차 중심차분을 쓰면 농도가 음수로 내려가 다음 스텝에 침강 플럭스가 부호를 뒤집는다. 전변분 감소 제한자나 MUSCL 재구성이 붙는 것은 정확도가 아니라 비음성 보장을 위해서다.
  • 형태학적 가속계수(MORFAC). 하상 변화가 너무 느려 실시간으로 100년을 돌릴 수 없다. 그래서 하상변동률에만 배수 ff 를 곱해 흐름 1년으로 하상 ff 년을 모의한다. Delft3D의 시그니처 기능인데, ff 를 너무 키우면 하상이 흐름 응답보다 빨라져 위 스케일 분리가 스스로 무너진다. 모델러가 자기 손으로 Fr→1\mathrm{Fr}\to1 상황을 만드는 셈이라 민감도 검토가 필수다.

8. 무엇이 나오는가[편집]

  • 사구(dune). 하상 기복이 저절로 자란다는 것은 선형 안정성 문제이며, 핵심은 위상 지연이다. 하상 전단응력의 최대점이 사구 마루와 정확히 일치하면 기복은 자라지도 줄지도 않는다. 흐름의 관성과 이송의 반응 지연 때문에 최대점이 마루보다 약간 상류로 밀리면 마루가 성장한다(케네디 1963, 엥글룬 1970 계열). 그래서 국소 평형 이송식만 쓰는 모델은 원리적으로 사구를 못 만든다 — 지연을 넣어야 나온다.
  • 하상 저하와 armoring. 댐을 세우면 하류로 내려오던 유사 공급이 끊긴다. 방류수는 여전히 이송능력을 갖고 있으므로 자기 하상을 파먹기 시작하고(“hungry water”), 잔 입자부터 골라 쓸어 간다. 남은 굵은 입자가 표층에 깔려 하상을 덮으면 장갑화(armoring)가 일어나 저하가 스스로 멈춘다. 댐 직하류에서 수 m 규모의 하상 저하가 관측되다가 수십 년 뒤 안정되는 전형적 궤적이 이것이다. 모델링하려면 단일입경으로는 안 되고 입도별 다층 하상(활성층·저장층) 부기가 필요하다.
  • 저수지 퇴사·교각 세굴·해안 표사. 각각 세류사, 국소 3차원 와류, 파랑-흐름 연성이 주역이라 1차원 하상 부기로는 접근이 안 된다. 교각 세굴은 아직도 경험식(HEC-18)이 실무 표준이고, 난류 해상 CFD로 대체하려는 시도는 검증 단계다.

9. 실무 코드[편집]

코드차원성격
HEC-RAS Sediment1D(2D 추가)준부정류 + 입도별 하상 부기. 미국 하천 실무 표준
Delft3D / D-Morphology2D·3D반 레인 계열, MORFAC. 하구·해안의 국룰
TELEMAC-GAIA2D·3D구 SISYPHE 후속. 유럽 오픈소스 계보
MIKE 21 ST2D해안·항만 상용
SRH-2D · CCHE2D2D미개간지·연구용 하천 모형

전부 공통적으로 입력 불확실성이 수치오차를 압도한다. 입도 분포, 상류 유사 공급량, 이송 공식 선택, 하상 활성층 두께 — 어느 하나도 잘 모르는 채로 들어간다. 그래서 이 분야의 검증 및 확인은 “정답과 얼마나 가까운가”보다 **“과거 관측 하상변동을 재현하도록 보정한 뒤 그 보정이 다른 사상에도 버티는가”**를 묻는 형태로 굳어졌다. 절대 예측이 아니라 대안 비교(준설을 여기 할까 저기 할까)에 쓰라는 지하수 유동 모델링의 교훈이 여기서도 그대로 반복된다.

10. 관련 문서[편집]

11. Footnotes[편집]

  1. 루스 수는 정확히 “떠 있으려는 힘 대 가라앉으려는 힘”의 비다. κu∗\kappa u_* 가 난류 연직 혼합의 대표 속도이므로, PP 가 1보다 크면 침강이 이기고 작으면 난류가 이긴다. 무차원수 하나로 세 이송양식을 다 가르는 것이 지나치게 편해 보이지만, 실제로는 경계값이 문헌마다 0.8/1.2/2.5로 다르게 적혀 있고 어느 표를 인용했느냐에 따라 보고서의 결론 문장이 바뀐다. 무차원수는 정직한데 그걸 쓰는 사람이 정직하지 않은 전형적 사례. ↩

  2. 이런 일이 벌어지는 이유는 단순하다. 원논문은 1948년 스톡홀름 IAHR 회보에 실린 취리히 수리연구소 보고이고, 교과서에 인용되는 과정에서 리플 팩터가 “실무상 1로 두어도 무방”이라는 각주와 함께 사라졌다. 그리고 각주는 다음 세대 교과서에서 또 사라졌다. 계수 8은 그렇게 60년을 살아남았다. 코드 안에서 이송량을 2배 부풀리는 상수는 대개 이런 식으로 태어난다. ↩

  3. 상승부와 하강부의 농도가 다르므로, 유량 하나에 유사량이 두 개 대응한다. 그래서 유량-유사량 관계식(sediment rating curve)을 회귀로 한 줄 그으면 연간 총유사량이 계통적으로 어긋난다. 관측점 산포도가 소시지 모양으로 벌어져 있는데 그 한가운데로 직선을 긋고 결정계수를 자랑하는 그래프를 보면, 그건 히스테리시스를 오차로 착각한 것이다. ↩