基于改进的HYDRUS-2D模型暗管排水过程模拟

刘 义,曾文治，敖 畅,雷国庆,伍靖伟,黄介生

(武汉大学水资源与水电工程科学国家重点实验室，湖北 武汉 430072)

土壤盐渍化是一个危害全球农业环境的难题，它不仅能够抑制作物的生长，降低作物的产量和品质，同时还能降低肥料的利用率，从而增加土壤养分流失的风险[1-2]。暗管排水通常被认为是治理土壤盐渍化的重要水利工程措施[3-4]。暗管的埋深、间距和外包材料等是影响暗管排水排盐效果的重要因素，国内外学者对此展开了一系列的研究[5-7]。张展羽等[8]在滨海地区模拟了暗管在不同埋深和不同间距下的排水排盐效果，确定了适宜在滨海地区埋设暗管的技术参数。Bonaiti等[9]对意大利东北部的地下灌溉和暗管排水的效果进行了研究，表明当管道之间的距离为8 m时，暗管排水效果最佳。王洪义等[10]在地下水埋深较浅的大庆盐碱地开展了不同埋深下的暗管排水降渍试验，结果表明暗管埋深0.8 m是当地暗管最佳工程布设参数。王振华等[11]研究了滴灌淋洗条件下的暗管排盐效果，发现当地下水埋深较大时，暗管间距为15 m排盐效果最佳。Mastrocicco等[12]的研究中也得到了类似的研究结果。Wang等[13]得出浅层暗管排水能够有效地降低耕作层土壤的含盐量，从而提高作物的产量。钱颖志等[14]研究了新疆地区田间的暗管排水试验，得出暗管上方土壤的脱盐率随着埋深的增加而增加。Tao等[15]通过在暗管周围添加具有较高渗透性能的过滤材料改进了暗管的排水能力。Filipovic等[16]研究证明了在暗管周围包裹砾石能够显著提高暗管的排水效果，主要表现在通过暗管排水能够显著减少地表积水和地表径流。然而，以上的这些试验通常需要耗费大量的人力物力并且周期一般较长，其结果通常是针对当地的土壤环境，所以缺乏普遍适用性。此外，由于涉及暗管布局的参数较多，很难对这些影响因素进行系统定量试验。因此，有必要通过数值模型快速、系统地研究暗管排水下的水盐运移过程，从而优化暗管的布局。

目前，考虑暗管布局的常用模型有HYDRUS-2D、DRAINMOD、SWAT等。相比之下，HYDRUS-2D具有灵活的边界条件可以较为准确地模拟二维非饱和土中的水盐运移过程。因此，许多研究将HYDRUS-2D应用于暗管排水的模拟。石培君等[17]运用HYDRUS-2D模拟了膜下滴灌条件下的暗管排水排盐效果，发现模拟结果与试验结果吻合较好。Castanheira等[18]通过HYDRUS-2D模拟暗管排水试验发现，利用暗管排水可以有效降低土壤盐碱化的风险。Adobea等[19]使用HYDRUS-2D通过对田间土壤含水率变化的模拟，确定了一年中可用于田间作业的工作天数。李显溦等[20-21]通过田间试验，基于土壤剖面水盐的率定结果测试了HYDRUS-2D在模拟暗管排水中的适用性。尽管上述研究中模型的模拟精度是可以接受的，但是他们往往仅从土壤剖面水盐含量的角度来率定模型，而忽略了排水和排盐过程，这可能会导致HYDRUS-2D在模拟暗管排水方面具有一定的误差。此外，以往的研究中暗管边界主要采用的是渗透边界，然而渗透边界的排水能力远大于实际暗管的排水能力，因此在HYDRUS-2D中仅利用渗透边界很难模拟暗管的排水、盐过程。所以有必要为暗管开发新的边界条件以提高HYDRUS-2D的模拟精度。

