黃土區裸露坡地徑流養分流失模型的建立與驗證*

邵凡凡，吳軍虎，李玉晨

（西安理工大學省部共建西北旱區生態水利國家重點實驗室，西安 710048）

溶質從表層土壤進入地表徑流的過程非常復雜。雨滴擊濺、土壤中的養分濃度以及對流、擴散作用和土壤顆粒吸附均會對這一過程產生影響[1-3]。Ahuja等[4]1981年通過在飽和土層不同深度位置放置一定量的32P，發現表層土壤的溶質進入徑流的可能性最高，隨著土層深度的增加，進入徑流的概率呈指數遞減。然后，Ahuja和Lehman[5]在1983年發現，觀察到的交換層深度較通過擬合模型和實測數據獲得的深度小得多。結合交換層理論和Rose土壤侵蝕模型，Gao等[6]在2003年建立了基于物理過程并考慮雨滴擊濺和擴散作用的溶質運移模型；該模型的所有參數均可通過實驗測量獲得。鑒于黃土高原地區養分流失的特點，王全九等[7]提出了一種新的方法來改進等效對流傳質模型。在此基礎上，Dong等[8]假設交換層被混合層代替，且交換率被雨滴誘導水分轉移率所代替。Yang等[9]結合質量守恒方程和降雨對表層土壤的剝蝕過程，建立了可預測黃土高原地表養分隨徑流流失的數學模型。但是，這些模型無法描述徑流發生之前土壤溶質濃度的變化，并且僅能通過擬合曲線來獲得初始土壤溶質濃度。Tong等[10]基于質量守恒方程和水平衡方程建立了一維兩層的溶質運移模型。該模型結合了入滲和擴散作用，并由不完全混合參數來描述。基于該模型，Tong等[11]使用集合卡爾曼濾波數據同化方法（EnKF）來校準參數并更新可溶性化學物質從土壤至地表徑流的轉移過程，并消除了實驗觀測數據的誤差。因此，基于混合層理論的模型由于具有明晰的物理意義而被廣泛用于預測斜坡上的溶質運移。

養分流失的模擬是在徑流模擬的基礎上進行的。通常使用 Saint-Venant方程（即連續性方程和動量方程）來描述地表徑流過程[12]。然而，由于Saint-Venant方程是高度非線性的，很難獲得解析解，這意味著僅能使用數值方法對其進行求解。但當忽略 Saint-Venant方程的加速度項時，可使用擴散波方程來對其進行簡化[13]。若同時忽略Saint-Venant方程的加速度和壓力項時，Saint-Venant方程可表示為運動波方程。Luce和Cundy[14]通過使用菲利普（Philip）入滲方程修改了運動波方程來預測超滲降雨條件下的產流過程。Yang等[15]通過假設水深與入滲率之間的線性關系簡化運動波模型中的水深項，并結合Philip入滲方程得到了運動波模型的近似解析解。該模型因其參數簡單易獲取而被廣泛使用于坡面徑流的模擬。在養分隨徑流流失的研究中，Gao等[16]在2004年提出了一個基于溶質守恒方程的模型，該模型考慮了雨滴飛濺和徑流沖刷作用。但是，該模型僅用于模擬積水條件下飽和土壤的養分流失過程，這與黃土區初始非飽和土壤條件下的流失過程存在較大差異。因此，本文以 Yang等[15]建立的坡面徑流近似解析解為基礎，進一步延伸于養分隨地表徑流流失過程的模擬中，并修改了Gao等[16]的模型以適應本文的測試條件。通過模型參數分析揭示了不同因素對養分流失的貢獻作用，提出了防止養分流失的有效措施。該研究可為防治農田退化和農業面源污染提供有力基礎。

1 材料與方法

1.1 理論與模型

1.1.1 坡面徑流運動過程 采用運動波模型來描述次降雨條件下的坡面水流流動過程[15]，其中超滲凈雨可用降雨強度與入滲率的差值來表示，見式（1）：

式中，h為徑流深，cm；t為徑流時間，min；q為單寬流量，cm2·min–1；x為坡面任一位置距離入流口的長度，cm；p為降雨強度，mm·h–1；i為土壤入滲率，cm·min–1。

由于坡面水深與入滲率之間存在關聯關系[9]，徑流水量為超滲凈雨所產生的，Yang等[15]用Philip公式表示降雨條件下的入滲過程，進一步求解了單寬流量和坡面水深，見式（2）和式（3）：

