1．はじめに
今回は『表面流モデル（kinematic runnoff model）』について紹介したいと思います．
表面流モデルは短期流出解析によく用いられ，大出水時の流出量をうまく再現できるそうです．
ちなみに表面流モデルは次式で表されます．
※式中の各記号については後で説明します．
\[
\begin{align}
&\dfrac{\partial h}{\partial t}+\dfrac{\partial q}{\partial x}=q _{*} &(連続式)\\
&h=K q ^{m} &(運動方程式)
\end{align}
\]
\[
\begin{align}
dq(x,t)&=\dfrac{\partial q(x,t)}{\partial x}dx+\dfrac{\partial q(x,t)}{\partial t}dt \\
\dfrac{dq(x,t)}{dt}&=\left(\dfrac{dx}{dt} \right) \dfrac{\partial q(x,t)}{\partial x}+\dfrac{\partial q(x,t)}{\partial t} \\
\end{align}
\]
\( \tfrac{dx}{dt} \) は\( x〜t \)平面上における擾乱の進行経路を表す曲線で特性曲線と呼ばれる．
また，この式は1階波動方程式といわれる．
1階の波動方程式をみたす解をkinematic waveという．
図に示すような斜面を流れる表面流の単位幅流量を\( q(x,t) \)，斜面上流端からの距離を\( x \)，降水量を \( q_{*}(t) \)として考える．
このとき，特性曲線\( \tfrac{dx}{dt} \)は次のグラフに示すような形状になる．
（縦軸，横軸の数値は便宜的なもの）
雨を球で表現し，降雨，表面流，特性曲線を動画で表すと以下のようになる．
定常降雨の場合，斜面下流端における流量は上昇域，定常域，下降域に分けられる．
上昇域：最初に斜面上流端に降った雨が斜面下流端に到達するまでの間
定常域：最初に斜面上流端に降った雨が斜面下流端に到達してから雨が止むまでの間
下降域：雨が止んでから流出が終了するまでの間
ただし，動画はわかりやすさを重視して，厳密性を欠いているため注意．
3.2．連続式
斜面を流れる雨水流のうち微小区間\( dx \)に注目する．
位置\( x \)の断面を通過する単位幅流量を\( q(x,t) \)とし，位置\( x \)，時刻\( t \)における水深を\( h(x,t) \)とする．
斜面幅\( B \)を1（単位幅），降水量を\( q_{*}(t) \)とする．
水路床勾配\( i \)が十分に小さいと仮定する（\( i \fallingdotseq sin \theta \fallingdotseq tan \theta \)）．
説明の都合上，位置\( x \)の上下流方向それぞれに\( \tfrac{1}{2}dx \)区間をとることにする．
テイラー展開より\( q(x+\tfrac{1}{2}dx,t) \)は次式になる．
\[
\begin{align}
q(x+\tfrac{1}{2}dx,t) &= q(x,t) + \dfrac{\partial q(x,t)}{\partial x}\dfrac{1}{2}dx + \dfrac{\partial ^{2} q(x,t)}{\partial x ^{2}} \left(\dfrac{1}{2}dx \right) ^{2} + \cdot \cdot \cdot \\
&\fallingdotseq q(x,t) + \dfrac{\partial q(x,t)}{\partial x}\dfrac{1}{2}dx
\end{align}
\]
同様に，テイラー展開より\( q(x-\tfrac{1}{2}dx,t) \)，\( h(x,t+dt) \)はそれぞれ次式になる．
\[
\begin{align}
&q(x-\tfrac{1}{2}dx,t) = q(x,t) - \dfrac{\partial q(x,t)}{\partial x}\dfrac{1}{2}dx \\
&h(x,t+dt) = h(x,t) + \dfrac{\partial h(x,t)}{\partial t}dt
\end{align}
\]
微小区間\( dx \)における時刻\( t \)，\( t+dt \)それぞれの水位の代表値が\( h(x,t) \)，\( h(x,t+dt) \)とすると，微小区間\( dx \)の水収支は
\[
\begin{align}
\left( h(x,t+dt) -h(x,t) \right) \cdot dx \cdot 1 &= q(x-\dfrac{1}{2}dx,t) \cdot dt -q(x+\dfrac{1}{2}dx,t) \cdot dt + q _{*}(t) \cdot dx \cdot 1 \cdot dt\\
\dfrac{\partial h(x,t)}{\partial t}dtdx &= -\dfrac{\partial q(x,t)}{\partial x}dxdt + q _{*}(t)dxdt \\
\dfrac{\partial h(x,t)}{\partial t} + \dfrac{\partial q(x,t)}{\partial x} &= q _{*}(t)
\end{align}
\]
これが表面流モデルの連続式である．
独立変数を書くと式が見えづらいので，省略して一般に次のように表記する．
\[
\dfrac{\partial h}{\partial t} + \dfrac{\partial q}{\partial x} = q _{*}
\]
3.3．運動方程式
水路床勾配\( i \)が十分に小さいと仮定する（\( i \fallingdotseq sin \theta \fallingdotseq tan \theta \)）．
微小区間\( dx \)の水の質量を\( m \)，微小体積を\( dV \)とし，この微小体積に作用する重力を\( W \)，摩擦力を\( F_{f} \)，全水圧を\( P \)とする．
時刻\( t \)に位置\( x \)を通過する流れの速度\( v(x,t) \)とすると，その全微分は次式になる．
\[
dv(x,t)=\dfrac{\partial v(x,t)}{\partial x}dx+\dfrac{\partial v(x,t)}{\partial t}dt\\
\]
よって，時刻\( t \)に位置\( x \)を通過する流れの加速度は次式になる．
\[
a(x,t)=\dfrac{dv(x,t)}{dt}=\dfrac{\partial v(x,t)}{\partial x}\dfrac{dx}{dt} + \dfrac{\partial v(x,t)}{\partial t}\\
\]
位置\( x \)時刻\( t \)における全水圧を\( P(x,t) \)とするとテイラー展開より
\[
\begin{align}
&P(x-\tfrac{1}{2}dx,t) = P(x,t) - \dfrac{\partial P(x,t)}{\partial x}\dfrac{1}{2}dx \\
&P(x+\tfrac{1}{2}dx,t) = P(x,t) + \dfrac{\partial P(x,t)}{\partial x}\dfrac{1}{2}dx \\
\end{align}
\]
水の密度を\( \rho \)，水路床と水の間に働く摩擦応力を\( \tau \)，潤辺を\( p_{w}(x,t) \)，重力加速度を\( g \)とすると，微小区間\( dx \)おける水の運動方程式は次式になる．
\[
\begin{align}
ma(x,t) &= W - F_{f} + P(x-\dfrac{1}{2}dx,t)-P(x+\tfrac{1}{2}dx,t) \\
\rho \cdot dV \left( \dfrac{\partial v(x,t)}{\partial x}\dfrac{dx}{dt} + \dfrac{\partial v(x,t)}{\partial t} \right) &= \rho \cdot dV \cdot g \cdot sin \theta - \tau \cdot dx \cdot p_{w}(x,t) + \dfrac{\partial P(x,t)}{\partial x}dx \\
\rho dV \left( \dfrac{\partial v(x,t)}{\partial x}\dfrac{dx}{dt} + \dfrac{\partial v(x,t)}{\partial t} \right) &= \rho dV g sin \theta - \tau dx p_{w}(x,t) + \rho g A(x,t) \dfrac{\partial h(x,t)}{\partial x}dx\\
\end{align}
\]
となる．
\( \tfrac{\partial P(x,t)}{\partial x} \)の導出方法は次の通りである．
水圧が静水圧分布すると仮定すると，水路床から図心までの高さ\( \overline{z}(x,t) \)は次式になる．
\[
\overline{z}(x,t) = \dfrac{1}{A(x,t)} \int _{0} ^{h(x,t)}z(x,t)f(z(x,t))dz(x,t)
\]
ちなみに式が見やすいように\( x \)と\( t \)を省略すると次式になる．
\[
\overline{z} = \dfrac{1}{A} \int _{0} ^{h}zf(z)dz
\]
\[
\begin{align}
\dfrac{\partial A}{\partial h} &= \dfrac{\partial }{\partial h} \int _{0} ^{h} f(z)dz = T \\
\end{align}
\]
\[
\begin{align}
\dfrac{\partial \overline{z}}{\partial h} &= \dfrac{\partial \overline{z}}{\partial A} \dfrac{\partial A}{\partial h} \\
&= \left(\dfrac{\partial }{\partial A} \left( \dfrac{1}{A} \right) \int _{0} ^{h} zf(z)dz + \dfrac{1}{A} \dfrac{\partial }{\partial A} \left( \int _{0} ^{h} zf(z)dz \right)\right) T\\
&= -\dfrac{1}{A ^ {2}} A T \overline{z} + \dfrac{T}{A} \dfrac{\partial }{\partial A} \left( \int _{0} ^{h} zf(z)dz \right) \\
\end{align}
\]
ここで，\( zf(z) \)の原始関数を\( G(z) \)とおくと，
\[
\begin{align}
\dfrac{\partial \overline{z}}{\partial h} &= -\dfrac{1}{A ^ {2}} A T \overline{z} + \dfrac{T}{A} \dfrac{\partial }{\partial A} \left( G[h]-G[0] \right) \\
&= -\dfrac{T}{A} \overline{z} + \dfrac{T}{A} \dfrac{\partial h}{\partial A} \dfrac{d \left( G[h]-G[0] \right)}{dh} \\
&= -\dfrac{T}{A} \overline{z} + \dfrac{T}{A} \dfrac{1}{T} hf(h)\\
\end{align}
\]
となる．
\( f(h) \)の原始関数を\( F(h) \)とおくと，
\[
\begin{align}
A &= \int _{0} ^{h}f(z)dz \\
&= F(h)-F(0) \\
\dfrac{\partial A}{\partial h} &= \dfrac{dF(h)}{dh} -\dfrac{dF(0)}{dh} \\
\dfrac{\partial A}{\partial h} &= f(h) = T \\
\end{align}
\]
\[
\begin{align}
\dfrac{\partial \overline{z}}{\partial h} &= -\dfrac{T}{A} \overline{z} + \dfrac{T}{A} \dfrac{1}{T} hf(h) \\
&= -\dfrac{T}{A} \overline{z} + \dfrac{T}{A} \dfrac{1}{T} hT \\
&= \dfrac{T}{A} (h-\overline{z}) \\
\end{align}
\]
\[
\begin{align}
\dfrac{\partial P(x,t)}{\partial x} &= \dfrac{\partial P}{\partial h} \dfrac{\partial h}{\partial x} \\
&= \dfrac{\partial }{\partial h} \left( \rho g (h- \overline{z})A \right) \dfrac{\partial h}{\partial x} \\
&= \left( \rho g A \left( 1-\dfrac{\partial \overline{z}}{\partial h}\right) + \rho g (h - \overline{z}) \dfrac{\partial A}{\partial h}\right) \dfrac{\partial h}{\partial x} \\
&= \left( \rho g A - \rho g A \dfrac{T}{A}(h - \overline{z}) + \rho g h T -\rho g \overline{z} T \right) \dfrac{\partial h}{\partial x} \\
&= \rho g A \dfrac{\partial h}{\partial x} \\
\end{align}
\]
潤辺\( p(x,t) \)は1である．
両辺を\( \rho dV g (=\rho A(x,t)dxg) \)で割ると，
\[
\begin{align}
\rho dV \left( \dfrac{\partial v(x,t)}{\partial x}\dfrac{dx}{dt} + \dfrac{\partial v(x,t)}{\partial t} \right) &= \rho dV g sin \theta - \tau dx p_{w}(x,t) + \rho g A(x,t) \dfrac{\partial h(x,t)}{\partial x}dx \\
\dfrac{1}{g} \left( \dfrac{\partial v(x,t)}{\partial x}\dfrac{dx}{dt} + \dfrac{\partial v(x,t)}{\partial t} \right) &= sin \theta - \dfrac{\tau dx p_{w}(x,t)}{\rho A(x,t) dx g} + \dfrac{\rho g A(x,t) dx}{\rho g A(x,t) dx} \dfrac{\partial h(x,t)}{\partial x} \\
\dfrac{1}{g} \left( v(x,t)\dfrac{\partial v(x,t)}{\partial x} + \dfrac{\partial v(x,t)}{\partial t} \right) &= i - \dfrac{\tau}{\rho h(x,t) g} + \dfrac{\partial h(x,t)}{\partial x} \\
\end{align}
\]
斜面勾配\( i \)と摩擦勾配\( S_{f} \)が卓越しており，他項が省略できるとすると，
\[
i = \dfrac{\tau}{\rho h(x,t) g} = S_{f} \\
\]
摩擦
表面流モデルのように非定常不等流の場合でも以下のように定常等流と同様に振舞うと一般に仮定される．
摩擦（層流）
以下に示す層流における摩擦力の式は定常・等流の条件下で成り立つ．
微小区間\( dx \)において，ある時刻\( t \)に注目し，定常・等流条件が成り立つと仮定して考える．
微小区間\( dx \)における河床からの高さ\( z \)における断面において，重力の流れ方向成分は次式になる．
\[
\rho (h-z)g dx sin \theta
\]
ニュートンの粘性方程式より，せん断力\( \tau (z) \)は次式になる．
\[
\tau (z) = \mu \dfrac{du}{dz}dx = \rho \nu \dfrac{du}{dz}dx
\]
ただし，\( u \)：水平流速，\( \mu \)：粘性係数，\( \nu \)：動粘性係数
河床からの高さ\( z \)における断面において，重力の流れ方向成分とせん断力がつりあう．
\[
\begin{align}
\rho (h-z)g dx sin \theta &= \rho \nu \dfrac{du}{dz}dx \\
i(h-z)g dz &= \nu du \\
\int _{0} ^{z} i(h-z)g dz &= \int _{u(0)} ^{u(z)}\nu du \\
\int _{0} ^{z} i(h-z)g dz &= \int _{u(0)} ^{u(z)}\nu du \\
i \left( hz-\dfrac{1}{2}z ^{2} \right) g &= \nu u(z) \\
u(z) &= \dfrac{hgi}{\nu}\left( z-\dfrac{z ^{2}}{2h} \right) \\
\end{align}
\]
平均流速\( \overline{u} \)を求める．
\[
\begin{align}
\overline{u} &= \dfrac{1}{A}\int _{A} u(z)dA \\
\overline{u} &= \dfrac{1}{A}\int _{0} ^{h} \dfrac{hgi}{\nu}\left( z-\dfrac{z ^{2}}{2h} \right) \cdot 1 \cdot dz \\
\overline{u} &= \dfrac{hgi}{\nu A} \left( \dfrac{1}{2} [z ^{2}] _{0} ^{h}-\dfrac{1}{6h}[z ^{3}] _{0} ^{h} \right) \\
\overline{u} &= \dfrac{hgi}{\nu A} \left( \dfrac{1}{2} h ^{2}-\dfrac{1}{6h}h ^{3} \right) \\
\overline{u} &= \dfrac{h ^{3} gi}{3 \nu A} \\
\end{align}
\]
\(i=S_{f} \)，\( A=1 \cdot h \)より
\[
\begin{align}
S_{f} &= \dfrac{3 \nu A}{h ^{3} g} \overline{u} \\
&= \dfrac{3 \nu A}{h ^{3} g} \dfrac{q}{A}\\
&= \dfrac{3 \nu q}{h ^{3} g}\\
\end{align}
\]
\( S_{f} = sin \theta \)より
\[
\begin{align}
S_{f} &= sin \theta = \dfrac{3 \nu q}{h ^{3} g} \\
h &= \left( \dfrac{3 \nu}{g sin \theta} \right) ^{\tfrac{1}{3}} q^{\tfrac{1}{3}} \\
\end{align}
\]
摩擦（乱流）
以下に示す乱流における摩擦力の式は等流の条件下で成り立つ．
微小区間\( dx \)において定等流条件が成り立つと仮定すると，動水勾配\( I= \)摩擦勾配\(S_{f}(=i=sin \theta) \)とみなすことができる．
平均流速を\( v \)，粗度係数を\( n \)とすると，マニングの公式より
\[
\begin{align}
v &= \dfrac{1}{n} R^{\tfrac{2}{3}} I^{\tfrac{1}{2}} \\
\dfrac{Q}{A} &= \dfrac{1}{n} h^{\tfrac{2}{3}} \left( sin \theta \right) ^{\tfrac{1}{2}} \\
\dfrac{q \cdot 1}{h \cdot 1} &= \dfrac{1}{n} h^{\tfrac{2}{3}} \sqrt{sin \theta} \\
h^{\tfrac{5}{3}} &= \dfrac{n}{\sqrt{sin \theta}} q \\
h &= \left( \dfrac{n}{\sqrt{sin \theta}} \right) ^{\tfrac{3}{5}} q^{\tfrac{3}{5}} \\
\end{align}
\]
層流と乱流の摩擦をまとめると次のように表すことができる．
\[
h = K q^{m}
\]
ただし，\( K \)と\( m \)は次式の通りである．
層流の場合：\( K = K_{L} = \left( \dfrac{3 \nu}{g sin \theta} \right) ^{\tfrac{1}{3}} \)，\( m=\dfrac{1}{3} \)
乱流の場合：\( K = K_{T} = \left( \dfrac{n}{\sqrt{sin \theta}} \right) ^{\tfrac{3}{5}} \)，\( m=\dfrac{3}{5} \)
これが，表面流モデルの運動方程式である．
連続式と運動方程式を独立変数を省略せずに書くと次のようになる．
\[
\begin{align}
\begin{cases}\dfrac{\partial h(x,t)}{\partial t}+\dfrac{\partial q(x,t)}{\partial x}=q _{*}(t) & (連続式)\\ \\ h(x,t)=K q ^{m}(x,t) & (運動方程式) \end{cases}
\end{align}
\]
3.4．単位幅流量と斜面上流端からの距離
運動方程式を\( t \)で偏微分する
\[
\begin{align}
h &= K q^{m} \\
\dfrac{\partial h}{\partial t} &= K \dfrac{\partial q^{m}}{\partial t} \\
&= K \dfrac{\partial q^{m}}{\partial q} \dfrac{\partial q}{\partial t} \\
&= mK q^{m-1} \dfrac{\partial q}{\partial t} \\
\end{align}
\]
これを連続式\( \dfrac{\partial h}{\partial t} + \dfrac{\partial q}{\partial x} = q _{*} \)に代入する．
\[
\begin{align}
\dfrac{\partial h}{\partial t} + \dfrac{\partial q}{\partial x} &= q _{*} \\
mK q^{m-1} \dfrac{\partial q}{\partial t} + \dfrac{\partial q}{\partial x} &= q _{*} \\
\dfrac{\partial q}{\partial t} + \left( \dfrac{q^{1-m}}{mK}\right) \dfrac{\partial q}{\partial x} &= \left( \dfrac{q^{1-m}}{mK}\right)q _{*} \\
\end{align}
\]
これを\( \dfrac{dq}{dt} =\left(\dfrac{dx}{dt} \right) \dfrac{\partial q}{\partial x}+\dfrac{\partial q}{\partial t} \)と比較すると
\[
\begin{cases}\dfrac{dx}{dt}=\dfrac{q^{1-m}}{mK}\\ \dfrac{dq}{dt}= \dfrac{q^{1-m}}{mK}q _{*} \end{cases}
\]
時刻\( t=\tau \)から出発する特性曲線について解くと，\( q(x,t) \)は次のようになる．
\[
\begin{align}
\dfrac{dq}{dt} &= \dfrac{q^{1-m}}{mK}q _{*} \\
mq^{m-1}dq &= \dfrac{q _{*} }{K}dt \\
\int _{q(0,\tau)} ^{q(x,t)} mq^{m-1}dq &= \int _{\tau} ^{t} \dfrac{q _{*} }{K}dt \\
q^{m}(x,t) - q^{m}(0,\tau) &= \int _{\tau} ^{t} \dfrac{q _{*} }{K}dt \\
\end{align}
\]
\( q(0,\tau) \)は時刻\( \tau \)の斜面上流端の単位幅流量のため\( q(0,\tau)=0 \)となるので，
\[
\begin{align}
q^{m}(x,t) - q^{m}(0,\tau) &= \int _{\tau} ^{t} \dfrac{q _{*} }{K}dt \\
q^{m}(x,t) &= \int _{\tau} ^{t} \dfrac{q _{*} }{K}dt \\
q(x,t) &= \left\{ \int _{\tau} ^{t} \dfrac{q _{*} }{K}dt \right\} ^{ \tfrac{1}{m}} \\
\end{align}
\]
\[
\begin{align}
\dfrac{dx}{dt} &= \dfrac{q^{1-m}}{mK} \\
dx &= \dfrac{1}{mK}q^{1-m}dt \\
\int _{0} ^{x} dx &= \int _{\tau} ^{t} \dfrac{1}{mK}q^{1-m}dt \\
\end{align}
\]
\( q(x,t) \)の\( t \)を便宜的に別のアルファベット\( s \)に置き換えると，
\[
\begin{align}
\int _{0} ^{x} dx &= \int _{\tau} ^{t} \dfrac{1}{mK}q^{1-m}dt \\
\int _{0} ^{x} dx &= \int _{\tau} ^{t} \dfrac{1}{mK} \left\{ \int _{\tau} ^{s} \dfrac{q _{*} }{K}dt' \right\} ^{ \tfrac{1-m}{m}} ds \\
x &= \dfrac{1}{m K ^\tfrac{1}{m}} \int _{\tau} ^{t} \left\{ \int _{\tau} ^{s} q _{*}dt' \right\} ^{ \tfrac{1-m}{m}} ds \\
\end{align}
\]
3.5．おわりに
本記事の内容をまとめると以下のようになります．
表面流モデルの連続式：
\[
\dfrac{\partial h}{\partial t}+\dfrac{\partial q}{\partial x}=q _{*}
\]
\[
h=K q ^{m}
\]
\[
q(x,t) = \left\{ \int _{\tau} ^{t} \dfrac{q _{*} }{K}dt \right\} ^{ \tfrac{1}{m}}
\]
\[
x = \dfrac{1}{m K ^\tfrac{1}{m}} \int _{\tau} ^{t} \left\{ \int _{\tau} ^{s} q _{*}dt' \right\} ^{ \tfrac{1-m}{m}} ds
\]
