王睿垠 魏永霞 張翼鵬 張雨鳳 馬德才

(1.東北農(nóng)業(yè)大學(xué)理學(xué)院， 哈爾濱 150030； 2.東北農(nóng)業(yè)大學(xué)水利與土木工程學(xué)院， 哈爾濱 150030；3.農(nóng)業(yè)部農(nóng)業(yè)水資源高效利用重點實驗室， 哈爾濱 150030； 4.中山大學(xué)中法核工程與技術(shù)學(xué)院， 珠海 519082)

0 引言

東北黑土區(qū)是世界三大黑土帶之一，也是國家商品糧生產(chǎn)基地，糧食產(chǎn)量對于保障國家安全具有重要意義。水分是農(nóng)業(yè)生產(chǎn)的重要因素，生物炭作為土壤改良劑對土壤水分運動的影響尚不明確。目前針對東北黑土施加生物炭的研究多為試驗研究，少有生物炭對土壤水分影響的數(shù)值研究。土壤是一個十分復(fù)雜的系統(tǒng)，是植物生長中水分和養(yǎng)分的主要來源，是由固、液、氣三相物質(zhì)組成的體系。生物炭是生物有機質(zhì)在無氧或缺氧情況下，經(jīng)高溫裂解的產(chǎn)物[1]。研究表明，在某些情況下，土壤中添加生物炭，可以修復(fù)土壤[2]并改變土壤的持水、導(dǎo)水特性[3-4]。土壤的導(dǎo)水率是土壤粒徑、黏土含量、有機質(zhì)含量以及粒子大小分布等的函數(shù)，是全面反映土壤導(dǎo)水特性的一個參量[5]。土壤具有較高的導(dǎo)水率時可以加快滲透和排水[6]，土壤快速排水有利于減少徑流的發(fā)生[7]，但滲透速度過快，也會減少水中營養(yǎng)物和農(nóng)用化學(xué)品充分溶解過濾的機會[8]。黑土是中國重要的土壤資源，由于黑土表土疏松、底土黏重以及人為不合理的耕作，使得東北黑土持水能力降低，農(nóng)田生產(chǎn)力下降。生物炭作為土壤改良劑，在改善土壤的持水和導(dǎo)水能力、提高作物的生產(chǎn)力方面與其他有機添加物相比，更加安全和持久[9]。研究證實了土壤添加生物炭可以影響作物的產(chǎn)量[10-12]，其物理機制多歸因于土壤孔隙度的增加、土壤容重的降低[13]、土壤保水導(dǎo)水能力的提高以及土壤可利用養(yǎng)分的增加[14]。研究表明，土壤結(jié)構(gòu)在時間上還會受到濕潤與干燥循環(huán)、生物活性以及農(nóng)業(yè)作業(yè)的影響[15]。

雖然添加生物炭對土壤水力學(xué)特性產(chǎn)生不同程度的影響，但有關(guān)預(yù)測生物炭定量改善土壤水分特性、生物炭決定水分分布的確切解的研究[16]鮮見報道。有些研究給出了土壤水動力的空間變化，主要是土壤粒子分布的不均勻性[17]。土壤中的水分運動，取決于固體顆粒的尺寸、形狀和分布[18]。土壤的多樣性以及水分運動的復(fù)雜性，使得試驗的數(shù)量和耗時巨大。本文推導(dǎo)添加生物炭土壤的毛細(xì)管模型[19]，并在此模型下利用有限元法或差分法對土壤進(jìn)行數(shù)值模擬。對比此模型與添加生物炭土壤土柱滲透的數(shù)值分析結(jié)果和試驗結(jié)果，以證實此模型預(yù)測不同條件下的土壤水分分布或者校驗試驗結(jié)果的正確性。

1 材料與方法

1.1 研究區(qū)概況

取土地點位于黑龍江省農(nóng)墾北安管理局紅星農(nóng)場實驗基地(48°10′E，127°1′N)。該地區(qū)為典型的東北黑土區(qū)，屬寒溫帶大陸性季風(fēng)氣候，冬季寒冷干燥，夏季高溫多雨，降雨時間集中，年降水量多介于450～650 mm，其中7—9月累積降水量約402 mm，約占年平均降水量的75%，因此旱澇災(zāi)害頻繁，水土流失嚴(yán)重。該區(qū)域的土壤多為壤土和粘壤土，以粗粉沙和粘粒最多，各占30%～40%，容重1.0～1.3 g/cm3，透水、持水、通氣性均較好，總孔隙度多在40%～60%，毛管孔隙可占20%～30%，通氣孔隙度約占20%。

1.2 試驗材料與設(shè)備

供試材料為玉米秸稈生物炭，購自遼寧金和福農(nóng)業(yè)開發(fā)有限公司，制備方式為450℃無氧條件下高溫裂解?；纠砘再|(zhì)：粒徑范圍1.5～2.5 mm，其中2 mm粒徑占60%以上。各指標(biāo)如下：全碳質(zhì)量分?jǐn)?shù)72.21%，全氮質(zhì)量分?jǐn)?shù)1.56%，全磷質(zhì)量分?jǐn)?shù)0.72%，全鉀質(zhì)量分?jǐn)?shù)1.64%，pH值為9.35，經(jīng)測量生物炭密度約596 kg/m3。