本研究旨在改进HYDRUS-2D模型在模拟暗管排水过程中的性能。为此，本文提出了虚拟土层(VSL)和实际开孔面积(AHA)两种方案来代替HYDRUS-2D暗管的原有渗透边界(PSB),并且利用试验数据对以上暗管边界在HYDRUS-2D的模拟精度进行了评价。

1 材料与方法

1.1 试验设计

试验在武汉大学水资源与水电工程科学国家重点试验室进行。试验设备是由钢板焊接而成的淋洗槽，其长度、宽度和高度分别为150、51 cm和95 cm(图 1)。为了控制水深并确保土壤在大水漫灌中得到持续的淋洗，试验前在淋洗槽的前面和两个侧面距顶部3.5 cm 处钻了6个溢流孔(φ=7.5 mm)。暗管埋设在淋洗槽的中心距顶部65 cm的位置，开孔率为3%。坡度为2‰，管道的外径为54 mm，内径为50 mm。在管道外包裹一层可渗透的无纺布，用来防止土壤颗粒堵塞暗管的排水孔。

试验所用土壤样品采集于中国新疆阿拉尔(40°33′58″N，79°49′32″E)，该地区近年来主要采用暗管排水作为土壤盐渍化的防治技术。土壤粒径分布通过激光粒度分析仪(Mastersizer 2000，Malvern仪器有限公司，英国)测得；土壤容重通过环刀法采集原状土壤测得(表1)。试验前将土壤自然风干，并用2 mm土壤筛除去杂质进行预处理。在室内试验中为了模拟田间土壤的层状分布，将填充至淋洗槽中的土壤分成四层。在填土过程中每隔5 cm将土壤压实1次，保证每次的土壤压实度相同。试验共分两次进行，在两次试验中分别向该淋洗槽中注入200 mm和95 mm的水。通过自记式流量计测得暗管的排水量，并且每30 min收集排水样品以测得电导率。在每次试验结束后用土钻采集10、15、20、25、30、35、40、45、50、55、60、65、70 cm和75 cm深度处的土壤样品，用来测量土壤剖面中水和盐的分布。土壤体积含水率用烘箱干燥法测得(105℃， 8 h)，土壤电导率用电导率仪(DDS-307，上海雷磁仪器厂)测得。

1.2 模型的基本原理

采用HYDRUS-2D模型模拟暗管排水条件下土壤水盐的动态变化。HYDRUS-2D可以基于有限元流动方程的数值解来模拟二维或三维轴对称流动。在前人的研究中HYDRUS-2D常被用于模拟土壤水分和溶质运移，并且模型的模拟结果在许多情况下均表现良好[22-23]。

在HYDRUS-2D模型中，假定土壤为具有各向异性结构的多孔介质，并且忽略空气和温度对土壤水和盐运移的影响。用二维饱和-非饱和Richards方程来描述HYDRUS-2D中的土壤水分运动:

(1)

表1 土壤的物理性质

图1 室内试验装置示意图Fig.1 Schematic diagram of laboratory experiment device

式中，θ为土壤体积含水率(cm3·cm-3)；K(θ)为非饱和土壤导水率(cm·h-1)；t为从开始入渗到排水结束的时间(h)；z为土壤深度(cm)。

模型采用土壤水力参数vanGenuchten模型来描述土壤水力特性：

(2)

(3)

模型中使用二维对流弥散方程来描述土壤盐分的运移:

(4)

式中，C为溶质浓度(g·cm-3)；t为从开始入渗到排水结束的时间(h)；qi为在xi方向的达西流速(cm·h-1)；Dij饱和-非饱和水动力弥散系数(cm2·h-1)，i,j分别表示x,z方向；xi,xj为空间坐标。

1.3 边界设置和初始条件

