郝玉潔 ,欒永霞,鄭西來, 3 * ,鄭天元,能 惠

(1.中國海洋大學(xué)環(huán)境科學(xué)與工程學(xué)院,山東 青島 266100;2.青島市城陽區(qū)農(nóng)業(yè)農(nóng)村局,山東 青島 266109;3.中國海洋大學(xué)海洋環(huán)境與生態(tài)教育部重點實驗室,山東 青島 266100)

引 言

野外監(jiān)測對于探究DON污染現(xiàn)狀具有重要意義,但在大尺度下針對土壤特性進行逐一實驗探究耗時費力,模型是進行相關(guān)研究的有力方式。Liang等[14]利用WHCNS_veg模擬了不同的農(nóng)業(yè)管理方式下溫室蔬菜系統(tǒng)中DON和硝態(tài)氮的遷移轉(zhuǎn)化過程。Tian等[15]利用DRAINMOD-N II實現(xiàn)了景觀生態(tài)系統(tǒng)中DON的運移過程的模擬。但這些模擬研究主要關(guān)注以點為尺度的均一土壤,未充分考慮土壤空間變異與DON淋失風(fēng)險的關(guān)系,不利于做出合理的管理決策,因此有必要將土壤空間變異納入DON淋失狀況評估中。

靈敏性分析可以量化輸入變量對模型輸出的影響程度。其中,局部敏感性分析方法通過評估單個輸入因子對模型的影響來進行參數(shù)的重要性排名。全局敏感性方法認(rèn)為模型輸出結(jié)果的改變是輸入因子的改變以及不同因子之間相互交互作用的結(jié)果,擴展傅立葉幅度敏感性檢驗法(Extended Fourier Amplitude Sensitivity Test,EFAST) 作為一種基于方差分析的全局敏感詞分析方法,能夠量化每個參數(shù)對模型輸出的無條件方差的貢獻,已廣泛應(yīng)用于不同研究[16-17]。綜合運用不同敏感性分析方法可以幫助我們識別影響DON淋失風(fēng)險的關(guān)鍵性土壤特性。

本研究擬建立土壤DON遷移轉(zhuǎn)化模型,探究水分滲漏和DON淋失對土壤初始狀態(tài)和土壤水力學(xué)參數(shù)的空間變異的響應(yīng)特征。同時,定量分析土壤空間變異對大尺度水分滲漏和DON淋失的影響,為評估DON淋失風(fēng)險提供依據(jù)。

1 土壤水分-溶質(zhì)運移模型構(gòu)建

在模型中構(gòu)建垂向高度為40 cm的均質(zhì)土壤剖面。設(shè)置灌水定額為50 mm,灌溉強度為10 mm h-1,灌溉時間為5 h,假定水分在土壤表面均勻輸入,故可將模擬情況簡化為水分和溶解性有機氮(DON)在土壤剖面上的一維垂向運動。選擇糞肥作為施加氮肥,施加量為 200 kg N ha-1,灌溉施肥過程設(shè)置為1/4W-1/2N-1/4W的模式[18]?？紤]到灌溉施肥后水分及DON再分布過程中均有可能發(fā)生水分滲漏和DON的淋失,模擬時長設(shè)置為自灌溉開始96 h。

1.1 土壤水分運動模型

COMSOL中運用Richards方程來描述一維非飽和帶水分運動,表達(dá)式為:

(1)

式中,h為土壤基質(zhì)勢,cm;z為垂向坐標(biāo),向上為正,cm;K(h)為非飽和導(dǎo)水率,cm d-1;運用van Genuchten模型描述水頭特征曲線。

初始條件可以表示為:

h0=H-z,-40≤z≤0 cm,t=0

(2)

式中,h0為土壤初始總水勢,cm;H為基質(zhì)勢,cm。

模擬區(qū)域上邊界設(shè)定為大氣邊界,大氣邊界的水流量取決于土壤入滲量,下邊界設(shè)定為自由排水邊界,邊界條件可以表示為:

(3)

(4)

式中,E(t)表示土壤表面入滲量,cm·d-1。

1.2 土壤DON遷移轉(zhuǎn)化模型

COMSOL中運用對流彌散方程描述非飽和帶溶質(zhì)的遷移轉(zhuǎn)化過程。DON隨水分向下遷移的過程中經(jīng)歷的反應(yīng)包括土壤釋放DON以及DON的礦化和吸附。因此,DON在土壤中的運移過程可以表示為:

(5)