脫水面半徑為9.8 cm的CR- 21G3型高速離心機，用于室內(nèi)測定土壤水分特征曲線；自制直徑8 cm、長為80 cm的有機玻璃筒，距供水端10、20、30、40、50、60、70 cm處設(shè)有取土孔，用于測量不同位置的土壤含水率。

1.3 試驗方法與內(nèi)容

土樣為田間試驗土樣。取土小區(qū)面積為20 m×5 m，坡度3°。設(shè)置5種生物炭用量添加比例，即不添加生物炭處理(CK)，以及生物炭添加量為25、50、75、100 t/hm2的處理。試驗小區(qū)采用相同的水肥管理，每個處理重復(fù)3次。2017年5月中旬將生物炭人工均勻混入耕層黑土(0～20 cm)，生物炭占土壤的體積比約為2%(C1)、4%(C2)、6%(C3)、8%(C4)，取樣時間為2017年9月25日。用環(huán)刀采樣器取對照土樣(CK)，取土深度為10～15 cm，環(huán)刀高5 cm、體積100 cm3，取土后將環(huán)刀用蒸餾水吸水48 h，至飽和后稱量。生物炭放置環(huán)刀中，用蒸餾水吸水并壓實48 h后稱量。用離心機法分別測定土壤和生物炭的水分特征曲線。離心機轉(zhuǎn)速設(shè)定為500～6 000 r/min，離心時間為100 min[20]，離心結(jié)束后稱量，土壤與生物炭各重復(fù)3次。

室內(nèi)試驗用定水頭法分別測定生物炭和5種處理的飽和導(dǎo)水率。利用馬氏瓶供水，維持水頭恒定為5 cm，每隔10 min記錄出流水量，出流水量保持穩(wěn)定時結(jié)束試驗，計算飽和導(dǎo)水率，各處理重復(fù)3次。

測定土柱的含水率。在5個處理小區(qū)取原狀土，取土深度10～15 cm，土樣經(jīng)風(fēng)干、破碎、篩分后，實地測量土壤干容重為1.1 g/cm3。將土樣按此容重每層3 cm分層回填入土柱，裝填好的土柱直徑8 cm、高80 cm。進(jìn)水端采用馬氏瓶供水，視進(jìn)水端的土壤為飽和，當(dāng)濕潤鋒到達(dá)土柱約3/4處，在各取土孔處迅速取土，測定各處的含水率。各處理重復(fù)3次。

1.4 數(shù)據(jù)處理方法

采用Excel 2010進(jìn)行數(shù)據(jù)整理，所用數(shù)據(jù)均為各重復(fù)測定的平均值，土壤水分特征曲線采用Retc軟件進(jìn)行擬合，采用python語言建模編程，采用通用的有限元分析軟件Abaqus 6.13進(jìn)行數(shù)值模擬計算。

2 理論分析

2.1 土壤水分特征曲線

目前學(xué)者大多采用VAN- GENUCHTEN(VG)模型給定土壤水分特征曲線(SWCC)。模型由VAN GENUCHTEN[21]在1980年提出，適用于描述低含水率時的土壤水分特征曲線。其表達(dá)式為

(1)

其中

m=1-1/n(n>1)

式中θr——土壤殘余含水率，cm3/cm3

θs——土壤飽和含水率，cm3/cm3

h——土壤吸力，cm

a、n、m——曲線擬合參數(shù)

θ——土壤含水率

本研究采用VAN- GENUCHTEN模型，利用試驗測得土壤含水率，擬合土壤飽和含水率和殘余含水率。

2.2 添加生物炭土壤的飽和導(dǎo)水率

2.2.1假設(shè)生物炭在土壤中呈層狀分布的導(dǎo)水率

穩(wěn)定水流分別穿過土壤和生物炭時，設(shè)K1為土壤飽和導(dǎo)水率，K2為生物炭飽和導(dǎo)水率。當(dāng)生物炭在土壤中呈層狀分布時，設(shè)生物炭與土壤的混合土層的總厚度為L，令混合土層中土壤體積比為α、生物炭體積比為β，則α+β=1，因此土壤層的厚度為αL，生物炭層厚度為βL，結(jié)構(gòu)如圖1a所示。設(shè)土壤層與生物炭層水力梯度分別為i1、i2，整個混合土層有效導(dǎo)水率為Ke，整個混合土層水頭損失應(yīng)等于兩層水頭損失之和，即

iL=αi1L+βi2L

(2)

式中i——總的水力梯度

JW=Kei=K1i1=K2i2

(3)

式中JW——水流通量

結(jié)合式(2)可以得到

(4)

圖1 生物炭分布示意圖Fig.1 Sketches of media in different directions

一般情況下，土壤體積比α遠(yuǎn)大于生物炭體積比β，另外，生物炭的飽和導(dǎo)水率K2也遠(yuǎn)大于草甸黑土的飽和導(dǎo)水率K1，這樣式(4)可以進(jìn)一步簡化為

(5)

