2. 厚真川小流域の崩壊斜面からの流出土砂の観測及び厚真川流域
における地震前後の土砂侵食リスクの定量的な評価
宮崎功太郎1・星野
剛2・山田 朋人3・泉典洋4・土田 宏一5
Kotaro MIYAZAKI, Tsuyoshi HOSHINO, Tomohito YAMADA,Norihiro IZUMI and Koichi TSUCHIDA
1

北海道大学大学院工学院河川流域工学研究室 修士1年
北海道大学大学院工学研究院土木工学部門河川流域工学研究室 博士研究員
3
北海道大学大学院工学研究院土木工学部門河川流域工学研究室 准教授
4北海道大学大学院工学研究院土木工学部門河川流域工学研究室 教授
5（株）福田水文センター
2

要旨

2018年9月6日に発生した北海道胆振東部地震により，胆振地方では多数の土
砂崩れが発生した．崩壊面積は明治以降で我が国最大となる48.4 km2にも及び，崩壊斜面から厚真川への土砂流出量の推定は喫緊の課題となっている．本研究の目的は1イベントの降雨による厚真川富里地区の一支川からの土砂流出量の観測による算出，および厚真川流域全体の経験式による年土砂侵食率[t・ha-1・
yr-1]（1 haあたりの1年間の土砂侵食量）の増加率の定量化である．土砂流出量は，観測対象流域の末端での流量，濁度の観測によって算出した．2019年から2020年の7つの降雨イベントを対象に調べた結果，総降雨量は同程度であっても，土砂流出量には大きな違いが生じることが明らかになった．また，地震前後の年間土砂侵食率を算出した結果，流域全体の総土砂侵食率は地震前と比べて8倍から29倍の範囲で増加することがわかり，特に崩壊斜面が多くの面積を占めるサブ流域においては増加が顕著であることが明らかとなった．地震による裸地化と地形変化それぞれの土砂侵食率への寄与度を調べた結果，土砂崩
壊による裸地化が土砂侵食率の増大に支配的であり，地形変化は各斜面内において侵食の増大と減少のどちらの影響ももたらしうることがわかった．さらに，厚真川流域の降雨の年変動を考慮したところ，地震後の年最小土砂侵食量
は，地震前の年最大土砂侵食率を上回ることが示された．これは地震後の厚真川流域では，地震前の年最大土砂侵食率以上の土砂侵食が毎年のように発生することを示唆する結果である．
《キーワード：厚真川；土砂崩壊；土砂流出；現地観測；年土砂侵食率；降雨変動》― 184 ―

2. 1. はじめに
2018年9月6日に発生した北海道胆振東部地震により，北海道厚真町，安平町，むかわ町を中心に大規模な土砂崩壊が生じた（図-1（c））．Kasai and Yamada（2019）によると，崩壊面積は我が国において明治以降最大となる48.4 km2であり，土砂崩壊の発生箇所は6,117箇所にも及ぶ1）．関口・佐藤（2006）によると，本地震以前の既往最大の崩壊面積を記録した新潟県中越地震における崩壊面積は11.2 km2であり，本地震ではその約4倍の崩壊面積となる広域の土砂崩壊が引き起こされた2）．
特に崩壊が集中した厚真川流域では，降雨による土砂流出リスクの増大が予想される．厚真川流域に形成された崩壊斜面は流域面積の13 %を占めており，降雨による厚真川への流出土砂量の増加が考えられる．厚真川本川では，河床への土砂の堆積による洪水流下能力の低下や，厚真川上流の厚真ダムと厚幌ダムへの土砂流入による貯水容量の低下が懸念される．このような被害リスクを考える上で，厚真川流域の流出土砂量の定量的な評価が必須となる．
本研究では，地震後の厚真川流域の降雨と土砂流出の関係性を把握するために，厚真川の小流域を対象とした水文，気象観測を行い，2019年と2020年の降雨イベントによる流出土砂量の算出を行い，得られた結果より降雨と土砂流出との関係性の定量化を試みた．また，厚真川流域全体の土砂侵食リスクの増加率を評価するために，経験的な土砂侵食式を用いた地震前後の土砂侵食の増加率の定量化を行った．

水たまり図-1 （a）厚真川流域の全体図と土砂崩壊，堆積地の分布図．丸い点線は（b）の観測対象流域を示す．
（b）観測対象流域と観測地点．赤い枠線は（c）の土砂崩壊斜面を示す．
（c）観測対象流域の土砂崩壊斜面（2019年7月25日に著者らにより撮影）

2. 2. 調査対象河川の特徴
2. 2. 1. 厚真川流域及び観測対象支川流域の概要
解析対象流域は北海道南西部に位置する厚真川流域である（図-1（a））．厚真川は延長が52.3 km，流域面積が382.9km2であり，夕張山地から太平洋へ南西に向かって流れ，上流に厚真ダムと厚幌ダムを擁する．標高の範囲は16 mから641 mである．本流域の主な土地利用は，上流域が広葉針葉混交林であり，下流の平野部が水田，畑地，市街地である．土質は道央の火山噴火，カルデラ形成に由来する火山灰質土であり，下層から支笏カルデラ（Spfa-1），恵庭岳（En），樽前山（Ta-a，Ta-b，Ta-c，Ta-d）の順に堆積している．小山内ら（2019）は，本地震で発生した土砂崩壊は主にTa-dを滑り面として崩落していると指摘している3）．
水文観測の対象流域は，厚真町富里地区の厚真川の支川流域である（図-1（b））．この支川流域の流域
面積は国土地理院が提供する5 m標高データから0.76 km2と算出された．また，国土地理院の衛星画像から判別した崩壊斜面の面積は流域面積の54%にあたる0.41 km2であった．これは厚真川流域全体での― 185 ―

崩壊面積率である13%と比べて大きく，対象流域の土砂侵食は地震による土砂崩壊の影響を強く受けていることが予想される．著者らのドローンによる観測から，本流域において二箇所の窪地の水たまりの存在を確認した．