式中，c為入滲率參數[15]；S為土壤吸滲率，cm·min–0.5；Δt=3S2/（16p2）。

式中，n為曼寧糙率系數，；S0為水力梯度，本文中坡度為15°，故S0為Sin 15°。

1.1.2 徑流養分流失過程 降雨條件下土壤表層養分在雨滴擊濺和水分入滲的作用下隨徑流遷移并在土壤中重新分配，因此土壤剖面的水和養分運移系統自上而下可分為3層：徑流積水層、養分交換層和交換層以下土壤，如圖1所示。

養分交換層是徑流積水層與土壤剖面交界面以下厚度較薄的土層。交換層中化學物質的傳輸主要受入滲、水動力彌散和雨滴飛濺侵蝕控制[16-17]。

式中，de為交換層深度，cm；Ce為交換層中溶質濃度，mg·L–1；Cw為徑流中溶質濃度，mg·L–1；is為徑流水進入交換層的入滲速率，cm·min–1；ix為交換層水分進入更深土層的入滲速率，cm·min–1；er為雨滴誘導水分轉移速率，cm·min–1；λCw為徑流層進入交換層的溶質濃度（0 ≤λ≤ 1，Gao等[6]研究表明計算模型對該參數不敏感，取λ=0），mg·L–1；J為較深土壤層與交換層內部的溶質擴散通量，mg·cm–2·min–1。

為了簡化計算過程，Gao等[16]對J進行近似求解如下：

式中，Ds為養分在土壤中的擴散性，cm2·min–1；Cs為更深土層的溶質濃度，mg·g–1；γ為土壤容重，g·cm–3；K為土壤吸附系數，mL·g–1；β=er/（αde）。

從降雨開始，將整個降雨過程劃分為3個階段。

第一階段：從降雨開始t0至交換層完全飽和tsa。在這一階段，土壤入滲率為降雨強度，土壤表層未產生徑流，故i=p，q=0，交換層完全飽和的時間tsa可以表示為：

式中，tsa為交換層完全飽和所需時間，min；θs為飽和含水率，cm3·cm–3；θ0為初始含水率，cm3·cm–3。

第二階段：從交換層完全飽和tsa至土壤表層出現徑流tp。在這個階段，徑流層溶質濃度Cw和雨滴誘導水分轉移速率er的取值為0，ix = p。交換層中溶質濃度見式（9）：

式中，C0為初始溶質濃度，mg·L–1。

當式（9）中t=tp時，開始產流時刻交換層中的溶質濃度見式（10）：

式中，A=ix/（αde）。

第三階段：從開始產流至降雨結束。這一過程中，徑流中的養分濃度遠低于交換層中的養分濃度，因此為了簡化計算過程，忽略了徑流養分在入滲作用下對交換層養分的微小補給作用，ix=0.01 cm·min–1。結合這一階段的起始產流時間，即：t=tp，可求解得到交換層中溶質濃度表示如下：

式中，B=（ix+er）/（αde）。

產流過程中，徑流中化學溶質的質量守恒關系可表示為：

結合式（1）和式（12），可得到：

為了簡化計算，忽略了入滲和擴散作用，式（13）可表示為：

對式（14）進行積分得到徑流中化學溶質的濃度，見式（15）：

徑流中化學溶質的濃度與徑流量的乘積便為徑流中化學溶質的流失速率，結合式（2）、式（3）和式（15），徑流中溶質流失速率可表示為：

式中，Mw為養分隨徑流流失的速率，mg·min–1。

1.2 試驗區概況

試驗于2019年5月在中國科學院水利部水土保持研究所長武黃土高原農業生態試驗站擔水溝流域野外模擬降雨小區（35°12′N，107°10′E）進行，試驗區平均海拔為1 200 m，氣候屬暖溫帶半濕潤大陸性季風氣候，年平均氣溫為 9.1℃，年平均降水量580 mm，地下水位50～80 m，無灌溉條件，屬典型的旱作雨養農業區。該流域內典型土壤為粉砂質壤土，母質為深厚的中壤質馬蘭黃土，具體土壤物理特性見表1，流域塬面面積占35%，梁坡占35.6%，溝谷占29.4%，各約占1/3；流域地貌屬典型的黃土高原溝壑區[18-20]。

表1 試驗區土壤物理化學特性Table 1 Physical and chemical properties of the soil tested

1.3 試驗方法

