流域產輸沙模型中考慮泥沙來源的動態粒徑計算方法

徐東坡，周祖昊，蔡靜雅，劉佳嘉，賈仰文，王 浩

（1.中國水利水電科學研究院 流域水循環模擬與調控國家重點實驗室，北京 100038；2.黃河勘測規劃設計研究院有限公司，河南 鄭州 450003；3.水利部黃河流域水治理與水安全重點實驗室（籌），河南 鄭州 450003）

目前已有大量學者開展了關于土壤侵蝕模型的研究，現有土壤侵蝕模型基本分為經驗模型、概念模型和物理模型［1］。USLE、RUSLE等經驗模型主要通過統計方法分析侵蝕量、輸沙率等與降雨量、徑流量的經驗關系，實現土壤侵蝕過程的模擬與預測，計算過程較簡單且對數據要求較低，但存在缺乏機理性、模擬精度低、推廣性有限等缺點［2］。LASCAM等概念模型對土壤侵蝕過程進行一定概化，包含一定的物理機制，能在一定程度上反映下墊面變化對土壤侵蝕過程的影響［3－4］，但模擬部分侵蝕過程時仍停留在經驗關系上，物理機制不完備。WEPP［5］、DYRIM［6］和WEP－SED［7］等物理模型是基于質量、動量和能量三大守恒原理建立土壤顆粒運動的物理方程，具有較完備的物理機制［8］。相較于其他土壤侵蝕模型，物理模型能充分反映降雨和下墊面等因素的時空異質性，因此該類模型被較多學者應用。采用物理模型計算泥沙沉速時往往基于Stokes公式［9］等，計算水流挾沙力時往往基于Govers公式［10］、Abrahams公式［11］、沙玉清公式［12］等，此類公式中泥沙粒徑是一個關鍵參數，其數值大小與準確性可以直接或間接地影響模型計算結果。許多學者采用泥沙中值粒徑這一特征值來表征泥沙粒徑分布特征［13－14］，通常采用固定的泥沙中值粒徑值在長時間序列上進行模型模擬計算，但根據大量研究結果可知，泥沙中值粒徑的年內分布特征與年際分布特征存在差異，其年內分布特征表現為汛期數值比非汛期數值小［15－17］，難以準確反映泥沙粒徑分布特征在不同時間尺度上的變化過程。

鑒于泥沙粒徑分布特征對模型模擬計算的重要性，本文以瑪曲水文站以上區域為研究對象，根據不同季節的泥沙來源不同這一特點，選用具有“坡面—溝壑—河道”三級匯流結構的WEP－SED分布式流域產輸沙模型，根據不同季節的坡面、谷溝、河道的水沙動力學條件計算泥沙中值粒徑，分析泥沙中值粒徑在時間尺度上的動態變化過程，以對流域產輸沙過程進行更詳盡的刻畫，提高產輸沙模型的模擬精度。

1 研究區概況與數據來源

1.1 研究區概況

黃河發源于青藏高原巴顏喀拉山脈北麓的約古宗列盆地，流經青海、四川、甘肅、寧夏、內蒙古、山西、陜西、河南、山東九省（區），于山東省墾利、利津兩縣之間流入渤海。黃河干流在瑪曲河段繞阿尼瑪卿山回轉180°，瑪曲河段有“天下黃河第一彎”之稱。選擇黃河干流瑪曲水文站以上區域為研究區域（見圖1），瑪曲水文站控制面積約8.6萬km2、河長約910 km，其中：黃河河源至達日段流經較平整的高原草甸區，該區域支流較多，擁有眾多的湖泊沼澤；達日至久治段河道兩岸多為較高且平緩的山坡，河谷逐漸縮窄；久治至瑪曲段流經眾多沼澤區域，該區域是河源區的第一大暴雨區，同時有支流黑河和白河匯入使區域整體水量大增，成為唐乃亥水文站以上區域的主要產洪區［18］。

圖1 瑪曲水文站控制區域

1.2 數據來源

流域產輸沙模型所需數據主要包括瑪曲水文站以上區域的DEM數據、水系拓撲信息、水文氣象數據和土地利用類型數據等。DEM數據為源自地理空間數據云網站的30 m精度的ASTER GDEM數據。水系拓撲信息為源自國家地理信息數據庫的1∶250 000空間分辨率的流域水系圖。日降雨量、日平均氣溫、日照時數、相對濕度以及風速等水文氣象數據源自中國氣象局和黃河水利委員會水文局。用于對模型模擬結果進行率定和驗證的實測數據源自《中華人民共和國水文年鑒》，主要包括1956—2018年瑪曲水文站實測的月平均流量和月平均輸沙率。土地利用類型數據源自中國多時期土地利用土地覆被遙感監測數據集（CNLUCC）［7］。河道初始泥沙中值粒徑由模型的輸入文件設定，數據源自《黃河泥沙公報》中多年平均泥沙中值粒徑，其在模型中作為可調參數進行模型率定。

