生物炭对东北草甸黑土水力特性影响的数值化研究

王睿垠 魏永霞 张翼鹏 张雨凤 马德才

(1.东北农业大学理学院， 哈尔滨 150030； 2.东北农业大学水利与土木工程学院， 哈尔滨 150030；3.农业部农业水资源高效利用重点实验室， 哈尔滨 150030； 4.中山大学中法核工程与技术学院， 珠海 519082)

0 引言

东北黑土区是世界三大黑土带之一，也是国家商品粮生产基地，粮食产量对于保障国家安全具有重要意义。水分是农业生产的重要因素，生物炭作为土壤改良剂对土壤水分运动的影响尚不明确。目前针对东北黑土施加生物炭的研究多为试验研究，少有生物炭对土壤水分影响的数值研究。土壤是一个十分复杂的系统，是植物生长中水分和养分的主要来源，是由固、液、气三相物质组成的体系。生物炭是生物有机质在无氧或缺氧情况下，经高温裂解的产物[1]。研究表明，在某些情况下，土壤中添加生物炭，可以修复土壤[2]并改变土壤的持水、导水特性[3-4]。土壤的导水率是土壤粒径、黏土含量、有机质含量以及粒子大小分布等的函数，是全面反映土壤导水特性的一个参量[5]。土壤具有较高的导水率时可以加快渗透和排水[6]，土壤快速排水有利于减少径流的发生[7]，但渗透速度过快，也会减少水中营养物和农用化学品充分溶解过滤的机会[8]。黑土是中国重要的土壤资源，由于黑土表土疏松、底土黏重以及人为不合理的耕作，使得东北黑土持水能力降低，农田生产力下降。生物炭作为土壤改良剂，在改善土壤的持水和导水能力、提高作物的生产力方面与其他有机添加物相比，更加安全和持久[9]。研究证实了土壤添加生物炭可以影响作物的产量[10-12]，其物理机制多归因于土壤孔隙度的增加、土壤容重的降低[13]、土壤保水导水能力的提高以及土壤可利用养分的增加[14]。研究表明，土壤结构在时间上还会受到湿润与干燥循环、生物活性以及农业作业的影响[15]。

虽然添加生物炭对土壤水力学特性产生不同程度的影响，但有关预测生物炭定量改善土壤水分特性、生物炭决定水分分布的确切解的研究[16]鲜见报道。有些研究给出了土壤水动力的空间变化，主要是土壤粒子分布的不均匀性[17]。土壤中的水分运动，取决于固体颗粒的尺寸、形状和分布[18]。土壤的多样性以及水分运动的复杂性，使得试验的数量和耗时巨大。本文推导添加生物炭土壤的毛细管模型[19]，并在此模型下利用有限元法或差分法对土壤进行数值模拟。对比此模型与添加生物炭土壤土柱渗透的数值分析结果和试验结果，以证实此模型预测不同条件下的土壤水分分布或者校验试验结果的正确性。

1 材料与方法

1.1 研究区概况

取土地点位于黑龙江省农垦北安管理局红星农场实验基地(48°10′E，127°1′N)。该地区为典型的东北黑土区，属寒温带大陆性季风气候，冬季寒冷干燥，夏季高温多雨，降雨时间集中，年降水量多介于450～650 mm，其中7—9月累积降水量约402 mm，约占年平均降水量的75%，因此旱涝灾害频繁，水土流失严重。该区域的土壤多为壤土和粘壤土，以粗粉沙和粘粒最多，各占30%～40%，容重1.0～1.3 g/cm3，透水、持水、通气性均较好，总孔隙度多在40%～60%，毛管孔隙可占20%～30%，通气孔隙度约占20%。

1.2 试验材料与设备

供试材料为玉米秸秆生物炭，购自辽宁金和福农业开发有限公司，制备方式为450℃无氧条件下高温裂解。基本理化性质：粒径范围1.5～2.5 mm，其中2 mm粒径占60%以上。各指标如下：全碳质量分数72.21%，全氮质量分数1.56%，全磷质量分数0.72%，全钾质量分数1.64%，pH值为9.35，经测量生物炭密度约596 kg/m3。