試驗用地為3年閑置坡耕地，模擬降雨試驗的小區設置尺寸為1.0 m×1.0 m，并根據當地典型坡耕地坡度和坡面侵蝕的臨界坡度，設置小區坡度為15°；使用針孔式人工模擬降雨裝置進行降雨試驗（圖2），其主要由主體支架（可調節高度）、底板布設有針孔的水槽（可根據降雨強度更換不同孔徑的針孔）和供水裝置 3部分組成，有效降雨面積為1.0 m2。經過測試：該套人工降雨器的平均雨滴直徑為2 mm，降雨均勻度在80%以上，雨滴終速符合天然降雨特征[3]。

為了消除土壤前期含水率對試驗結果的影響，每次開始降雨24 h前以25 mm·h–1的降雨強度在試驗小區進行預降雨，直至開始產流時停止降雨。開始正式降雨試驗前測定小區內土壤表層 0～20 cm剖面的初始含水率，采用烘干法測定質量含水率為0.11± 0.003 7 g·g–1（ 即 體 積 含 水 率 為 0.15±0.005 cm3·cm–3）時開始試驗。為了提高土壤初始養分濃度值，使養分在土壤中均勻分布，預降雨結束后，在小區土壤表面均勻噴灑氯化銨和硝酸鉀混合溶液，其具體操作方式為：將預先配置好的5.0 g·L–1的氯化銨溶液和10.0 g·L–1的硝酸鉀溶液各取1 L混合均勻（為消除溶液噴施次序對其分布的影響），并用壓力噴壺在每個小區分別定量噴灑2 L氯化銨和硝酸鉀的混合溶液，為盡可能減小噴壺壓力對表層土壤造成的壓實作用，將噴嘴調節至霧化度最強位置處，噴嘴霧化半徑為5 cm，采用左右往復的方式將混合溶液均勻噴灑在小區表層，使其在2 cm范圍內均勻分布。并在小區坡面上、中、下三個部位分別取土樣來測定表層 2 cm土壤中的溶質濃度作為模型計算所需的初始土壤養分濃度（表2）。根據研究區暴雨實測資料及降雨分級標準[20]，設計30、45、60、75、90 mm·h–1的 5種降雨強度，設計總降雨歷時為120 min，按照0～10 min之間，每隔2 min承接1次徑流，10～120 min之間，每隔5 min承接1次徑流的頻率用量杯承接出口處徑流，并用量筒進行精確測量，通過沉淀過濾除去徑流中的泥沙，用50 mL的塑料瓶收集徑流水樣并存放于實驗室冰箱中，用全自動高通量間斷分析儀（SmartChem450，AMS Allinace公司，意大利）測定徑流和土壤中的養分濃度。

1.4 模型基本參數

通過現場測定和文獻查閱等方式獲取了模型計算所需參數[20-21]（表2）：

表2 模型基本參數Table 2 Basic parameters of the calculation model

1.5 數據處理

所有試驗實測數據均為 3次重復試驗的平均值，使用Matlab 2015b進行參數求解和模型模擬；使用SPSS 24.0進行數據分析，使用Origin 2018進行圖表繪制和函數擬合。

2 結果與討論

2.1 產流過程分析及模擬

2.1.1 產流過程分析 土壤表層的養分通常會隨著地表徑流而流失。因此，探究產流規律是模型準確預測養分隨徑流流失過程的基礎。在5種降雨強度下，開始產流的時間點分別為20.5、8.5、4.8、3.0和 1.8 min（圖3）；90 mm·h?1較 30 mm·h?1提前19 min產流；說明隨降雨強度的增大，開始產流時間開始顯著縮短。降雨強度與產流時間的關系可用冪函數來描述，R2=0.997 6。單寬流量開始產流后快速增大，而后進入穩定產流階段（圖4）。這是由于表層土體中的黏粒分散堵塞了土壤的孔隙，并伴隨著雨滴的飛濺使表層土壤變得密實，降低了土壤的入滲能力[13]。5種降雨強度下實測單寬流量的標準差分別為 0.01～0.16、0.01～0.16、0.10～0.41、0.13～0.71和0.23～0.69。這可能是由于土壤的非均質性、蟻穴和植物根系對小區土壤入滲過程的影響，以及雨滴飛濺和徑流發育過程中微地形的形成，可能導致徑流滯后。進一步分析發現，90 mm·h?1在穩定產流階段的單寬流量分別較其他降雨強度依次增加6.3倍、2.7倍、1.6倍和1.2倍，這表明降雨強度的增加顯著增大了坡面徑流率。

