魯舒心，萬愉快，孫伯顏，農(nóng)睿，朱磊

（1.寧夏大學(xué)土木與水利工程學(xué)院，銀川 750021；2.寧夏回族自治區(qū)黃河水聯(lián)網(wǎng)數(shù)字治水重點實驗室，銀川 750021）

0 引 言

降雨徑流從坡面流下的過程中會發(fā)生下滲[1]，與此同時，溶質(zhì)也會隨著水分運動在地表與土壤間進行遷移轉(zhuǎn)化[2]。由水分運動驅(qū)動的溶質(zhì)遷移往往是坡地養(yǎng)分流失的主要原因，不僅降低了耕地肥力造成經(jīng)濟損失，還會導(dǎo)致農(nóng)業(yè)面源污染，產(chǎn)生河湖等地表地下水體污染和富營養(yǎng)化的風(fēng)險[3,4]。在水資源短缺，水污染尤其是面源污染問題日益嚴(yán)峻的當(dāng)下，研究地表水土壤水水分運動與溶質(zhì)遷移規(guī)律，對防治和解決水資源與水生態(tài)問題具有重要意義。

目前，有關(guān)地表水土壤水運動與溶質(zhì)運移過程的研究，大多是通過野外現(xiàn)場監(jiān)測和室內(nèi)模擬試驗相結(jié)合的方式進行。野外監(jiān)測能夠通過現(xiàn)場試驗直觀地理解水分和溶質(zhì)遷移過程[5]，獲取較為真實的數(shù)據(jù)資料，為數(shù)值模型的開發(fā)奠定基礎(chǔ)。室內(nèi)試驗的數(shù)據(jù)采集和監(jiān)測更加方便準(zhǔn)確，可以人為控制實驗條件并具有可重復(fù)性。近年來，隨著計算性能的提高和數(shù)值求解算法的改進，構(gòu)建數(shù)值模擬模型成為定量化研究地表水土壤水與溶質(zhì)遷移的一個重要手段[6]。相比于野外監(jiān)測和室內(nèi)試驗，數(shù)值模擬克服了試驗過程人為因素多、監(jiān)測手段受限等不足，不僅可以獲取整個研究區(qū)內(nèi)水分和溶質(zhì)濃度的空間分布狀態(tài)以及指定位置處的連續(xù)變化過程[7]，還可以高效地改變水力與介質(zhì)參數(shù)和邊界條件，對不同場景中的水分和溶質(zhì)運移做出預(yù)測和解析，已經(jīng)被國內(nèi)外眾多學(xué)者廣泛運用。

在應(yīng)用數(shù)學(xué)模型方法對地表水與土壤水進行模擬研究時，以往的研究中常將兩者分開單獨進行模擬，忽略了兩者間的相互作用[8]。而地表水土壤水耦合模擬則是從系統(tǒng)分析的角度考慮水分在地表和土壤兩部分中的運移轉(zhuǎn)化，選用一定的耦合方法對表述兩者的圣維南方程和連續(xù)性方程進行耦合。根據(jù)耦合方法的不同，可以將地表水土壤水水分耦合模型分為松散耦合模型、順序迭代耦合模型、完全耦合模型等三類[9,10]。松散耦合模型[11-13]通常先求解地表模塊，再求解地下水流，兩模塊間利用交換通量進行連接。此方法雖然計算效率較高，但由于采用了過多簡化計算，且未考慮土壤水對地表水的反饋作用，使其模擬精度較低并可能導(dǎo)致質(zhì)量平衡錯誤[10]，難以表述地表和土壤存在強交互作用時的水分運動。順序迭代耦合模型[14-16]通過將交換量作為源匯項，交替求解地表和地下模塊，實現(xiàn)了土壤水與地表水的雙向反饋，提高了模擬精度，但計算的復(fù)雜度增加，使得計算效率降低。完全耦合模型[1]則是用全局隱式格式將地表水土壤水方程進行統(tǒng)一離散并同步求解，能夠同時解決飽和模塊、非飽和模塊和界面邊界條件問題[17]，模擬精度得到進一步提高，但由于計算的復(fù)雜性，國內(nèi)外研究較少。在水量耦合研究的基礎(chǔ)上，加入溶質(zhì)運移方程，從物理機制上對地表地下水分運動與溶質(zhì)遷移進行模擬的全耦合模型研究則更鮮見。