为了改进HYDRUS-2D在暗管排水中的应用，考虑了3种暗管边界的概化方式(图2)。第一种方案，设置暗管边界为原始的渗透边界(PSB)作为其它两种方案的对照组。第二种方案，在暗管周围1cm范围内添加虚拟土层(VSL)，其作用是模拟无纺布和波纹管的结合界面。第三种方案，考虑到暗道的实际开孔率非常小，将暗管边界设置为无通量和自由渗透相交替的边界，用来模拟暗管的实际开孔面积(AHA)。因为试验中使用的淋洗槽的周边和底部是不透水的，所以将模拟区域的左右边界和下边界设置为无通量边界。在本试验中，由于上边界的水头随时间变化是逐渐减小的，所以根据实测水头数据，上边界可简化为变水头边界，当上边界无积水存在时则变为大气边界。

模拟区域内各层土壤的初始含水率和含盐量在填槽前进行测量，假设土壤含水率和含盐量在土壤剖面中呈线性分布，具体见表2。

1.4 模型的率定与验证

本次数值模拟以“小时(h)”为时间单位，模拟周期分别为90 h和60 h。其中第一次试验用于模型的率定，第二次试验用于模型的验证。输入的土壤参数由HYDRUS-2D内嵌的人工神经网络计算工具(Rosetta)根据测量的每层土壤的容重和粒径分布来计算。在参数的率定方面，首先对暗管的排水和排盐过程进行率定，然后对土壤剖面的最终含盐量和含水率进行率定。在模拟过程中，首先提取暗管5、10、15、20、25、30、35、40、45、50、55、60、65、70 h和75 h的累积排水量，以及10、20、30、40、50、60 h和70 h的累积排盐量来率定暗管的排水和排盐过程。然后，在10、15、20、25、30、35、40、45、50、55、60、65、70 cm和75 cm的深度处设置观察点，并提取观测点土壤剖面的最终含盐量和含水率。最后，通过PEST的多目标校正法将实测值与模拟值进行对比。当模拟值与实测值足够接近时，则选取率定的参数值作为模型最终输入参数值。上述过程由FORTRAN编程完成。率定后的土壤水力参数如表3所示。在参数率定之后，将第二个试验用于模型的验证。

1.5 评价指标

均方根误差(RMSE)和决定系数(R2)这两个指标用于评估模型的模拟精度：

(5)

(6)

表2 土壤含水率和含盐量的初始值

表3 三种边界率定的土壤水力参数

2 结果与分析

2.1 模型模拟结果

3种暗管边界的模拟结果如图3和图4所示。对于VSL和AHA来说，在模拟初期，暗管出流量的模拟值均小于暗管出流量的实测值，而在模拟后期，暗管出流量的模拟值则略大于实测值。但是，总体来看两种模型的排水历程与实测排水过程较为一致。然而，在PSB中则出现了相反的模拟结果。尽管PSB的最终排水量与实测值一致，但是排水过程的时长却远小于实测排水过程。对比模拟值与实测值，VSL的均方根误差RMSE为9.48 L，小于AHA和PSB中相对应的值，并且VSL的决定系数(R2)为0.930，在3个模型中最高，因此，VSL模拟排水过程时精度最高。排盐与排水过程相类似(图3b和图4b)，PSB的排盐过程均比其它2个模型结束得早。此外，在3个模型中，VSL在排盐过程的率定和验证阶段的RMSE最小而且R2最高，这表明VSL在对暗管排盐过程的模拟中表现最佳。

图3c、3d和图4c、4d显示了3种边界对土壤剖面的最终含盐量和土壤含水率的模拟结果。在模型的率定与验证的过程中，模型均略微低估了0～40 cm

图3 PSB、VSL和AHA在模型率定阶段的表现Fig.3 The performance of PSB, VSL and AHA during model calibration stage

图4 VSL和AHA在模型验证阶段的表现Fig.4 The performance of VSL and AHA during model validation stage