2.1.2 產流過程模擬 通過將已知參數S、p（表2）和實測單寬流量代入式（2）中推求入滲率參數c，采用R2、均方根誤差（RMSE）和納什效率系數（NSE）對模擬結果進行評價。從表3可看出，入滲率參數c隨著降雨強度的增大呈減小趨勢，且分布在0.003 1～0.006 0之間；R2均在0.89以上，隨著降雨強度的增大，RMSE也隨之增大，取值分布在0.406～1.052之間，NSE均大于0.397，而當降雨強度大于等于 60 mm·h–1時，NSE 則進一步增大至0.783以上；說明降雨強度越大，模型計算值和實測值的匹配度也隨之提高。指數函數可很好地擬合參數c與降雨強度之間的關系（圖5），決定系數R2為0.976 5，表達式為：

表3 入滲率參數c的最佳擬合值Table 3 Optimal fitting values of c，R2，RMSE and NSE relative to rainfall intensity

圖6分別顯示了5種降雨強度下的單寬流量模擬過程。可以看出，產流模型能夠較好地模擬地表徑流過程，且隨著降雨強度的增大，模擬趨勢變得更加準確。在產流初期，計算值的上升趨勢均慢于實測值；在穩定產流階段，30 和45 mm·h?1下的模擬值均大于實測值，而在 60、75、90 mm·h?1下，實測值與計算值的匹配程度較好。這可能是由于降雨強度較小時，雨滴動能的濺蝕增加了前期土壤表面粗糙度；雨滴擊濺形成的微地形和洼地攔截部分地面徑流，從而削減了連續徑流的沖刷作用。這表明產流模型可準確模擬大于等于 60 mm·h?1的產流過程。此外，雨滴的擊濺使表層土壤被壓實，容重增大，土壤表層形成密封層，降低了土壤的入滲能力。但該模型未考慮地表土壤容重和孔隙率的變化，導致模型計算的土壤入滲能力明顯大于實測值。

2.2 養分流失過程分析及模擬

2.2.1 養分隨徑流流失過程分析 養分從土壤至徑流的傳輸是通過雨滴擊濺作用和徑流溶解作用來完成的[16]。不同降雨強度下硝態氮和銨態氮隨時間的流失過程可用單峰形式來描述；即：徑流初期養分流失速率迅速增大，到達峰值后開始減少，最后進入穩定流失階段的趨勢（圖7）。這是由于土壤水和土壤顆粒表面吸附的養分在雨滴擊濺作用下進入徑流所引起的；隨著土壤結皮厚度和徑流深度的增大，土壤表層形成“堅實的保護殼”削弱了雨滴動能，延緩土壤水和徑流的交換作用，使得進入徑流的土壤水和溶解態氮顯著減少。此外，徑流中氮素濃度的降低也是由于隨著降雨過程的推移，表層土壤中氮含量的逐漸減少所引起的。同時，降雨強度、入滲能力、養分濃度和徑流率對養分的峰值流失速率及其發生時間均有一定的影響。簡言之，硝態氮和銨態氮的峰值損失率隨降雨強度的增大而增大。以硝態氮流失過程為例，當降雨強度為 30 mm·h–1時，硝態氮流失速率在 25 min左右達到峰值5.74 mg·min–1，而當降雨強度為 45、60、75 和90 mm·h–1時，分別在 13、8、6和 5 min達到硝態氮流失速率的峰值：35.21、121.3、280.4和468.4 mg·min–1。因此，降雨強度對硝態氮的峰值流失速率具有較大貢獻。通過對比養分的峰值流失速率出現時間和穩定產流時間，各降雨強度下養分流失速率峰現時間分別為25、13、8、6、5 min，而穩定產流時間分別為40、28、18、16、14 min；由此看出，硝態氮流失速率的峰值出現時間要早于徑流速率達到穩定階段所需的時間，這可能是隨著產流時間的推移，交換層土壤中硝態氮濃度的不斷減小和徑流量的增大共同作用所造成的。5種降雨強度下硝態氮損失率的標準誤差分別分布在0.03～0.53、0.05～5.13、0.03～6.9、0.1～21.0 和 1.1～31.7。