當(dāng)生物炭在土壤中層狀分布時，混合土壤的飽和導(dǎo)水率可以近似地認(rèn)為是土粒的飽和導(dǎo)水率除以土粒在土壤中的比例。

2.2.2假設(shè)生物炭在土壤中呈條狀分布的導(dǎo)水率

當(dāng)生物炭在土壤中條狀分布時，土壤區(qū)的面積視為αA，生物炭區(qū)的面積為βA，結(jié)構(gòu)如圖1b所示。整個混合土壤截面的有效導(dǎo)水率為Ke、水力梯度為i，整個混合土層的水流通量應(yīng)等于兩區(qū)域水流通量之和，即

KeAi=K1αAi+K2βAi

(6)

故

Ke=αK1+βK2

(7)

當(dāng)生物炭在土壤中呈條狀分布時，混合土壤的飽和導(dǎo)水率可以近似地用式(7)估計。

由于生物炭在土壤中是隨機均勻分布的，雖然式(5)和式(7)均不能給出飽和導(dǎo)水率的準(zhǔn)確值，但通過這2個公式大致能夠估算出飽和導(dǎo)水率。

2.3 毛細(xì)管模型推導(dǎo)添加生物炭土壤的導(dǎo)水率

前人以試驗為基礎(chǔ)，提出很多測定土壤水分運動參數(shù)的方法。由于一些模型參數(shù)較難獲得[22]，因此限制了水分運動的數(shù)值模擬計算。本文利用相對容易獲得的土壤水分特征曲線，推導(dǎo)了添加生物炭土壤的毛細(xì)管理論，給出了土壤水分特征曲線和非飽和導(dǎo)水率的關(guān)系。

當(dāng)水與毛細(xì)管管壁完全浸潤時，水在毛細(xì)管內(nèi)上升的最大高度為

(8)

式中ρW——水的密度

H——水在毛細(xì)管內(nèi)上升的最大高度

g——重力加速度

R——毛細(xì)管的曲率半徑

σ——水的表面張力系數(shù)

假定添加生物炭的土壤可以用不同形態(tài)的毛細(xì)管束來描述，每個毛細(xì)管束的長度為Le，土壤的實際長度為L，土壤的橫截面積為A，土壤兩端的水頭差為ΔH。根據(jù)泊肅葉公式，對于曲率半徑為Rj的毛細(xì)管，其流量為

(9)

式中μ——液體粘滯系數(shù)

當(dāng)土壤中的所有毛細(xì)管充滿水時，水流通量為

(10)

其中

式中Nj——半徑為Rj的毛細(xì)管的數(shù)量

P——不同半徑毛細(xì)管的數(shù)量

nj——單位面積內(nèi)半徑為Rj的毛細(xì)管數(shù)量

由于土壤的基質(zhì)吸力與土壤的毛細(xì)管孔徑有關(guān)，因此土壤水分特征曲線也描述了毛細(xì)管模型的尺寸分布規(guī)律。通過未添加生物炭土壤(CK)的水分特征曲線以及生物炭的水分特征曲線，能夠給出添加生物炭后土壤的非飽和導(dǎo)水率的范圍。

根據(jù)土壤的含水率將土壤水分特征曲線畫成M個相同寬度的區(qū)間，如圖2所示，每個區(qū)間寬度為Δθ，每一個區(qū)域?qū)?yīng)一個基質(zhì)吸力hj(即與該吸力對應(yīng)的所有毛細(xì)管都具有相同孔徑)和土壤含水率θs-jΔθ，即

hj=h(θs-jΔθ)

(11)

圖2 計算非飽和導(dǎo)水率的毛細(xì)管模型Fig.2 Capillary model for calculating unsaturated hydraulic conductivity

另外，由于生物炭和土壤顆粒共同在一起，所以同一位置處吸力相同，即在hj吸力下，根據(jù)生物炭的水分特征曲線θW(h)得到生物炭中的含水率差值為

ΔθWj=θW(hj-1)-θW(hj)

(12)

θ(h0)=θWS，θWS為生物炭的飽和含水率。根據(jù)式(8)，基質(zhì)吸力hj對應(yīng)的毛細(xì)管孔徑為

(13)

對于土壤和生物炭，對應(yīng)的j區(qū)的單位面積內(nèi)毛細(xì)管數(shù)量分別為

(14)

(15)

式中nk——j區(qū)生物炭內(nèi)單位面積的毛細(xì)管數(shù)量

ΔθWk——生物炭的含水率區(qū)間寬度

Rk——生物炭內(nèi)毛細(xì)管半徑

當(dāng)h=hj時，所有孔徑小于Rj的毛細(xì)管充滿了水，提供了土壤的水流通量，并且這些毛細(xì)管是由土壤顆粒和生物炭顆粒共同提供的。另外，一般認(rèn)為[23]生物炭的導(dǎo)水率遠(yuǎn)大于土壤的導(dǎo)水率，意味著生物炭中毛細(xì)管的密度大于土壤。如果假定生物炭中多于土粒中毛細(xì)管的倍數(shù)可以由飽和導(dǎo)水率的比值近似給出，再同時考慮式(10)，得到混合土壤的水流通量表達(dá)式為