2. 2. 2. 観測手法，土砂流出量の算出
対象流域の水文観測および気象観測を行った．水文観測地点は厚真川本川と観測対象支川との合流部のコンクリート矩形水路である（図-2（a））．観測諸量は水位，水深，濁度，降雨量であり，それぞれ超音波水位計，濁度・図-2 （a）
気象観測地点の転倒ます雨量計．水文観
測地点から東へ500mに設置（b）
水文観測地点の超音波水位計，濁度・水深計，インターバルカメラ（2020年7月16日設置時撮影）

水深計（赤外後方散乱濁度計，
圧力式水深計），
転倒マス雨量計を設置して観測した（図-2（b））．なお濁度計と水深計は一つの測器に内蔵されている．降雨と水位の観測インターバルは5分に設定し，濁度と水深は5分毎に1秒間隔で10回計測を行った．超音波水位計の観測値は，観測開始時の水位からの相対的な水位の変化であるため，観測開始時の水深（2019年，2020年ともに5 cm）を現地で計測し，相対的な水位の変化を水深に変換した．
観測は2019年7月25日から11月13日，2020年7月16日から11月20日にかけて行われた．濁度計は2019年
7月25日から8月14日まで河床に，8月14日から9月13日まで河道中央の水面上方に設置されていたが，それぞれ出水時に流れてきた土砂にセンサーが覆われる，流下物がひっかかるなど計測に不適な状態になった．そのため2019年9月13日以降は左岸水面上方に設置し，2020年は全期間にわたり右岸水面上方に設置した．本観測では濁度と水深は平水時には計測されず，濁度計が浸かる高水時のみ計測される（図-3）．また，2020年観測期間中に濁度・水深計の電源が喪失したため，2020年8月11日から9月2日にかけて濁度，水深が欠測となっている．このように圧力式水深計は時間的に連続した水深を観測できていな
いため，本解析では，超音波水位計の水位を水深に換算した値を主に用いた．また，超音波水位計にもデータの欠損期間があるため，その期間の解析には高水時の圧力式水深計の水深を用いた．
観測された水深から流量を算出するために，河道近傍に30フレームレート（30枚/秒）のトレイルカ
メラ（昼間は可視カメラ，夜間は近赤外線カメラ）を設置し，2019年10月4日の出水時の表面流速[m/s]を推定した．推定手法はSpace-Time Image Velocimetry（STIV）法4）である．STIV法は主として一次元（主流方向）の平均流速を求める非接触型流速測定法であり，水表面の流れ方向に想定した検査線上の（a）平水時（2020年9月2日）
図-3（b）出水時（2020年9月5日）

水文観測サイトのインターバルカメラ画像― 186 ―

輝度分布の時間変化を時空間プロットし，そこに現れる縞パターンの勾配から表面流速を求めるものである．なお断面平均流速は表面流速に浮子観測（表面浮子）の更生係数である0.85を乗じて算出し5），
断面平均流速から流量Qを算出した．このようにして求められた流量Qに対応する水深hからh-√Q関係を作成し，連続的な流量Qを算出した．2019年10月4日イベントのピーク流量は5分間平均で0.20 m3/sと全イベントの中で比較的大きかったため，h-√Q関係を作成するためのイベントとして適当と判断した
（図-4
（a））．
赤外後方散乱濁度計の濁度は水中の浮遊物質からの赤外線の反射率で測定されているため，土砂量を求める際にキャリブレーションを行う必要がある．そこで反射率濁度-浮遊物質濃度関係を得るために2019年10月4日，2020年10月23日の出水時に採水を実施し，その分析を行った．10月4日イベントのピー
ク浮遊物質濃度は5分間平均で33.8 g/m3と全イベントの中で比較的大きかったため，反射率濁度-浮遊物
質濃度関係の作成ためのイベントとして適当である．得られた反射率濁度-浮遊物質濃度関係から連続的な浮遊物質濃度を算出した（図-4（b））．
以上述べた流量，浮遊物質濃度を乗じることで流砂量が得られ，さらに時間積分することによりイベントごとの総土砂流出量を算出した．

図-4 （a）2019年10月4日および2020年10月23日出水イベント時のh-√Qの散布図および線形近似
2019年10月4日および2020年10月23日出水イベント時の反射率濁度-浮遊物質濃度の散布図（b）
および線形近似

2. 3. 観測によるイベントごとの土砂流出量の算出結果
観測期間における観測結果の日別値を図-5に示す．観測期間中の最大総降雨イベントは2019年9月23日であり，総降雨量は63.8 mmである．水位は，降雨イベントの前と比べて降雨イベント後の方が小さく，あるいは大きくなるケースが多数見られることから，河床が洗掘される，あるいは流出土砂が堆積していることが示唆される．観測期間中最大の流出土砂量を記録した降雨イベントは2019年10月4日の降雨イベントであり，総降雨量44.8 mmに対し，流出土砂量は26.3 tであった（図-6）．
観測期間における主要な降雨イベントの流量-浮遊物質濃度関係を図-7（a）に示す．流量が同じでも浮遊物質濃度には大きな差が生じており，この差はイベント毎に異なるだけでなく，同じイベント中においても生じている．このことから，土砂流出量を求めるモデルの構築にあたり，流量だけでなく，浮遊物質濃度の変化も考慮しなければならないことが示唆される．また，イベント中の流出土砂量の時間変化は流量ピークよりも浮遊物質濃度のピークが早い初期高濃度型と濃度ピークが遅い後期高濃度型の2つに分類される．図-7（a）において，右回りのループを描いているイベントが初期高濃度型であり，― 187 ―

図-5

日別観測値（水位，水深，濁度，降水量）（a）2019年（b）2020年．
破線で囲われているイベントを対象に解析を行った．

図-6

2019年10月4日の降雨イベント時の観測値（水深，水位，濁度，降雨量）