脱水面半径为9.8 cm的CR- 21G3型高速离心机，用于室内测定土壤水分特征曲线；自制直径8 cm、长为80 cm的有机玻璃筒，距供水端10、20、30、40、50、60、70 cm处设有取土孔，用于测量不同位置的土壤含水率。

1.3 试验方法与内容

土样为田间试验土样。取土小区面积为20 m×5 m，坡度3°。设置5种生物炭用量添加比例，即不添加生物炭处理(CK)，以及生物炭添加量为25、50、75、100 t/hm2的处理。试验小区采用相同的水肥管理，每个处理重复3次。2017年5月中旬将生物炭人工均匀混入耕层黑土(0～20 cm)，生物炭占土壤的体积比约为2%(C1)、4%(C2)、6%(C3)、8%(C4)，取样时间为2017年9月25日。用环刀采样器取对照土样(CK)，取土深度为10～15 cm，环刀高5 cm、体积100 cm3，取土后将环刀用蒸馏水吸水48 h，至饱和后称量。生物炭放置环刀中，用蒸馏水吸水并压实48 h后称量。用离心机法分别测定土壤和生物炭的水分特征曲线。离心机转速设定为500～6 000 r/min，离心时间为100 min[20]，离心结束后称量，土壤与生物炭各重复3次。

室内试验用定水头法分别测定生物炭和5种处理的饱和导水率。利用马氏瓶供水，维持水头恒定为5 cm，每隔10 min记录出流水量，出流水量保持稳定时结束试验，计算饱和导水率，各处理重复3次。

测定土柱的含水率。在5个处理小区取原状土，取土深度10～15 cm，土样经风干、破碎、筛分后，实地测量土壤干容重为1.1 g/cm3。将土样按此容重每层3 cm分层回填入土柱，装填好的土柱直径8 cm、高80 cm。进水端采用马氏瓶供水，视进水端的土壤为饱和，当湿润锋到达土柱约3/4处，在各取土孔处迅速取土，测定各处的含水率。各处理重复3次。

1.4 数据处理方法

采用Excel 2010进行数据整理，所用数据均为各重复测定的平均值，土壤水分特征曲线采用Retc软件进行拟合，采用python语言建模编程，采用通用的有限元分析软件Abaqus 6.13进行数值模拟计算。

2 理论分析

2.1 土壤水分特征曲线

目前学者大多采用VAN- GENUCHTEN(VG)模型给定土壤水分特征曲线(SWCC)。模型由VAN GENUCHTEN[21]在1980年提出，适用于描述低含水率时的土壤水分特征曲线。其表达式为

(1)

其中

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

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

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

h——土壤吸力，cm

a、n、m——曲线拟合参数

θ——土壤含水率

本研究采用VAN- GENUCHTEN模型，利用试验测得土壤含水率，拟合土壤饱和含水率和残余含水率。

2.2 添加生物炭土壤的饱和导水率

2.2.1假设生物炭在土壤中呈层状分布的导水率

稳定水流分别穿过土壤和生物炭时，设K1为土壤饱和导水率，K2为生物炭饱和导水率。当生物炭在土壤中呈层状分布时，设生物炭与土壤的混合土层的总厚度为L，令混合土层中土壤体积比为α、生物炭体积比为β，则α+β=1，因此土壤层的厚度为αL，生物炭层厚度为βL，结构如图1a所示。设土壤层与生物炭层水力梯度分别为i1、i2，整个混合土层有效导水率为Ke，整个混合土层水头损失应等于两层水头损失之和，即

iL=αi1L+βi2L

(2)

式中i——总的水力梯度

JW=Kei=K1i1=K2i2

(3)

式中JW——水流通量

结合式(2)可以得到

(4)

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

一般情况下，土壤体积比α远大于生物炭体积比β，另外，生物炭的饱和导水率K2也远大于草甸黑土的饱和导水率K1，这样式(4)可以进一步简化为