2 研究方法

2.1 WEP－SED模型原理

為詳盡刻畫黃河流域產輸沙過程，WEP－SED模型將WEP－L模型中“坡面—河道”匯流結構改進為“坡面—溝壑—河道”三級匯流結構［19］，并分別進行各環節的產輸沙量計算（見圖2）。其中：坡面侵蝕分為雨滴濺蝕、薄層水流侵蝕以及坡面細溝侵蝕，谷溝侵蝕主要有重力侵蝕和沖刷侵蝕等，河道侵蝕主要是沖刷侵蝕。谷溝、河道產輸沙模擬主要基于不平衡輸沙理論，值得注意的是，模型使用不平衡輸沙理論時谷溝和河道的側向來沙存在一定差異，谷溝的側向來沙包含當前子流域部分坡面來沙，河道的側向來沙除有部分坡面來沙外，還有谷溝匯入河道的沙量。

圖2 模型模擬產輸沙過程示意

在模型模擬計算谷溝環節和河道環節的產輸沙量時，水流挾沙力和泥沙沉速等參數的計算極其重要，并且這2個參數與泥沙粒徑大小具有一定關系。水流挾沙力公式如下［20］：

式中：T為水流挾沙力，kg／m3；C為水流含沙量，kg／m3；γs、γm分別為泥沙、渾水的容重，N／m3；K、K0分別為渾水、清水的卡門常數；g為重力加速度，m／s2；R為水力半徑，m；ω為泥沙沉速，m／s；h為水深，m；v為流速，m／s；d50為泥沙中值粒徑，mm，即粒徑小于d50的泥沙質量占全部泥沙質量的50%。

采用Stokes公式［9］計算泥沙沉速，公式為

式中：D為泥沙沉降粒徑，mm，在模型計算時等同于泥沙中值粒徑d50；ρs為泥沙密度，g／cm3；ρw為清水密度，g／cm3；ν為水的運動黏滯系數，cm2／s。

2.2 泥沙動態粒徑計算方法

基于WEP－SED模型的基本原理，針對模型中谷溝環節和河道環節的產輸沙量計算，把不同泥沙來源的沙量作為權重對泥沙中值粒徑進行迭代計算，從而模擬泥沙中值粒徑隨沖淤變化的動態變化過程，計算步驟如下。

步驟一：基于WEP－SED模型采用初始泥沙中值粒徑分別進行坡面、谷溝、河道環節的初次產輸沙量計算。

步驟二：根據步驟一計算的谷溝環節的產輸沙量判斷其是否發生沖刷，若發生沖刷，則以進入谷溝的坡面產輸沙量和谷溝環節自身的重力侵蝕量和沖刷量作為權重計算谷溝環節的綜合泥沙中值粒徑；若不發生沖刷，則選取的權重不包括谷溝環節的沖刷量。谷溝環節的綜合泥沙中值粒徑計算公式如下：

式中：D谷溝綜合m為谷溝環節的綜合泥沙中值粒徑，第一次計算時其為谷溝環節的初始泥沙中值粒徑，迭代計算時其為上次計算得到的谷溝綜合泥沙中值粒徑（D谷溝綜合m－1）；Q坡面－谷溝為進入谷溝環節的坡面產輸沙量；D坡面為坡面環節的泥沙中值粒徑；Q谷溝為谷溝環節的重力侵蝕量與沖刷量之和，若谷溝環節不發生沖刷，則Q谷溝只包含谷溝環節的重力侵蝕量。

步驟三：根據步驟二計算的谷溝環節的綜合泥沙中值粒徑，第二次進行谷溝環節的產輸沙量計算；再重復步驟二第二次進行谷溝環節的綜合泥沙中值粒徑計算，然后以該結果第三次進行谷溝環節的產輸沙量計算；最終得到谷溝環節的產輸沙量和綜合泥沙中值粒徑。

步驟四：根據步驟一計算的河道環節的產輸沙量判斷其是否發生沖刷，若發生沖刷，則以進入河道的坡面產輸沙量、谷溝環節最終的產輸沙量、上游河道的產輸沙量以及本級河道的沖刷量作為權重計算河道環節的綜合泥沙中值粒徑；若不發生沖刷，則選取的權重不包括本級河道的沖刷量。

