劉 昭，徐燕星，聶小飛，胡 皓，鄭海金

(1.江西省土壤侵蝕與防治重點(diǎn)實(shí)驗(yàn)室，江西 南昌 330029；2.江西省水土保持科學(xué)研究院，江西 南昌 330029； 3.江西水利職業(yè)學(xué)院，江西 南昌 330013)

坡耕地作為南方紅壤區(qū)的重要農(nóng)業(yè)生產(chǎn)資源被高強(qiáng)度開(kāi)發(fā)利用，加之南方紅壤區(qū)易出現(xiàn)集中強(qiáng)降雨情況，因此極易發(fā)生水土流失，造成大量氮素流失，而該區(qū)域特殊的地形和土壤理化性質(zhì)，使得驅(qū)動(dòng)土壤氮素運(yùn)移的地表-土壤水文過(guò)程和土壤侵蝕過(guò)程具有十分復(fù)雜的特征。地表徑流和土壤水運(yùn)動(dòng)過(guò)程相互作用特征明顯，氮素遷移隨地表徑流、泥沙及土壤水文過(guò)程在時(shí)間和空間上呈動(dòng)態(tài)變化，是地表徑流、土壤水運(yùn)動(dòng)和泥沙輸移等多種過(guò)程耦合作用的結(jié)果。目前國(guó)內(nèi)外對(duì)坡地氮素流失的研究主要集中在徑流流失方面，對(duì)地表徑流氮素遷移過(guò)程的影響因素和相關(guān)的調(diào)控措施進(jìn)行了深入的分析[1-2]，但對(duì)氮素流失與水文泥沙過(guò)程(尤其是土壤水文過(guò)程)的相互關(guān)系尚未深入研究。因此，構(gòu)建氮素隨水文泥沙過(guò)程遷移的數(shù)學(xué)模型，將有助于防治坡地氮素流失和減少水體污染。

土壤水流、地表徑流和泥沙輸移具有明確的動(dòng)力學(xué)機(jī)制，目前已有大量模型考慮了這些動(dòng)力學(xué)過(guò)程。在土壤水流運(yùn)動(dòng)方面，以達(dá)西定律和質(zhì)量守恒定律推導(dǎo)的Richards方程具有堅(jiān)實(shí)的物理基礎(chǔ)，對(duì)于水分通量計(jì)算具有較大的優(yōu)勢(shì)[3]，在各種水文模型中得到了廣泛應(yīng)用，如HYDRUS[4]、SWAP[5]等，但是這些模型以土壤水流模擬為主，很少詳細(xì)考慮地表水沙動(dòng)力學(xué)過(guò)程。由經(jīng)典的Navier-Stokes流體運(yùn)動(dòng)方程，可以得到在坡面地表徑流模擬中廣泛使用的二維淺水方程、一維Saint Venant方程，以及更為簡(jiǎn)化的運(yùn)動(dòng)波近似方程；對(duì)于泥沙輸移模擬，一般結(jié)合地表水動(dòng)力學(xué)模擬結(jié)果，采用泥沙輸移方程求解。目前，也有許多基于這類(lèi)地表水沙動(dòng)力學(xué)方程的模型，如國(guó)外的WEPP模型[6]、EROSEM模型[7]等，國(guó)內(nèi)的雷廷武模型[8]、龍滿生模型[9]等。但是上述這類(lèi)地表水沙模型往往忽略土壤水分過(guò)程或以經(jīng)驗(yàn)?zāi)Ｐ瓦M(jìn)行描述。因此，坡地土壤水沙動(dòng)力學(xué)模型構(gòu)建的關(guān)鍵點(diǎn)之一是如何耦合土壤水流、地表徑流和泥沙輸移的動(dòng)力學(xué)過(guò)程。

本研究利用江西水土保持生態(tài)科技園大型土壤水分滲漏裝置的地表徑流、泥沙、壤中流及其中的氮素含量實(shí)測(cè)數(shù)據(jù)，構(gòu)建氮素輸出與水文泥沙過(guò)程的數(shù)學(xué)關(guān)系，耦合土壤水分運(yùn)移模型HYDRUS-2D和坡面土壤侵蝕模型WEPP，構(gòu)建可以動(dòng)態(tài)模擬水文泥沙過(guò)程的水沙動(dòng)力學(xué)模型，最終形成一套同時(shí)可定量描述地表徑流、泥沙、壤中流、氮素輸出過(guò)程的綜合模型，以期為紅壤坡地水土流失預(yù)測(cè)和農(nóng)業(yè)面源污染防治提供計(jì)算工具。

