流域產(chǎn)輸沙模型中考慮泥沙來源的動(dòng)態(tài)粒徑計(jì)算方法

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

（1.中國水利水電科學(xué)研究院 流域水循環(huán)模擬與調(diào)控國家重點(diǎn)實(shí)驗(yàn)室，北京 100038；2.黃河勘測(cè)規(guī)劃設(shè)計(jì)研究院有限公司，河南 鄭州 450003；3.水利部黃河流域水治理與水安全重點(diǎn)實(shí)驗(yàn)室（籌），河南 鄭州 450003）

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

鑒于泥沙粒徑分布特征對(duì)模型模擬計(jì)算的重要性，本文以瑪曲水文站以上區(qū)域?yàn)檠芯繉?duì)象，根據(jù)不同季節(jié)的泥沙來源不同這一特點(diǎn)，選用具有“坡面—溝壑—河道”三級(jí)匯流結(jié)構(gòu)的WEP－SED分布式流域產(chǎn)輸沙模型，根據(jù)不同季節(jié)的坡面、谷溝、河道的水沙動(dòng)力學(xué)條件計(jì)算泥沙中值粒徑，分析泥沙中值粒徑在時(shí)間尺度上的動(dòng)態(tài)變化過程，以對(duì)流域產(chǎn)輸沙過程進(jìn)行更詳盡的刻畫，提高產(chǎn)輸沙模型的模擬精度。

1 研究區(qū)概況與數(shù)據(jù)來源

1.1 研究區(qū)概況

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

圖1 瑪曲水文站控制區(qū)域

1.2 數(shù)據(jù)來源

流域產(chǎn)輸沙模型所需數(shù)據(jù)主要包括瑪曲水文站以上區(qū)域的DEM數(shù)據(jù)、水系拓?fù)湫畔?、水文氣象?shù)據(jù)和土地利用類型數(shù)據(jù)等。DEM數(shù)據(jù)為源自地理空間數(shù)據(jù)云網(wǎng)站的30 m精度的ASTER GDEM數(shù)據(jù)。水系拓?fù)湫畔樵醋試业乩硇畔?shù)據(jù)庫的1∶250 000空間分辨率的流域水系圖。日降雨量、日平均氣溫、日照時(shí)數(shù)、相對(duì)濕度以及風(fēng)速等水文氣象數(shù)據(jù)源自中國氣象局和黃河水利委員會(huì)水文局。用于對(duì)模型模擬結(jié)果進(jìn)行率定和驗(yàn)證的實(shí)測(cè)數(shù)據(jù)源自《中華人民共和國水文年鑒》，主要包括1956—2018年瑪曲水文站實(shí)測(cè)的月平均流量和月平均輸沙率。土地利用類型數(shù)據(jù)源自中國多時(shí)期土地利用土地覆被遙感監(jiān)測(cè)數(shù)據(jù)集（CNLUCC）［7］。河道初始泥沙中值粒徑由模型的輸入文件設(shè)定，數(shù)據(jù)源自《黃河泥沙公報(bào)》中多年平均泥沙中值粒徑，其在模型中作為可調(diào)參數(shù)進(jìn)行模型率定。

2 研究方法

2.1 WEP－SED模型原理

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

圖2 模型模擬產(chǎn)輸沙過程示意

在模型模擬計(jì)算谷溝環(huán)節(jié)和河道環(huán)節(jié)的產(chǎn)輸沙量時(shí)，水流挾沙力和泥沙沉速等參數(shù)的計(jì)算極其重要，并且這2個(gè)參數(shù)與泥沙粒徑大小具有一定關(guān)系。水流挾沙力公式如下［20］：

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

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

式中：D為泥沙沉降粒徑，mm，在模型計(jì)算時(shí)等同于泥沙中值粒徑d50；ρs為泥沙密度，g／cm3；ρw為清水密度，g／cm3；ν為水的運(yùn)動(dòng)黏滯系數(shù)，cm2／s。

2.2 泥沙動(dòng)態(tài)粒徑計(jì)算方法

基于WEP－SED模型的基本原理，針對(duì)模型中谷溝環(huán)節(jié)和河道環(huán)節(jié)的產(chǎn)輸沙量計(jì)算，把不同泥沙來源的沙量作為權(quán)重對(duì)泥沙中值粒徑進(jìn)行迭代計(jì)算，從而模擬泥沙中值粒徑隨沖淤變化的動(dòng)態(tài)變化過程，計(jì)算步驟如下。

步驟一：基于WEP－SED模型采用初始泥沙中值粒徑分別進(jìn)行坡面、谷溝、河道環(huán)節(jié)的初次產(chǎn)輸沙量計(jì)算。

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

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