本級河道上游來沙的綜合泥沙中值粒徑計算公式如下：

式中：Q上游i為當前子流域第i個上游河道的來沙量；D上游綜合i為當前子流域第i個上游河道的綜合泥沙中值粒徑，若無上游河道，則D上游綜合＝0。

河道環節的綜合泥沙中值粒徑計算公式如下：

式中：D河道綜合m為本級河道的綜合泥沙中值粒徑，第一次計算時其為河道的初始泥沙中值粒徑，迭代計算時其為上次計算得到的河道綜合泥沙中值粒徑（D河道綜合m－1）；Q坡面－河道為本級河道的坡面來沙量；Q谷溝－河道為本級河道的谷溝來沙量；Q上游為上游河道的來沙量；D上游綜合為上游河道的綜合泥沙中值粒徑；Q河道為本級河道的沖刷量，若河道不發生沖刷，則Q河道＝0。

步驟五：根據步驟四計算的河道環節的綜合泥沙中值粒徑第二次計算河道環節的產輸沙量；再重復步驟四第二次計算河道環節的綜合泥沙中值粒徑，以該結果第三次計算河道環節的產輸沙量；最終得到河道環節的產輸沙量和綜合泥沙中值粒徑。

3 結果與分析

利用WEP－SED模型對黃河河源至瑪曲水文站的水沙運移過程進行模擬，通過計算不同來源的泥沙中值粒徑分析其對流域產輸沙過程模擬結果的影響，分別采用納什系數（NSE）、相對誤差（RE）對模擬結果進行評價。

3.1 動態泥沙中值粒徑對模型模擬輸沙率的影響

基于WEP－SED模型分析分別采用泥沙固定粒徑計算方法和不同來源的泥沙動態粒徑計算方法對瑪曲水文站逐月輸沙率的影響，模擬結果分別見圖3和圖4。與采用泥沙固定粒徑計算方法的模擬結果相比，采用泥沙動態粒徑計算方法的逐月輸沙率模擬值尤其汛期的輸沙率模擬值明顯增大。為量化模型采用泥沙動態粒徑計算方法對逐月輸沙率模擬結果的改進程度，分別計算2種模擬結果的NSE值和RE值。采用泥沙動態粒徑計算方法時模型模擬結果的NSE值和RE值分別為0.706和0.04%，采用泥沙固定粒徑計算方法時模型模擬結果的NSE值和RE值分別為0.601和－1.09%，表明采用動態泥沙中值粒徑可提高WEPSED模型的模擬精度。

圖3 采用泥沙動態粒徑時瑪曲水文站1956—2018年逐月輸沙率模擬結果

圖4 采用泥沙固定粒徑時瑪曲水文站1956—2018年逐月輸沙率模擬結果

3.2 泥沙中值粒徑的動態變化過程

為明確WEP－SED模型模擬汛期輸沙率提升的根本原因，分析采用泥沙動態粒徑計算方法時模型模擬的泥沙中值粒徑逐日變化過程，并對采用泥沙動態粒徑計算方法的泥沙中值粒徑逐月變化過程進行統計，其結果見圖5。可以看出，采用泥沙動態粒徑計算方法實現了對泥沙中值粒徑動態變化過程的模擬。由圖5（b）可知，采用泥沙動態粒徑計算方法時模型模擬的汛期泥沙中值粒徑小于非汛期泥沙中值粒徑，因此模型模擬結果體現了汛期更易發生沖刷過程，實現了汛期輸沙率模擬精度的提高，該結論與馬勇等［15］、全棟［16］、孫維婷等［17］、薛博文等［21］、那巍等［22］研究黃河流域實測輸沙率的年內變化規律呈現出一致性。

圖5 河道泥沙中值粒徑年內變化情況

4 結 論

本文以黃河干流瑪曲水文站以上區域為研究對象，基于WEP－SED模型，考慮不同季節泥沙來源不同的特點進行泥沙中值粒徑的計算，實現了對泥沙中值粒徑動態變化過程的模擬。與采用泥沙固定粒徑計算方法相比，采用泥沙動態粒徑計算方法使河源至瑪曲段1956—2018年逐月輸沙率的模擬精度有所提升，逐月輸沙率模擬結果的NSE值從0.601提升至0.706，并且使得模型模擬的泥沙中值粒徑呈現汛期小、非汛期大的動態變化規律，模型模擬的泥沙中值粒徑逐月變化規律與相關研究結果及實測資料一致。