1 材料與方法

1.1 試驗(yàn)裝置及數(shù)據(jù)采集

試驗(yàn)場(chǎng)地設(shè)置在江西水土保持生態(tài)科技園，園區(qū)氣候、地貌、土壤等自然條件在南方紅壤區(qū)均具有代表性。園區(qū)屬于亞熱帶季風(fēng)氣候區(qū)，多年平均降水量1 469 mm，多年平均氣溫16.7 ℃；地貌為淺丘崗地，坡度5°～25°，海拔30～100 m；土壤主要為第四紀(jì)紅黏土發(fā)育的紅壤，呈酸性至微酸性。試驗(yàn)裝置為大型土壤水分滲漏裝置(裸地)，如圖1所示。小區(qū)坡度為14°，水平投影長(zhǎng)15 m、寬5 m。試驗(yàn)裝置的周?chē)暗装逵娩摻罨炷翝仓?，地板上設(shè)砂礫反濾層，坡腳修筑梯形鋼筋混凝土擋土墻，擋土墻高出地表30 cm以阻止水分進(jìn)出小區(qū)，從而形成一個(gè)封閉排水式的土壤入滲裝置。小區(qū)土壤為第四紀(jì)紅黏土發(fā)育的紅壤，土層深105 cm，剖面為三層，其中0～30 cm為Ah層(淋溶層)，60 cm以下為Bsv層(淀積層)，30～60 cm為Bs層(過(guò)渡層)。坡底設(shè)置地表徑流和30、60、105 cm深壤中流的收集裝置。

圖1 土壤水分滲漏裝置示意

2015年5月22日在小區(qū)施尿素300 kg/hm2，施肥后進(jìn)行逐場(chǎng)次降雨的徑流、泥沙和氮素輸出試驗(yàn)觀測(cè)。降雨量采用試驗(yàn)區(qū)內(nèi)設(shè)置的虹吸式自記雨量計(jì)和全自動(dòng)氣象采集系統(tǒng)進(jìn)行監(jiān)測(cè)，可以獲得每次降雨的降雨量和降雨歷時(shí)；地表徑流量通過(guò)地表徑流池池壁水尺讀數(shù)計(jì)算；不同層次壤中流量通過(guò)自記水位計(jì)數(shù)據(jù)讀出；產(chǎn)流結(jié)束后，從地表徑流池取樣并分析徑流中泥沙和總氮含量，從不同層次壤中流收集裝置中取等量水樣充分混勻后分析其總氮含量，泥沙量測(cè)定采用稱重法進(jìn)行，總氮的測(cè)定采用堿性過(guò)硫酸鉀消解紫外分光光度法。至2016年4月22日，共收集到21場(chǎng)次降雨的地表徑流數(shù)據(jù)，由于地表徑流產(chǎn)流歷時(shí)短，因此記錄數(shù)據(jù)為每場(chǎng)降雨發(fā)生的累積地表徑流量；而壤中流持續(xù)時(shí)間很長(zhǎng)(尤其是底層壤中流)，故壤中流數(shù)據(jù)為日記錄數(shù)據(jù)，共337組。因人工采樣受限于時(shí)間和設(shè)備等因素，故泥沙量、地表氮素輸出量和壤中氮素輸出量共有10、12、14場(chǎng)次降雨的數(shù)據(jù)。

1.2 模擬評(píng)價(jià)相關(guān)說(shuō)明

本研究將2015年5月22日至2016年1月31日間的數(shù)據(jù)用于經(jīng)驗(yàn)?zāi)Ｐ蜆?gòu)建及機(jī)理性模型參數(shù)校正，這段時(shí)期稱為模型校正期，期間地表徑流量、壤中流量、泥沙量、地表氮素輸出量和壤中氮素輸出量數(shù)據(jù)分別有16、255、6、8和9組；而將2016年2月1日至2016年4月22日間的數(shù)據(jù)用于模型驗(yàn)證，這段時(shí)期稱為模型驗(yàn)證期，期間地表徑流量、壤中流量、泥沙量、地表氮素輸出量和壤中氮素輸出量分別有5、82、4、4和5組樣本數(shù)據(jù)。