(5)

当生物炭在土壤中层状分布时，混合土壤的饱和导水率可以近似地认为是土粒的饱和导水率除以土粒在土壤中的比例。

2.2.2假设生物炭在土壤中呈条状分布的导水率

当生物炭在土壤中条状分布时，土壤区的面积视为αA，生物炭区的面积为βA，结构如图1b所示。整个混合土壤截面的有效导水率为Ke、水力梯度为i，整个混合土层的水流通量应等于两区域水流通量之和，即

KeAi=K1αAi+K2βAi

(6)

故

Ke=αK1+βK2

(7)

当生物炭在土壤中呈条状分布时，混合土壤的饱和导水率可以近似地用式(7)估计。

由于生物炭在土壤中是随机均匀分布的，虽然式(5)和式(7)均不能给出饱和导水率的准确值，但通过这2个公式大致能够估算出饱和导水率。

2.3 毛细管模型推导添加生物炭土壤的导水率

前人以试验为基础，提出很多测定土壤水分运动参数的方法。由于一些模型参数较难获得[22]，因此限制了水分运动的数值模拟计算。本文利用相对容易获得的土壤水分特征曲线，推导了添加生物炭土壤的毛细管理论，给出了土壤水分特征曲线和非饱和导水率的关系。

当水与毛细管管壁完全浸润时，水在毛细管内上升的最大高度为

(8)

式中ρW——水的密度

H——水在毛细管内上升的最大高度

g——重力加速度

R——毛细管的曲率半径

σ——水的表面张力系数

假定添加生物炭的土壤可以用不同形态的毛细管束来描述，每个毛细管束的长度为Le，土壤的实际长度为L，土壤的横截面积为A，土壤两端的水头差为ΔH。根据泊肃叶公式，对于曲率半径为Rj的毛细管，其流量为

(9)

式中μ——液体粘滞系数

当土壤中的所有毛细管充满水时，水流通量为

(10)

其中

式中Nj——半径为Rj的毛细管的数量

P——不同半径毛细管的数量

nj——单位面积内半径为Rj的毛细管数量

由于土壤的基质吸力与土壤的毛细管孔径有关，因此土壤水分特征曲线也描述了毛细管模型的尺寸分布规律。通过未添加生物炭土壤(CK)的水分特征曲线以及生物炭的水分特征曲线，能够给出添加生物炭后土壤的非饱和导水率的范围。

根据土壤的含水率将土壤水分特征曲线画成M个相同宽度的区间，如图2所示，每个区间宽度为Δθ，每一个区域对应一个基质吸力hj(即与该吸力对应的所有毛细管都具有相同孔径)和土壤含水率θs-jΔθ，即

hj=h(θs-jΔθ)

(11)

图2 计算非饱和导水率的毛细管模型Fig.2 Capillary model for calculating unsaturated hydraulic conductivity

另外，由于生物炭和土壤颗粒共同在一起，所以同一位置处吸力相同，即在hj吸力下，根据生物炭的水分特征曲线θW(h)得到生物炭中的含水率差值为

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

(12)

θ(h0)=θWS，θWS为生物炭的饱和含水率。根据式(8)，基质吸力hj对应的毛细管孔径为

(13)

对于土壤和生物炭，对应的j区的单位面积内毛细管数量分别为

(14)

(15)

式中nk——j区生物炭内单位面积的毛细管数量

ΔθWk——生物炭的含水率区间宽度

Rk——生物炭内毛细管半径

当h=hj时，所有孔径小于Rj的毛细管充满了水，提供了土壤的水流通量，并且这些毛细管是由土壤颗粒和生物炭颗粒共同提供的。另外，一般认为[23]生物炭的导水率远大于土壤的导水率，意味着生物炭中毛细管的密度大于土壤。如果假定生物炭中多于土粒中毛细管的倍数可以由饱和导水率的比值近似给出，再同时考虑式(10)，得到混合土壤的水流通量表达式为