左回りのループを描いているイベントが後期高濃度型である．初期高濃度型イベントは2019年9月23日
イベントのみであり，それ以外のイベントは後期高濃度型であった．2019年9月23日イベントは本観測期間中，最大総降水量を記録したイベントであること，また解析対象とされた7イベント中最初のイベントであることから，総降水量が流出土砂の時間変化に影響を及ぼすことや観測流域の土砂流出形態が9月23日イベントを境に不可逆的に変化した可能性がある．
各イベントで算出された総流出土砂量は0.2 tから26.3 tであった．流域面積で除すると0.3 g/m2から
34.6 g/m2となる．各イベントにおいて算出された総流出土砂量と総降雨量，平均降雨強度，1時間あたり最大降雨強度との散布図をそれぞれ図-7（b），（c），（d）に示す．総降水量と1時間あたり最大降雨強度は総流出土砂量と強い相関（相関係数:0.85，0.9）が見られるが，平均降雨強度と総土砂流出量の間の相関は弱い（相関係数:0.16）．土砂浸食は主に雨滴の衝撃による土砂粒子の剥離と，水流による浸食によって発生する．これら侵食の激しさは総降雨量，降雨強度だけでなく，降雨の時空間分布や流域の初期土壌水分量からも影響を受ける．したがって土砂流出量を求めるモデルの構築にあたり，降雨の時
空間分布および初期土壌水分量といった情報も考慮しなければならない可能性がある．

― 188 ―

図-7 （a）
流量-浮遊物質濃度のプロット．
初期高濃度型（右回りループ）のイベントは2019年9月23日のみであり，それ以外のイベン
トは後期高濃度型（左回りループ）のイベント（b）各イベントの総降雨量と総土砂流出量の散布図（c）各イベントの平均降雨強度と総土砂流出量の散布図（d）各イベントの1時間あたり最大降雨強度と総土砂流出量の散布図

2. 4. 経験的な土砂侵食式を用いた土砂侵食率の変化の定量化
2. 4. 1. 経験的な土砂浸食式（USLE）
降雨による土 砂 侵 食 を 推 定 す る モ デ ル には概念，物理，経験に基づく多くのモデルが存 在す
る6）, 7）, 8）, 9）, 10）．その中で，経験的なモデルであるUniversal Soil Loss Equation（USLE）は広く用いられる．
USLEは Wischmeier and Smith（1977）が提案した，農地における降雨による年間土壌侵食率を算定するために開発された経験式であり，5つの係数（地形，作物，降雨，土質，人為的な保全）で構成される11）．さらにWischmeier and Smith（1977）は，作物係数を樹冠，林床の植生，腐食層からなる補助係数で構成することを提案し，USLEの対象領域を山地森林へ拡張させた．また，USLEの係数は一般化された経験的な物理式に基づいて決定されるため，米国のみならず，各国の山地森林で適用されている．一方で，USLEはその制限の解消や精度の向上のために，現在にかけて改良が進められている12）,13）．
USLEにおいて土砂侵食率は降雨および対象地域の地形，植生，土質，人為的な保全から経験的に決― 189 ―

定される係数の乗算によって算出される．

𝐴𝐴 � �� 𝐶𝐶 𝑅𝑅 𝐾𝐾 𝑃𝑃（1）

ここで，A:年間土砂侵食率[t･ha−1･yr−1]，LS:地形係数[-]，C:作物係数[-]，R:降雨係数[MJ・mm・ha−1・
h−1・yr−1]，K:土壌係数[t・h・MJ−1・mm−1]，P:保全係数[-]である．USLEには制限が存在する．本式の
降雨-土壌侵食関係は米国における傾斜枠試験によって経験的に得られている．傾斜枠試験は傾斜した農地を囲う枠を設置し，降雨によって流出した土壌を枠下部に堆積させることで土壌侵食量を計測する．
降雨による土砂侵食の形式は雨滴の衝撃による剥離と水流による侵食に大別される．傾斜枠試験を実施した農地の流下方向の基準斜面長は22.13 mであり，これは水流による侵食の発生には不十分な距離である．加えて，一般的な流域における，土砂侵食が発生した地点から河川までの距離に比べて非常に短い．したがってUSLEは雨滴の衝撃による土砂侵食のみを対象としており，水流による土砂侵食を表現できない．また，河川のある地点における上流域からの土砂流出量も算出することはできない．しかし，そのような制限の解消，および精度の向上のために現在にかけて改良が進められており，本研究では改良された係数を用いる．
地形係数LSは斜面長係数Lと勾配係数Sの積である．斜面長は斜面において土砂侵食を伴う流れが発生してから，斜面の勾配が十分に減少して堆積するまでの距離と定義される．斜面長係数は斜面長が22.13 mである基準斜面の土砂侵食量に対する，ある斜面長かつ斜面長以外の条件が基準斜面と同じ斜面の土砂侵食量の比である．同様に勾配係数Sは勾配が9%である基準斜面の土砂侵食量に対する，ある
勾配かつそれ以外の条件が同じ斜面の土砂侵食量の比である．
斜面長の計測は対象地域が複雑な山岳地形である場合困難である．それを解消し，またガリ侵食よりも小さい水流による侵食であるリル侵食を考慮するために，Desmet and Govers（1996）は斜面長の代わりに上流域面積Accを用いた斜面長係数Lを提案した14）．

���𝐴𝐴𝑐𝑐𝑐𝑐 � 𝐷𝐷� ���� � 𝐴𝐴𝑐𝑐𝑐𝑐 ���𝐷𝐷��� �sin 𝜃𝜃 � cos 𝜃𝜃�22.13�

（2）Acc:対象グリッドの上流域に位置するグリッド群の面積[m2]，D:標高データの単位グリッド長さ[m]，θ:対象グリッドの勾配，m:斜面長のべき指数である．また，22.13[m]は傾斜枠試験の基準斜面長である．
斜面長のべき指数mは雨滴の衝撃による侵食に対するリル侵食の比βによる算出法がFoster etal.（1977）によって提案されている15）．

��𝛽𝛽
1 � 𝛽𝛽（3）βは対象グリッドの勾配による算出が McCool et al.（1989）によって提案されている16）．

���sin 𝜃𝜃 /0.0896�
3.0�sin 𝜃𝜃��.� � 0.56（4）勾配係数SはMcCool et al.（1987）が山地の急勾配斜面に対応した勾配係数Sを提案している17）．

� � 10.8 sin � � 0.03 �� � 9%�

� � 16.8 sin � � 0.50 �� � 9%�

― 190 ―

（5）（6）

