沼液穴灌入渗影响因素分析及数值模拟研究

郑 健，周鹏杰，张彦宁，殷李高，张恩继

(1.兰州理工大学西部能源与环境研究中心，甘肃 兰州 730050；2.甘肃省生物质能与太阳能互补供能系统重点实验室，甘肃 兰州 730050；3.兰州理工大学能源与动力工程学院，甘肃 兰州 730050；4.陕西省土地工程建设集团，陕西 西安 710075)

沼液是禽畜粪便等有机物经过厌氧发酵后的残留液体，其总固体含量小于1%[1]，含有水溶性及多种养分的速效肥料，保留了90%以上的发酵原料中的氮、磷和钾，且氮素结构较原污水更优，施用可以显著改善土壤的理化性质，提高土壤的保水保肥能力，缓解土壤的板结酸化[2]。与清水和常规施肥相比，施用沼液会降低土壤细菌值，增加微生物数量并提高土壤酶活性[3]，同时还能促进作物生长，提高农产品产量和品质[4-6]，降低病虫害发生率[7]，沼液浸种还可提高种子发芽率和生长水平[8-10]。目前关于沼液施用的研究主要集中在作物产量、品质及土壤肥力等方面。

沼液穴灌技术是从设施农业生产中总结出的一种节肥技术，目前该技术的实施主要根据经验，与其相关的试验研究和理论分析尚少。土壤水分运动的分析是研究农田灌溉技术的基础，准确掌握穴灌的土壤水分运动变化规律对于促进穴灌技术的发展具有极其重要的意义。国内外众多学者对多种灌溉方式下的土壤水分分布问题做了研究，但从文献内容来看，多集中在滴灌[11-12]、沟灌[13]、膜孔灌[14-15]等方面。而将沼液与穴灌技术结合，关于沼液灌溉入渗规律和灌溉技术参数等方面的研究鲜有报道。

基于此，本文通过设置不同沼液配比、穴孔直径和土壤容重，研究这些因素对沼液穴灌入渗的影响，并运用COMSOL Multiphysics软件模拟分析沼液穴灌入渗湿润锋运移及土壤含水率分布规律，以期为水—沼液一体化灌溉技术提供依据，使沼液穴灌技术建立在更为科学合理的基础上。

1 材料与方法

1.1 试验材料及装置

试验土壤选取甘肃省兰州市周边正常使用的温室大棚内表层0～40 cm的土壤，土壤经过自然风干、粉碎后，过2 mm孔径土壤筛备用，土壤颗粒级配见表1。按照中国土壤分类标准，其土壤类型为粉质壤土，田间持水量为22.23%(质量含水率)。沼液取自发酵物总固体质量分数(TS)为8%，以牛粪和麦秸为主要原料，正常使用的户用沼气池，静置2个月，待其理化性质稳定后，用纱布过滤掉沼液中较大的颗粒备用。

试验装置主要由土箱和马氏瓶组成，考虑到穴孔入渗断面的对称性，采用实际穴灌土壤湿润体的1/12作为研究对象[16]，土箱为30°扇形有机玻璃装置，高和宽分别为50 cm和40 cm，马氏瓶用来提供恒定水头，容量520 mL。试验装置如图1所示。

1.2 试验设计及方法

试验设置土壤容重为1.35 g·cm-3和1.45 g·cm-3；依据前期预试验结果设置3个沼液与水的配比(体积比)，即1∶4、1∶6、1∶8和1个纯水处理(对照)，3个穴孔直径(3 cm、5 cm和7 cm)；设定孔深为5 cm。每组试验做3次，共计72组试验，取平均值作为试验结果，试验方案见表2。

采用室内土箱试验，土体采用分层填装(每层5 cm)，以设定土壤容重计算每层装土质量，采用精度0.01 g电子秤称量；最后5 cm填土时，在土箱30°角处预留试验设计尺寸的穴孔；将按比例配比好的沼液装入马氏瓶，用乳胶软管分别将马氏瓶出水孔和排气孔与装完土的试验土箱连接，为保证恒定水头，将排气管固定在入渗孔与土体表面平齐的位置。为方便记录试验数据，在土箱外壁绘制间距5 cm的网格线。然后同时打开马氏瓶的出水孔和排气孔阀门，当沼液进入穴孔时开始计时，随着试验的进行，当湿润锋运移到与外壁水平向(H)和垂向(V)的网格线相切时，记录一次入渗时间和马氏瓶的读数，并用记号笔在土箱外壁绘制湿润锋运移曲线。当入渗量达到520 ml时或入渗时间约12 h时同时关闭马氏瓶的进水孔和排气孔阀门，待穴孔内沼液渗完，试验结束。为避免土壤水分再分布造成土壤含水率的变化，在试验结束后立即沿着土槽中轴线水平方向和垂直方向每3 cm的交叉位置用土钻取土一次，土样编号并放入烘箱内烘干(105℃，8 h)后测定土壤质量含水率。