式中,c為土壤溶液中DON的濃度,mg N cm-3;s為吸附在土壤顆粒上的DON的含量,mg N kg1;ρ為土壤密度,g cm-3;t為時間,d;q為土壤水通量,cm d-1;D為水動力彌散系數(shù),cm2d-1;μw,1和μs,1指液相中和固相中DON的礦化速率常數(shù),d-1;μw,2和μs,2指液相中和固相中DON從土壤中的釋放速率常數(shù),d-1。

DON從土壤中的釋放這一過程的模擬借鑒了DAISY模型中有機質(zhì)周轉(zhuǎn)過程的設(shè)置[19],DAISY模型指出,DON的形成是土壤有機質(zhì)和微生物共同作用的結(jié)果,DON的釋放速率可以表示為:

μw,2=μs,2=fDON(DSMB1+DSMB2+kSOM2)

(6)

式中,SMB指土壤微生物庫;SOM指土壤有機質(zhì)庫,每個主庫可以劃分為兩個部分:快速分解的庫(SMB1和SMB2)和緩慢分解的庫(SOM1和SOM2)。DOM指溶解性有機質(zhì)庫;fDON指從DOM到DON的分配系數(shù);D指微生物死亡速率常數(shù),d-1;k指分解速率常數(shù),d-1,需要根據(jù)實際的黏粒含量,溫度和體積含水量來進行修正,校準(zhǔn)系數(shù)分別為0.7、0.9和1。

DON的吸附可以用Freundlich吸附模型來描述[20-21]:

s=KFcN

(7)

式中,KF指Freundlich模型中吸附常數(shù),mg1-NL-Nkg-1;N為吸附線性指數(shù),L mg-1。

假設(shè)土壤DON的初始濃度呈線性分布(以大沽河流域2019年部分區(qū)域氮素含量為參考[22],設(shè)置土壤剖面底層DON累積量較多,此時剖面底部DON淋失量較大),設(shè)置變化速率為0.000 5,土壤剖面上不同深度處土壤溶液中DON的初始濃度可以表示為:

c(z,0)=c0-0.0005×(40+z), -40≤z≤0 cm,t=0

(8)

式中,c0為土壤溶液中DON的初始濃度(z = -40 cm),mg N cm-3;

模擬區(qū)域上邊界允許溶質(zhì)輸入,輸入量取決于施肥強度,下邊界為自由排泄邊界,邊界條件可以表示為:

(9)

(10)

式中,ca為肥料溶液中DON的濃度,mg N cm-3。

1.3 模型參數(shù)

1.3.1 土壤水力學(xué)參數(shù)

土壤水分特征曲線的模擬運用了van Genuchten模型,這個模型中包含5個獨立的水力學(xué)參數(shù):θr,土壤殘余體積含水率;θs,土壤飽和體積含水率;α和n為經(jīng)驗擬合參數(shù);Ks為飽和導(dǎo)水率,L T-1,利用Rosetta人工神經(jīng)網(wǎng)絡(luò)可以求得不同顆粒組成和容重土壤的水力學(xué)參數(shù)。為了探究土壤空間變異對土壤水分滲漏和DON淋失情況的影響,設(shè)定變異系數(shù)CV<0.1為弱變異、0.1

表1 不同變異程度土壤物理性質(zhì)統(tǒng)計特征值

表2 不同變異程度土壤水力學(xué)參數(shù)統(tǒng)計特征值

1.3.2 溶質(zhì)運移參數(shù)

縱向彌散系數(shù)取為DL=10 cm,DON在自由水中的分子擴散系數(shù)取為D=1.2 cm2d-1[14]。DON的礦化速率為μw,1=μs,1=0.02 d-1[14]。DON從土壤中的釋放過程涉及到的相關(guān)參數(shù)設(shè)置如下:DSMB1=1.85×10-4d-1,DSMB2=0.02 d-1,kSOM2=1.40×10-4d-1[19]。以大沽河流域土壤為研究對象,該區(qū)域非飽和帶土壤溶液碳氮比為26.3:2.79[22],假定在整個模擬過程中土壤溶液的碳氮比是恒定不變的,根據(jù)從DOM到DOC的分配系數(shù)(fDOC=0.05)和碳氮比計算可以得到從DOM到DON的分配系數(shù)值[19]。Freundlich吸附模型參數(shù)取值為KF= 2.619,N=1.788[21]。

1.4 敏感性分析方案

1.4.1 局部敏感性分析

探究不同初始濃度(c0取值為0.08、0.16、0.24、0.32、0.40 mg N cm-3)和初始基質(zhì)勢(H取值為-1、-2、-3、-4、-5 m)下土壤累積水分滲漏量和累積DON淋失量的影響,其它參數(shù)取均值,將待分析參數(shù)值輸入COMSOL軟件進行模擬。