θ< 9%ならば（5）式，θ> 9%ならば（6）式を適用する．
作物係数Cは植生が一切ない裸地面の土砂侵食量と植生で覆われた斜面の土砂侵食量の比である（以降，被覆係数と記述する）．農地では主に作付けされている農作物の種類毎に，調査で得られた経験的な被覆係数が与えられる．本研究では実地調査によって得られた森林被覆と被覆係数の対応関係を用いる．国内の屋外における半年以上かつ植生と裸地面との詳細な比較が行われている被覆係数Cの調査がKitahara et al.（2000）によってまとめられている18）．Kitahara et al.（2000）は天然林，人工林にかかわらず被覆係数は0.005から0.02の範囲の値をとり，林床がない場合を除き，0.01が適切であると指摘している．裸地の被覆係数は植生が存在しないため，1.0としている．
降雨係数は雨滴による表面土層の土砂の剥離に伴う侵食を代表する指標である．一般的に降雨強度と雨滴の粒径は正比例し，雨滴の粒径とその終端速度も正比例する．したがって降雨強度は雨滴の終端速度を表す変数となることから，Brown and Foster（1987）は次式を提案した19）．

𝑒𝑒� � 0.29�1 � 0.72𝑒𝑒𝑒𝑒𝑒𝑒��0.05𝑖𝑖� ��（7）

er:雨滴の運動エネルギー量[MJ・ha−1・mm−1]，im:1降雨イベントの平均降雨強度[mm/h]である．これ
らから各降雨イベントの総降雨侵食指数EI30及びその合算である年間降雨係数Ryearが算出される．

� � �� 𝑣𝑣�

（8）𝑅𝑅���� � � �� 𝐼𝐼����

（9）

�

���

k, j:各年に起こった降雨イベントの数，vr:1降雨イベントの総降雨量[mm]，I30:各降雨イベントの最大
30分間降雨強度[mm/h]，EI30:各降雨イベントの総降雨侵食指標[MJ･ha−1･mm−1]である．対象となる土
砂侵食性降雨イベントは総降雨量12.7 mm以上か最大15分間降雨量6.35 mm以上の降雨イベントと定義されている．土砂侵食性降雨イベントとみなされない降雨イベントを計算の対象外とすることで計算コストは大幅に低減される．降雨イベントは無降雨継続時間が6時間で区分される．また，降雨係数では降雪による土砂侵食は考慮されていない．降雪による土砂侵食は着雪よりも融雪の影響が大きいことが考えられるため，考慮する場合，USLE以外の融雪-土砂侵食モデルを適用しなければならない．
土壌係数Kは土壌の由来物質，粒径，透水性，有機物の含有量に基づいたモノグラフ法によって算出される．谷山（2003）は日本の農地の代表的な土質に対応した土壌係数をまとめている20）．
保全係数Pは農地においては畝や耕起などの地表面の形状効果を示す係数である．山地森林においては治山砂防工事による山腹工や法面保護がこれにあたる．厚真川流域において，実際には砂防工事が入っているが本研究では地震による影響を調べるためにこれらを考慮せず，保全係数Pは1.0とする．

2. 4. 2. 解析に用いたデータ解析に用いたデータは厚真川流域の地震前後のDigital Elevation Model（DEM）と斜面崩壊および堆積
分布図，土壌分布図，植生分布図，降雨量データである．厚真川流域のDEMは土砂崩壊の地形的特徴の把握，USLEの地形係数の算出のために用いた．Benavidez（2018）はDEMの解像度が粗いと，USLEの土砂侵食率は過小評価される傾向にあると指摘している21）．したがって，地震前の標高データは国土地理院が提供している，厚真川流域において使用可能な最小解像度である10 m解像度のDEMを用いた．
また，地震後の標高データは北海道開発局から提供された2018年9月の0.5 m解像度のLaser Profiler（LP）
データを10 mに内挿して用いた．
土砂崩壊，堆積の判別のために国土地理院が提供している83 m解像度の斜面崩壊・堆積分布図を用い― 191 ―

た．さらにUSLEの降雨係数，土壌係数，作物係数の算出のために気象庁が提供しているAutomatedMeteorological Data Acquisition System（AMeDAS）10分間降雨量データ，農林水産省が提供している20
万分の1日本土壌図，環境省が提供している2.5万分の1現存植生図を用いた．
また，解析結果が妥当であるか検証するために地球観測光学衛星Sentinel-2の可視光域バンド画像（撮影日:2018年10月20日）を合成することで厚真川流域のカラー合成画像（コンポジット画像）を作成した．

2. 5. 土砂崩壊の特徴の解析結果
以降の解析は崩壊斜面が集中している厚真川上中流域を対象領域として行う．土砂崩壊は主に解析対象流域の中流域から下流域の山地で発生しており，崩壊斜面，堆積面積は解析対象流域の13%であった．
厚真川流域の各支川の崩壊斜面率を比較するために厚真川流域を平均3.2 km2，69個のサブ流域に分割
し，サブ流域ごとの崩壊斜面率を算出した．なお，崩壊面積率は土砂崩壊・堆積部の面積を全面積で除した値とした．これ以降，土砂崩壊面と堆積部をまとめて裸地と記述する．算出結果を図-8（a）に示す．
解析対象流域において崩壊面積率が20%を超えるサブ流域は69流域中16流域存在し，広範囲で土砂崩壊が生じたことがわかる．崩壊斜面の面積率が最も大きいサブ流域は下流に位置する東和川流域で41%である．
土砂崩壊の地形的な要因を把握するために地震前の標高，勾配ごとに裸地の頻度を調べた．勾配は標
高データからある斜面の最急勾配値を代表値として算出した．結果を図-8に示す．図-8（b）,（c）は解析
対象流域の地震前の標高と勾配分布である．分布図の中央の白抜き部分のうち，左下が厚幌ダム，右上が厚真ダムである．図-8（d）,（e）はある範囲ごとの解析対象流域の裸地と非裸地の頻度分布である．図

図-8 （a）解析対象流域におけるサブ流域ごとの崩壊斜面，

堆積部の面積率.灰色の線で囲われているサブ流域は東和川流域と日高幌内川流域，（b）地震前の標高分布図，（c）地震前の勾配分布図，（d）標高の範囲ごとの裸地面積の相対頻度分布グラフ，（e）勾配の範囲ごとの裸地面積の相対頻度グラフ― 192 ―