表1 土壤颗粒级配

1.马氏瓶； 2.进水管； 3.排气管； 4.穴孔；5.土槽； 6.湿润锋； 7.取样点； 8.土体1. Markov bottle; 2. Inlet pipe; 3. Venting pipe; 4. The hole;5. Trough; 6. Wetting front; 7. Sampling point; 8. Soil.图1 试验装置及取样示意图Fig.1 Schematic diagram of the test equipment and sampling

处理Treatment土壤容重/( g·cm-3)Soil bulk density沼液配比Biogas slurry:water ratio穴孔直径/cmHole diameter处理Treatment土壤容重/( g·cm-3)Soil bulk density沼液配比Biogas slurry:water ratio穴孔直径/cmHole diameterT11.351∶43T131.451∶43T21.351∶45T141.451∶45T31.351∶47T151.451∶47T41.351∶63T161.451∶63T51.351∶65T171.451∶65T61.351∶67T181.451∶67T71.351∶83T191.451∶83T81.351∶85T201.451∶85T91.351∶87T211.451∶87T101.35纯水 Water only3T221.45纯水 Water only3T111.35纯水 Water only5T231.45纯水 Water only5T121.35纯水 Water only7T241.45纯水 Water only7

1.3 基于COMSOL软件模型建立

本文采用多物理场耦合软件COMSOL Multiphysics进行水、沼液渗流模拟分析。COMSOL以有限元为基础，通过求解偏微分方程(单场)或者偏微分方程组(多场)来实现真实物理现象的仿真[17]，具有大规模的计算能力。求解器是基于C++程序采用最新的计算技术编写，它的求解器类型有直接、迭代、分离、多重格点；求解类型有稳态、瞬态、参数等；处理器支持多核与集群。根据需要，用户可以通过MATLAB GUI界面选择单场求解或者多物理场求解。内嵌丰富的CAD建模工具的有限元软件COMSOL，可以直接进行二维以及三维建模。COMSOL在网格剖分方面具有强大的功能，它的网格生成器可以将腔体划分为三角形或者四面体的网格单元，通过仿真系统的交互界面用户可以根据自己的需要选择剖分的精细程度。

图2 直径3 cm穴孔模型及网格划分图Fig.2 Cavity model and mesh map of 3 cm hole

穴孔灌的土壤水分运动属于三维点源入渗条件下的土壤水分运动，具有轴对称性。在此次模拟分析中，假定土壤均匀且各向同性，同时将穴孔灌的水分运动简化为轴对称条件下的二维土壤入渗问题进行处理，本次模拟将沼液与水配比后的液体视为混杂了颗粒物的浑水，模拟对象统一定为边长40 cm的正方形轴对称模型，在对称轴顶点处留有不同尺寸的穴孔，见图2。

土壤水分运动控制方程可表达为[18]:

(1)

式中，θ为土壤体积含水率；h为土壤负压水头，cm；Kh为土壤非饱和导水率，cm·min-1；t为时间，min；x为径向坐标，cm，x向右为正；z为纵向坐标，cm。

1.4 模型参数

根据COMSOL模型模拟需要，对土壤水分特征曲线参数进行拟合，土壤水分特性曲线θh可用van Genuchten模型来表示，即

(2)

式中，θr为残余土壤体积含水率；θs为饱和土壤体积含水率；α、n、m为土壤水分特征曲线的形状参数，m=1-1/n，n>1。

根据表1土壤级配表，采用ROSETTA软件拟合了两种容重土壤的水力参数，结果见表3，其中Ks为土壤饱和导水率。

沼液中含有有机悬浮颗粒，且粘性较水高，在入渗过程中会影响沼液入渗[18]。因此，本文在模拟中沼液入渗的饱和导水率采用变量设置。为保证变量设置的合理性，本文采用课题组实测数据拟合后获得的函数[19]。函数为幂函数的形式：

Kst=Ks1·t-α

(3)

式中，Kst为t时刻的入渗速率(cm·min-1)；t为入渗历时(min)；Ks1为第一时间单位末的入渗速率，根据试验测得；α为经验入渗指数，通过Origin软件对试验数据进行拟合获得。

1.5 定解条件的确定

1.5.1 初始条件 假定土壤初始含水量在所研究的区域内分布均匀，初始条件可表示为:

hx,z,t=h00≥z≥Z,0≤x≤X，t=0

(4)

式中，h0为土壤初始负压水头(cm)；Z、X分别指研究区域垂直方向、水平方向最大距离(cm)。

1.5.2 边界条件 参照图3可以看出：上边界AB，试验过程中不考虑蒸发，水流通量为零，则