(16)

其中

q=Kcs/Kss

式中q——生物炭与土壤饱和导水率的比值

Kcs——生物炭饱和导水率

Kss——土壤饱和导水率

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

(17)

式中hk——区间k内毛细管的吸力

N——生物炭水分特征曲线划分的最大区间数

(18)

其中

l=L/Le

式中l——毛细管弯曲度

式(18)与达西定律相比，得到混合土壤的非饱和导水率为

(19)

当毛细管全部充满水，就得到混合土壤的饱和导水率

(20)

结合式(19)和式(20)得到混合土壤的非饱和导水率

(21)

基于土壤和生物炭的水分特征曲线，再根据添加生物炭的比例以及土壤和生物炭的饱和导水率，可以用式(21)近似给出添加不同量生物炭后土壤的非饱和导水率。

3 结果与分析

3.1 土壤水分特征曲线

用离心机法测定的土壤和生物炭的水分特征曲线如图3所示。

图3 土壤与生物炭的水分特征曲线Fig.3 Water characteristic curves of soil and biochar

采用VAN- GENUCHTEN模型，利用Retc软件对离心机法的实测数据进行拟合，拟合的饱和含水率为0.418 1 cm3/cm3，残余含水率0.009 8 cm3/cm3。由图3可以看出，当土壤吸力小于30 cm时，生物炭与土壤几乎平行，当土壤吸力介于30～2 200 cm，土壤的SWCC位于生物炭上方,当吸力大于2 200 cm时，生物炭的SWCC位于土壤上方，即生物炭在中高吸力段的含水率大于土壤，因此添加生物炭后土壤持水能力的增加，主要依靠的是生物炭在中高吸力段持水能力的提升。也有文献表明，生物炭的添加减小了土壤容重、增加了土壤孔隙度[24]，使得土壤持水能力增强。

3.2 毛细管模型计算土壤的非饱和导水率

由于毛细管模型综合考虑了生物炭和土壤的性能，因此根据毛细管模型公式(21)以及测定的土壤和生物炭的水分特征曲线，可以计算不同添加量生物炭土壤的非饱和导水率。采用图3数据计算毛细管模型的数据如表1所示。计算中采用的生物炭和土壤的饱和导水率实际测量值见表2。由表1可以看出，[K(θs-jΔθ)/Ks]1/3项可以近似看作有效饱和度。在低吸力段，有效饱和度随着生物炭含量的增加从0.710 446增加到0.739 290，有效饱和度的增高可看作土壤导水性能得到提高[20]；在35.82～1 590 cm吸力段，有效饱和度随生物炭含量的增加而减小，土壤导水性能降低。导水性能的增加是由于生物炭在低吸力下的高导水率，而高吸力下导水性能的降低，是由于生物炭在高吸力下的低导水率。

3.3 土壤饱和导水率测量结果

表2为实际测量的CK、C1、C2、C3、C4的饱和导水率，以及采用式(5)和式(7)计算的估计值。

表1 采用毛细管模型计算不同添加量生物炭土壤的非饱和导水率Tab.1 Calculated unsaturated hydraulic conductivity by using capillary model

表2 土壤饱和导水率的测量值和估计值Tab.2 Tested and estimated values of soil saturatedhydraulic conductivity mm/h

由表2可以看出，C1、C2、C3、C4饱和导水率的实际测量值，均在层状分布和条状分布的估计值之间，但非常接近于层状分布，与层状分布相比的最大偏差为8.9%。分析其原因，可能是当土壤中的生物炭施加量较少时，不足以让生物炭粒子连接起来，所以应该看作层状。另外，试验数据表明[16]，添加生物炭土壤中的水分流动，多为水在生物炭和土粒之间的空隙流动。由于层状分布的估计值没有考虑生物炭和土粒之间的空隙，当生物炭添加后，空隙间的水分流动会对饱和导水率产生影响，因而导致实际测量的饱和导水率比层状估算值大。