國內(nèi)外對地表地下水溶質(zhì)耦合模型的研究多基于現(xiàn)有的商業(yè)軟件及其二次開發(fā)，如GMS，MODFLOW，Hydrus-2d等軟件及各軟件中的模塊組合以建立耦合模型：付強等[18]利用GMS-FEMWATER構(gòu)建了河流與地下水水流溶質(zhì)運移耦合模型，李躍鵬等[19]將GMS-Borehole模塊和MODFLOW模型結(jié)合，用以模擬淺層地下水水流，通過Hydrus-2D，RT3D，MT3DMS模型，建立了飽和-非飽和帶污染物聯(lián)合模擬模型。陳志慧等[20]通過構(gòu)建SWAT-MODFLOW-MT3DMS耦合模型對沙潁河流域地下地表水硝酸鹽通量過程做出了定量解析。此外，一些學(xué)者提出的部分水文模型也開發(fā)了溶質(zhì)耦合模擬功能，如將地表和地下水互為邊界條件結(jié)合菲克定律進行溶質(zhì)運移耦合的PAWS模型[21]；通過溶質(zhì)梯度計算地表地下交界處溶質(zhì)通量的WASH123D模型[22]；以及耦合了地表模型FLUXNoah、PIHM模型、生物地球化學(xué)模型（RT）的BioRT-Flux-PIHM模型等[23]。上述模型在求解溶質(zhì)運移問題時往往將地下水和地表水單獨處理，通過較為松散的結(jié)果將兩者進行耦合，導(dǎo)致地表與地下水相互影響的區(qū)域計算結(jié)果會有較大的偏差，進而影響溶質(zhì)運移的模擬精度。為提高交界處的水通量和溶質(zhì)交換量的計算準(zhǔn)確度，有學(xué)者提出了雙層結(jié)點法[24,25]，與公共節(jié)點法[26]相比，雙層結(jié)點法不依賴于兩個系統(tǒng)間的水頭連續(xù)性假設(shè)，而是通過達西流關(guān)系來描述地表與地下水體之間的交換量[27]，提供了更多的靈活性，在描述地表-地下界面處的水頭連續(xù)性更為準(zhǔn)確，計算效率更高[28]。

本文針對地表水土壤水水分運動與溶質(zhì)運移完全耦合模型研究相對較少的現(xiàn)狀，基于二維擴散波方程和傳統(tǒng)的三維Richards方程，采用雙層結(jié)點法進行耦合，使用二維和三維對流彌散方程分別描述地表水和土壤水中的溶質(zhì)運移，根據(jù)地表水和土壤水中的溶質(zhì)濃度與兩者間的交換水量計算溶質(zhì)的交換量，從而構(gòu)建地表水土壤水水分運動與溶質(zhì)運移全耦合模型。根據(jù)已發(fā)表文獻中不同算例的實驗結(jié)果，驗證本研究所提出模型的合理性和準(zhǔn)確性并分析降雨條件下地表水土壤水運動與溶質(zhì)遷移規(guī)律。

1 全耦合模型的構(gòu)建

1.1 地表水土壤水運動耦合模型

1.1.1 土壤水運動數(shù)學(xué)模型

使用經(jīng)典的變飽和土壤水分運動方程(Richards方程)表述地下區(qū)域的土壤水分運動[1,29]：

式中：θs為土壤飽和含水率，cm3/cm3；Sw為飽和度，其值大小等于土壤體積含水率θ和θs的比值；?為土壤水分負壓，cm；Qs為地表以下的源匯項，s-1；Qos為土壤水地表水交換量，s-1；xi(i=1，2，3)為三維笛卡爾空間坐標(biāo)，對于垂直截面的平面流動，x1=x表示水平坐標(biāo)，x1=z表示垂直坐標(biāo)，取向下為負；t表示時間，s；Kij為土壤飽和導(dǎo)水系數(shù)張量，cm/s；Krw為相對非飽和導(dǎo)水率。

1.1.2 地表水運動數(shù)學(xué)模型

采用二維擴散波方程描述地表徑流運動[30]：

