1. 引言

黄土地区厚度巨大的土层具有储存水分的良好基础，其中的土壤水可为植被生长提供水分和养料 [1] [2] ，而适度埋深的土壤水分零通量面是维持土壤水分持续供给的必要条件之一。但近50年来，气候变化加上工农业开发，黄土高原区域尤其是灌区内土壤水分的变化日趋剧烈，黄土灌区土壤水分动态研究的需求也日趋迫切 [3] [4] 。目前学者们在土壤水分时空变化与水量平衡，植被承载力等方面取得很多的研究成果，为黄土灌区土壤水分亏缺的解决提供了大量科研依据 [5] 。但进一步的深层次研究需探明水循环中包括深层渗漏和蒸散发量等水文变化量，也就是水分零通量面的时空变化情况，而这类变量在土层深厚黄土地区的动态监测较为困难 [6] ，尤其是供水频繁的灌区，因此相应的研究较少。而使用数值方法，依据少量监测数据获得接近实际状态的模拟结果可以有效解决问题，Hydrus就是其中一类常用的模拟方法 [7] [8] 。

宝鸡峡灌区位于陕西省关中西部，地处大陆性半干旱气候区，年均降雨量约600 mm左右，年均蒸发量约1100 mm。灌区总灌溉面积1.943 × 10⁵ hm²。按自然地形和工程布局，分为塬上和塬下两大灌域，塬上灌域也称黄土塬灌区，塬下灌域称为渭河阶地灌区。黄土塬灌区的地表水源主要来自于渭河干流林家村渠首，灌溉面积1.241 × 10⁵ hm²，塬下灌区由渭河干流魏家堡自流引水，灌溉面积7.020 × 10⁴ hm²。经过50多年运行，灌区基本上形成了地表水为主、地下水为辅的灌溉格局，用水模式在大型黄土灌区中极具代表性。灌区主要作物以玉米、小麦、油菜、果树等为主 [9] ，是全省粮、油、果的生产和供应基地，灌区的发展对全省的粮食安全和关中经济社会发展具有重要支撑作用。因此研究该区黄土中水分零通量面的运动特征十分必要。

2. 研究方法

Hydrus2D是一款用于土壤水分运移和热运移建模的软件，基本原理是建立数学模型来描述土壤中水分和热量的运动和变化。通过离散化土壤系统的物理特性，将系统建模为一系列方程，然后求解这些方程得到系统各部分的水、热分布情况。Hydrus2D的水分运移模型基于Richards非线性偏微分方程，方程描述了水分在土壤中的运动速度与土壤水势梯度间的关系。在Hydrus2D中，Richards方程被离散化为一组有限差分方程，求解后可得到土壤中的水分含量、水势以及水分的流动情况 [10] 。

本文使用Hydrus2D进行水分运移模拟，通过参数的设置，构建地质模型和多层黄土水流入渗模型，利用有限元生成模块进行长方形区域网格剖分，设置初始水头压力条件和边界条件，通过求解水流控制方程等方程组，不断迭代，直到所有设置的时间步长都计算出结果，得到研究区域网格上每个节点的水头值、土壤含水量等数据，据此构建多层黄土中水分入渗和土壤水分零通量面位置变化模型，进行数值模拟。

3. 情景设定和地层分析

3.1. 情景设定

通过分析资料确定模拟地层结构，然后使用Hydrus2D Water Flow模块构建多层结构黄土入渗模型(假定制取长10 m的土柱)，供水条件设置为人工灌溉和天然降雨状态，土壤水力参数设置见表2。

人工灌溉供水条件的强度通过分析灌区玉米和冬小麦的灌水定额及灌溉期进行设定，本文设置次灌溉定额为60 m³/亩，灌溉期为30天，每天灌溉8小时，灌溉强度0.03749 cm/h。而天然降雨强度受季节等因素条件影响较大，本文设定天然降雨强度为0.0458 cm/h，降雨历时为4小时，降雨间隔时间8小时。据陕西省气象局资料，宝鸡市年平均蒸发量约为1200 mm，平均每天蒸发量约为3.3 mm，综合考虑取蒸发强度为0.02083 cm/h。

边界条件的设置：顶边界为大气边界，底边界设置为自由排水边界，两侧边界假定无水分交换作用，设置为零通量边界。最后通过分析模拟结果，总结零通量面位置的变化规律及其上下运动。

3.2. 地层结构分析

通过比对宝鸡峡塬上灌区钻孔剖面图，从东至西选取三个代表钻孔：84号钻孔(武功)、86号钻孔(乾县)和93号钻孔(礼泉) (图1)。钻孔分析发现，宝鸡峡塬上灌区的地质结构自上而下大都包含有上更新统黄土、中更新统上部黄土、中更新统下部黄土以及下伏砂砾卵石、砂、粘质砂土等岩性。按照量算的地层比例，确定钻孔各地层的虚拟厚度，结果见表1，各地层水力参数见表2。