探究不同土壤水力學(xué)參數(shù)對累積水分滲漏量和累積DON 淋失量的影響,其它變量取固定值(c0取值為0.16 mg N cm-3,H取值為-3 m),待分析參數(shù)波動值為±10%、±20%,將待分析參數(shù)值輸入COMSOL軟件進行模擬,與均值參數(shù)組模擬結(jié)果進行對照。

1.4.2 全局敏感性分析

擴展傅立葉幅度敏感性檢驗法(EFAST)結(jié)合了Sobol’法和傅立葉幅度敏感性檢驗法(FAST)的優(yōu)點,是一種基于方差分解的全局敏感性分析方法。它的基本思想如下:

模型的總方差可以分解為:

(11)

式中,V表示模型總方差;Vi表示xi的自身波動引發(fā)的方差;Vij表示xi通過xj作用貢獻的方差;Vijm表示xi通過xj、xm作用貢獻的方差;V12…k表示xi通過其余k-1個變量作用貢獻的方差。

變量xi的一階敏感性指數(shù)表示xi對模型輸出總方差的直接貢獻率,可以表示為:

(12)

變量xi的全局敏感性指數(shù)表示xi對模型輸出總方差的直接貢獻率和變量之間的相互作用對模型輸出總方差的貢獻率之和,可以表示為:

(13)

式中,-i表示除變量xi之外的其它參數(shù),V-i表示不包括xi的其它變量的波動引發(fā)的方差之和。

利用蒙特卡洛法對5個土壤水力學(xué)參數(shù)進行隨機采樣,每種變異程度采樣445次(EFAST法認(rèn)為有效采樣次數(shù)需要大于參數(shù)個數(shù)×65),將待分析參數(shù)值輸入COMSOL軟件進行求解,利用MATLAB調(diào)動COMSOL軟件進行模擬大批量的模擬,利用EFAST方法對模擬結(jié)果進行分析。

2 結(jié)果與分析

2.1 初始狀態(tài)的影響

設(shè)置不同初始濃度(c0取值為0.08～0.40 mg N cm-3)和初始基質(zhì)勢(H取值為-1～-5 m)值,對水力學(xué)參數(shù)取均值,得到不同條件下土壤剖面上DON分布、含水率分布(T=96 h)情況以及底部累積水分滲漏量、累積DON淋失量(圖1)。由于c0僅與土壤剖面上DON濃度分布和底部淋失通量有關(guān),所以沒有分析c0與土壤剖面上含水率分布和底部水分滲漏量的關(guān)系。

圖1 不同初始DON濃度(c0)條件下土壤剖面上DON濃度分布和底部累積DON淋失量

如圖1所示,當(dāng)T=96 h時,DON濃度在整個剖面上呈下降趨勢,大部分DON累積在土壤表層,0～10 cm土壤層中DON含量占整個剖面總量的32.70%～50.51%。隨著c0的增加,土壤剖面各個深度的DON濃度都明顯增大,累積DON淋失量也顯著增加,說明土壤中存儲的DON是DON淋失的重要來源。由圖2可知,H對土壤剖面上體積含水率以及累積土壤水分滲漏影響較小。當(dāng)T=96 h時,土壤體積含水率隨深度的下降而降低。隨著H的增加,土壤剖面上各個深度含水率小幅增加,土壤累積水分滲漏量隨之增加。由圖3可知,對于不同的H值,土壤剖面上各個深度DON含量沒有明顯差異,但是累積DON淋失通量隨H的增加而增加。這是因為土壤基質(zhì)勢較高時,土壤易發(fā)生水分滲漏和溶質(zhì)淋失[24],DON隨水分大量淋失,土壤含水量降低,但是土壤中DON總量也降低,所以土壤剖面上DON濃度大小沒有明顯變化。

圖2 不同土壤初始基質(zhì)勢條件下土壤剖面上含水率分布和底部累積水分滲漏量

圖3 不同土壤初始基質(zhì)勢條件下土壤剖面上DON濃度分布和底部累積DON淋失量

2.2 土壤水力學(xué)參數(shù)的影響

2.2.1 局部敏感性分析

對五個土壤水力學(xué)參數(shù)(θr,θs,α,n,Ks)進行局部敏感性分析,參數(shù)波動值為±10%、±20%,將其它變量取固定值(c0取值為0.16 mg N cm-3,H取值為-3 m),可以得到不同條件下土壤累積水分滲漏量和累積DON淋失量,與均值參數(shù)組模擬結(jié)果進行對比,可以得到單個土壤水力學(xué)參數(shù)值發(fā)生波動時土壤累積水分滲漏量和累積DON淋失量的變化幅度(圖4)。為了規(guī)避不合理的模擬結(jié)果,忽略當(dāng)n<1.1的情形[25]。