式中：h0為地上水面高程，cm；h0=z0+d0，z0是地表高程，cm；d0是地表水深，cm；Qo是地表水的源匯項，cm/s；kox和koy分別為x、y方向上的地表水力傳導(dǎo)系數(shù)，cm/s。

1.1.3 地表水與土壤水運動的耦合

朱磊等[1]采用雙層結(jié)點法對土壤水運動模型和地表水運動模型進行全耦合，并通過對比分析模型試驗結(jié)果驗證了該方法的合理性和可靠性，因此本文采用其提出的全耦合模型模擬地表徑流和土壤水分運動，其關(guān)系式表示為：

兩者之間水量交換的驅(qū)動因素主要為水力梯度，交換水量Qos為[1]：

式中：kso為垂直方向上的土壤飽和滲透系數(shù)張量，cm/s；hs為土壤水水頭，cm；h0為地表水水頭，cm；lso為耦合長度，cm；對于kr，當(dāng)hs＞h0時，kr為地表的相對滲透率，當(dāng)hs＜h0時，kr為土壤相對非飽和導(dǎo)水率。

1.1.4 模型參數(shù)的處理

采用van Genuchten模型來表示土壤水分特征曲線和飽和非飽和土壤水力傳導(dǎo)率[31]：

式中：Ks為飽和水力傳導(dǎo)度，cm/s；Se為有效飽和度，無量綱；h為土水勢，cm；α為進氣吸力倒數(shù)，cm-1；a和b是與土壤水分特征曲線相關(guān)的經(jīng)驗常數(shù)；l為孔隙連通性參數(shù)，對大多數(shù)土壤常取0.5；θr為殘余體積含水率，cm3/cm3。

對于地表徑流運動參數(shù)，kox和koy可以由曼寧方程推導(dǎo)求出：

式中：nx、ny分別為x、y方向上的曼寧糙率系數(shù)，s/cm1/3。

1.2 地表水土壤水運動與溶質(zhì)運移全耦合模型

采用經(jīng)典的對流彌散方程描述溶質(zhì)在二維地表水與三維地下水中的運移，根據(jù)達西通量和溶質(zhì)濃度得到地表水與土壤水的溶質(zhì)交換量，在水分耦合模型的基礎(chǔ)上建立水分運動與溶質(zhì)運移耦合模型，兩者采用相同的空間和時間離散。

土壤水中的溶質(zhì)運移方程[32]：

地表水中的溶質(zhì)運移方程[32]：

地表水土壤水之間的溶質(zhì)交換量：

式中：q和qo分別是土壤水滲流通量和地表水滲流通量，cm/s；C和Co分別為土壤水和地表水中溶質(zhì)的濃度，mg/cm3；D和Do分別為土壤水和地表水的水動力彌散張量，cm2/s；Ωo、Ωs和Ωos分別為地表水溶質(zhì)的源匯項、土壤水溶質(zhì)的源匯項、地表水和土壤水之間溶質(zhì)的交換通量，mg/(L·s)；R和Ro分別為溶質(zhì)在土壤水和地表水中的阻滯系數(shù)，無量綱；?代表梯度算子，是在三維空間各方向上的全微分，?ˉ代表平面上的二維梯度算子；CL+1為上游節(jié)點處的溶質(zhì)濃度，mg/L，當(dāng)水分由地表向地下補給時，CL+1是地表水中的溶質(zhì)濃度，當(dāng)水分由地下向地表補給時，CL+1是土壤水中的溶質(zhì)濃度。

水動力彌散張量D和Do由下式計算：

式中：αl和αt分別為縱向和橫向彌散系數(shù)，cm；|q|是達西通量的大小，τ是基質(zhì)的彎曲度，無量綱；Dfree是自由溶液擴散系數(shù)，cm2/s；I為單位張量。

對于溶質(zhì)阻滯系數(shù)R和Ro，可由下式給出[33]：

式中：ρb為土壤容重，g/cm3；K′為Freundlich吸附等溫線的平衡分布系數(shù)。

2 模型驗證