深度处的土壤脱盐率，但40～80 cm深度处土壤剖面含盐量的模拟值与测量值吻合较好。3种模型的含盐量模拟精度的RMSE和R2范围分别为0.757～1.075 g·kg-1和0.964～0.995，这表明模型在模拟土壤的盐分分布时具有较高的精度。对于土壤含水量，3种模型的RMSE均很小，说明3种模型对土壤剖面水分分布的模拟精度较高。总的来说，与PSB和AHA相比，VSL在模拟排水过程、排盐过程和土壤剖面水盐方面表现最好。

2.2 虚拟土层相关参数的确定

由前述所得，VSL是改善HYDRUS-2D用于模拟暗管排水的最佳方案。其中，VSL的饱和渗透系数(Ks)和虚拟土层的圆环宽度(RW)是影响模型模拟精度的两个关键参数。为此，我们分析了暗管排水过程、排盐过程、土壤剖面含盐率和含水率对参数Ks和RW的敏感性。当参数Ks和RW变化50%分别导致暗管排水过程和排盐过程变化32.04%～35.82%，24.63%～27.91%；而土壤剖面盐分分布和水分分布变化范围为0.11%～8.93%。这些结果表明暗管排水过程和排盐过程对参数Ks和RW有较高的敏感性，土壤剖面盐分和水分对这两个参数的敏感性较低。可见，Ks和RW的变化对排水、排盐过程影响最大。为了定量表征它们之间的关系，将虚拟土层环宽(RW)分别设置为0.5、1、1.5、2 cm和2.5 cm，然后用PEST率定相应的Ks。发现在不同的RW情景下，通过调整Ks，该模型能够很好地模拟土壤剖面的盐分和水分以及暗管的排水和排盐过程。将所得的RW与其对应的Ks进行线性拟合(图5)，拟合效果较好，关系如下:

Ks=0.019RW+0.0058R2=0.99

(7)

式中，Ks为虚拟土层的饱和渗透系数(cm·h-1);RW为虚拟土层的宽度(cm)。

该拟合公式的R2为0.99，接近于1，可见该拟合公式能较好地反映虚拟土层的饱和渗透系数与其圆环宽度之间的定量关系。

图5 虚拟土层的饱和导水率(Ks)与环宽(RW)之间的关系Fig.5 Relationship between saturated hydraulic conductivity(Ks) and ring width (RW) of virtual soil layer

3 讨 论

在数值模拟中，边界设置是建模的关键，它关系到模型的模拟精度[24]。目前，HYDRUS-2D中还没有专门用于反映暗管排水性能的边界。在以往的研究中，研究者们主要应用渗透边界(PSB)来模拟暗管排水[17, 20, 25]。然而，当周围土壤饱和时，在重力水的作用下，渗透边界会立即排水。渗透边界的排水面积大于试验中所用暗管的开孔面积，所以PSB的排水和排盐过程可以较快完成，这导致模型模拟的排水过程比实测过程要短。因此，在HYDRUS-2D中仅应用渗透边界只能模拟暗管排水排盐过程的总趋势，不能模拟暗管的整个排水过程的动态变化量，这与李显溦等[20]的研究结果相近。考虑暗管的实际排水面积、物理结构和开孔率较小(3%)等特性，提出了2种方式(VSL和AHA)来概化暗管的边界，以改进HYDRUS-2D在模拟暗管排水时的精度。