4 添加生物炭土壤渗水过程的数值分析

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

为研究相同渗透时间后，不同生物炭添加量土样中的水分分布，首先要对添加生物炭的土壤进行数值模拟。由于实际土壤的复杂性，计算的模型拟采用二维形式。添加生物炭的土壤被分成土壤区和生物炭区。计算时，这两个区具有各自的非饱和导水率和初始水分特征。

由于生物炭是随机分布在土壤中，采用Python语言编程的方法给定土壤截面。将生物炭随机分布在土壤截面中的步骤如下：首先建立计算土壤的面积区域[25]。面积选择过大，生物炭颗粒的数量随之变大，导致计算量过大难以实现。因此，计算的面积选取10 mm×157 mm。其次确定土样中生物炭的个数，生物炭的个数通过添加生物炭的比例来校对。本研究中，添加生物炭的体积比约为2%、4%、6%、8%。根据添加比例和试验所用生物炭的平均粒径给出生物炭的个数(生物炭粒径范围1.5～2.5 mm，其中2 mm粒径约占60%以上，故取2 mm)。随机产生生物炭所在位置的圆心坐标，画圆代表生物炭。校验是否有圆相交，如果有，删除后再补充相应数量的圆。最后，对生成的所有生物炭颗粒求面积，核算添加的比例和实际掺炭比例是否相等。完成上述工作后赋予土壤区和生物炭区不同的水分运动参数。图4为最终生成的随机添加生物炭的土壤模型。

图4 随机分布生物炭颗粒的土壤截面 Fig.4 Profiles of biochar particles distributed randomly in soil

4.2 有限元采用的模型

采用通用的有限元分析软件Abaqus进行数值计算。 对于饱和土，渗透系数可以视为常量，但对于液体与气体并存的非饱和土，饱和度的大小直接影响渗透的结果，因此渗透系数K是饱和度S的函数。Abaqus采用BROOKS等[26]提出的渗透系数和饱和度的关系。利用土壤水分特征曲线给出的3个参数，即土壤的进气值、有效饱和度以及空隙尺寸分布指标λ，渗透系数函数与有效饱和度的关系为(下文中出现的饱和度均指有效饱和度)

(22)

其中

式中Sr——有效饱和度

δ——经验常数,取3

(ua-uw)——基质吸力

(ua-uw)b——进气值

4.3 实际土样与毛细管模型土样的数值计算

为研究不同土样经过相同渗透时间后的水分分布，数值模拟部分首先计算添加生物炭土壤的饱和度，其次采用前面推导的毛细管模型等效的非饱和导水率计算掺炭土壤的饱和度，初始给定的土壤吸力是相同的。

当区域为对照土样(CK)时，此时α=1，采用的吸力和对应的饱和度如表1所示，实际上这是采用单一土壤的毛细管模型给出土壤的非饱和导水率，饱和导水率Ks采用试验测得的值1.152 mm/h。边界条件为：左边界初始吸力为零表示进水，其他边界为不透水边界。全部区域的初始吸力和饱和度分别为200 cm、0.152 7，渗透时间5 000 s，计算出CK的饱和度随位置的分布规律如图5a所示。

当区域为C4时，区域如图4c所示，采用的吸力和饱和度如表1所示，土壤区为α=1，生物炭区为α=0。土壤和生物炭的饱和导水率分别为1.152、22.068 mm/h，边界条件仍为左边界初始吸力为零，土壤的初始吸力和饱和度为200 cm、0.152 7，生物炭的初始吸力和饱和度为为200 cm、0.060 7，渗透时间5 000 s，计算出的C4的饱和度随位置的分布规律如图5b所示。

当区域采用C4的毛细管模型C4*时，采用的吸力和饱和度如表1所示。饱和导水率为试验测得的1.371 6 mm/h，边界条件仍为左边界初始吸力为零，C4*初始吸力和饱和度为200 cm、0.086 5，渗透时间5 000 s，计算后C4*的饱和度随位置的分布规律如图5c所示。图5表明，毛细管模型很好地描述了施加生物炭后，土壤饱和度随位置的分布规律。