(5)

表3 土壤水力参数

图3 穴孔入渗模型Fig.3 Hole infiltration model

右边界BC和下边界DC，由于所选研究区域较大，可认为水分无法到达该边界，边界条件描述为：

右边界BC，

hX,z,t=h0t>0

(6)

下边界DC，

hx,Z,t=h0t>0

(7)

边界OD为对称轴,水流通量为零，

(8)

边界ab、Ab为穴孔中入渗面，与入渗液体接触，边界ab，

hx,z,t=5cm (z=-5 cm)

(9)

边界Ab，

hx0,z,t=hi

(10)

式中，x0为穴孔半径，hi为作用于孔壁各点的水头，cm，其数值等于孔壁各点到水面的距离。

2 结果及分析

2.1 孔径对入渗的影响

图4中黑色实线表示湿润锋的实际运移情况，不同颜色色块代表湿润锋不同时段运移的模拟结果。从T7、T8、T9处理的湿润锋运移情况可以看出，在垂向入渗距离达到15 cm时，所需入渗时间分别为214.8、130.4、83.2 min，T8相对于T7、T9相对于T8的入渗耗时衰减率分别为39.3%和36.2%，即孔径越大，水分入渗速率越快，且在相同时间内，孔径越大，湿润体的体积也就越大。这主要是由于水分在土壤入渗过程中初始阶段主要受到土壤基质势作用，在入渗的后期，土壤基质势作用逐渐减弱，而重力势作用逐渐加强，随穴孔直径的增大，穴孔内的沼液体积增加，导致重力场的变化，因而对入渗过程的作用愈加显著，使入渗速率随着穴孔直径的增大而增大。为进一步分析孔径对入渗的影响，本文以实际垂向入渗距离达到15 cm时的湿润锋为研究对象，分析实际水平入渗距离与垂直入渗距离的比值变化规律，结果表明T7、T8和T9处理的比值分别为1.007、0.952和0.945，T10、T11和T12处理的比值分别为1.050、0.950和0.902，T19、T20和T21处理的比值分别为1.061、0.938和0.948，T22、T23和T24处理的比值分别为1.127、0.959和0.886，由比值变化可以看出孔径越大，比值越小，即湿润体纵向长度随孔径增大而逐渐大于横向长度(穴孔尺寸不计)。初步分析产生的原因是：随着孔径的增大，穴孔中入渗面面积也随之变大，在孔深不变的情况下，即孔壁入渗面面积不变，孔底入渗面面积变大，相应的过水通道变多，使垂向入渗距离逐渐大于水平入渗距离。

图4 不同处理湿润锋模拟值与实测值对比图Fig.4 Comparison of simulated and measured values of wetting front under different treatments

2.2 土壤容重对土壤水分入渗的影响分析

为说明土壤容重对入渗的影响，选取T7与T19处理、T8与T20处理、T9与T21处理进行分析，从图4中可以看出，T7与T19处理、T8与T20处理、T9与T21处理实际垂向入渗距离达到15 cm处所用时间分别为214.8 min与335.3 min、130.4 min与154.8 min、83.2 min与94.7 min，可以看出土壤容重1.45 g·cm-3的入渗速率明显低于1.35 g·cm-3，即土壤容重大的入渗速率慢，这一试验现象与很多学者的研究成果相符合[20-23]。张振华等[24]的研究指出，在入渗过程中，土壤容重对入渗率的影响具有稳定性和一致性，表现为对应于同一时刻，容重越大入渗率越低。分析形成的主要原因是随着土壤容重的增大，土壤的团粒结构丧失，土壤变得紧密，从而土体内的大孔隙比例降低，导致土壤的基质吸附作用和毛管作用降低，使得入渗速率也随之而降低。图5中(a)、(b)、(c)和(d)分别表示1∶4(沼液∶水)、1∶6、1∶8和纯水CK对照条件下相同孔径在不同土壤容重条件下的累积入渗量对比，从图5中可以看出，累积入渗量在相同时间内土壤容重呈1.45 g·cm-3<1.35 g·cm-3。例如图5(c)中累积入渗量呈T7>T19处理，T8>T20处理，T9>T21处理。

2.3 沼液配比对入渗的影响

入渗液体浓度的差异会对土壤水分入渗产生影响[25]。图4中T14、T17和T20处理表示孔径5 cm，土壤容重1.45 g·cm-3，沼液水配比分别是1∶4、1∶6和1∶8的处理，从图4湿润锋运移规律可以看出，当垂向入渗距离达到15 cm时，分别耗时244.1、210.0、154.8 min，即入渗速率随着沼液配比度的升高而减慢。从图6累积入渗量的变化规律可以看出，随着水沼液混合液中水体积的升高，累积入渗量逐渐增大，并且混合液的入渗量衰减速度明显快于纯水，初步分析产生上述规律的主要原因为沼液溶液中微小的悬浮颗粒含量随着沼液配比的增加而增多，在入渗过程中，悬浮颗粒堵塞了土壤中的孔隙[26]，降低了土壤的基质吸附作用和毛管作用，从而降低了沼液的入渗速率。

