平原区饱和-非饱和土壤水运动模型及数值算法研究

陈景波，王船海，杜世鹏，杨 海，张梦菲，朱立国

(1.河海大学水文水资源与水利工程科学国家重点实验室，江苏南京210098；2.浙江省水利水电勘测设计院，浙江杭州310002；3.江苏省水文水资源勘测局无锡分局，江苏无锡214000)



平原区饱和-非饱和土壤水运动模型及数值算法研究

陈景波1，王船海1，杜世鹏2，杨 海1，张梦菲1，朱立国3

(1.河海大学水文水资源与水利工程科学国家重点实验室，江苏南京210098；2.浙江省水利水电勘测设计院，浙江杭州310002；3.江苏省水文水资源勘测局无锡分局，江苏无锡214000)

以土壤水运动理论为基础，构建了研究区三维饱和-非饱和土壤水运动模型，通过网格剖分模拟研究区各个位置土壤的体积含水率以及变换过程，并采用多个算例验证模型的稳定性和准确性。计算结果可知，土壤水分及地下水水位变化过程符合土壤水运动原理和日常经验，重现了特定土壤水运动的基本过程，为探究平原区水循环过程、水分转化关系以及进行大尺度平原区水文预报提供理论基础。

土壤水；水循环；饱和-非饱和流；数值模拟

0 引 言

平原区土壤水是水文循环中的关键纽带，运移交换频繁[1- 3]。现阶段研究土壤水运动主要通过耦合现有的地表水和地下水模型来模拟计算，连接方式多以地下水模型进行三角形或矩形差分构建基本单元[4]。在计算过程中，两模型之间实际上是分离的，地表水和地下水模型分别独立进行计算，模拟实际情况下的产流和地下水动态变化，模型之间只是通过下渗量的传递来进行简单、松散的耦合[5]，如SWATMOD[6]模型、MODFLOW-DAFLOW[7]模型等。这就导致模型不能准确地模拟土壤水，特别是平原区饱和-非饱和带土壤水的运动过程，对耦合模型要重点解决的诸如潜水蒸发、土壤水变化及产流量等问题不能很好地解决。程勤波等[8]在可变坡土槽进行了室内人工降水入渗试验，以此来分析MODFLOW-UZF1模型模拟降水入渗补给过程的精度。结果表明，模型不能准确地模拟入渗过程中土壤含水量的渐变过程，忽略了非饱和带土壤水分运移中基质势的影响，在地下水浅埋地区模拟效果较差。

本文基于三维理查德方程推导并构建了适用于平原区三维饱和-非饱和土壤水流动模型，并在太湖流域的金坛市建立野外实测试验基地，以此来观测分析和研究地势平坦、地下水埋深浅的特殊条件下的土壤水运动和浅层地下水水位的变化，以及土壤水分对降水和蒸发的响应。模型以土水势为研究变量，首先对研究区域进行网格剖分，利用有限差分的格点法将研究对象离散化，进而模拟平原区浅层土壤含水率以及非饱和区与饱和区的变换过程。在模型求解过程中，运用“边循环”的概念和矩阵标识法，进一步提高了模型的精确度和计算效率。

1 模型的求解

采用有限差分法进行求解。求解过程中，运用“边循环”的概念和矩阵标识法达到快速准确求解模型的目的。采用有限差分方法将渗流区域划分网格，并采用格点法，对方程中的左端项采用中心差分，右端时间项采用向前差分，方程离散如下

式中，K为x，y，z方向上的导水率；H为总水头；Sf为土壤空隙中含水百分比；us为给水度；C为比水容量；h为压力水头；W为源汇项的水量。

方程联系着格点(i,j,k)在t时刻与t+Δt时刻的水头，同时还联系着与这个格点相邻的6个格点的水头。对每个未知水头的格点都列出上述方程，当格点靠近边界时用给定水头或给定流量的边界条件，这样可以组成一个线性代数方程组。其中，方程左端的水头全部采用t+Δt时刻的值即采用全隐式解法；将上述方程化简为m×x=y方程的格式，并与其周围相邻网格的水头关系方程式进行联立，构成如MX=Y形式的矩阵以方便求解。系数矩阵M有以下特点：每个网格最多会产生7个非零元素，对于宽广计算区域，网格很多，方程阶数很高，对这样的高稀疏矩阵，计算的效率就显得尤为重要。本文采用矩阵标识法求解方程，使网格编码具有任意性，避免非零元素的存储和计算，又不用考虑收敛问题，较好的提高了模型的精确度和计算效率。

2 模型验证

2.1 垂向一维降雨入渗模拟