计算出的CK、C4中的土壤、C4中的生物炭以及C4*的饱和度随位置的分布规律如图6a所示。

图6 渗透5 000 s后的土壤饱和度曲线Fig.6 Saturation of soil after infiltration of 5 000 s

由图6a可以看到，C4*的饱和度分布更类似于C4中生物炭的饱和度分布形式，这是由于采用毛细管模型(式(21))的推导结果。由于生物炭中毛细管的密度远多于土壤，与同体积的土壤相比，能够吸收更多的水分，因此对土壤的饱和度影响更大。式(21)中的q值决定了C4*饱和度曲线的形状。图6a中，C4*的饱和度曲线在水平距离上略微超前于C4中生物炭的饱和度曲线，说明计算C4*饱和度时，采用测量的饱和导水率偏大。这可能是由于C4的计算中，没有考虑生物炭和土壤的空隙，而这是影响饱和导水率的重要因素。

设定渗透时间5 000 s，按上述方法计算C1、C2、C1*、C2*的饱和度分布规律。C1*、C2*饱和导水率为测得的1.195 2、1.249 2 mm/h。毛细管模型采用的土壤初始吸力均为200 cm，C1*和C2*的饱和度分别为0.114 4和0.099 9。图6b为渗透时间5 000 s后土样的饱和度随位置的分布曲线。C1*、C2*为毛细管模型模拟的土样(由于图线过于密集，C3没有画出)。

由图6可以看到，水分渗透一定时间后，C1*、C2*、C4*曲线均落在CK上方，即相同位置上添加生物炭土壤的饱和度都比CK大，这说明添加生物炭可以提高土壤的导水率。并且导水率随着生物炭添加量的增大而增大。C4*曲线落在CK下方，是由于高吸力下生物炭导水性能的降低造成的。毛细管模型综合考虑了生物炭和土壤的特性，能较好地反映添加生物炭后土壤的导水特性。

4.4 土柱的数值计算和试验数据的比较

本试验采用的土柱尺度较大，因此在有限元计算时，如果考虑实际的生物炭粒子，划分的网格必须很密，这样会导致计算量过大，所以一般无法对土柱直接计算，这时可以采用毛细管模型近似计算。计算时，由于C1*、C2*、C4*的饱和含水率和残余含水率各不相同，因此C1*、C2*、C4*的初始饱和度是不同的。为简化计算，给定C1*、C2*、C4*的初始饱和度，与试验测得CK的初始饱和度相同。

为验证毛细管模型的正确性，设定水分渗透时间为24 h，实际测量了C1、C2、C4在80 cm土柱中的含水率，采用CK拟合的饱和含水率和残余含水率值计算饱和度。计算的初始数据如下：实际测定的CK初始饱和度0.294，与之对应的土壤的吸力约为68.6 cm。毛细管模型在此饱和度下的吸力分别为51.9 cm(C1*)、43.4 cm(C2*)、40.4 cm(C4*)。图7为水分渗透24 h后，饱和度随土柱位置变化关系图。

图7 渗透24 h后土柱饱和度的计算值和测量值Fig.7 Calculated and measured values of soil saturation after infiltration of 24 h

由图7的计算结果看出，土柱中水分渗透24 h后，对于同一生物炭添加比例，毛细管模型的饱和度与实际测量值相比略低。这是由模型中没有考虑实际土壤中生物炭粒子和土壤之间的空隙造成的。生物炭体积比为2%土壤的误差范围为4.27%～7.24%，4%土壤的误差范围为3.06%～8.6%，8%土壤的误差范围为2.90%～10.11%。可见，随着生物炭体积比的增大，数值计算结果与毛细管模型的计算结果偏差有所增加，这种偏差是因为在毛细管模型假设中，生物炭中的不同孔径的毛细管数量是均匀分布的，总是用一个恒定的饱和导水率的比值近似给出，这与实际是有偏差的。由图7结果看出，在距供水端较远处，即土壤吸力较大处，土壤的实际导水率大于毛细管模型的值，说明生物炭中的细小毛细管的数量多于估计值。因此，采用毛细管模型的数值计算，对于中等及以下生物炭施用量，在土壤水力特征的定量分析上更加准确。毛细管模型可以简化土壤的水力特性，可以用来对实际的大块土地进行数值分析。