2.4 入渗模拟模型有效性验证

2.4.1 土壤含水率模拟结果验证 对24个处理进行模拟后共得到466个观测点质量含水率数据，与实测值对比分析后，相对误差在0～5%之间的点有290个，占62.2%，5%～10%之间的点有162个，占34.8%，10%～15%之间的有9个，占1.9%，大于15%的有5个，占1.1%。图7和图8分别显示了相同土壤容重在距离OD轴6 cm和9 cm处取样点含水率实测值与模拟值对比，从图中可以看出模拟与实测值吻合程度较好，纯水处理的模拟结果精确性要高于沼液；整个模拟结果中，误差较大处大多出现在第一个和最后一个取样点。实测值与模拟值的对比分析结果表明所建立的模型能较好反映穴孔自由入渗的土壤含水率分布情况。

图5 不同处理累积入渗量曲线Fig.5 Cumulative infiltration curve of different treatments

2.4.2 湿润锋模拟结果验证 图4中，实线表示不同时刻实际湿润锋位置，不同颜色区域表示不同时刻模拟湿润锋位置，从图中可以看出：(1)各时刻湿润体形状的模拟结果与实测结果一致性较好，湿润锋的模拟值与实测值可以较好的吻合；(2)由于沼液中颗粒物的存在，入渗速度衰减较快，模拟精度比纯水要低；(3)在入渗后期实际湿润锋运移速度快于模拟湿润锋运移速度。

2.4.3 影响模拟的因素分析 模拟结果表明，土壤含水率和湿润锋的模拟与实测结果存在一定误差，初步分析是受到土壤条件、沼液配比、试验装置和边界条件设置的影响。主要体现在：(1)实际试验中土壤采用分层装填，可能存在局部的不均匀性；(2)模拟中对上边界条件的简化以及采用的土壤水分运动控制方程时忽略了温度与土壤气体对水分运动的影响，同时试验土箱的侧壁为光滑的有机玻璃，与土体接触性较差，试验过程中可能存在一定程度的优先流也造成了实测数值与模拟值之间的差异；(3)湿润锋边界附近含水率变化较大，用土钻取样较容易产生误差，并且湿润锋附近土体含水率值较小，测量容易产生偏差；(4)水和沼液的混合液在实际入渗过程中，土壤水分运动参数是动态变化的，本文采用实测数据的拟合方程代入，对模拟精度也造成了一定的影响。

图6 不同沼液配比累积入渗量Fig.6 Relationship of accumulated infiltration volume withinfiltration time under different biogas slurry ratios

注：M表示实测值，S表示模拟值，下同。Note: M represents the measured value, S represents the simulated value, the same below.图7 1.35 g·cm-3含水率模拟值与实测值对比图Fig.7 Comparison of the simulated values of soil water content with the measured valuesunder soil bulk density of 1.35 g·cm-3

图8 1.45 g·cm-3含水率模拟值与实测值对比图Fig.8 Comparison of the simulated values of soil water content with the measured valuesunder soil bulk density of 1.45 g·cm-3

3 结 论

1)水-沼液一体穴孔自由入渗能力受穴孔孔径、土壤容重和沼液配比的影响较为显著。累积入渗量和土壤湿润体随穴孔直径的增加而增大，随沼液配比的增加而减小；水平入渗距离与垂向入渗距离的比值随穴孔孔径增大而减小，湿润体形状随穴孔直径的增加逐渐由水平轴大于垂直轴的椭球体趋向于水平轴小于垂直轴；1.45 g·cm-3土壤容重条件下的入渗速率和累积入渗量明显低于1.35 g·cm-3土壤容重，表明随土壤容重增大，相同时间内入渗速率、累积入渗量和湿润体面积呈减小趋势；随沼液配比的增大，入渗速率和累积入渗量均随之减小。在实际操作过程中，当土壤容重较大时可通过增加穴孔直径和降低沼液水配比来增加土壤入渗量。

2)采用COMSOL Multiphysics建立穴孔入渗模型，并对室内试验条件下的穴孔入渗试验进行了数值模拟研究。各处理条件下湿润锋运移距离的实测值与模拟值吻合较好，均表现出初期运移速度较快，随试验时间的延长渐渐减缓，并趋于稳定。土壤含水率的实测值与模拟值对比来看，相对误差<10%的占97%，说明所建立的模型合理，利用COMSOL Multiphysics可以准确描述水-沼液一体穴孔自由入渗过程中的水分运动。