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

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

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

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

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

1 材料與方法

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

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

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

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

1.2 模擬評價相關(guān)說明

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

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

(1)

(2)

(3)

(4)

(5)

(6)

2 模型構(gòu)建

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

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

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

(7)

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

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

(8)

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

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

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

表1 氮素輸出經(jīng)驗?zāi)Ｐ驮u價

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

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

(9)

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

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

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

(10)

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

泥沙輸移方程可表示為

(11)

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

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

圖4 綜合模型框架

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

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

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

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

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

表3 參數(shù)修正前后模擬評價

表4 各土層飽和水力傳導度校正前后比較

4 結(jié) 論

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

圖5 模型驗證期動態(tài)模擬結(jié)果

數(shù)據(jù)類別樣本數(shù)量評價指標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ù)學模型，并可定量描述壤中流和地表徑流耦合下的土壤氮素遷移輸出過程。這將為紅壤坡地水土流失預測和農(nóng)業(yè)面源污染防治提供計算工具。受限于時間，模型研究中還需要進一步加強和完善的有以下方面：一是今后應(yīng)加強觀測，減少漏測現(xiàn)象，提高觀測精度(如土壤水分、泥沙量等)，因為觀測數(shù)據(jù)的質(zhì)量在很大程度上決定了模型及參數(shù)的準確程度；二是模型自身缺陷還需修補，例如HYDRUS-2D模型遇強降雨時不收斂的問題，可以通過修改其源代碼數(shù)值格式進行修補。