選取已發(fā)表的野外坡面降雨產(chǎn)流試驗與室內(nèi)土槽坡面氮素遷移試驗對本文構(gòu)建的模型進行驗證。利用野外坡面降雨產(chǎn)流試驗對水分運動耦合模型進行檢驗；坡面氮素遷移試驗中包含硝態(tài)氮遷移過程，用來驗證本研究中的水分運動與溶質(zhì)運移耦合模型。通過比較地表徑流量和徑流溶質(zhì)濃度的觀測值和模擬值，來進一步評價模型效果。采用均方根誤差（root mean square of error，RMSE)和平均相對誤差（averaged relative error，ARE）作為判別指標(biāo)，RMSE和ARE的數(shù)學(xué)計算公式為：

式中：Swsim,i是通過實測參數(shù)并優(yōu)化參數(shù)進行建模得到的模擬值；Swobs,i是試驗觀測值；N為樣本總數(shù)；RMSE和ARE的值越小代表擬合效果越好。當(dāng)觀測值和模擬值相同時，RMSE和ARE的值為0。

2.1 水分運動耦合模型的驗證

算例1選自XING等[34]的野外邊坡人工降雨試驗，實驗設(shè)置如圖1所示，裸露的斜坡寬5 m，長10 m，試驗期間沒有作物生長。實驗前在試驗地塊上施加一定強度的降雨，直到產(chǎn)生徑流，以此對坡面進行預(yù)濕處理，保證各處理的初始含水率一致。試驗過程中分別以75 mm/h和50 mm/h的雨強在土壤容重為1.45 g/cm3的不同坡度坡面進行試驗，每次模擬降雨持續(xù)50 min，地表徑流在坡面出流口處收集。詳細實驗步驟參見文獻[34]。

圖1 實驗設(shè)置示意Fig.1 Schematic diagram of the experimental setup

根據(jù)實驗設(shè)置，本模型中將模擬區(qū)上部邊界設(shè)置為降水輸入邊界，四周為滲流面邊界，底部為透水邊界，地表出流口為臨界深度排水邊界。在模型離散中，空間步長在坡面和垂直坡面方向均為10 cm，垂向上劃分為10層。為在表層水流強交互作用處加密網(wǎng)格，設(shè)置上面5層網(wǎng)格尺寸為2 cm，下面5層網(wǎng)格尺寸為4 cm，模擬區(qū)域共生成56 661個節(jié)點。迭代計算時，設(shè)定初始時間步長0.001 s，最小時間步長乘數(shù)為0.5，最大時間步長乘數(shù)為5.0，最大時間步長為30 s，飽和度Sw迭代精度為0.05。模型計算所采用的主要參數(shù)可分為實測和經(jīng)驗參數(shù)兩部分。實測參數(shù)直接從參考文獻[34]中選取，參考文獻中未提供的參數(shù)主要有van Genuchten模型中的參數(shù)，依據(jù)試驗區(qū)土壤質(zhì)地選取經(jīng)驗值并經(jīng)過率定得到。模型主要輸入?yún)?shù)如表1所示。

表1 模型主要參數(shù)Tab.1 Main parameters of the model

選取研究中的3種坡度（5°、10°和15°）和兩種降雨強度（50 mm/h和75 mm/h）共6種試驗情景，以其中50 mm/h降雨強度下3個不同坡度坡面的試驗數(shù)據(jù)對本模型參數(shù)進行率定，75 mm/h雨強的3組試驗數(shù)據(jù)用于驗證。將整個模擬期內(nèi)坡面出口處徑流單寬流量的模擬值與實測值進行對比，兩種降雨強度下不同坡度的產(chǎn)流過程分別如圖2和圖3所示。從圖2和圖3中可以看出，在幾種實驗情景下產(chǎn)流后的單位徑流流量均急劇增加，后隨時間推移逐漸趨于穩(wěn)定。通常情況下，降雨初期土壤入滲能力受土壤水勢影響較大。隨著降雨入滲使土壤水分含量增大，土壤水勢減小，使得土體滲透率快速減小，而后徑流量隨著滲透率的逐漸穩(wěn)定也趨于平穩(wěn)。

圖2 雨強為50 mm下不同坡度坡面的單寬流量Fig.2 Unit width flux of different slopes under 50 mm rainfall intensity

圖3 雨強為75 mm下不同坡度坡面的單寬流量Fig.3 Unit width flux of different slopes under 75 mm rainfall intensity