-8（d）より，標高50-100 mの流域面積は26.5 km2であり，解析対象流域の11.7%と比較的小さい一方で裸地面積率は25%と最も高い．下流域である標高0–200 mの範囲には裸地の総面積の84%が集中している．つまり厚真川流域では下流域において土砂崩壊，堆積リスクが高いことが示唆される．図-8（e）より，各勾配範囲の間に裸地面積率の差はほとんど見られず，20-30°の範囲に裸地の総面積の39%が分布している．比較的勾配が小さい0-10°の範囲は主に堆積が発生していると考えられる．したがって，
急勾配地形では土砂崩壊率は高いが，急勾配地形の割合は小さいため全体としては斜面の裸地化による土砂侵食リスクは勾配の大小に関係せず存在することが示唆される．

2. 6. 厚真川流域を対象としたUSLEの適用結果
2. 6. 1. USLEの係数の算出結果
地形係数LSの算出結果を図-9に示す．Renard et al.（1997）はUSLEの適用が不適当である地形の条件を2つ設け，条件に当てはまるグリッドの地形を計算の対象外とした22）．ひとつは河道である．USLEは
河道の土砂輸送を表現することができないので上流域面積Accが5,000 m2以上のグリッドを河道とみな

図-9

 形係数LSおよび地震前と比べた地震後の地形係数の倍率の分布図
地
LSの色は（a）の解析対象流域内の平均値から1/2標準偏差値ごとに分類されている．
（a）地震前の地形係数（b）地震後の地形係数（c）地震前と比べた地震後の倍率（d）（a）の黒枠の拡大図（e）（b）の黒枠の拡大図（f）（c）の黒枠の拡大図
（g）
（a）
 （b）
（c）の黒枠の部分の拡大衛星写真．衛星写真は地球観測光学衛星Sentinel-2の可視光域バンド画像を重ねて作成したカラー合成衛星画像（コンポジット画像）（撮影日:2018年10月20日）― 193 ―

し，計算の対象外とした．もうひとつは十分に勾配が緩やかな平地である．平地では侵食ではなく堆積が発生するので，斜面勾配5%未満のグリッドは平地とみなし，計算の対象外とした．加えてPanagos etal.（2015）は山地森林の勾配が50%以上の斜面に脆い土壌はほとんど存在し得ず，雨滴の衝撃による土砂侵食は発生しないため勾配係数の算出で無視されるべきであると指摘している23）．しかし，解析対象
流域の急勾配斜面の大半を占める崩壊斜面の表層土質は比較的脆いTa-dであるため,本研究では勾配の上限は設けない．図-9（a）,（b）はそれぞれ地震前と地震後の地形係数の分布図である．図-9（c）は地震前に対する地震後の地形係数の倍率である．図-9（d）～（g）は上部の枠で囲われた部分の拡大図である．図-9（f）から峰の部分では地形係数が大きく増加していることがわかる．これは地震後，崩壊斜面の上方境界の地形が急勾配に変化し，勾配係数が増加したことに起因する．反対に崩壊斜面の下部，および支川周辺は地形係数が減少している．減少の要因は勾配と斜面長の2つの要因に分けられる．前者
は土砂の堆積による緩勾配斜面の形成，後者は上流域面積の減少である．上流域面積の減少とは斜面における水流が集中しうる谷地形が土砂堆積によって埋められて平滑化することを示す．さらに流域の全グリッドの地形係数を積算した値を地震前と地震後で比較すると地震後は0.9倍と変化はわずかである．
後述する2.6.2.節の（10）式をLとSについて適用したところ減少に支配的な係数は斜面長係数であることが示された．また，流域全体では勾配係数は増加傾向，斜面長係数は減少傾向にあることが示された．
図-10は地震前と地震後の地形係数の相対頻度グラフである．図-10（b）はLSafの下位99%区間における相対頻度グラフである．地形係数の相対頻度は10未満の範囲では増加傾向にある一方で，12以上では減少傾向にある．これは地震後の地形は一様に土砂侵食が激しくなる地形に変化したのではなく，土砂侵食が発生しにくい地形に変化した部分もあることを示している．

図-10（a）地形係数LSの範囲ごとの面積の相対頻度グラフ（b）（a）の地震後の地形係数の下位99%区間における相対頻度グラフ
被覆係数Cの分布図を図-11に示す．土地利用の判別には2.5万分の1現存植生図を適用し，崩壊斜面の判
別には斜面崩壊・堆積分布図を適用した．解析対象流域の主な土地利用は森林や裸地，市街地，伐採跡地である．森林の被覆係数はKitahara et al.（2000）が示した0.01，裸地には1.0を適用した．市街地では土砂侵食が発生しないため，解析の対象外とした．同様に伐採跡地は面積率が解析対象流域の2.9%と小さく，流出土砂量への寄与度が小さいため解析の対象外とした．図-11（c）,（e）を比べると航空写真で判別できる裸地面と裸地面の被覆係数は一致している．図-11（b）より，厚真ダム上流部では崩壊斜面がほとんど分布しない一方で，厚幌ダムの上流域では斜面崩壊が発生しており，土砂流入のリスクは比較的大きい可能性がある．流域の全グリッドの被覆係数の平均値については地震前と比べて地震後は12倍増加している．

― 194 ―

図-11

被覆係数Cの分布図（a）地震前（b）地震後（c）（a）の黒枠部の拡大図（d）（b）の黒枠部の拡大図（e）（a）
 （b）の黒枠の部分の拡大衛星写真．衛星写真は地球観測光学衛星Sentinel-2の可視光域バンド