(16)

其中

q=Kcs/Kss

式中q——生物炭與土壤飽和導(dǎo)水率的比值

Kcs——生物炭飽和導(dǎo)水率

Kss——土壤飽和導(dǎo)水率

再由式(13)～(15)，得到

(17)

式中hk——區(qū)間k內(nèi)毛細(xì)管的吸力

N——生物炭水分特征曲線劃分的最大區(qū)間數(shù)

(18)

其中

l=L/Le

式中l(wèi)——毛細(xì)管彎曲度

式(18)與達(dá)西定律相比，得到混合土壤的非飽和導(dǎo)水率為

(19)

當(dāng)毛細(xì)管全部充滿水，就得到混合土壤的飽和導(dǎo)水率

(20)

結(jié)合式(19)和式(20)得到混合土壤的非飽和導(dǎo)水率

(21)

基于土壤和生物炭的水分特征曲線，再根據(jù)添加生物炭的比例以及土壤和生物炭的飽和導(dǎo)水率，可以用式(21)近似給出添加不同量生物炭后土壤的非飽和導(dǎo)水率。

3 結(jié)果與分析

3.1 土壤水分特征曲線

用離心機法測定的土壤和生物炭的水分特征曲線如圖3所示。

圖3 土壤與生物炭的水分特征曲線Fig.3 Water characteristic curves of soil and biochar

采用VAN- GENUCHTEN模型，利用Retc軟件對離心機法的實測數(shù)據(jù)進(jìn)行擬合，擬合的飽和含水率為0.418 1 cm3/cm3，殘余含水率0.009 8 cm3/cm3。由圖3可以看出，當(dāng)土壤吸力小于30 cm時，生物炭與土壤幾乎平行，當(dāng)土壤吸力介于30～2 200 cm，土壤的SWCC位于生物炭上方,當(dāng)吸力大于2 200 cm時，生物炭的SWCC位于土壤上方，即生物炭在中高吸力段的含水率大于土壤，因此添加生物炭后土壤持水能力的增加，主要依靠的是生物炭在中高吸力段持水能力的提升。也有文獻(xiàn)表明，生物炭的添加減小了土壤容重、增加了土壤孔隙度[24]，使得土壤持水能力增強。

3.2 毛細(xì)管模型計算土壤的非飽和導(dǎo)水率

由于毛細(xì)管模型綜合考慮了生物炭和土壤的性能，因此根據(jù)毛細(xì)管模型公式(21)以及測定的土壤和生物炭的水分特征曲線，可以計算不同添加量生物炭土壤的非飽和導(dǎo)水率。采用圖3數(shù)據(jù)計算毛細(xì)管模型的數(shù)據(jù)如表1所示。計算中采用的生物炭和土壤的飽和導(dǎo)水率實際測量值見表2。由表1可以看出，[K(θs-jΔθ)/Ks]1/3項可以近似看作有效飽和度。在低吸力段，有效飽和度隨著生物炭含量的增加從0.710 446增加到0.739 290，有效飽和度的增高可看作土壤導(dǎo)水性能得到提高[20]；在35.82～1 590 cm吸力段，有效飽和度隨生物炭含量的增加而減小，土壤導(dǎo)水性能降低。導(dǎo)水性能的增加是由于生物炭在低吸力下的高導(dǎo)水率，而高吸力下導(dǎo)水性能的降低，是由于生物炭在高吸力下的低導(dǎo)水率。

3.3 土壤飽和導(dǎo)水率測量結(jié)果

表2為實際測量的CK、C1、C2、C3、C4的飽和導(dǎo)水率，以及采用式(5)和式(7)計算的估計值。

表1 采用毛細(xì)管模型計算不同添加量生物炭土壤的非飽和導(dǎo)水率Tab.1 Calculated unsaturated hydraulic conductivity by using capillary model

表2 土壤飽和導(dǎo)水率的測量值和估計值Tab.2 Tested and estimated values of soil saturatedhydraulic conductivity mm/h

由表2可以看出，C1、C2、C3、C4飽和導(dǎo)水率的實際測量值，均在層狀分布和條狀分布的估計值之間，但非常接近于層狀分布，與層狀分布相比的最大偏差為8.9%。分析其原因，可能是當(dāng)土壤中的生物炭施加量較少時，不足以讓生物炭粒子連接起來，所以應(yīng)該看作層狀。另外，試驗數(shù)據(jù)表明[16]，添加生物炭土壤中的水分流動，多為水在生物炭和土粒之間的空隙流動。由于層狀分布的估計值沒有考慮生物炭和土粒之間的空隙，當(dāng)生物炭添加后，空隙間的水分流動會對飽和導(dǎo)水率產(chǎn)生影響，因而導(dǎo)致實際測量的飽和導(dǎo)水率比層狀估算值大。

4 添加生物炭土壤滲水過程的數(shù)值分析

4.1 添加生物炭后的土壤模型

為研究相同滲透時間后，不同生物炭添加量土樣中的水分分布，首先要對添加生物炭的土壤進(jìn)行數(shù)值模擬。由于實際土壤的復(fù)雜性，計算的模型擬采用二維形式。添加生物炭的土壤被分成土壤區(qū)和生物炭區(qū)。計算時，這兩個區(qū)具有各自的非飽和導(dǎo)水率和初始水分特征。