5 讨论

本研究针对东北黑土区的草甸黑土，推导了施加生物炭后土壤的毛细管模型，利用此模型，研究了一定渗透时间后土壤中的水分分布，选择饱和度等参数将添加不同量生物炭的土壤进行数值模拟，最后利用土柱的实际测定结果对数值模拟结果进行验证，以此开展生物炭对土壤水分影响的数值化定量化研究。

从土壤水分特征曲线看出，添加生物炭后土壤的持水能力增加，由于生物炭在中高吸力段的含水率大于参照组土壤，因此添加生物炭后土壤持水能力的增加，主要依靠的是生物炭在中高吸力段持水能力的提升。毛细管模型表明，施加生物炭后土壤性能的改善，是由于低吸力段生物炭的良好导水特性。利用毛细管模型计算的土壤非饱和导水率，依赖于土壤的水分特征曲线和饱和导水率，为得到精确的数值计算结果，必须精确测量土壤水分特征曲线。

土柱一端加水渗透24 h后，数值模拟的土壤水分分布值比实际测量值略低，主要是由于毛细管模型中生物炭多于土壤的毛细管数量，总是用一个恒定的值近似给出，这与实际有偏差。另外是由于模型中生物炭与土壤的空隙分布不能与实际情况完全吻合造成的。

虽然直接测定的结果更接近真实值，但测量过程繁琐、消耗时间长，因此采用毛细管模型的有限元数值计算优势明显。另外，对于只能通过试验研究解决的问题，如一定范围的生物炭施加量、不同的施加面积、地势情况(如坡耕地等)、不同的渗透时间等情况，数值模拟能够更好更快地给出近似答案。采用毛细管模型还可以对试验数据进行校验，找出试验中有偏差的数据，有助于提高试验结果的准确性，从而为研究土壤的水力特性提供更多的方法与保障。

6 结论

(1)试验测得的水分特征曲线表明，当土壤吸力大于2 200 cm时，生物炭的残余含水率高于对照组，因此添加生物炭后土壤持水能力的增加，主要依靠的是生物炭在中高吸力段持水能力的提升。

(2)通过饱和导水率试验数据和理论估计值的比较可以看出，实际测量的饱和导水率与生物炭层状分布的饱和导水率非常接近，最大偏差为8.9%，说明添加生物炭土壤更应考虑为层状分布。试验数据和理论数据的差别可能是生物炭和土壤之间的空隙增加造成的。

(3)毛细管模型结果显示，添加生物炭土壤导水率的提高，主要是由生物炭的导水特性决定的，添加的生物炭主要改变了土壤在低吸力下的导水率。

(4)数值分析显示，施用生物炭能够提高东北草甸黑土的导水率。添加生物炭的比例越高，导水率增加越大。

(5)生物炭中的不同孔径的毛细管数量不是均匀分布的，孔径越小，毛细管的数量越多。

(6)毛细管模型简单有效，生物炭体积比小于8%时，数值分析的结果与试验数据吻合良好，说明该模型可以为土壤水分运动参数的数值计算提供比较准确的方法。

1 何绪生, 耿增超, 佘雕, 等. 生物炭生产与农用的意义及国内外动态[J]. 农业工程学报, 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 魏永霞, 刘志凯, 冯鼎锐, 等. 生物炭对草甸黑土物理性质及雨后水分动态变化的影响[J/OL]. 农业机械学报, 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 勾芒芒, 屈忠义, 杨晓, 等. 生物炭对砂壤土节水保肥及番茄产量的影响研究[J/OL]. 农业机械学报, 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]. 河海大学学报, 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.