画像を重ねて作成したカラー合成衛星画像（コンポジット画像）（撮影日:2018年10月20日）
胆振地方東部の降雨係数Rはアメダスの厚真観測所の1994年から2020年までの4月から11月の10分間
降雨量から算出された．一般的に降雨係数は対象地域の全ての過去降雨イベントの降雨侵食エネルギーの平年値を用いるが，本研究では降雨の年変動に伴う土砂侵食率の年変動を評価するために各年の年間の降雨係数Ryearを算出した．算出結果を図-12に示す．平年値は1,055.2 MJ・ha−1・mm−1であり，地震後の2019年，2020年の降雨係数はそれぞれ945.2 MJ・ha−1・mm−1，415.3 MJ・ha−1・mm−1でありいずれも平
年値よりも小さかった．そのため将来，平年値を上回る降雨が発生した場合，現在よりも激しい土砂侵
食が起こる可能性がある．図-12より最小，最大年間侵食性降雨量を記録した年は1997年，2008年である一方，最小，最大降雨係数を記録した年は2004年，2000年である．このように降雨係数はイベントベースで算出されるため年侵食降雨量との相関はあるが，最大，最小値を記録した年は一致しない．
土壌係数は谷山（2003）によってまとめられた土壌係数を適用した．本流域の主な表層土壌は産総研の地質図よりTa-d，褐色森林土，灰色低地土であった．Ta-dは淡色黒ボク土に該当するため0.032とした．
褐色森林土と灰色低地土は粒度の情報がなかったため，土質係数はそれぞれの土質の平均値をとり0.028，0.047とした．

― 195 ―

図-12 各年の年降雨係数Ryearと年間侵食性降雨量（侵食性降雨:総降雨量12.7 mm以上か15分間降雨量6.4 mm以上の降雨イベント）

2. 6. 2. 年間土砂侵食率の算出結果
USLEによって算出された地震前後の平年土砂侵食率とその分布図を図-13（a）～（c），相対頻度分布図
を図-14に示し，表1にまとめた．地震前後の平年土砂侵食率は地震前後の地形に1994年から2020年までの年間降雨係数の平均値をそれぞれ適用して算出した（表-1，Case1）．なお単位は[t・ha-1・yr-1]から1
グリッド（100 m2）あたりの年間土砂侵食率[t・grid-1・yr-1]に換算している．
表-1Case1より土砂侵食率の流域全体の積算値は地震前が1.06×107 t/yr，地震後が16.00×107 t/yrであり，地震前に対する地震後の土砂浸食率は15倍に増加している．図-13（b）と図-1（a）を比較すると土砂侵食率は裸地面で顕著に増加していることが確認された．また，図14-（b）より地震後，地震前には存在しなかった1.5 t・grid-1・yr-1以上の土砂侵食が発生しうる急勾配，裸地斜面が形成されていることが確認された．観測対象流域における年土砂侵食率の積算値は地震前が3.02 t/yr，地震後が161.02 t/yr
であり，50倍増加している．厚真川流域全体よりも増加率が大きいことは流域全体の崩壊面積率が13%であることに対し，観測対象流域の面積率が54 %であることに起因すると考えられる．
表-1 Kitahara et al.が示した森林における被覆係数Cの値ごとに算出された年間土砂侵食率の流域積算値― 196 ―

図-13年間土砂侵食率の分布図.
色は（a）の平均値からの1/2標準偏差値ごとに分類されている（a）地震前 （b）地震後 （c）地震前と比べた地震後の倍率
サブ流域ごとの年間土砂侵食率の積算値.
色は（e）の平均値からの1/2標準偏差値で分類されている（d）地震前 （e）地震後 （g）地震前と比べた地震後の倍率

厚真川流域のサブ流域の土砂侵食リスクを比較するために図-8（a）と同じサブ流域ごとに地震前，地震後の平年土砂侵食率の積算値を算出した．結果を図-13（d）,（e）,（f）に示す．図-13（e）と図-8（a）を比較すると崩壊面積率が高いサブ流域ほど，土砂侵食率は大きい．地震後，日高幌内川では大規模な深層崩壊が発生しているため特に土砂侵食が激しくなっている．図-13（f）より土砂侵食率が20倍より大きいサブ流域は69
流域中，19流域である．図-8（a）より裸地の面積率が41%と解析対象流域内のサブ流域の中で最も大きい東
和川の増加率は49倍である．また，厚真ダム上流域において土砂崩壊がほとんど発生していないためサブ
流域では土砂侵食率はほとんど増加していない．一方で厚幌ダム上流域においては土砂侵食率が増加しているサブ流域が存在する．そのため，厚真ダムよりも厚幌ダムの土砂流入リスクが大きいことが示された．

図-14 （a）年間土砂侵食率の範囲ごとの面積の相対頻度グラフ（b）（a）のAafの99%区間における面積の相対頻度グラフ― 197 ―

地震による土砂侵食率の増加について支配的な係数の特定を行った．あるグリッドの土砂侵食率は地形係数LSと被覆係数Cの変化に伴い増加した．LS（またはC）以外の係数を固定し，地震前LS（C）と地震後LS（C）を適用して土砂浸食率を比較し，LS（C）の影響を定量化する感度分析を行うことは有用である．しかし現実にはLSとC両方が変化しているため土砂浸食率はLSとCの共変効果を受けている．
共変効果を考慮するために地形係数LSと被覆係数Cの積の地震前後の変化量を求め，それを地形変化の影響項，被覆変化の影響項，共変変化項の3つに分離して地形と被覆変化の影響を比較する．

ߜሺ‫ܵܮܥ‬ሻ� ��� 𝐿𝐿𝐿𝐿�� � ��� 𝐿𝐿𝐿𝐿��（10）� �𝐶𝐶�� � �𝐶𝐶��𝐿𝐿𝐿𝐿�� � �𝐿𝐿𝐿𝐿� � 𝐶𝐶�� 𝐿𝐿𝐿𝐿��ൌ Ɂ‫ܥܵܮ‬௕௙ ൅ ‫ܵܮ‬௕௙ Ɂ‫ ܥ‬൅ Ɂ‫ܵܮ‬Ɂ‫ܥ‬

下付き文字bf, afはそれぞれ地震前，地震後を示し，δは地震前後の変化量であることを示す．第1項は地形係数，第2項は被覆係数の変化による効果を示し，第3項は地形係数と被覆係数の両者の共変効果による影響を示す．各項の分布図を図-15に示す．各項のグリッドごとに絶対値をとったときの流域積算

図-15

式（10）から算出された各項の分布図（a） δ𝐿𝐿𝐿𝐿𝐿𝐿��

（b） 𝐿𝐿𝐿𝐿�� δ𝐶𝐶（c） δ𝐿𝐿𝐿𝐿𝐿𝐿𝐿（d）（a）の黒枠部分の拡大図（e）（b）の黒枠部分の拡大図（f）（c）の黒枠部分の拡大図
（g）
（a）
 （b）