地表徑流過程同時受降雨強度和坡度的影響。由表2可以看出，模擬后期出流趨于穩(wěn)定時，在75 mm/h的降雨強度下，坡度為5°、10°、15°的坡面出口處單寬流量分別為1.596 cm2/s、1.697 cm2/s、1.741 cm2/s。這是由于在徑流過程中，當(dāng)斜坡坡度增大時，徑流沿坡面方向的重力分力隨之增加，徑流速度的提高縮短了入滲時間，降水滲入到土體的部分減少，從而產(chǎn)流量不斷增大，因此同等降雨強度下，徑流流量與坡面坡度呈現(xiàn)正相關(guān)關(guān)系。在50 mm/h的降雨強度下也可以觀察到同樣現(xiàn)象。數(shù)值模擬結(jié)果也與此前研究一致。在坡度為5°的坡面上，50 mm/h和75 mm/h的降雨強度下的坡面出口處單寬流量分別為1.001 cm2/s、1.596 cm2/s。即坡度一定的情況下，當(dāng)雨強增大時，地表徑流量隨之增加，同時徑流過程能夠更快地趨于平穩(wěn)。在對10°和15°坡面的模擬中也得到同樣結(jié)果。這主要是由于降雨過程中土壤孔隙逐漸被水分填滿，土體下滲強度在減小后趨于穩(wěn)定，使得多余的水量積于地表形成徑流。一般來說，隨著降雨強度的增加，滲透率更容易達到穩(wěn)定狀態(tài)。

表2 不同雨強和坡度下的徑流開始時間、穩(wěn)定單寬流量及單寬流量實測值與模擬值誤差分析Tab.2 Runoff start time，stable unit width flux and error analysis of measured and simulated values of unit width flux under different rain intensities and slopes

坡面產(chǎn)流速度也與降雨強度和坡度相關(guān)。由表2可知，在50 mm/h的降雨強度下，坡度為5°、10°、15°的坡面徑流開始時間分別為97.4、66.8、64.3 s，說明降雨強度一定時，坡面產(chǎn)流速度隨坡面坡度的增加而趨于加快。在75 mm/h的降雨條件下亦有此規(guī)律。這是由于降雨條件相同時，隨著坡面坡度的增大，降雨的入滲損失減小，同時重力沿坡向的分力使水分在坡向的運動速度更大，土壤表面的水流匯集速度更快，從而使徑流形成的時間更短。此外，對比相同坡度坡面在不同降雨強度下的產(chǎn)流時間，可以發(fā)現(xiàn)降雨強度越大時，坡面產(chǎn)流時間越短。一般情況下，隨著降雨強度增大，土壤孔隙能夠更快被填滿，也提高了水分在地表的積累速度，從而使得產(chǎn)流時間更早。

不同降雨強度和坡度下模型模擬計算的評價結(jié)果列于表2。坡面單寬流量的RMSE值在0.058 cm2/s到0.147 cm2/s之間，ARE均在14.8%以內(nèi)，且水分運動耦合模型的模擬結(jié)果與該算例實驗數(shù)據(jù)擬合情況良好（圖2和圖3），驗證了本文方法和程序的正確性及適用性。

2.2 水分運動與溶質(zhì)運移全耦合模型的驗證

算例2來自朱焱等[35]2017年的室內(nèi)非飽和坡面降雨徑流和氮素遷移土槽試驗，試驗包括2個長1.5 m，高0.6 m寬0.3 m的鋼制土槽，坡角分別為8.85°和7.55°，編號分別為1號和2號。實驗設(shè)置示意見圖4。兩個土槽試驗的進行時間分別為3 260和3 200 s。試驗期間，兩土槽降雨強度分別為0.00 181、0.00 145 cm/s。詳細試驗步驟參見文獻[35]。采用本文構(gòu)建的水分運動與溶質(zhì)運移耦合模型對該試驗的地表徑流過程和徑流中硝態(tài)氮濃度進行模擬。

圖4 土槽試驗設(shè)置示意圖（1號）Fig.4 Schematic diagram of the soil tank experimental setup（No.1）