圖4 土壤水力學(xué)參數(shù)波動與累積水分滲漏量和累積DON淋失量變幅的關(guān)系

由圖4(a)可知,五個水力學(xué)參數(shù)對土壤累積水分滲漏量影響程度從大到小依次為:α>n>Ks>θs>θr,其中n、α、Ks的影響較為明顯,θr與水分滲漏過程量大小無明顯關(guān)系。此外,累積水分滲漏量與n、Ks呈正相關(guān),與α、θs呈負(fù)相關(guān)。隨著α的增加,累積水分滲漏量顯著下降,原因是α與土壤進氣吸力是倒數(shù)關(guān)系,α的增加說明土壤進氣吸力變小,土壤導(dǎo)水能力減弱[26],進而降低了土壤水分向深層的滲漏。隨著n的增加,累積水分滲漏量下降,這是因為n值與土壤介質(zhì)中孔隙分布有關(guān),隨著n值的增加,土壤孔隙中細(xì)小孔隙數(shù)量變少,土壤粘滯力降低,土壤導(dǎo)水能力增強[27]。隨著Ks的增加,累積水分滲漏量呈線性增大的趨勢。由圖4(b)可知,五個水力學(xué)參數(shù)對累積DON淋失量影響程度從大到小依次為:α>n>Ks>θs>θr,與對累積土壤水分滲漏量的影響程度排名相同,且參數(shù)波動造成的累積土壤水分滲漏量和累積DON淋失量的變化幅度也基本一致,說明大量水分的滲漏可能伴隨著大量DON的淋失,累積DON淋失量對土壤水力學(xué)參數(shù)中n、α、Ks的波動更為敏感。

2.2.2 全局敏感性分析

利用蒙特卡洛法對5個土壤水力學(xué)參數(shù)進行隨機采樣,每種變異程度產(chǎn)生445個模擬案例,將其它變量取固定值(c0取值為0.16 mg N cm-3,H取值為-3 m),可以得到不同條件下土壤累積水分滲漏量和累積DON淋失量,利用EFAST方法對模擬結(jié)果進行分析,計算得到五個土壤參數(shù)對兩者的敏感性指數(shù)(圖5、6)。如圖5,6所示,五個水力學(xué)參數(shù)對累積水分滲漏量和DON累積淋失量的影響程度排序與局部敏感性分析結(jié)果不同,是因為土壤顆粒組成的變化對五個參數(shù)值的影響程度不同,即對于某一區(qū)域,五個水力學(xué)參數(shù)的變異程度是不一致的(表1、2),對水力學(xué)參數(shù)設(shè)定相同的變幅進行局部敏感性分析與實際情況會有偏差, 所以有必要首先確定五個參數(shù)值的取值范圍和分布類型,得到符合實際分布的參數(shù)序列,進而進行敏感性分析得到的結(jié)果更為準(zhǔn)確。

由圖5(a)可知,土壤在弱變異條件下,五個水力學(xué)參數(shù)對累積土壤水分滲漏量影響程度依次為:θs>Ks>α>n>θr,由圖1可知,θs和Ks對土壤累積水分滲漏量的影響最明顯,對累積水分滲漏量的變化的方差貢獻分別達(dá)到82.20%和82.13%,這是因為θs與土壤剖面可能存儲的水分總量密切相關(guān),θs較大時,同一體積含水率相應(yīng)的負(fù)壓值升高,土壤滲透性降低,水分易滯留在土壤孔隙中[26]。Ks與水分在剖面上的運移速度密切相關(guān),Ks的空間變異程度(CV=0.354)明顯大于其它參數(shù),是造成土壤水分滲漏量不確定性的主要參數(shù)之一。由圖5(b)可知,五個水力學(xué)參數(shù)對累積DON淋失量影響程度依次為Ks>α>n>θr>θs,Ks對累積DON淋失量全局敏感性指數(shù)為0.83(即可解釋累積DON淋失量變化方差的83.00%),對DON淋失有較大影響。一階敏感性指數(shù)的差異說明Ks、α的影響主要通過與其它參數(shù)的交互實現(xiàn),而n、θr、θs的影響主要是直接作用引起的。