（c）の黒枠部分の衛星写真．衛星写真は地球観測光学衛星Sentinel-2の可視光域バンド画像を重ねて作成したカラー合成衛星画像（コンポジット画像）（撮影日:2018年10月20日）― 198 ―

図-16 各
 年の降雨係数と地震前後の降雨係数以外の係数から算出された年土砂浸食率の相対頻度グラフ(年土砂浸食率は流域積算値であり，地震前の年土砂浸食率の標準偏差ごとに階級分けされている.)
値は第1項から0.16×106，3.8×106，2.4×106であったことから，主な構成成分は第2項である．したがって地震による土砂侵食率の変化は被覆係数の変化が支配的であることが示唆される．衛星写真で確認された裸地と第2項の値が大きい地点は合致している（図-15（e）,（g））ことから，土砂侵食率は裸地で増加していることが改めて確認された．また，図-15（f）より共変効果によって河道周辺で土砂侵食率は増加，減少している．共変効果による土砂侵食率の増加，減少は裸地化による被覆係数の増加と谷地形の堆積による平滑化または緩勾配化による地形係数の減少の差で決定されることが示唆される．
被覆係数が支配的であることから，土砂侵食リスクを評価する上で被覆係数の値に幅をもたせた土砂
侵食率の算出が重要と考えられる．森林の被覆係数の値にKitahara et al.（2000）が示した最小，最大値である0.005，0.02を適用して地震前後の土砂侵食率を算出した結果が表-1のCase2，Case3である．
Case2.Case3より地震前からの土砂侵食率の増加率は8倍〜29倍の範囲となった．
厚真川流域における土砂侵食リスクの評価は降雨による年変動を考慮することも重要である．そのため，1994年から2020年までの年降雨係数Ryearを用いて地震前後の年土砂侵食率をそれぞれ26通り，計52通り算出した．図-16はその算出結果である．年土砂侵食率は流域積算値であり，地震前の年土砂侵食率の標準偏差0.05×108 t/yrごとに階級分けされている．相対頻度は地震前と地震後それぞれにおいて，母数を26通りとして算出した．地震後の最小年土砂侵食率は地震前の最大年土砂侵食率を超過している．これは年降雨係数の倍率が平均値に対して0.3倍から1.9倍である（図-12）ことに比べ，地震前後の
年土砂侵食率の増加率が最小であっても8倍（表-1Case 3）であり，降雨変動よりも地震の方が年土砂侵食率への影響がはるかに大きいことに起因する．したがって地震後の厚真川流域では地震前の年最大
土砂侵食率を上回る土砂侵食が毎年発生することが示唆される．また，年土砂侵食率の標準偏差は地震
前が0.05×108 t/yrであることに対し，地震後が0.77×108 t/yrであった．よって，厚真川流域における最小から最大の土砂侵食リスクの幅は大きく広がったことが示された．この結果より，厚真川流域における土砂侵食リスクについての議論は平均値に基づくのではなく，極値に基づいてなされる方が適当であることも示唆される．このように現在の厚真川流域の土砂侵食ポテンシャルは地震前に比べて非常に高い状態にあり，さらに極端気象現象が増加すると予想されることを鑑みると，厚真川流域の土砂侵食リ
スクも今後増加し続けていくと考えられる．

― 199 ―

2. 7. まとめ
本研究では厚真川の一支川流域からの流出土砂量の観測と厚真川流域全体のUSLEによる土砂侵食リ
スク増加の定量化を行った．観測期間中の本流域からのイベントごとの総土砂流出量は0.2 tから26.3 tであった．本観測地点の2019年と2020年の日最大降雨量は63.8 mmである一方，アメダス厚真観測点における既往最大の日降雨量は187 mmである．このことから観測期間中に発生した以上の土砂流出が発生する可能性が示唆される．また，本観測では濁度計が水面に浸かっている間の浮遊している土砂のみ
観測しているため，低水時の流出土砂および掃流砂を含めると実際の総流出土砂量は算出値以上であると考えられる．また，本流域の土砂流出の関係性はイベントごとで大きく異なる．総流出土砂量は総降
雨量とはある程度の相関がみられるが，平均降雨強度とは線形関係ではないことが示唆された．加えて，流出土砂濃度の時間変化が観測期間中最大の降雨イベントのみ異なったことは降雨と土砂流出の関係を把握する上で重要だと考えられる．
将来，総流出土砂量および降雨と土砂流出の関係は降雨の特徴と流域の特徴によって変動すると考えられる．前者は総降雨量や平均降雨強度，降雨の時空間分布などが，後者は崩壊斜面の植生の回復，土層の水平，鉛直分布などがあげられる．今後，土砂流出のモニタリングを続け，流出土砂量の推移と降雨，流域特性の関係から土砂流出機構および土砂流出リスクの将来変化の推定を試みる．
厚真川流域全体の土砂侵食率の増加の定量的な評価およびその支配要因の特定を行った．その結果，厚真川流域全体で土砂侵食率は8倍から29倍増加しており，サブ流域では最大49倍増加していることが示唆された．また，地震による土砂侵食リスクの増加は降雨の年変動に伴う土砂侵食率の変動よりも大きく，増加の支配要因は植生の裸地化であることが示された．
厚真川下流域は上流域に比べて土砂侵食率の増加率が大きく，特に東和川流域の土砂侵食率の増加率
は厚真川流域のサブ流域の中で最も大きい．河道への土砂堆積を考える上でサブ流域における侵食土砂の下流（市街地）への輸送過程を考慮する必要があるが，USLEでは水流による土砂輸送を表現できない．
したがって今後は現地観測データと力学モデルの適用により，厚真川流域の降雨-土砂流出モデルを構築する必要がある．
地震による植生，地形の変化の内，土砂侵食率の変化に支配的な要因は植生の変化である．そのため，今後の土砂侵食率の変動を把握するために長期にわたる植生のモニタリングがなされるべきである．
Lin et al（2006）は地震後の増加した土砂侵食率は6年後には植生の変化に伴い地震前とほとんど同じ量にまで減少していると報告している24）．しかし，植生遷移は10から102年スケールで遷移し，また植生遷
移と地形の侵食は相互に干渉し合う．よって，土砂侵食率は今後，非線形に回復する可能性がある．し
たがって，今後厚真川流域において土砂侵食率を推定する場合，植生のモニタリングと現地観測に基づ
く土砂侵食と植生繁茂の相互干渉を考慮したモデルを構築する必要がある．
地震後の年最小土砂侵食率は地震前の年最大土砂侵食率を上回ると試算された．したがって，地震後の厚真川流域では地震前の年最大土砂侵食率以上の土砂侵食が毎年発生することが示唆された．本研究では過去の降雨変動のみを対象としたが，将来極端降雨が増加すると予想されていることを踏まえると，
侵食性降雨イベントの増加が示唆される．したがって，将来気候モデルを用いた将来極端降雨により発生しうる最大の土砂侵食率の予測を含めて降雨と地震の要因比較を行うべきである．
謝辞
本研究に用いられたデータは北海道開発局から提供を受けた．現地への立ち入り及び水文観測施設の設置については北海道および厚真町に多大な協力をいただいた．ここに謝意を表する．