4. 模型建立

根据模拟设计的构想，根据84号、86号、93号钻孔资料，分别设定人工灌溉和天然降雨状态两种边界条件，建立6个模型。

4.1. Water Flow模型原理和假设

4.1.1. 模型原理

Water Flow模型主要用于描述土壤中水分的运移过程。该模型基于Richards基本方程：

$\frac{\partial \theta }{\partial t}=\frac{\partial }{\partial x_i}\left[K\left(K_{ij}^A\frac{\partial h}{\partial x_j}+K_{iz}^A\right)\right]-S$ (1)

θ为体积含水量(cm³/cm³)；h为压力水头(cm)；K为饱和导水率(cm/d)； $K_{ij}^A$ 为无量纲张量K^A的第i个主分量在j方向上的分量；源汇项S(d⁻¹)用于解释根系吸水 [11] 。

模拟过程中，Water Flow模型将土壤划分为一系列离散的网格单元，通过求解Richards方程来描述土壤中的水分运移过程，并考虑土壤温度、孔隙度、压缩性等因素的影响。

4.1.2. 模型假设

Water Flow模型的假设包括：多孔介质土壤，水分运移受毛管力、重力和惯性力影响。土壤孔隙度、渗透率、吸附性能、压缩性、温度等参数会随时间和空间变化。土壤水分是连续的，无间隙或气泡。土壤中的水分流速与水分梯度成正比，且满足Darcy定律 [12] 。土壤毛管力是由土壤颗粒表面张力和水–气界面张力间的相互作用所产生。土壤水分运移也受温度影响，升温会降低土壤孔隙度和水分吸附性能，从而影响水分运移速度。

4.2. 模型参数的确定和构建

4.2.1. 几何参数确定

6个模型均为长100 cm，深1000 cm，有结构的有限元网格简单几何的二维流动运移的垂直面。

4.2.2. 水流参数确定

6个模型均在“主要过程(Main Processes)”对话框中选用WaterFlow(水流)模块；在“时间信息(Time Information)”对话框设定：模拟起始时间为0 h，模拟结束时间为72 h (三天)，初始时间步长0.002 h，最小时间步长0.001 h，最大时间步长1 h。可变时间边界数设置为72。

天然降雨供水条件模型3个，“输出信息(Output Information)”参数设置为模型每运行10个时间步长输出一个数据，输出详细数据数组总数为360 (即输出结果时间间隔为0.2 h)。其可变时间边界参数设置为：降雨强度0.0458 cm/h，蒸发强度0.02083 cm/h。间隔8 h降雨一次，每次降雨历时4 h，降雨时忽略蒸发。

人工灌溉供水条件模型3个，同样运行10个时间步长输出一个数据，输出详细数据数组总数为108 (即输出结果时间间隔为0.667 h)。可变时间边界参数设置为：灌溉强度为0.03749 cm/h，蒸发强度为0.02083 cm/h，每天灌溉8 h。

6个模型的迭代标准(Iteration Criteria)：van Genuchten-Mualem水力模型，不考虑滞后性。

根据表2的各地层颗粒比，使用神经网络预测模块进行土壤参数预测。在区域属性窗口下，输入土壤分层参数。

6个模型的区域离散化x (长)方向取11个节点，z (深)方向取1001个节点。

4.2.3. 初始条件和边界条件的设置

6个模型的初始压力水头均使用Hydrus软件默认的0~100 cm值。顶边界为大气边界，底边界设置为自由排水边界，两侧边界假定无水分交换作用，设置为零通量边界。

4.3. 模型运行

所有设置完毕后运行模型。模型输出结果包括：压力水头图、土壤含水量图、入渗速率图、速度矢量图、观测点结果及边界通量等。图2、图3分别为天然降雨下93号钻孔土柱12.20 h (上层50 cm内)的入渗速率图、速度矢量图。

5. 模拟结果分析

在Hydrus中，若速度矢量箭头图中的某个点的速度箭头为零，则这些点连成的面就是零通量面。水份在零通量面这个速度场中的特殊区域内被认为是静止的，不会向任何方向移动 [13] [14] 。

根据零通量面理论，蒸发状态下出现的水分零通量面为土壤水势曲线的极大值，又称发散型零通量面。与之相对土壤水势曲线的极小值出现的水分零通量面又称收敛型零通量面。在速度矢量图内，若某个面上方的速度矢量箭头向上，下方的速度矢量箭头向下，则该面为发散型零通量面；若某个面上方的速度矢量箭头向下，下方的速度矢量箭头向上，则该面为收敛型零通量面 [15] 。