由于生物炭是隨機分布在土壤中，采用Python語言編程的方法給定土壤截面。將生物炭隨機分布在土壤截面中的步驟如下：首先建立計算土壤的面積區(qū)域[25]。面積選擇過大，生物炭顆粒的數(shù)量隨之變大，導(dǎo)致計算量過大難以實現(xiàn)。因此，計算的面積選取10 mm×157 mm。其次確定土樣中生物炭的個數(shù)，生物炭的個數(shù)通過添加生物炭的比例來校對。本研究中，添加生物炭的體積比約為2%、4%、6%、8%。根據(jù)添加比例和試驗所用生物炭的平均粒徑給出生物炭的個數(shù)(生物炭粒徑范圍1.5～2.5 mm，其中2 mm粒徑約占60%以上，故取2 mm)。隨機產(chǎn)生生物炭所在位置的圓心坐標(biāo)，畫圓代表生物炭。校驗是否有圓相交，如果有，刪除后再補充相應(yīng)數(shù)量的圓。最后，對生成的所有生物炭顆粒求面積，核算添加的比例和實際摻炭比例是否相等。完成上述工作后賦予土壤區(qū)和生物炭區(qū)不同的水分運動參數(shù)。圖4為最終生成的隨機添加生物炭的土壤模型。

圖4 隨機分布生物炭顆粒的土壤截面 Fig.4 Profiles of biochar particles distributed randomly in soil

4.2 有限元采用的模型

采用通用的有限元分析軟件Abaqus進(jìn)行數(shù)值計算。 對于飽和土，滲透系數(shù)可以視為常量，但對于液體與氣體并存的非飽和土，飽和度的大小直接影響滲透的結(jié)果，因此滲透系數(shù)K是飽和度S的函數(shù)。Abaqus采用BROOKS等[26]提出的滲透系數(shù)和飽和度的關(guān)系。利用土壤水分特征曲線給出的3個參數(shù)，即土壤的進(jìn)氣值、有效飽和度以及空隙尺寸分布指標(biāo)λ，滲透系數(shù)函數(shù)與有效飽和度的關(guān)系為(下文中出現(xiàn)的飽和度均指有效飽和度)

(22)

其中

式中Sr——有效飽和度

δ——經(jīng)驗常數(shù),取3

(ua-uw)——基質(zhì)吸力

(ua-uw)b——進(jìn)氣值

4.3 實際土樣與毛細(xì)管模型土樣的數(shù)值計算

為研究不同土樣經(jīng)過相同滲透時間后的水分分布，數(shù)值模擬部分首先計算添加生物炭土壤的飽和度，其次采用前面推導(dǎo)的毛細(xì)管模型等效的非飽和導(dǎo)水率計算摻炭土壤的飽和度，初始給定的土壤吸力是相同的。

當(dāng)區(qū)域為對照土樣(CK)時，此時α=1，采用的吸力和對應(yīng)的飽和度如表1所示，實際上這是采用單一土壤的毛細(xì)管模型給出土壤的非飽和導(dǎo)水率，飽和導(dǎo)水率Ks采用試驗測得的值1.152 mm/h。邊界條件為：左邊界初始吸力為零表示進(jìn)水，其他邊界為不透水邊界。全部區(qū)域的初始吸力和飽和度分別為200 cm、0.152 7，滲透時間5 000 s，計算出CK的飽和度隨位置的分布規(guī)律如圖5a所示。

當(dāng)區(qū)域為C4時，區(qū)域如圖4c所示，采用的吸力和飽和度如表1所示，土壤區(qū)為α=1，生物炭區(qū)為α=0。土壤和生物炭的飽和導(dǎo)水率分別為1.152、22.068 mm/h，邊界條件仍為左邊界初始吸力為零，土壤的初始吸力和飽和度為200 cm、0.152 7，生物炭的初始吸力和飽和度為為200 cm、0.060 7，滲透時間5 000 s，計算出的C4的飽和度隨位置的分布規(guī)律如圖5b所示。

當(dāng)區(qū)域采用C4的毛細(xì)管模型C4*時，采用的吸力和飽和度如表1所示。飽和導(dǎo)水率為試驗測得的1.371 6 mm/h，邊界條件仍為左邊界初始吸力為零，C4*初始吸力和飽和度為200 cm、0.086 5，滲透時間5 000 s，計算后C4*的飽和度隨位置的分布規(guī)律如圖5c所示。圖5表明，毛細(xì)管模型很好地描述了施加生物炭后，土壤飽和度隨位置的分布規(guī)律。

計算出的CK、C4中的土壤、C4中的生物炭以及C4*的飽和度隨位置的分布規(guī)律如圖6a所示。

圖6 滲透5 000 s后的土壤飽和度曲線Fig.6 Saturation of soil after infiltration of 5 000 s