1 m×1 m×2.3 m的垂向土柱，垂向上剖分成230层，dx=dy=1 m、dz=0.01 m。土柱的初始压力水头为-0.8 m，上边界给4 mm/h的降雨，下边界、四周边界没有水量交换。模型计算结果见图1。

图1 不同降雨历时土壤水分剖面

从图1可知，降雨开始后，在积水点之前，降雨全部下渗到土柱中，湿润峰面下移；当地表含水率达到相对饱和，湿润峰面全部进入土体中，即t=3 h时，土体开始积水入渗，此后湿润峰继续下移；t=13 h时，入渗率开始出现拐点，入渗率明显下降，此后当湿润峰全部通过底部出溢面时，入渗率达到饱和点，土柱开始饱和入渗。

2.2 降雨模拟

2.2.1 单网格降雨入渗

2.5 m×2.5 m×1 m的土柱剖分成25×25×100的网格，即dx=dy=0.1 m、dz=0.01 m。选上边界中任意网格，给予连续10 mm/h 的18 h降雨，下边界、四周边界都没有水量交换，各层的土壤水力参数设为相同，设土柱的初始压力水头为-5 m。z方向的饱和导水率为0.003 44 cm/s，x、y方向的饱和导水率是z方向的1/4。40 h降雨计算结果见图2。

图2 土壤含水率等值线

从图2可知，降雨落到网格后，向下和四周扩散。当x、y方向上的导水率值相等时，由中心网格流入周边网格的流速、水量都是一致的。图2a中，有降雨网格的土壤含水率由0.257上升到0.295，降雨影响到的外围网格的土壤含水率上升到0.260；图2b中，有降雨网格的土壤含水率由0.257上升到0.268，降雨影响到的外围网格的土壤含水率由0.257上升到0.258。对比可知，降雨落到网格上后，10 cm处受降雨影响的面积大约为1 m2，20 cm处大约为0.49 m2，由上到下，受影响面积逐渐减小。对比模型初始和计算完成后的水量与降雨入渗及产流量，模型计算符合水量平衡要求。

2.2.2 多网格降雨入渗

运用上述土柱条件，给上层网格中相互对称的4个角落网格连续10 mm/h的降雨。计算结果见图3。从图3可知，降雨分布在平面对称的4个网格中。在网格土壤特性一致的前提下，土壤水的运动也保持一致，土壤含水率的值相同，且4个角上的等值线保持对称，降雨所在网格下20 cm处的含水率为0.445，降雨影响到的最大区域为0.9 m，该处的土壤含水率为0.426。

图3 土壤含水率等值线

2.3 蒸发模拟

土柱信息同上。在最上层的中心网格给1 mm/h的蒸发，其他网格没有降雨和蒸发。计算结果见图4。从图4可知，当中间网格由于蒸发导致土壤含水率降低时，四周网格的水势平衡被破坏，土壤水在土水势的作用下，由土水势的高处向低处运动。因此，中心网格(有蒸发)的土壤含水率减小幅度最大，由0.258减小到0.253，由中心网格扩散，四周的网格土壤含水率减小幅度逐步减小。

图4 土壤含水率等值线

2.4 浅层地下水水位模拟

平原区地下水埋深较浅，一次降雨对潜水位影响较大。在此，模拟降雨条件下的浅层地下水水位选取3 m×3 m×2 m的一个大土柱，土柱上部高程一致，将土柱剖分成30×30×200的网格，即dx=dy=0.1 m、dz=0.01 m，相当于分成了90个2 m的小土柱，在整个大土柱的左上角的0.5 m×0.5 m的区域内给148 mm/h的降雨。地下埋深3 m处的不同历时x方向上浅层地下水水位随时间变化见图5。从图5可知，土柱含水率在时间上成上升趋势，但在x方向上从左向右有减小的趋势。

图5 x方向上浅层地下水水位随时间变化

40 h土柱地下水水位渲染见图6。从图6可知，降雨发生前，大土柱的地下水水位保持稳定，且各网格的地下水水位持平。降雨发生在某一片区域后，经土壤水垂直下渗和水平运移，降雨所在网格的小土柱的地下水水位率先上升，由降雨中心向外，各小土柱的地下水水位依次增大。而离有降雨网格最远的对角网格的地下水水位保持不变，从地下水水位最高处到最低处，有规律的下降，下降曲线应是光滑的。由于模型计算误差和网格划分不够密，导致图形存在锯齿。

图6 40 h土柱地下水水位渲染

不同地下水水位土壤水分剖面见图7。从图7可知，地下水水位上升过程中，土壤水分剖面曲线先是迅速由扁平变成尖瘦；在地下水水位由0.695 m上升到0.855 m的过程中，土壤水分剖面变动很小。

图7 不同地下水水位土壤水分剖面

2.4 实际应用