步驟三：根據(jù)步驟二計(jì)算的谷溝環(huán)節(jié)的綜合泥沙中值粒徑，第二次進(jìn)行谷溝環(huán)節(jié)的產(chǎn)輸沙量計(jì)算；再重復(fù)步驟二第二次進(jìn)行谷溝環(huán)節(jié)的綜合泥沙中值粒徑計(jì)算，然后以該結(jié)果第三次進(jìn)行谷溝環(huán)節(jié)的產(chǎn)輸沙量計(jì)算；最終得到谷溝環(huán)節(jié)的產(chǎn)輸沙量和綜合泥沙中值粒徑。

步驟四：根據(jù)步驟一計(jì)算的河道環(huán)節(jié)的產(chǎn)輸沙量判斷其是否發(fā)生沖刷，若發(fā)生沖刷，則以進(jìn)入河道的坡面產(chǎn)輸沙量、谷溝環(huán)節(jié)最終的產(chǎn)輸沙量、上游河道的產(chǎn)輸沙量以及本級(jí)河道的沖刷量作為權(quán)重計(jì)算河道環(huán)節(jié)的綜合泥沙中值粒徑；若不發(fā)生沖刷，則選取的權(quán)重不包括本級(jí)河道的沖刷量。

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

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

河道環(huán)節(jié)的綜合泥沙中值粒徑計(jì)算公式如下：

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

步驟五：根據(jù)步驟四計(jì)算的河道環(huán)節(jié)的綜合泥沙中值粒徑第二次計(jì)算河道環(huán)節(jié)的產(chǎn)輸沙量；再重復(fù)步驟四第二次計(jì)算河道環(huán)節(jié)的綜合泥沙中值粒徑，以該結(jié)果第三次計(jì)算河道環(huán)節(jié)的產(chǎn)輸沙量；最終得到河道環(huán)節(jié)的產(chǎn)輸沙量和綜合泥沙中值粒徑。

3 結(jié)果與分析

利用WEP－SED模型對(duì)黃河河源至瑪曲水文站的水沙運(yùn)移過程進(jìn)行模擬，通過計(jì)算不同來源的泥沙中值粒徑分析其對(duì)流域產(chǎn)輸沙過程模擬結(jié)果的影響，分別采用納什系數(shù)（NSE）、相對(duì)誤差（RE）對(duì)模擬結(jié)果進(jìn)行評(píng)價(jià)。

3.1 動(dòng)態(tài)泥沙中值粒徑對(duì)模型模擬輸沙率的影響

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

圖3 采用泥沙動(dòng)態(tài)粒徑時(shí)瑪曲水文站1956—2018年逐月輸沙率模擬結(jié)果

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

3.2 泥沙中值粒徑的動(dòng)態(tài)變化過程

為明確WEP－SED模型模擬汛期輸沙率提升的根本原因，分析采用泥沙動(dòng)態(tài)粒徑計(jì)算方法時(shí)模型模擬的泥沙中值粒徑逐日變化過程，并對(duì)采用泥沙動(dòng)態(tài)粒徑計(jì)算方法的泥沙中值粒徑逐月變化過程進(jìn)行統(tǒng)計(jì)，其結(jié)果見圖5?？梢钥闯?，采用泥沙動(dòng)態(tài)粒徑計(jì)算方法實(shí)現(xiàn)了對(duì)泥沙中值粒徑動(dòng)態(tài)變化過程的模擬。由圖5（b）可知，采用泥沙動(dòng)態(tài)粒徑計(jì)算方法時(shí)模型模擬的汛期泥沙中值粒徑小于非汛期泥沙中值粒徑，因此模型模擬結(jié)果體現(xiàn)了汛期更易發(fā)生沖刷過程，實(shí)現(xiàn)了汛期輸沙率模擬精度的提高，該結(jié)論與馬勇等［15］、全棟［16］、孫維婷等［17］、薛博文等［21］、那巍等［22］研究黃河流域?qū)崪y(cè)輸沙率的年內(nèi)變化規(guī)律呈現(xiàn)出一致性。

圖5 河道泥沙中值粒徑年內(nèi)變化情況

4 結(jié) 論

本文以黃河干流瑪曲水文站以上區(qū)域?yàn)檠芯繉?duì)象，基于WEP－SED模型，考慮不同季節(jié)泥沙來源不同的特點(diǎn)進(jìn)行泥沙中值粒徑的計(jì)算，實(shí)現(xiàn)了對(duì)泥沙中值粒徑動(dòng)態(tài)變化過程的模擬。與采用泥沙固定粒徑計(jì)算方法相比，采用泥沙動(dòng)態(tài)粒徑計(jì)算方法使河源至瑪曲段1956—2018年逐月輸沙率的模擬精度有所提升，逐月輸沙率模擬結(jié)果的NSE值從0.601提升至0.706，并且使得模型模擬的泥沙中值粒徑呈現(xiàn)汛期小、非汛期大的動(dòng)態(tài)變化規(guī)律，模型模擬的泥沙中值粒徑逐月變化規(guī)律與相關(guān)研究結(jié)果及實(shí)測(cè)資料一致。