由圖6a可以看到，C4*的飽和度分布更類似于C4中生物炭的飽和度分布形式，這是由于采用毛細(xì)管模型(式(21))的推導(dǎo)結(jié)果。由于生物炭中毛細(xì)管的密度遠(yuǎn)多于土壤，與同體積的土壤相比，能夠吸收更多的水分，因此對土壤的飽和度影響更大。式(21)中的q值決定了C4*飽和度曲線的形狀。圖6a中，C4*的飽和度曲線在水平距離上略微超前于C4中生物炭的飽和度曲線，說明計算C4*飽和度時，采用測量的飽和導(dǎo)水率偏大。這可能是由于C4的計算中，沒有考慮生物炭和土壤的空隙，而這是影響飽和導(dǎo)水率的重要因素。

設(shè)定滲透時間5 000 s，按上述方法計算C1、C2、C1*、C2*的飽和度分布規(guī)律。C1*、C2*飽和導(dǎo)水率為測得的1.195 2、1.249 2 mm/h。毛細(xì)管模型采用的土壤初始吸力均為200 cm，C1*和C2*的飽和度分別為0.114 4和0.099 9。圖6b為滲透時間5 000 s后土樣的飽和度隨位置的分布曲線。C1*、C2*為毛細(xì)管模型模擬的土樣(由于圖線過于密集，C3沒有畫出)。

由圖6可以看到，水分滲透一定時間后，C1*、C2*、C4*曲線均落在CK上方，即相同位置上添加生物炭土壤的飽和度都比CK大，這說明添加生物炭可以提高土壤的導(dǎo)水率。并且導(dǎo)水率隨著生物炭添加量的增大而增大。C4*曲線落在CK下方，是由于高吸力下生物炭導(dǎo)水性能的降低造成的。毛細(xì)管模型綜合考慮了生物炭和土壤的特性，能較好地反映添加生物炭后土壤的導(dǎo)水特性。

4.4 土柱的數(shù)值計算和試驗數(shù)據(jù)的比較

本試驗采用的土柱尺度較大，因此在有限元計算時，如果考慮實際的生物炭粒子，劃分的網(wǎng)格必須很密，這樣會導(dǎo)致計算量過大，所以一般無法對土柱直接計算，這時可以采用毛細(xì)管模型近似計算。計算時，由于C1*、C2*、C4*的飽和含水率和殘余含水率各不相同，因此C1*、C2*、C4*的初始飽和度是不同的。為簡化計算，給定C1*、C2*、C4*的初始飽和度，與試驗測得CK的初始飽和度相同。

為驗證毛細(xì)管模型的正確性，設(shè)定水分滲透時間為24 h，實際測量了C1、C2、C4在80 cm土柱中的含水率，采用CK擬合的飽和含水率和殘余含水率值計算飽和度。計算的初始數(shù)據(jù)如下：實際測定的CK初始飽和度0.294，與之對應(yīng)的土壤的吸力約為68.6 cm。毛細(xì)管模型在此飽和度下的吸力分別為51.9 cm(C1*)、43.4 cm(C2*)、40.4 cm(C4*)。圖7為水分滲透24 h后，飽和度隨土柱位置變化關(guān)系圖。

圖7 滲透24 h后土柱飽和度的計算值和測量值Fig.7 Calculated and measured values of soil saturation after infiltration of 24 h

由圖7的計算結(jié)果看出，土柱中水分滲透24 h后，對于同一生物炭添加比例，毛細(xì)管模型的飽和度與實際測量值相比略低。這是由模型中沒有考慮實際土壤中生物炭粒子和土壤之間的空隙造成的。生物炭體積比為2%土壤的誤差范圍為4.27%～7.24%，4%土壤的誤差范圍為3.06%～8.6%，8%土壤的誤差范圍為2.90%～10.11%?？梢姡S著生物炭體積比的增大，數(shù)值計算結(jié)果與毛細(xì)管模型的計算結(jié)果偏差有所增加，這種偏差是因為在毛細(xì)管模型假設(shè)中，生物炭中的不同孔徑的毛細(xì)管數(shù)量是均勻分布的，總是用一個恒定的飽和導(dǎo)水率的比值近似給出，這與實際是有偏差的。由圖7結(jié)果看出，在距供水端較遠(yuǎn)處，即土壤吸力較大處，土壤的實際導(dǎo)水率大于毛細(xì)管模型的值，說明生物炭中的細(xì)小毛細(xì)管的數(shù)量多于估計值。因此，采用毛細(xì)管模型的數(shù)值計算，對于中等及以下生物炭施用量，在土壤水力特征的定量分析上更加準(zhǔn)確。毛細(xì)管模型可以簡化土壤的水力特性，可以用來對實際的大塊土地進(jìn)行數(shù)值分析。

5 討論

本研究針對東北黑土區(qū)的草甸黑土，推導(dǎo)了施加生物炭后土壤的毛細(xì)管模型，利用此模型，研究了一定滲透時間后土壤中的水分分布，選擇飽和度等參數(shù)將添加不同量生物炭的土壤進(jìn)行數(shù)值模擬，最后利用土柱的實際測定結(jié)果對數(shù)值模擬結(jié)果進(jìn)行驗證，以此開展生物炭對土壤水分影響的數(shù)值化定量化研究。