― 200 ―

引用文献
1）Kasai, Mio, and Takashi Yamada.“Topographic effects on frequency-size distribution of landslides triggeredby the Hokkaido Eastern Iburi Earthquake in 2018.,”Earth, Planets and Space 71. 1（2019）: 1-12.
2）関口 辰夫,佐藤 浩,新潟県中越地震における斜面崩壊の特徴と分布,日本地すべり学会誌,2006,43 巻,3号,p. 142-154
3）小山内信智ら,“平成 30 年北海道胆振東部地震による土砂災害.”砂防学会誌 71.5（2019）: 54-65
4）藤田一郎:河川表面流速の画像計測ソフトウェア KU-STIV 開発の背景とその応用,建設工学研究所
論文報告集第 57 号[報告],2015 年 11 月
5）国土交通省:河川砂防技術基準調査編,第2章第 4 節
6）Hafzullah Aksoy, M. Levent Kavvas, A review of hillslope and watershed scale erosion and sedimenttransport models,CATENA, Volume 64, Issues 2–3,（2005）:247-27.
7）R.A. Young, C.A. Onstad, D.D. Bosch, W.P.Anderson, AGNPS: a nonpoint-source pollution model forevaluating agricultural watersheds, Journal of Soil and Water Conservation（1989）, pp. 168-173,（March–
April）Viney, Neil R., and Murugesu Sivapalan.“A conceptual model of sediment transport: application to
the Avon River Basin in Western Australia.”Hydrological Processes 13.5（1999）: 727-743.
8）Beasley, D. B., L. F. Huggins, and ampEJ Monke.“ANSWERS: A model for watershed planning.”
Transactions of the ASAE 23.4（1980）: 938-0944.
9）Nearing, Mark A., et al.“A process-based soil erosion model for USDA-Water Erosion Prediction Projecttechnology.”Transactions of the ASAE 32.5（1989）: 1587-1593.
10）Morgan, R. P. C., et al.“The European Soil Erosion Model（EUROSEM）: a dynamic approach forpredicting sediment transport from fields and small catchments.”Earth Surface Processes and Landforms:
The Journal of the British Geomorphological Group 23.6（1998）: 527-544.
11） Wischmeier, Walter H., and Dwight David Smith. Predicting rainfall erosion losses: a guide to conservation
planning. No. 537. Department of Agriculture, Science and Education Administration, 1978.
12）Williams, J.R., 1975. Sediment-yield prediction with Universal Equation using runoff energy factor, present
and prospective technology for predicting sediment yield and sources. ARS-S-40. Brooksville, FL: US
Department of Agriculture, Agricultural Research Service, 244–252.
13）Renard, Kenneth G. Predicting soil erosion by water: a guide to conservation planning with the Revised
Universal Soil Loss Equation（RUSLE）. United States Government Printing, 1997
14）Desmet, P. J. J., and Gerard Govers.“A GIS procedure for automatically calculating the USLE LS factor on
topographically complex landscape units.”Journal of soil and water conservation 51.5（1996）: 427-433.
15）Foster, G. R., L. D. Meyer, and C. A. Onstad.“A runoff erosivity factor and variable slope length exponents
for soil loss estimates.”Transactions of the ASAE 20.4（1977）: 683-0687.
16）McCool, Don K., et al.“Revised slope length factor for the Universal Soil Loss Equation.” Transactions of
the ASAE 32.5（1989）: 1571-1576.
17）McCool, D. K., et al.“Revised slope steepness factor for the Universal Soil Loss Equation.”Transactions of
the ASAE 30.5（1987）: 1387-1396.
18）Hikaru Kitahara, Yoichi Okura, Toshiaki Sammori & Akiko Kawanami（2000）Application of Universal Soil
Loss Equation（USLE）to Mountainous Forests in Japan, Journal of Forest Research, 5:4, 231-236, DOI:
10.1007/BF02767115
19）Brown, L. C., and G. R. Foster.“Storm erosivity using idealized intensity distributions.” Transactions of the
ASAE 30.2（1987）: 379-0386.

― 201 ―

20）谷山一郎.“農耕地からの表面流去水の発生に係わる土壌要因の解明とMIの作成.” 農林水産業および農林水産物貿易と資源 環境に関する総合研究（2003）: 149-152.
21）Benavidez, R. A.: Understanding the effect of changing land use on floods and soil erosion in the Cagayan deOro catchment,（Un- published doctoral dissertation）, Victoria University of Welling- ton, New Zealand,
2018.
22）Renard, Kenneth G. Predicting soil erosion by water: a guide to conservation planning with the Revised
Universal Soil Loss Equation（RUSLE）. United States Government Printing, 1997.
23）Panagos, Panos & Borrelli, Pasquale & Meusburger, K..（2015）. A New European Slope Length and Steepness
Factor（LS-Factor）for Modeling Soil Erosion by Water. Geosciences. 2015. 117-126.
24）Lin, Wen-Tzu, Chao-Yuan Lin, and Wen-Chieh Chou.“Assessment of vegetation recovery and soil erosion
at landslides caused by a catastrophic earthquake: a case study in Central Taiwan.” EcologicalEngineering 28.1（2006）: 79-89.

― 202 ―