AHA与VSL均能准确模拟土壤剖面中的水分和盐分分布(图3c、3d和 图4c、4d)。然而在模拟排水量与排盐量方面，AHA和VSL的最终模拟值均高于实际值。这是由于试验中所用各层土壤的初始含水率和含盐量可能存在分布不均的现象，而模拟中采用的土壤初始含水率和含盐量在各层中是均匀分布的，导致了土壤初始水量与盐量的模拟值大于实际值[26-27]。从而在相同灌水量的条件下，两种模型的最终排水量与排盐量大于实测值。在本研究中，为了表征暗管的实际开孔率和实际排水面积，我们设置了AHA方案。与实测值相比，AHA的排水和排盐过程具有一定的滞后性。其中的原因是，与PSB相比,尽管AHA的排水面积与暗管的实际排水面积相同，但是AHA所应用的部分渗透边界的导水率仍与PSB相同且大于实际情况，所以较高的导水率仍会导致暗管的排水速率大于实测值[26]。为了降低AHA的排水速率，上层土壤的Ks被调整为0.12 cm·h-1(表3)，这比下层土壤(1.07～2.07 cm·h-1)小得多。由于饱和导水率较小，上层土壤的透水量减少，下层土壤的饱和速率也随之降低[28]。暗管周围土层的饱和速率降低，所以AHA中的渗透边界的排水时间延后，整体上呈现出一定的滞后性。在模型的验证阶段，AHA的暗管出流时间比模型的率定阶段要早，这是由于第二次试验的土壤各层的初始含水量大于第一次试验，造成了暗管周围土壤饱和更快。然而就VSL而言，暗管的排水和排盐的过程与实测过程几乎一致。与AHA和PSB的模拟精度相比,VSL的R2最大，RMSE最小，可见添加虚拟土层后模型的精度有了很大的改进。这是因为尽管渗透边界的排水能力较强，但是虚拟土层的饱和导水率Ks本身较低，流经虚拟土层的水流的速度相对较小[29]，所以最终渗透边界流出的水量较小，接近实际流出水量。可见，可以通过减少虚拟土层的饱和导水率来削弱渗透边界排水能力，无需对上层土壤的饱和导水率Ks做出较大的改动。从这点来看，VSL更能反映暗管排水过程的实际情况。

本文的研究结果还表明虚拟土层的圆环宽度(RW)与Ks有很好的线性关系。这是由于虚拟土层本身的Ks较低，随着土层的扩展，渗透边界附近的水流流入量会越来越小[27]，应该增大暗管周围虚拟土层的饱和渗透系数(Ks)，使得模拟值与实测值基本相符。虚拟土层区别于实际土壤，虚拟土层是用来表征暗管性能的特殊“土层”，本研究中的虚拟土层所选用的参数能够较为准确地模拟暗管的真实出流情况。然而在不同的试验或二维模拟区域中，虚拟土层的参数和环宽往往是不同的。为了简化虚拟土层参数的选取，公式(7)仅考虑了环宽改变时对饱和导水率的影响并建立了线性关系，利用此公式能够快速确定不同模拟情景中的虚拟土层参数的取值，不同情景中的模拟结果也验证该公式的准确性。此外，在建模阶段，虚拟土层的环宽(RW)不宜超过埋设暗管的土层厚度，否则可能会对土壤的物理性质和虚拟土层饱和导水率的选取产生较大的影响，并最终导致模型的模拟精度较低。

4 结 论

1)根据暗管的物理结构、实际排水面积和排水速率，提出了两种策略(VSL和AHA)来概化HYDRUS-2D中的暗管边界，并通过2次室内暗管排水试验的实测值对模拟的排水过程、排盐过程、土壤剖面水分和盐分分布的精度进行了评估。模拟值与实测值具有很好的一致性。

2)在3种模型的各项模拟结果中，VSL在对暗管排水和排盐过程的模拟最为准确。VSL能够很好地反映暗管的排水规律和水盐运移情况。因此，VSL可以对干旱地区应用暗管排盐技术的理论研究提供帮助。

3)基于VSL的模拟试验，发现虚拟土层的环宽(RW)与饱和导水率(Ks)之间也存在定量的线性关系。在不同的模拟情景中，应用此公式可以确定虚拟土层中的相关模型参数。

在本文的研究中，仅考虑了土壤剖面盐分和水分分布的静态值来率定模型而忽略了过程中的动态变化，并且边界研究仅限于暗管边界。下一步的研究中需要考虑土壤水盐的动态变化和上边界的概化形式以进一步提高模型的模拟精度。