研究中為了判別模擬值與觀測(cè)值之間的擬合精度，需要引入統(tǒng)計(jì)學(xué)指標(biāo)。常用于判別擬合精度的統(tǒng)計(jì)學(xué)指標(biāo)有決定系數(shù)R2(R2∈[0,1])、均方根誤差RMSE(RMSE∈[0,∞))、偏差BIAS(BIAS∈[0,∞))、一致性指標(biāo)IA(IA∈[0,1])等。其中：BIAS表示的是平均誤差與觀測(cè)均值之比，與R2、IA都是相對(duì)指標(biāo)，無(wú)量綱；RMSE的量綱與統(tǒng)計(jì)對(duì)象相同。各指標(biāo)計(jì)算公式為

(1)

(2)

(3)

(4)

(5)

(6)

2 模型構(gòu)建

2.1 氮素輸出經(jīng)驗(yàn)?zāi)Ｐ?/h3>

如圖2所示，分析模型校正期有效觀測(cè)數(shù)據(jù)發(fā)現(xiàn)，壤中氮素輸出量與壤中流量顯著相關(guān)(相關(guān)系數(shù)為0.95)。在實(shí)際中，若無(wú)壤中流產(chǎn)生，則氮素不會(huì)從土壤中輸出，因此本研究采用無(wú)截距的二次函數(shù)表示壤中氮素輸出量與壤中流量的關(guān)系，即

圖2 壤中氮素輸出經(jīng)驗(yàn)?zāi)Ｐ?/p>

(7)

式中：Ng為壤中氮素輸出量，g；Qg為壤中流量，mm。

相關(guān)性分析表明，地表氮素輸出量與泥沙量的相關(guān)系數(shù)達(dá)到0.92，而地表氮素輸出量與地表徑流量的相關(guān)系數(shù)為0.78，可見(jiàn)地表氮素輸出量與土壤侵蝕關(guān)系更為密切。通過(guò)分析，本研究采用冪函數(shù)表示地表氮素輸出量與泥沙量的關(guān)系，即

(8)

式中：Ns為地表氮素輸出量，g；Me為泥沙量，kg。

圖3 地表氮素輸出經(jīng)驗(yàn)?zāi)Ｐ?/p>

由表1可知，式(7)和式(8)均具有較高的模擬精度。

表1 氮素輸出經(jīng)驗(yàn)?zāi)Ｐ驮u(píng)價(jià)

2.2 壤中流和地表徑流耦合下氮素輸出模型構(gòu)建

Richards方程由達(dá)西定律和質(zhì)量守恒定律推導(dǎo)而來(lái)，對(duì)于水分通量計(jì)算具有較大的優(yōu)勢(shì)，可用于描述飽和-非飽和土壤水流動(dòng)力學(xué)過(guò)程。本研究假設(shè)垂直于坡向上的水流過(guò)程是相同的，則土壤水流動(dòng)力學(xué)過(guò)程以二維Richards方程表示，即

(9)

式中：t為時(shí)間，h為土壤壓力水頭；K(h)為水力傳導(dǎo)度；C(h)為容水度；β為參數(shù)，在土壤非飽和時(shí)為0，飽和時(shí)為1；μs為彈性釋水系數(shù)；x為坡向水平投影方向坐標(biāo)；z為距離地表的深度，向下為正。

式(9)的定解條件可根據(jù)實(shí)際情況進(jìn)行設(shè)置，下邊界根據(jù)地下水位情況設(shè)為定水頭邊界或者根據(jù)不透水層位置設(shè)為0通量邊界；在降雨條件下，上邊界可設(shè)為給定通量邊界，即降雨量與地表徑流量之差。以混合型Richards方程為基礎(chǔ)，采取有限單元法進(jìn)行時(shí)空離散，采取Picard迭代方法進(jìn)行數(shù)值求解，能夠模擬土壤水分、壓力和邊界通量動(dòng)態(tài)過(guò)程，具有較高的模擬精度，在農(nóng)田水利、水文水資源及地下水污染防治方面得到了廣泛應(yīng)用，因此本研究選取HYDRUS-2D模型對(duì)式(9)進(jìn)行求解。