根據(jù)試驗建立三維模型，初始狀態(tài)下土壤表層含有一定濃度的硝態(tài)氮，試驗過程中上界面有指定強度和溶質(zhì)濃度的降雨輸入，因此將模擬區(qū)上部邊界設(shè)置為降水邊界和指定質(zhì)量通量邊界；鋼制土槽四周不透水，因此側(cè)向均為隔水邊界；底部由鐵砂網(wǎng)和石英砂，將其概化為透水邊界；地表出流邊緣設(shè)置為臨界深度排水邊界，使邊界處的深度等于臨界深度。采用等間距矩形網(wǎng)格對模擬區(qū)域進行空間離散，設(shè)置坡面方向網(wǎng)格尺寸為5 cm，垂直坡面方向尺寸為10 cm，垂向上剖分成20層。因土壤表層中含有一定濃度溶質(zhì)，為在水分與溶質(zhì)強交互作用處加密單元格，設(shè)置上面10層網(wǎng)格尺寸為0.5 cm，下面10層網(wǎng)格尺寸為3.5 cm，模擬區(qū)域共生成7 497個節(jié)點。兩土槽模擬時長分別為3 260和3 200 s。在進行迭代計算時，設(shè)定初始時間步長0.001 s，最大時間步長乘數(shù)為5.0，最小時間步長乘數(shù)為0.2，最大時間步長為20 s，飽和度Sw迭代精度為0.01，溶質(zhì)濃度求解迭代精度為1.0×10-10。溶質(zhì)運移相關(guān)參數(shù)直接從文獻中選取，因此利用試驗的徑流數(shù)據(jù)對模型水力參數(shù)進行率定后，即可根據(jù)徑流中溶質(zhì)濃度數(shù)據(jù)對水分與溶質(zhì)耦合運動模型做出驗證。試驗條件與模型主要輸入?yún)?shù)如表3所示。

表3 土槽試驗?zāi)Ｐ洼斎雲(yún)?shù)Tab.3 Input parameters of simulated soil tank test model

兩土槽坡面出口位置的地表徑流單寬流量、徑流中硝態(tài)氮濃度的試驗數(shù)據(jù)與模擬計算結(jié)果的擬合情況如圖5所示。整體上看，針對坡面徑流出口處單寬流量的模擬，在初始產(chǎn)流階段與實測情況有所偏差。如圖5（a）所示，1號土槽在初始產(chǎn)流階段流量實測值波動較大，而數(shù)值模擬結(jié)果變化則較為平緩，如圖5（c）所示，2號土槽在模擬前期的實測流量與模型模擬結(jié)果偏差較大，多點實測數(shù)據(jù)均大于模型模擬值。造成該現(xiàn)象的原因可能是降雨初始階段非飽和坡面的入滲情況較為復(fù)雜，試驗中的土槽坡面也并非完全均質(zhì)平整，致使該階段模擬結(jié)果產(chǎn)生偏差。隨著試驗的進行，坡面徑流過程逐漸趨于穩(wěn)定，后期實測流量與模型模擬結(jié)果的擬合情況相對較好。1號土槽單寬徑流量計算的RMSE值為0.583 cm2/s，ARE值為11.8%；2號土槽單寬徑流量計算的RMSE值為0.833 cm2/s，ARE值為21.5%。造成2號土槽徑流計算誤差較大的原因主要是由于該土槽的實測數(shù)據(jù)點分布較為離散，波動較大，而利用模型計算的徑流與徑流溶質(zhì)濃度數(shù)值模擬結(jié)果均為平滑曲線，相對平緩。整體而言，本研究模型的徑流模擬計算結(jié)果良好。

圖5 兩土槽坡面徑流量、徑流中硝態(tài)氮濃度的實測與模擬對比Fig.5 Comparison between the experimentally measured and model simulated values of runoff and nitrate nitrogen concentration on the slopes of the two soil trenches