试验区选址在常州市金坛朱林镇红旗圩村。该地区地形平坦、坡度小、河网交织、湖塘遍布。土壤多为第四系冲、洪积松散沉积和湖相沉积的粉粘土，颗粒细微，渗透性差。该地区地下水年平均埋深1 m左右，地下水水位对降雨输入的响应迅速。该区域能较好地反映太湖流域平原河网区的基本水文特征。

降雨资料由Campbell TE525MM系列翻斗式雨量筒测定，精度为0.1 mm；土壤含水率采用TDR传感器观测。应用上述模型进行模拟，对比模拟计算含水率和实测含水率。10、20 cm处土壤含水率计算和实测对比见图8。

从图8可以看出，模型较好地模拟了平原区饱和-非饱和带土壤含水率的变化情况，特别是在表层10 cm处，计算拟合数据与实际观测数据的相关性R为0.64，个别场次的降雨模拟误差低于10%，能够反映土壤水的整体变化趋势。相对于实测值，模拟结果整体较低，需对模型进行改进。

3 结 语

本文以土壤水动力学理论为基础，辅以考虑土壤蒸发、降雨土壤入渗等水文循环

图8 土壤含水率计算和实测对比

过程模型，构建了三维饱和-非饱和土壤水运动模型，运用有限差分法推导出三维饱和-非饱和流方程的差分方程。运用“边循环”和矩阵标识法，缩短模型计算时间，求解速度大大提高。实际应用表明，该模型符合土壤水运动原理和日常经验，重现了特定土壤水运动的基本过程，为探究平原区水循环过程、水分转化关系以及进行大尺度平原区水文预报提供理论基础。

[1]DEMETRIOU C， PUNTHAKEY J F. Evaluating sustainable groundwater management options using the MIKE SHE integrated hydrogeological modelling package[J]. Environmental Modelling & Software， 1998， 14(2- 3)： 129- 140．

[2]刘新仁， 杨海舰. 土壤水动力学在平原水文模拟中的应用[J]. 河海大学学报， 1989， 17(4)： 9- 15．

[3]张蔚榛. 地下水与土壤水动力学[M]. 北京： 中国水利水电出版社， 1996．

[4]陈喜， 程勤波， 张志才. 饱和-非饱和水流数值模拟[M]. 北京： 中国水利水电出版社， 2011．

[5]雷志栋， 杨诗秀， 谢森传. 土壤水动力学[M]. 北京： 清华大学出版社， 1988．

[6]张雪刚， 毛媛媛， 董家瑞， 等. SWAT模型与MODFLOW模型的耦合计算及应用[J]. 水资源保护， 2010， 26(3)： 49- 52．

[7]XU X， HUANG G， ZHAN H， et al. Integration of SWAP and MODFLOW- 2000 for modeling groundwater dynamics in shallow water table areas[J]. Journal of Hydrology， 2012， 412(1)： 170- 181．

[8]程勤波， 陈喜， 赵玲玲， 等. 饱和与非饱和带土壤水动力耦合模拟及入渗试验[J]. 河海大学学报， 2008， 37(3)： 284- 289．

(责任编辑 杨 健)

Study on Flow Model and Numerical Simulation of Unsaturated and Saturated Soil Water in Plain Area

CHEN Jingbo1, WANG Chuanhai1, DU Shipeng2, YANG Hai1, ZHANG Mengfei1, ZHU Liguo3

(1. State Key Laboratory of Hydrology-Water Resources and Hydraulic Engineering, Hohai University, Nanjing 210098, Jiangsu, China; 2. Water Resources and Hydropower Survey and Design Institute of Zhejiang Province, Hangzhou 310002, Zhejiang, China; 3. Wuxi Branch, Hydrology and Water Resources Survey Bureau of Jiangsu Province, Wuxi 214000, Jiangsu, China)

Based on the theory of soil water movement, a 3D model is constructed to simulate the movement of unsaturated and saturated water in plain area. The finite difference method is used to simulate soil moisture content and transformation process in each location of study area. Some case examples are used to verify the stability and accuracy of model. The calculation results show that the process of soil water moisture and groundwater level change is reasonable in theory, the simulation results of actual soil water situation is relatively reasonable, which can reflect the hydrological cycle in plain area. The results provide theoretical basis for exploring water cycle process, water transformation relationship and large-scale hydrological forecasting in plain area．

soil water; water cycle; unsaturated and saturated flow; numerical simulation

2016- 04- 06

江苏省水利科技项目(2015007)；国家水体污染控制与治理科技重大专项(2014ZX07101- 011)

陈景波(1990—)，男，河北邯郸人，硕士研究生，研究方向为流域水文模拟及预报．

S152.7(253)

A

0559- 9342(2016)09- 0013- 04