對(duì)于徑流泥沙過(guò)程的模擬，可參考美國(guó)農(nóng)業(yè)部開(kāi)發(fā)的WEPP模型。其單次降雨模塊以一維運(yùn)動(dòng)波方程模擬地表徑流，描述為

(10)

式中：l為某點(diǎn)沿下坡方向的距離，h為地表水深，r為凈雨強(qiáng)，i為入滲速率，θ為坡面與水平面夾角，q為地表徑流通量。

泥沙輸移方程可表示為

(11)

式中：G為輸沙量；DL為從相鄰坡面流入的泥沙量；DF為水流對(duì)溝道的剝蝕量或水流中泥沙的沉積量。對(duì)于裸土而言，當(dāng)確定了土壤相關(guān)參數(shù)、土壤含水量和氣象因素后，即可借助WEPP模型對(duì)式(10)和(11)進(jìn)行求解。

為保證上邊界通量和土壤水分計(jì)算的一致性，本研究所采取的耦合方法如圖4所示，即：將HYDRUS-2D模型土壤水分模擬結(jié)果作為WEPP模型的輸入條件，然后以WEPP模型計(jì)算的徑流量分割降雨量(扣除蒸發(fā)量后)，從而得到HYDRUS-2D模型的上邊界通量；以FORTRAN 95程序編寫(xiě)WEPP模型和HYDRUS-2D模型的輸入輸出文件交互接口，并控制程序運(yùn)行，模擬出壤中流量和侵蝕泥沙量后，根據(jù)式(7)和式(8)即可得到壤中和地表氮素的輸出量。

圖4 綜合模型框架

3 參數(shù)校正和模型檢驗(yàn)

采用圖4所示模型進(jìn)行模擬，還需要確定機(jī)理性模型的參數(shù)。對(duì)于模型參數(shù)初值，一部分采用實(shí)測(cè)法獲取，其余則通過(guò)模型自帶的估算方法(如HYDRUS-2D中的土壤水動(dòng)力參數(shù)、WEPP中的土壤侵蝕參數(shù)等)獲取。其中實(shí)測(cè)參數(shù)見(jiàn)表2。

表2 實(shí)測(cè)的模型參數(shù)

研究中發(fā)現(xiàn)將上述參數(shù)初值直接用于模擬，會(huì)使得地表徑流、壤中流和泥沙的模擬值與實(shí)測(cè)值相差較大，因此需要進(jìn)一步校正參數(shù)。通過(guò)敏感性分析可知，對(duì)地表徑流量、壤中流量和泥沙量影響最大的為飽和水力傳導(dǎo)度(Ks)，其余參數(shù)影響相對(duì)較小，因此主要是對(duì)各層飽和水力傳導(dǎo)度進(jìn)行了校正。如前所述，模型校正期選為2015年5月22日至2016年1月31日，參數(shù)校正前后模擬值與觀測(cè)值的比較結(jié)果見(jiàn)表3，校正后的參數(shù)模擬精準(zhǔn)度得到了很大程度的提高(R2和IA增加，BIAS和RMSE減少)。以校正后的參數(shù)Ks進(jìn)行模擬，壤中流的模擬精度最高，泥沙量的模擬精度相對(duì)較差。值得注意的是，各土層飽和水力傳導(dǎo)度實(shí)測(cè)值與校正值可相差1.4～4.7倍(見(jiàn)表4)，這是導(dǎo)致模擬值與實(shí)測(cè)值相差較大的根本原因。