兩土槽徑流中硝態(tài)氮濃度的變化過程如圖5（b）和圖5（d）所示，由于降雨與土壤中的硝態(tài)氮存在較大濃度差，表層土壤中的硝態(tài)氮在對流擴散作用的驅(qū)動下隨徑流不斷流失，使得徑流中硝態(tài)氮濃度呈現(xiàn)先急劇下降后趨于穩(wěn)定的變化。模擬結(jié)果與試驗值的偏差在模擬初期較為明顯，多數(shù)模擬值均小于試驗結(jié)果，可能是受徑流計算偏差的影響。當(dāng)1號土槽穩(wěn)定出流且徑流中硝態(tài)氮濃度變化相對平穩(wěn)時，在1 285 s處有一實測的徑流硝態(tài)氮濃度數(shù)據(jù)與整體偏離較多，為更準(zhǔn)確評價模型的模擬效果，選擇在誤差計算時除去該點數(shù)據(jù)。1號土槽徑流中硝態(tài)氮濃度計算的RMSE值為1.334 mg/L，ARE值為24.7%，2號土槽徑流中硝態(tài)氮計算的RMSE值為0.248 mg/L，ARE值為12.1%。由圖4-14-2可知，1號土槽在初始產(chǎn)流階段流量實測值波動較大，且徑流硝態(tài)氮濃度模擬值與實測值的誤差主要是在模擬初期產(chǎn)生，故1號土槽硝態(tài)氮濃度計算誤差較大的原因可能是受模擬初期徑流計算偏差的影響?？傮w而言，模型模擬計算結(jié)果與試驗數(shù)據(jù)擬合較好。綜合以上分析，本文構(gòu)建的水分運動與溶質(zhì)運移耦合模型能夠準(zhǔn)確模擬地表徑流過程及徑流中硝態(tài)氮的濃度變化過程。

2.3 參數(shù)敏感性分析

在模型的構(gòu)建過程中，通過引入耦合長度lso這一參數(shù)來實現(xiàn)水分和溶質(zhì)在地表與土壤兩者間運動遷移的耦合，因此選擇分析該參數(shù)的敏感性。使用的敏感性分析方法是：通過調(diào)整模型中參數(shù)取值，比較模擬變量計算結(jié)果的變化情況，分析過程中采用控制變量法。具體來說，本文在最佳耦合長度lso的基礎(chǔ)上分別上調(diào)（+）和下調(diào)（-）10%后進行模型計算，觀察此參數(shù)對地表流量和徑流中溶質(zhì)濃度的影響。

算例1中降雨強度設(shè)置為50和75 mm/h，坡度設(shè)置為5°、10°和15°。算例1中情景1-6的降雨強度和坡度分別為50 mm/h和5°、50 mm/h和10°、50 mm/h和15°、75 mm/h和5°、75 mm/h和10°、75 mm/h和15°，計算結(jié)果如表4所示。由表4可知，不同的耦合長度設(shè)置對地表徑流和徑流氮素濃度模擬的影響不同。隨著耦合長度的增加，模擬的地表單寬流量和徑流硝態(tài)氮濃度均有不同程度的改變，但總體上看，兩者的變化均與耦合長度參數(shù)呈正相關(guān)關(guān)系。從表4中還可以看出，在對算例1的6種實驗情景的模擬中，耦合長度lso變化10%對地表出流量的改變量均在0.5%以內(nèi)；在對算例2兩個土槽試驗的模擬中，耦合長度lso變化10%對地表出流量和徑流硝態(tài)氮濃度模擬結(jié)果的改變均小于5%，因此可認為模型中該參數(shù)不敏感。

表4 參數(shù)敏感性分析Tab.4 Sensitivity analysis of parameters

3 結(jié) 論

本文采用二維擴散波方程和傳統(tǒng)的三維Richards方程描述地面徑流和土壤水分運動，選用雙層結(jié)點法對兩者進行耦合，根據(jù)達西定律和地表與土壤水中的溶質(zhì)濃度，計算地面徑流和土壤水分的溶質(zhì)交換量，從而構(gòu)建地表水土壤水運動與溶質(zhì)運移全耦合數(shù)值模型。野外邊坡降雨產(chǎn)流試驗觀測數(shù)據(jù)與模型模擬結(jié)果的對比顯示，本模型水分計算的平均相對誤差小于14.8%，均方根誤差小于0.147 cm2/s。室內(nèi)坡面徑流與氮素遷移試驗監(jiān)測數(shù)據(jù)與模型模擬結(jié)果的對比顯示，本模型水分計算的平均相對誤差小于21.5%，均方根誤差小于0.833 cm2/s，硝態(tài)氮濃度計算的平均相對誤差小于24.7%，均方根誤差小于1.332 mg/L。說明本文構(gòu)建的全耦合模型對地表水土壤水水分運動與溶質(zhì)運移過程的模擬準(zhǔn)確可靠，模型對地表和地下不同模塊間的耦合關(guān)系表述合理。