從土壤水分特征曲線看出，添加生物炭后土壤的持水能力增加，由于生物炭在中高吸力段的含水率大于參照組土壤，因此添加生物炭后土壤持水能力的增加，主要依靠的是生物炭在中高吸力段持水能力的提升。毛細(xì)管模型表明，施加生物炭后土壤性能的改善，是由于低吸力段生物炭的良好導(dǎo)水特性。利用毛細(xì)管模型計算的土壤非飽和導(dǎo)水率，依賴于土壤的水分特征曲線和飽和導(dǎo)水率，為得到精確的數(shù)值計算結(jié)果，必須精確測量土壤水分特征曲線。

土柱一端加水滲透24 h后，數(shù)值模擬的土壤水分分布值比實際測量值略低，主要是由于毛細(xì)管模型中生物炭多于土壤的毛細(xì)管數(shù)量，總是用一個恒定的值近似給出，這與實際有偏差。另外是由于模型中生物炭與土壤的空隙分布不能與實際情況完全吻合造成的。

雖然直接測定的結(jié)果更接近真實值，但測量過程繁瑣、消耗時間長，因此采用毛細(xì)管模型的有限元數(shù)值計算優(yōu)勢明顯。另外，對于只能通過試驗研究解決的問題，如一定范圍的生物炭施加量、不同的施加面積、地勢情況(如坡耕地等)、不同的滲透時間等情況，數(shù)值模擬能夠更好更快地給出近似答案。采用毛細(xì)管模型還可以對試驗數(shù)據(jù)進(jìn)行校驗，找出試驗中有偏差的數(shù)據(jù)，有助于提高試驗結(jié)果的準(zhǔn)確性，從而為研究土壤的水力特性提供更多的方法與保障。

6 結(jié)論

(1)試驗測得的水分特征曲線表明，當(dāng)土壤吸力大于2 200 cm時，生物炭的殘余含水率高于對照組，因此添加生物炭后土壤持水能力的增加，主要依靠的是生物炭在中高吸力段持水能力的提升。

(2)通過飽和導(dǎo)水率試驗數(shù)據(jù)和理論估計值的比較可以看出，實際測量的飽和導(dǎo)水率與生物炭層狀分布的飽和導(dǎo)水率非常接近，最大偏差為8.9%，說明添加生物炭土壤更應(yīng)考慮為層狀分布。試驗數(shù)據(jù)和理論數(shù)據(jù)的差別可能是生物炭和土壤之間的空隙增加造成的。

(3)毛細(xì)管模型結(jié)果顯示，添加生物炭土壤導(dǎo)水率的提高，主要是由生物炭的導(dǎo)水特性決定的，添加的生物炭主要改變了土壤在低吸力下的導(dǎo)水率。

(4)數(shù)值分析顯示，施用生物炭能夠提高東北草甸黑土的導(dǎo)水率。添加生物炭的比例越高，導(dǎo)水率增加越大。

(5)生物炭中的不同孔徑的毛細(xì)管數(shù)量不是均勻分布的，孔徑越小，毛細(xì)管的數(shù)量越多。

(6)毛細(xì)管模型簡單有效，生物炭體積比小于8%時，數(shù)值分析的結(jié)果與試驗數(shù)據(jù)吻合良好，說明該模型可以為土壤水分運動參數(shù)的數(shù)值計算提供比較準(zhǔn)確的方法。

1 何緒生, 耿增超, 佘雕, 等. 生物炭生產(chǎn)與農(nóng)用的意義及國內(nèi)外動態(tài)[J]. 農(nóng)業(yè)工程學(xué)報, 2011,27(2):1-7．

HE Xusheng, GENG Zengchao, SHE Diao, et al. Implications of production and agricultural utilization of biochar and its international dynamic[J]. Transactions of the CSAE, 2011, 27(2):1-7. (in Chinese)

2 LEHMANN J, JOSEPH S. Biochar for environmental management: an introduction [J]. Science and Technology, 2009, 25(1):15801-15811.

3 GITHINJI L. Effect of biochar application rate on soil physical and hydraulic properties of a sandy loam[J]. Archives of Agronomy and Soil Science, 2013:1-14.

4 STREUBEL J D, COLLINS H P, GARCIA P M, et al. Influence of contrasting biochar types on five soils at increasing rates of application[J]. Soil Biology & Biochemistry,2011,75:1402-1413.

5 WEST L T, ABREW M A, BISHOP J P. Saturated hydraulic conductivity of soils in the Southern Piedmont of Georgia, USA: field evaluation and relation to horizon and landscape properties[J]. CATENA, 2008, 73(2): 174-179.

6 ABEL S, PETERS A, TRINKS S, et al. Impact of biochar and hydrochar addition on water retention and water repellency of sandy soil[J]. Geoderma, 2013, 202-203:183-191.

7 魏永霞, 劉志凱, 馮鼎銳, 等. 生物炭對草甸黑土物理性質(zhì)及雨后水分動態(tài)變化的影響[J/OL]. 農(nóng)業(yè)機械學(xué)報, 2016, 47(8): 201-207. http:∥www.j-csam.org/jcsam/ch/reader/view_abstract.aspx?flag=1&file_no=20160825&journal_id=jcsam.DOI:10.6041/j.issn.1000-1298.2016.08.025.