進(jìn)行參數(shù)校正后，即可對(duì)模型驗(yàn)證期(2016年2月1日至2016年4月22日)的地表徑流量、 壤中流量、泥沙量、地表氮素輸出量、壤中氮素輸出量進(jìn)行模擬。如圖5所示，模擬結(jié)果不僅可以較好地匹配地表徑流量、壤中流量、泥沙量、地表氮素輸出量、壤中氮素輸出量等觀測(cè)值，還可以估計(jì)出未能實(shí)測(cè)到的數(shù)據(jù)。精度評(píng)價(jià)指標(biāo)計(jì)算結(jié)果(表5)表明，壤中流量和壤中氮素輸出量的模擬效果(決定系數(shù)達(dá)到0.84以上，一致性指標(biāo)達(dá)到0.96以上)要強(qiáng)于地表徑流量、泥沙量和地表氮素輸出量(決定系數(shù)為0.41～0.68，一致性指標(biāo)為0.58～0.87)。分析模擬結(jié)果(圖5)發(fā)現(xiàn)，在降雨產(chǎn)流期壤中流的氮素輸出量顯著高于地表徑流的氮素輸出量，且在未降雨期間氮素隨著壤中流持續(xù)輸出，說(shuō)明壤中流是紅壤坡地氮素的主要流失途徑，這與前人研究結(jié)論一致[1]。此外，壤中流產(chǎn)流高峰較降雨峰值延后，且一次降雨的壤中流會(huì)持續(xù)數(shù)天，說(shuō)明土壤起到了一定的緩沖作用。

表3 參數(shù)修正前后模擬評(píng)價(jià)

表4 各土層飽和水力傳導(dǎo)度校正前后比較

4 結(jié) 論

以實(shí)測(cè)數(shù)據(jù)分析為基礎(chǔ)，建立了氮素輸出經(jīng)驗(yàn)?zāi)Ｐ?，并將該?jīng)驗(yàn)?zāi)Ｐ团c土壤水運(yùn)動(dòng)模型(HYDRUS-2D)和坡面土壤侵蝕模型(WEPP)耦合，形成了一套模擬徑流-泥沙-壤中流-氮素輸出的綜合模型。以實(shí)測(cè)數(shù)據(jù)對(duì)模型參數(shù)進(jìn)行校正并進(jìn)行了模型檢驗(yàn)，校正后的參數(shù)對(duì)徑流量、泥沙量、壤中流量、氮素輸出量均具有較好的模擬效果(決定系數(shù)在0.41以上)，其中壤中流量和壤中氮素輸出量模擬效果更好， 決定系數(shù)分別達(dá)0.88和0.84。分析模擬結(jié)果表明，在降雨產(chǎn)流期壤中氮素輸出量顯著高于地表氮素輸出量，且在未降雨期間氮素隨著壤中流持續(xù)輸出，進(jìn)一步印證了壤中流是紅壤坡地氮素的主要流失途徑。此外，壤中流產(chǎn)流高峰較降雨峰值延后，且一次降雨的壤中流會(huì)持續(xù)數(shù)天，說(shuō)明土壤起到了一定的緩沖作用。

圖5 模型驗(yàn)證期動(dòng)態(tài)模擬結(jié)果

數(shù)據(jù)類(lèi)別樣本數(shù)量評(píng)價(jià)指標(biāo)R2RMSEBIASIA地表徑流量50.5140.9750.6130.585泥沙量40.6850.0060.3870.869地表氮素輸出量40.4100.0710.5720.648壤中流量820.8840.7650.2690.972壤中氮素輸出量50.8400.0640.1440.962

本研究已基本建立了具有一定物理意義、適合紅壤坡地產(chǎn)流的數(shù)學(xué)模型，并可定量描述壤中流和地表徑流耦合下的土壤氮素遷移輸出過(guò)程。這將為紅壤坡地水土流失預(yù)測(cè)和農(nóng)業(yè)面源污染防治提供計(jì)算工具。受限于時(shí)間，模型研究中還需要進(jìn)一步加強(qiáng)和完善的有以下方面：一是今后應(yīng)加強(qiáng)觀測(cè)，減少漏測(cè)現(xiàn)象，提高觀測(cè)精度(如土壤水分、泥沙量等)，因?yàn)橛^測(cè)數(shù)據(jù)的質(zhì)量在很大程度上決定了模型及參數(shù)的準(zhǔn)確程度；二是模型自身缺陷還需修補(bǔ)，例如HYDRUS-2D模型遇強(qiáng)降雨時(shí)不收斂的問(wèn)題，可以通過(guò)修改其源代碼數(shù)值格式進(jìn)行修補(bǔ)。