由圖6(a)可知,五個水力學(xué)參數(shù)對中等變異土壤累積水分滲漏量影響程度與弱變異條件下順序明顯不同,五個水力學(xué)參數(shù)對累積土壤水分滲漏量影響程度依次為:Ks>α>n>θr>θs,Ks的全局敏感性指數(shù)為0.76,對土壤累積水分滲漏量影響程度最大,α和n對累積土壤水分滲漏量影響也較為明顯。θs的全局敏感性指數(shù)為0.08,影響程度大幅度降低。由圖6(b)可知,五個水力學(xué)參數(shù)對累積DON淋失量影響程度依次為Ks>α>n>θr>θs,累積DON淋失量的變化主要由Ks和α造成。由表2可知,中等變異粘土Ks的空間變異程度(CV=0.704)明顯大于其它參數(shù),這加劇了其對土壤水分滲漏量和DON淋失量不確定性的影響。

圖6 中等變異土壤參數(shù)對土壤累積水分滲漏量和DON累積淋失量的全局和一階敏感性指數(shù)

2.3 不同變異程度下區(qū)域土壤水分滲漏和DON淋失量

利用蒙特卡洛法對5個土壤水力學(xué)參數(shù)進行隨機采樣,每種變異程度產(chǎn)生445個模擬案例,將其它變量取均值(c0取值為0.16 mg N cm-3,H取值為-3 m),可以得到不同條件下累積DON淋失量和土壤水分滲漏量,將445個案例進行排列組合,可用于分析大尺度下地區(qū)土壤水氮淋失情況。由表3可知,兩種變異程度下得到的累積水分滲漏量和累積DON淋失量變異系數(shù)都達(dá)到了中等變異程度(0.1

表3 不同變異程度累積土壤水分滲漏量和DON累積淋失量統(tǒng)計特征值

圖7、圖8給出了土壤在弱變異和中等變異條件下累積水分滲漏量和累積DON淋失量的頻數(shù)直方圖,縱坐標(biāo)表示該區(qū)間的模擬值的數(shù)量占比,對頻數(shù)直方圖進行擬合,可以看出,累積水分滲漏量和累積DON淋失量出現(xiàn)的頻率基本服從對數(shù)正態(tài)分布,擬合的得到的相關(guān)參數(shù)如表4所示。由圖7、圖8可知,在弱變異條件下,累積土壤水分滲漏量和累積DON淋失量范圍是1.009～8.647 mm和2.165～12.164 kg N ha-1,在中等變異條件下,兩者范圍是0.300～25.294 mm和1.316～34.438 kg N ha-1,隨著變異程度的升高,兩者的變化范圍都明顯增大。對于大尺度地區(qū),土壤空間變異性越大,土壤參數(shù)空間分布不均勻性增大,出現(xiàn)大規(guī)模水氮淋失的可能性越大。

圖7 弱變異累積土壤水分滲漏量和累積DON淋失量的頻率直方圖

圖8 中等變異土壤累積水分滲漏量和累積DON淋失量的頻率直方圖

表4 土壤累積水分滲漏量和累積DON淋失量頻率直方圖擬合參數(shù)值

3 結(jié)論

本研究擬建立土壤溶解性有機氮遷移轉(zhuǎn)化模型,探究了土壤初始狀態(tài)和土壤水力學(xué)參數(shù)的空間變異對土壤累積水分滲漏量和累積DON淋失量的影響,同時,定量分析土壤空間變異對大尺度水分滲漏和溶解性有機氮淋失的影響。主要結(jié)論如下:

(1) 土壤DON初始含量和土壤初始基質(zhì)勢的增加促使土壤剖面上DON含量和底部DON淋失量顯著增加,其中土壤DON初始含量的影響程度大于土壤初始基質(zhì)勢。

(2) 土壤累積水分滲漏量和累積DON淋失量與n、Ks值呈正相關(guān),與α、θs值呈負(fù)相關(guān),與θr的大小無明顯關(guān)系。

(3) 土壤在弱變異條件下,五個水力學(xué)參數(shù)對累積水分滲漏量和累積DON淋失量影響程度順序分別為:θs>Ks>α>n>θr和Ks>α>n>θr>θs;土壤在中等變異條件下,五個水力參數(shù)對累積水分滲漏量和累積DON淋失量影響程度順序分別為:Ks>α>n>θr>θs和Ks>α>n>θr>θs。

(4) 對于大尺度地區(qū),土壤在弱變異條件下,大部分區(qū)域土壤水分滲漏和DON淋失程度一致,隨著變異度增加,土壤水分和DON空間變異性增大,局部區(qū)域更可能出現(xiàn)大規(guī)模的水氮淋失。因此,在大尺度下評估土壤水分滲漏和DON淋失風(fēng)險時,應(yīng)當(dāng)關(guān)注土壤空間變異的影響。