通过观察72 h的6个模型速度矢量图模拟结果，记录零通量面出现的时间和位置，制成Excel表格，并根据整理好的时间位置数据绘图。图4、图5、图6、图7、图8、图9分别为人工灌溉下84号钻孔土柱零通量面分布图、86号钻孔土柱零通量面分布图、93号钻孔土柱零通量面分布图以及天然降雨条件下84号钻孔土柱零通量面分布图、86号钻孔土柱零通量面分布图和93号钻孔土柱零通量面分布图。

5.1. 人工灌溉条件下黄土灌区土层中零通量面的分布特征

如图4，84号钻孔土柱的发散型零通量面在蒸发状态下出现，且其随着蒸发时间的延长不断下移，厚度逐渐变大，持续16小时的下移极限在距地表34 cm处。在灌溉状态下，发散型零通量面开始上升且厚度不断变大，收敛型零通量面不断下移，极限值在距地表17 cm处。两种类型的零通量面均在供水开始一小时内消失。

如图5，86号钻孔土柱的发散型零通量面在蒸发状态下出现，且随着蒸发时间的延长不断下移，且厚度逐渐变大，持续16小时的下移极限在距地表34 cm处。灌溉状态下，发散型零通量面开始上升且厚度不断变大，收敛型零通量面不断下降，极限值在距地表17 cm处。两种类型的零通量面均在供水开始一小时内消失。与人工灌溉下84号钻孔土柱零通量面变化相近。

如图6，93号钻孔土柱的发散型零通量面在蒸发状态下出现，且随着蒸发时间的延长不断下移，厚度逐渐变大，持续16小时的下移极限在距地表34 cm处。在灌溉状态下，发散型零通量面开始上升，厚度增大，收敛型零通量面不断下降，极限值在距地表17 cm处。两种类型的零通量面均在供水开始一小时内消失，与人工灌溉下84、86号钻孔土柱零通量面变化相近。

5.2. 天然降雨条件下黄土灌区土层中零通量面的分布特征

天然降雨条件下，84号钻孔土柱的发散型零通量面在蒸发状态下出现(图7)，且随蒸发时间的延长不断下移，下移速度在蒸发开始的3小时后趋于稳定，持续8小时的下移极限在距地表20 cm处。在降雨状态下，发散型零通量面厚度变大，收敛型零通量面快速下移，极限值在距地表7 cm处。两种零通量面均在降雨开始半小时内消失。

86号钻孔土柱的发散型零通量面在蒸发状态下出现(图8)，且随着蒸发时间的延长不断下移，同时零通量面的厚度不断增大，持续8小时的下移极限在距地表30 cm处。在降雨状态下，发散型零通量面厚度变大且有上移的趋势，收敛型零通量面快速下移，极限值在距地表7 cm处。两种类型的零通量面均在供水开始半小时内消失。从图8可看出天然降雨条件下地层结构的变化对发散型零通量面变化影响较大。

图9中，93号钻孔土柱的发散型零通量面在蒸发状态下出现，且随蒸发时间的延长不断下移，下移速度在蒸发开始的3小时后趋于稳定，持续8小时的下移极限在距地表20 cm处。在降雨状态下，发散型零通量面厚度变大，收敛型零通量面快速下移，且极限值在距地表9 cm处。两种类型的零通量面均在供水开始半小时内消失，与天然降雨下84号钻孔土柱零通量面变化相近。

6. 结论

1) 在人工灌溉供水条件下发散型零通量面多出现在距土壤表面20~30 cm间，收敛型零通量面多出现在距土壤表面10~20 cm间。而天然降雨条件下发散型零通量面多出现在距土壤表面10~20 cm之间，收敛型零通量面比较稳定，出现在距土壤表面5~7 cm处。

2) 灌区主要作物玉米和冬小麦的根系层主要在15~30 cm处。当发散型零通量面移动到根系层时，可导致土壤含水量降低，作物生长受限。特别是在玉米和冬小麦的拔节期、抽穗期等关键生长期，水分供应不足可造成作物减产甚至死亡。本文中，发散型零通量面的变化表明，持续6~12小时的蒸发后，作物根系层土壤水分含量将不足，需要人工灌溉补充。

3) 对比不同供水条件下零通量面分布发现：人工灌溉条件下的发散型零通量面和收敛型零通量面下移的极限值比天然降雨条件下的要深。由86号钻孔土柱(两层黄土)模拟结果可知地层结构对天然降雨条件的发散型零通量面变化影响较大，而在人工灌溉条件(连续蒸发时间长)下地层结构对零通量面的变化影响不大。

4) 供水时段的合理分配对于维持土壤水分含量稳定(零通量面下移极限)有较大的影响。

基金项目

国家自然科学基金项目“变化环境下水文非平稳强波动序列随机建模与预测方法研究”(52079110)。

参考文献