WEI Yongxia, LIU Zhikai, FENG Dingrui, et al. Influences of the biochar on physical properties of the meadow black soil and dynamic changes of soil water after the individual rainfall[J/OL]. Transactions of the Chinese Society for Agricultural Machinery, 2016, 47(8):201-207.(in Chinese)

8 LI Y, GAO R, YANG R, et al. Using a simple soil column method to evaluate soil phosphorus leaching risk[J]. CLEAN-Soil, Air, Water, 2013, 41(11):1100-1107.

9 ZIMMERMAN A R. Abiotic and microbial oxidation of laboratory-produced black carbon (biochar)[J]. Environmental Science & Technology, 2010, 44(4): 1295-1301.

10 CHAN K Y, VAN ZWIETEN L, MESZAROS L, et al. Agronomic values of greenwaste biochar as a soil amendment[J]. Australian Journal of Soil Research,2007, 45:629-634.

11 FENG L, LI G T, LIN Q M, et al. Crop yield and soil properties in the first 3 years after biochar application to a Calcareous soil[J]. Journal of Integrative Agriculture,2014, 13(3): 525-532.

12 勾芒芒, 屈忠義, 楊曉, 等. 生物炭對砂壤土節(jié)水保肥及番茄產(chǎn)量的影響研究[J/OL]. 農(nóng)業(yè)機械學(xué)報, 2014, 45(1): 137-142.http:∥www.j-csam.org/jcsam/ch/reader/view_abstract.aspx?flag=1&file_no=20140122&journal_id=jcsam.DOI:10.6041/j.issn.1000-1298.2014.01.022.

GOU Mangmang, QU Zhongyi, YANG Xiao, et al. Study on the effects of biochar on saving water, preserving fertility and tomato yield[J/OL]. Transactions of the Chinese Society for Agricultural Machinery, 2014, 45(1):137-142. (in Chinese)

13 AKHTAR S S, LI G, ANDERSEN M N, et al. Biochar enhances yield and quality of tomato under reduced irrigation[J]. Agricultural Water Management, 2014, 138:37-44.

14 ASAI H, SAMSON B K, STEPHAN H M, et al. Biochar amendment techniques for upland rice production in Northern Laos 1. Soil physical properties, leaf SPAD and grain yield[J]. Field Crops Research, 2009, 111(1): 81-84.

15 BODNER G, SCHOLL P, KAUL H P. Field quantification of wetting-drying cycles to predict temporal changes of soil pore size distribution[J]. Soil & Tillage Research, 2013, 133(5):1-9.

16 LIU Z, DUGAN B, MASIELLO C A, et al. Biochar particle size, shape, and porosity act together to influence soil water properties.[J]. Plos One, 2017, 12(6):e0179079.

17 ALLETTO L, COQUET Y, VACHIER P, et al. Hydraulic conductivity, immobile water content, and exchange coefficient in three soil profiles[J]. Soil Science Society of America Journal, 2006, 70(4):1272-1280.

18 SHE D, GAO X, WANG P, et al. Comparison of soil hydraulic properties with different levels of soil salinity and sodicity[J]. Arabian Journal of Geosciences, 2015, 8(8):5351-5360.

19 ROGER E J, RICHARD G, DEXTER A R, et al. Integration of soil structure variations with time and space into models for crop management: a review[J]. Agronomy for Sustainable Development, 2009, 29(1):135-142.

20 尚熳廷,馮杰,劉佩貴,等. SWCC 測定時吸力計算公式與最佳離心時間的探討[J]. 河海大學(xué)學(xué)報, 2009, 37(1): 12-15.

SHANG Manting, FENG Jie, LIU Peigui, et al. On formulas for soil suction and optimum centrifugal time by use of SWCC[J]. Journal of Hohai University, 2009, 37(1): 12-15. (in Chinese)

21 VAN GENUCHTEN M T H. A closed-form equation for predicting the hydraulic conductivity of unsaturated soils[J]. Soil Science Society of America Journal, 1980, 44(5):892-898.

22 BARNES R T, GALLAGHER M E, MASIELLO C A, et al. Biochar-induced changes in soil hydraulic conductivity and dissolved nutrient fluxes constrained by laboratory experiments[J]. Plos One, 2014, 9(9):e108340.

23 YARGICOGLU E N, SADASIVAM B Y, REDDY K R, et al. Physical and chemical characterization of waste wood derived biochars[J]. Waste Management, 2015, 36:256-268.

24 MELO C M T, HOLANDA N M, MADARI B E, et al. Biochar increases plant-available water in a sandy loam soil under an aerobic rice crop system[J].Solid Earth,2014,5:939-952.

25 WONGSTO A, LI S. Micromechanical FE analysis of UD fibre-reinforced composites with fibres distributed at random over the transverse cross-section[J]. Composites Part A: Applied Science & Manufacturing, 2005, 36(9):1246-1266.

26 BROOKS R H, COREY A T. Hydraulic properties of porous media and their relation to drain design[J].Transactions of the ASAE, 1964,7(1):26-28.