2.2.2 養分流失過程模擬 交換層深度de和雨滴誘導水分轉移速率er是養分流失模型中的兩個重要參數。由于受室外實驗條件的限制，交換層深度de很難通過實地測量得到，因此需借助模型擬合實測的養分流失速率來反推交換層深度de。研究[4-5]發現交換層的深度de在2～3 mm的范圍內。Tong等[11]指出，交換層深度隨入滲率的增加而減小。有研究[21-23]指出，交換層深度隨著初始含水量的增大而增大。關于雨滴誘導水分轉移速率er，Gao等[16]在2004年提出了適用于初始飽和土壤的雨滴誘導水分轉移率er的計算方法，由于黃土區坡耕地在降雨前為非飽和土壤，因此，借助 Matlab非線性擬合的方法，將式（3）所計算出的坡面徑流深度代入式（16）來計算養分流失速率，并結合硝態氮流失速率的實測值推求出了式（16）中的交換層深度de和雨滴誘導水分轉移速率er（表4），并進一步模擬了銨態氮流失過程。可以看出：de和er均隨著降雨強度的增大而增大，其分別從0.68增至1.32、從0.006增至0.023。這與前述單寬流量隨著降雨強度的增大而增加是一致的[24]。硝態氮和銨態氮流失速率的R2值分別分布于 0.834～0.922和 0.800～0.921之間，RMSE值分別分布在1.188～58.50和0.974～58.37之間，NSE值分別分布在0.653～0.881和0.546～0.775之間。從圖8中可以看出，計算出的曲線可很好地模擬養分流失過程。當降雨強度為30 mm·h–1時，初始增大階段的測量值與計算值之間的差異較大，但隨著降雨強度的增大，差異逐漸減小。而在養分流失的穩定階段，實測數據大于計算值。這可能是由于本文建立的養分流失模型近似求解了對流彌散項所造成的，這使得穩定減小階段養分流失速率衰減得過快[9]。同時模型忽略了徑流層養分對交換層的微弱補給作用，然而在降雨開始時可能存在從交換層至徑流層的擴散過程[16]。從交換層完全飽和（tsa）到坡面開始產流（tp）的時間段內，模型假設交換層中的養分隨入滲水向土壤深處遷移的速率大小即為對應的降雨強度，從而使得計算出的開始產流時土壤表層的濃度Ce（tp）小于理論值。

表4 不同降雨強度下參數de和er的最佳擬合值Table 4 Optimal fitting values of de，er，R2，RMSE and NSE relative to rainfall intensity

以上分析表明本文建立的模型可以很好地模擬裸露坡面徑流和養分隨徑流遷移過程。但該模型未用于模擬不同坡度、坡長和土壤初始含水量條件下的養分流失過程，本研究所獲得的參數是否具有普遍適用性，需要在以后的研究中加以驗證。可以預見的是，坡長的增大將顯著增加徑流量和泥沙量，初始含水量的增大將提前產流時間并增加養分的峰值流失速率[25，26]。同時，土壤中植物根系的生長和土壤生物活性可能形成連通的土壤孔隙結構，導致優先流的出現，這將對模型的模擬精度產生較大影響。泥沙顆粒中通常吸附有大量養分，但本研究的模型并未將泥沙考慮在內，使得模型并不能完整模擬徑流過程所帶走的養分總量。簡言之，該近似解析模型充分考慮了非飽和土壤水分入滲對交換層中養分運移過程的影響，因此，該模型可用于預測干旱和半干旱氣候條件下裸露坡耕地的養分流失過程。但是，徑流過程的精確計算是進行養分流失模擬的基礎，應根據土壤質地、養分類型和雨水中養分濃度選擇合適的入滲公式和溶質吸附系數。

3 結 論

本研究以交換層理論為基礎，根據黃土區降雨量少，降雨前土壤通常為非飽和狀態，其產流需要較長時間的實際情況對降雨過程進行劃分，建立了基于坡面徑流養分遷移理論的機理模型，并通過 5個降雨強度的模擬降雨試驗對模型進行了驗證。3組實測重復試驗間的標準誤差均較小，試驗結果具有可靠性。模型驗證的結果表明，本文建立的養分流失近似解析模型能夠準確描述不同降雨強度下的坡面流和養分流失特征（R2> 0.8，NSE > 0.347）。參數c（入滲率參數）、de（交換層深度）、er（雨滴誘導水分轉移速率）均隨降雨強度的增大而增大。養分流失模型對交換層的深度de較雨滴誘導水分轉移速率er更敏感，de可顯著影響可交換溶質的量。因此，在施肥過程中應采取一些措施，如施肥后覆蓋坡面土壤、暴雨前避免施肥等，以達到防止土壤貧瘠化并控制農業面源污染的目的。