1、基于全耦合的地表径流与土壤水分运动数值模拟朱 磊 11,2，田军仓 21,2，孙骁磊 31,2(1.宁夏大学土木与水利工程学院,宁夏银川 750021；2.旱区现代农业水资源高效利用教育部工程研究中心，宁夏银川 750021)摘 要：针对降雨径流从坡面流下的过程中会发生下渗，导致土壤水非饱和带含水率增大这一动力学过程，本文从物理机制上对土壤水和地表水进行耦合，将二维平面地表水模型叠置在土壤水模型的顶部，对土壤水、地表水模型进行相同的空间和时间离散，在模型的计算过程中通过达西流关系对两者之间的水量交换进行计算（双层结点法耦合）或整合离散方程的整体法进行耦合。通过对两种耦合方法的比较，和与前人的实
2、验结果对比，本文模型与耦合方法能够准确模拟和预测地表径流与土壤水分运动过程。研究结果可为分析地表水流与饱和非饱和带水分与溶质耦合机理提供理论支持。关键词：地表径流；土壤水；水量交换；全耦合中图分类号：P64.2 文献标志码：A 文章编号：1001-6791前 言从降雨到径流的形成是一个复杂的动态过程，对地表产流情况的研究，必然要考虑土壤的入渗问题。在降雨和灌溉期，一部分随水流入渗到非饱和带和地下水中，当土壤饱和或来水量大于土壤表层的入渗能力时，地表将会形成多余的水量，或积于地表，或形成径流。由于集流面坡度、坡位、坡形及处理方式的不同，使地表坡面水量转化的规律变得比较复杂。地表所接受的雨水在向土
3、壤入渗过程中存在多维运动，且径流沿坡面不断累积 1。目前，国内外关于地表水与土壤水系统之间的互动关系的研究，多采用野外实地监测与室内试验相结合的方法。应用数学模型方法模拟研究时，由于观测资料的缺失、条件限制或所关注问题的侧重点不同，模型构建侧重于水循环系统的某些水文过程，而对其他水文过程采取简化处理2。而对于地表水流土壤水动力学耦合通常是将土壤水运动和地表水体流动作为两个不同的系统，分别采用连续性方程和圣维南方程进行描述，模拟区域水文过程时，两个系统在一定程度上实现过程耦合，包括外耦合、迭代耦合、完全耦合等 3 种模式 3-5。外耦合模式虽然计算效率较高，但因未考虑土壤水对地表水的反馈作用，致
4、使模拟精度较低。迭代耦合模式的计算效率相对较低，但模拟精度得到改善 6-8。完全耦合模式则是同步计算地表 土壤水动力学模型，是一种完全意义的耦合，对于各个系统以及内部边界条件共同进行求解。即在土壤介质和地表水，以及边界条件中的数值方程同时进行求解。具有震荡敏感性低和质量守恒误差小等优点 9，但由于计算的复杂性，国内外研究仍较少。本文针对完全耦合仍缺乏系统研究的现状，根据地表径流物理机制与饱和非饱和土壤水分运动理论，在一定降雨过程条件下，分别采用驱动为两者之间水头梯度的双层结点法和整合离散收稿日期： 网络出版日期：网络出版地址：基金项目：水利部公益性行业科研专项“西北渠井结合灌区配置标准与调控技
5、术研究” ，国家自然科学基金（51369025）作者简介：朱磊（1980） ，男， （汉族） ，宁夏中宁人，副教授，博士，主要从事土壤水分运动研究。Email:方程的整体法，两种完全耦合的方法，建立与土壤物理参数、灌排方式相关的地表径流与饱和非饱和土壤水运动全耦合模型，根据前人的实验结果，利用本文所提出的模型，在保证降雨与地表水、饱和非饱和土壤水水量平衡的情况下，精确描述了地表径流产流与土壤中水分运动过程并对两种完全耦合方法进行了对比分析。1 地表径流与土壤水全耦合运动模型1.1 土壤水运动数学模型 土壤水运动方程应用三维饱和非饱和土壤水分运动方程（Richards 方程）描述，方程的基本形式
6、为 10-11：()()swijrwGSjzSKQtxx（1） 式中， 为土壤饱和含水率 ，为饱和度，其值 ， 为土壤体积含水率 ,s3Lws3L为土壤基质势 ， 为根系吸水项或其它源汇项 ， 为空间坐标，当研究LQ31LTi(,23)x垂直截面的平面流时， 为横坐标， 为垂直坐标，且取向上为正， 是时间 ， 为1x3xztTijK土壤饱和渗透系数张量 ， 为相对非饱和水力传导率 。 为地表水土壤水交换量TrwKGSQ。31LT1.2 地表径流运动数学模型 采用二维扩散波方程描述地表径流运动 11-15。（2）()()ssssxsySsGShhhdkdkdt式中， 为水面高程 ， ， 为地面高
7、程 ； 为水深 ； 为源汇项 1LT，ShLSSzsZLLQ和 分别是 x 和 y 方向上的地表传导系数 1T， 为地表水土壤水交换量 。 sxky GS1.3 方程参数的处理 土壤水分特征曲线和非饱和土壤水力传导率采用 van Geneuchten 模型来表示 10-11：0()1srhrmnhs（3） （4）()0()srKh其中：（5）21/21/()mreeS其中：（6）res（7）1-,mn其中 为残余体积含水率 ， 为饱和含水率 ， 为饱和水力传导度 ，r3Ls3LsK1LT为有效饱和度 ， 、 和 是经验常数。eSn对于地表水模型参数， 和 可分别由下式求出 11-13sxky，
8、 （8）2/31/2sxsdnh2/31/2sysdknh式中的 和 分别是 x 和 y 方向上的曼宁糙率系数 ； 为边界上的法向地表传导系xny 1/3LTnk数 。1LT1.4 边界条件的处理对于土壤水渗出面边界，对程序模拟时间内有可能成为渗出面的所有结点予以事先设定。在每次迭代期间，潜在渗出面的饱和部分处理成水头为 0 的已知水头边界，而非饱和部分则处理成流量 为 0 的已知流量边界。两个部分的长度在迭代的过程中不断地调整，直到饱和部分流量Q的计算值（由式（9）计算）和非饱和部分水头的计算值都为负值为止，此时表明水流仅通过渗出面边界的饱和部分离开计算区域 10。（9）irwiQKk其中
9、为渗出面边界结点水头值， 为渗ii出面边界某结点流量值，渗出面边界流量 为所有该边界结点流量之和。对于地表出流边界，采用临界深度排水边界，单宽流量 为 11,13-14：sQ（10） 3ssgd其中 为重力加速度。2 地表径流与土壤水全耦合方法本文采用的两种完全耦合方法基本思路为（图 1）：首先将二维平面地表水模型重叠在土壤水模型的顶部，对土壤水、地表水模型进行相同的空间和时间离散，表层的地表水模型结点与土壤水模型顶部结点具有完全一致的空间坐标，每个地表水结点与相应的土壤水结点进行水力联系。在模型的计算过程中通过达西流关系对两者之间的水量交换进行计算（双层结点法耦合）或整合离散方程的整体法进行
10、耦合 11-15。图 1 地 表 水 与 土 壤 水 耦 合 示 意 图 Fig.1 Spatial discretization of the surface flow system and its conection to the subsurface 2.1 双层结点法耦合双层结点法是通过土壤水与地表水间的水量交换对两模型进行全耦合。其关系描述如下 11-15：（11）()()()swijrwGSjssssxsySsSzSKQt xhhhdkdkdt 其驱动主要为两者之间水头梯度，交换量 为： GS（12）o()rzexchkKQl式中， 为 方向土壤饱和渗透系数张量 ； 和 分别为土壤
11、水和地表水的水头L ；zKZ1LTo为耦合长度，本次计算中采用 0.1；计算 时采用迎风算法，当 时， 为地表介质相对exchl rkohrk渗透率，计算公式如下（其中， 为细沟贮存量，由于地表存在细微的沟状起伏，水必需先填ohd满这些起伏才能开始产生径流，填满这些起伏的水量就是细沟贮存量） ，当 时， 为土壤相ohrk对非饱和水力传导率 。2(1)001sohsohsrdssohdk（13）2.2 整体法由土壤水流运动离散方程（14）1 11*()()()L LLewsIIJIJIGSIIeBVQVqNdt 式中： 为结点 的体积影响因子 3， 表示与结点 相连接的结点序列， 为基质区结IV
12、I I IJ点 与结点 之间水势差， 表示采用迎风加权算法时，结点 与结点 之间的相对非饱和水力IJIJIJ传导率- ， 为系数矩阵， 是指在边界结点 所代表的附近边界上通过的流动通量 ，当Iq LT流速的方向为从模拟区域向外时，该流量值为正，否则为负，对于内部结点 ， ， 表示I*0qeB研究区域边界 ， 为单元 中相应的基函数。LINe可知，如果结点 为土壤与地表共有结点，该结点上包含两区的双重性质，若假定 ，IsIh则对于结点 的土壤水流动离散方程中， 中包含有两区的水量交换，而对于地表水I I()QVL+1IGS流的离散方程：（15）1 1111(/2)() ()sILLLLsIsI
13、sIJsIJsIsGSIsjAh QdAt式中： 为结点 的面积影响因子 ， 表示与结点地表结点 相连接的结点序列，sIA2sI为地表结点 与结点 之间水头差， 表示采用迎风加权算法时，结点 与结点 之间的地表sIJJsIJIJ相对传导率 ， 为系数矩阵。sIJ式（15）中 即为两区水量的交换量，可表示为：1LsGSIdQA（16）1111(/2)()sILLLsIsIssI sIJsIJsIjAh QAt将式（16）代入式（14）得整体法的两区耦合方程：（17）1 11*1(1/2)()() sIL LLewsIIJIJIIeBLsIsI sIJsIJsIjVVqNdtAh QAt 上式就是
14、采用了整体法的两区耦合离散方程 16-17。3 模型验证与耦合方法比较本文采用国外已有室内试验结果来检验模型的准确性。实验设置如图2所示，T1T13为土壤水势监测点，出流口1与2分别收集地下与地表的出流量，实验过程为先将土壤完全饱和后放置24小时，后进行第一次降雨过程，持续时间45分钟，降雨强度0.69cm/min，第一场降雨完成后放置60分钟，再进行第二次降雨过程，持续时长30分钟，降雨强度0.69cm/min。具体步骤见参考文献18 18， 因模型中采用的各参数公式与文献相同，且为了模型验证的准确性，直接选用参考文献18中参数，如表 1所示： 0 30 60 90 120150180030
15、6090120 出 流 口 2出 流 口 1模 拟 降 雨 T12T13T1T10T9T8T7T6T5T4T3T2T1D CBA高度(cm)长 度 (cm)图 2 实验设置示意图 4Fig.2 The schematic diagram of the experimental setup4表 1 模型主要参数表Table 1 The main parameters of the model孔隙率（ %）饱和水力传导度 K(cm/min)残余体积含水率 r土壤饱和含水率 s(1/cm)曼宁系数 n(min/cm-1/3)33 1.28 0.024 0.33 5.00 0.05 0.0012243
16、.1 观测点水势响应图 3 为沿坡度方向三个位置上，孔隙水压力实验值与模拟值的比较，表 2 为实验值与模拟值误差分析。在三个不同的深度位置上进行测量，在地表面（地面以下 5cm 处；T1，T5 和 T12） ，中间层土壤剖面（地面以下 32cm 处；T2，T6 和 T12）以及靠近底层（底层以上 1cm；T3，T7和 T13） ，各点位置如图 2 所示。由图 3 及表 2 可以看出，本文所构建的模型与实验实测结果匹配良好，根据参考文献18 ，土壤在实验进行前完全饱和后放置 24h，由于边界渗流排水，在模拟降雨开始前沙子整体干燥，由模型计算得到压力水头数据与实验数据相符，在实验开始前沙子整体干燥
17、，尽管在土层的底部，下边界（下游）附近，存在一个小区域（T13）仍保持湿润。从系统完全饱和并自由排水 46.0 h 以后，坡底的压力水头在靠近边坡出口（T13）处近似为零；相反，坡底的压力水头靠近上/上游边界（T3 ）为-36cm。在地表附近，实验土壤更加干燥，压力水头从-41cm（坡口附近）变化到-56cm （坡顶） 。总体而言，整个系统，在实验砂土上坡和近表面区域的初始水分分布不均匀。实验与数值模拟结果均显示，沿坡度方向，在同一深度土层不同位置处，压力水头随着时间的变化有着非常相似的趋势。在地表以下 5cm 的位置(T1，T5 和 T11)，降雨开始 2 分钟以后压力水头迅速增大，一直持续
18、到第 20 分钟。压力水头跳跃式的快速增长反映了由降雨渗透，穿过地表，进入初始干沙的前期湿润锋。尽管湿润锋的到来明显提高了局部压力水头，但仍是负压，表明土壤水饱和度接近但仍低于 100%。这种部分饱和现象持续了 20 分钟；在此期间，没有更多的水积累在浅层而是通过渗透到达深层。之后随着坡底积水量增大，湿润锋从下游边界推进朝向上游边界运动，直到坡顶达到饱和。与其他测量点不同，接近土壤剖面的底部测量点，压力水头从初始值一直增加到稳态值（图3） 。压力水头在这些地区的持续上升表明，流入水量超过了流出的水量，提高了当地的土壤含水量，压力也跟着上升。由于该系统的下界面是不透水层，所以水从土壤剖面的底部，
19、以平行于坡基面的方向由渗出面边界流出。与文献18相同，在第一次降雨水头值升高的模拟结果中，相比于实验观察预测水势响应还是有明显的滞后现象（图3） 。这种滞后的原因是由于Lisse效应的存在 19，不透水边界，最初干燥的土壤和地表土壤水分的快速交换导致空气滞留的效果明显，而在土壤含水率相对较大的第二次降雨期，压力水头的滞后响应明显减少(图3) 。0 50 10 150-4-200实 验 值双 层 结 点 法 整 体 法压力水头(cm)时 间 (min)T10 50 10 150-4-20020实 验 值双 层 结 点 法 整 体 法压力水头(cm)时 间 (min)T20 50 10 150-2
20、002040实 验 值双 层 结 点 法 整 体 法压力水头(cm)时 间 (min)T3 0 50 10 150-40-200实 验 值双 层 结 点 法 整 体 法压力水头(cm)时 间 (min)T50 50 10 150-4-20020实 验 值双 层 结 点 法 整 体 法压力水头(cm)时 间 (min)T60 50 10 15002040实 验 值双 层 结 点 法 整 体 法压力水头(cm)时 间 (min)T70 50 10 150-40-200实 验 值双 层 结 点 法 整 体 法压力水头(cm)时 间 (min)T110 50 10 150-200实 验 值双 层 结
21、点 法 整 体 法压力水头(cm)时 间 (min)T120 50 10 150510实 验 值双 层 结 点 法 整 体 法压力水头(cm)时 间 (min)T13图 3 孔隙水压力实验值与模拟值的比较Fig.3 The water pressure value comparison between experimental and the simulation表 2 孔隙水压力实验值与模拟值误差分析Table 2 Error analysis of water pressure between experimental and simulation valueT1 T2 T3 T5 T6
22、T7 T11 T12 T13标准差（cm ） 2.01 2.07 2.18 2.03 4.33 1.62 1.50 2.43 0.21双层结点法平均误差值(cm) 1.60 1.65 1.74 1.63 3.46 1.30 1.20 1.94 0.18标准差（cm ） 2.49 1.73 0.35 3.79 3.25 2.53 3.54 2.78 0.72整体法平均误差值(cm) 1.99 1.38 0.28 3.02 2.59 2.02 2.82 2.22 0.58图 4 为沿坡基面部署 1（T13） ，2（T10） ，3（T7）和 4（T3）点压力水头值随时间变化过程。可以看出，虽然与实验
23、室实验相比数值模拟（双层结点法）的结果存在明显的滞后现象，但也显示出前面所述湿润锋从下游边界推进朝向上游边界运动过程。5 10 15 20 25 30 35-40-200204060压力水头/cm时 间 /min12 3 4 零 压 力 水 头实 验 值15 20 25 30 35 40-40-200204060压力水头(cm)时 间 (min)12 3 4 零 压 力 水 头本 文 模 型 (双 层 结 点 法 )图 4 沿坡底各点水头值变化Fig.4 The head value change along the bottom （a） （b）双层结点法（c） （d）整体法图 5 降雨 15
24、 分钟与 27 分钟两种耦合方法流速与含水率的分布Fig.5 Rainfall for 15 minutes (a), (c) and 27 min (b), (d) two coupling method of flow velocity and water content distribution图 5 为两种耦合方法计算出降雨 15 分钟（a） 、 （c）与 27 分钟（b） 、 （d）流速与含水率的分布，提供了解释土壤水流和饱和区发展的关系。从图中可以看出，两种方法计算出的流速与含水率几乎相同，可见两种耦合方法在计算精度上差别不大。在降雨的初级阶段（15 分钟） ，水流运动主要以垂直方
25、向为主。只有在下边界坡底附近出现较小的饱和区域。在 27 分钟，由于饱和区形成和扩大，靠近上部饱和区，水流被分流在上坡方向，饱和区中的水向上坡移动，上坡水流速（达西速度）比在下坡方向较小。3.2 地表与渗出面边界排水过程图 6 为出流过程实测与模拟值比较，可以看出本文模型模拟结果较好，这也说明本文模型在边界条件处理上的正确性。这里需要说明的是在整体法耦合的时候，地表坡底边界结点既是渗出面边界又是临界深度排水边界，在处理渗出面边界时，潜在渗出面的饱和部分处理成水头为 0 的已知水头边界，但是在整体法中，假定 ，于是 ，而此时地表由于已积水，故 ，IshsI0h sIh整个地表出流流量过程计算出现
26、错误，为了解决整体法在这种情况下的缺陷，本文将地表坡底边界结点只设置为临界深度排水边界，出流量计算结果与实验值对比表明这种方法对于模型计算结果影响较小。0 50 10 1505010150 实 验 值双 层 结 点 法 整 体 法流量/(cm3/in)时 间 /min（A）地表出流流量过程0 50 10 150501015020250 实 验 值双 层 结 点 法 整 体 法流量/(cm3/in)时 间 /min（B）渗出面边界流量过程 0 50 10 15010203040 实 验 值双 层 结 点 法 整 体 法流量/(cm3/in)时 间 /min（C）总流量过程图 6 出流过程实测与模
27、拟值比较Fig.6 The comparison between flow measurement and simulation value表 3 出流过程实测与模拟值误差分析Table 3 Error analysis of water flow between measurement and simulation value地表出流流量过程 渗出面边界流量过程 总流量过程标准差（m 3/min) 8.32 84.15 44.96双层结点法 平均误差值（m 3/min)6.64 67.14 35.88标准差（m 3/min) 8.36 66.39 39.56整体法 平均误差值（m 3/min)6.67 52.97 31.573.3 双层结点法参数敏感性分析由于在双层结点法耦合中，引入了参数耦合长度 ，因此需对此参数进行参数敏感性分析，exchl本文将调整后的参数取值加减 10%以上（表 4） ，进行了参数敏感性分析，重点考虑了这些参数对压力水头、出流量的影响。从表 4 可以看出，耦合长度 对压力水头与出流量计算结果影响较exchl小，均在 5之内，模型对此参数不敏感。表 4 参数敏感性分析
