饱和-非饱和耦合模拟文献摘录

平原区饱和-非饱和土壤水运动模型及数值算法研究_陈景波（2016）

0引言

        平原区土壤水是水文循环中的关键纽带，运移交换频繁[1-3]。

现阶段研究土壤水运动主要通过耦合现有的地表水和地下水模型来模拟计算，连接方式多以地下水模型进行三角形或矩形差分构建基本单元[4]。

在计算过程中，两模型之间实际上是分离的，地表水和地下水模型分别独立进行计算，模拟实际情况下的产流和地下水动态变化，模型之间只是通过下渗量的传递来进行简单、松散的耦合[5]，如SWATMOD[6]模型、MODFLOW-DAFLOW[7]模型等。

这就导致模型不能准确地模拟土壤水，特别是平原区饱和-非饱和带土壤水的运动过程，对耦合模型要重点解决的诸如潜水蒸发、土壤水变化及产流量等问题不能很好地解决。

程勤波等[8]在可变坡土槽进行了室内人工降水入渗试验，以此来分析MODFLOW-UZF1模型模拟降水入渗补给过程的精度。

结果表明，模型不能准确地模拟入渗过程中土壤含水量的渐变过程，忽略了非饱和带土壤水分运移中基质势的影响，在地下水浅埋地区模拟效果较差。

本文基于三维理查德方程推导并构建了适用于平原区三维饱和-非饱和土壤水流动模型，并在太湖流域的金坛市建立野外实测试验基地，以此来观测分析和研究地势平坦、地下水埋深浅的特殊条件下的土壤水运动和浅层地下水水位的变化，以及土壤水分对降水和蒸发的响应。

模型以土水势为研究变量，首先对研究区域进行网格剖分，利用有限差分的格点法将研究对象离散化，进而模拟平原区浅层土壤含水率以及非饱和区与饱和区的变换过程。在模型求解过程中，运用“边循环”的概念和矩阵标识法，进一步提高了模型的精确度和计算效率。

参考文献：

[1]DEMETRIOU C,PUNTHAKEY J F.Evaluating sustainable groundw-ater management options using the MIKE SHE integrated hydrogeologi-cal 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 MOD-FLOW-2000 for modeling groundwater dynamics in shallow water tableareas[J].Journal of Hydrology,2012,412(1):170-181.

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

饱和—非饱和渗流三维数学模型及数值方法_（曹 渊，2013）

0 引言

        传统的地下水渗流模拟主要包括饱和流模拟和非饱和流模拟。

        无论是在饱和流模拟(以潜水面为上边界)还是在非饱和流模拟(以潜水面为下边界)中，都不可避免地遇到潜水面上下波动给数学建模带来的较大困难。

        实际工作中，很难将潜水面作为边界条件预先确定下来。

        将饱和流和非饱和流统一起来模拟，则可以把潜水面置于整个研究区域之中，不再作为边界条件，从而避免了边界面上下波动所带来的建模难题。

        尽管饱和-非饱和渗流模拟的难度较两者单独模拟大了很多，但研究这一问题仍具有重要的实用意义。

        目前国内外已经开展了饱和-非饱和渗流的数学建模研究，但主要采用压力水头的表示形式[1-4]。该形式的模型在数值求解时会有较大的数值振荡和数值弥散。

本文首先建立了以总水头表示的统一形式的饱和-非饱和渗流三维控制方程，给出了主要耦合参数的计算方法，简化了控制方程的性质，使得数值求解变得容易；其次，选取了合适的有限差分格式，确定了有效的数值求解方法；最后，通过实际算例验证了模型和解法。

参考文献

[1] Neuman S P,Narasimhan T N,Witherspoon P A.Ap-plication of Mixed Explicit-implicit Finite ElementMethod to Nonlinear Diffusion-type Problems[A].Finite Elements in Water Resources[C].London:Pen-tech Press,1977.1153-1186.

[2] 朱军,刘光廷,陆述远.饱和非饱和三维多孔介质非稳定渗流分析[J].武汉大学学报(工学版),2001,34(3):5-8.(Zhu J,Liu G T,Lu S Y.Saturated-unsatu-rated unsteady seepage analysis of porous medium[J].Engineering Journal of Wuhan University,2001,34(3):5-8.(in Chinese))

[3] 朱岳明,龚道勇.三维饱和非饱和渗流场求解及其逸出面边界条件处理[J].水科学进展,2003,14(1):67-71.(Zhu Y M,Gong D Y.Solution to 3-D unsteadysaturated-unsaturated seepage problem and accuratetreatment of saturated and unsaturated exit surfacesof seepage[J].Advances in Water Science,2013,14(1):67-71.(in Chinese))

[4] 薛禹群.地下水动力学原理[M].北京:地质出版社,1986.(Xue Y Q.Groundwater Dynamics[M].Bei-jing:Geological Publishing House,1986.(in Chi-nese))

[5] 杨金忠,蔡树英,王旭升.地下水运动数学模型[M].北京:科学出版社,2009.(Yang J Z,Cai S Y,Wang XS.Groundwater Flow Models[M].Beijing:SciencePress,2009.(in Chinese))

[6] 王洪涛.多孔介质污染物迁移动力学[M].北京:高等教育出版社,2008.(Wang H T.Dynamics of FluidFlow and Contaminant Transport in Porous Media[M].Beijing:Higher Education Press,2008.(in Chi-nese))

[7] 薛禹群,谢春红.地下水数值模拟[M].北京:科学出版社,2007.(Xue Y Q,Xie C H.Numerical Simulationfor Groundwater[M].Beijing:Science Press,2007.(inChinese))

[8] 陆金甫,关治.偏微分方程数值解法[M].北京:清华大学出版社,2009.(Lu J F,Guan Z.Numerical Solu-tion of Partial Differential Equations[M].Beijing:Ts-inghua University Press,2009.(in Chinese))

区域尺度饱和-非饱和水分及溶质运移集成模型研究与应用 

1.1  研究背景及意义

        非饱和带是大气水、地表水同地下水进行水分及能量交换的枢纽，是陆地植物赖以生存的源泉，也是地表污染物进入地下水的通道，在整个水文循环体系中具有十分重要的意义。

从区域范围来看，非饱和带仅仅是薄薄的一层（厚度通常在几米至几十米之间）覆在饱和含水层之上，其空间分布尺度远远大于自身垂向厚度。

但是，在此极薄的非饱和带中，具有非常复杂的物理、化学和生物过程，且变化极为频繁，是地下水文系统中最为活跃的部分。

因此，对于区域地下水资源和浅层地下水环境保护而言，非饱和带中的土壤水分及溶质运移是关键的控制要素。同样，非饱和带的状况也不容乐观。

土壤水分运动规律极其复杂，因此合理描述土壤水分运动规律及其与地下水耦合运动规律、高效准确地评估土壤水及溶质对地下水及溶质的影响具有重要的实际意义。

1.2  国内外研究现状

        数值模拟方法是目前研究土壤水地下水水分与溶质运移问题的主要工具和手段。

对于饱和带地下水分及溶质运移的模拟方法，现阶段已经较为成熟，以有限差分法和有限元法应用最为广泛（王军进等，2018）。

对于非饱和带土壤水分及溶质运移的模拟，主要为基于Richards方程和对流弥散方程的数值方法和各种简化方法。

非饱和带土壤水分及溶质运移过程非常复杂，近几十年来虽然有大量的研究成果，但学者们仍在试图构建更加健全、高效的数值模拟方法。

对于饱和-非饱和水分及溶质运移耦合模型，其理论与应用尚处于不断发展完善的过程中。

区域尺度土壤水盐运移规律的应用研究对于土壤盐碱化等问题具有重要的意义，尚有大量问题有待解决。

1.2.1  非饱和带水分及溶质运移模型研究进展

        土壤水分运动的理论起源于Darcy定律。

Richards (1931)将Darcy定律和能量概念应用于非饱和土壤水的流动，推导出多孔介质中水流运动的基本方程Richards方程，是土壤水运动定量研究的基础。

土壤溶质运移依赖于土壤水分运动，在土壤水分运动方程的基础上，许多科学家通过系统的研究，提出了对流-弥散方程(Advection Dispersion Equation，简称ADE) (Nielson and Biggar, 1961; Biggar and Nielson, 1962; Gardner, 1965; Bresler, 1973; Bear, 1972)。

以求解Richards方程和ADE方程为基础的数值方法是目前研究小尺度(实验室尺度或田间尺度)土壤水分及溶质运移规律的主要方法。

        原始的Richards方程存在两个变量，以水头h来表示势能，以含水量θ来表示土壤含水量，因此称为混合型Richards方程。

对混合型Richards方程的求解，需要引入本构关系表达式，即土壤水分特征曲线，将Richards方程中的两个未知量进行统一，最终得到以水头h为唯一自变量的水头型Richards方程或以含水量θ为唯一自变量的含水量型Richards方程 (查元源, 2014; 杨金忠等, 2016)。

        水头型Richards方程因其对各种土壤类型和水流条件良好的适应能力而得到广泛的应用，但是该类型的方程容易引起较大的总体质量误差。

较多学者提出了不同的改进数值格式并取得了良好的效果 (Celia et al., 1990)。

采用水头型Richards方程的模型有HYDRUS (Šimůnek et al., 2008)、SWAP (Kroes et al., 2017)、DAISY (Hansen et al., 2012)、MIKE-SHE (Graham and Butts, 2005)等。

含水量型Richards方程具有无条件质量守恒的优点，使得该类型方程在模拟非饱和土壤水分运动时具有很大的优势。

但是该类型方程不能描述土壤饱和的状况。

针对混合型Richards方程，存在“主变量转换技术”(Kirkland et al., 1992; Forsyth et al., 1995)、“方程切换技术”(Zeng et al., 2018)等求解方法。但是目前为止，以水头型Richards方程的应用最为广泛。

        Richards方程具有强烈的非线性性，输入参数多且测量困难。

此外，非饱和带在大气边界的作用下，含水量变化剧烈，使得Richards方程的数值解法对时空步长有严格的要求，计算成本较高。

因此，现阶段Richards方程模型较多用于小尺度的模拟，对区域问题的研究则较少。在区域尺度的模型中，基于均衡原理的简化模型应用更为广泛。这类模型也称为概念模型，其基本原理为采用一个个相互连通的存储单元来描述各水文过程，采用基于经验或半经验的物理公式描述各单元中的水分运动过程 (Devia et al., 2015; Li et al., 2015)。概念模型具有较高的计算效率与数据适应能力，模型稳定性良好。

        对于土壤水分运动过程的描述，概念模型一般采用“tipping-bucket”的方法。该方法将土壤视为一个个垂向相互离散的串联的“水桶”，当上部的“水桶”装满了水，就会向下倾泻。Ranatunga et al. (2008)根据模型的复杂程度，将基于物理机理的描述非饱和带土壤水分运动过程的概念模型分为两类，分别为单层水均衡模型和多层水均衡模型，如图1.1所示。单层模型仅采用一个“水桶”来描述非饱和带土壤水分运动，该“水桶”中的水量随着降雨而增多，随着蒸散发作用而减少。这类模型一般不考虑作物，因此并不区分蒸发与蒸腾作用。当代表土壤或者根系层的“水桶”处于饱和状态时，多余的入渗水将由地表径流排出。

        单层模型一般不考虑水平向和垂向的土壤变异性。此类模型产生较早，物理概念清晰，但是结构过于简单，采用这种策略的水均衡模型有WATBAL (Yates, 1996)、AWBM等，这类模型常常用于预测地下水补给和简单计算浅根作物根系层土壤水分动态 (Eilers et al., 2007)。多层模型可以考虑垂向土壤的非均质性与不同层含水量的变化，通常情况下考虑作物根系并将蒸散发划分为土壤蒸发与作物蒸腾。此外，多层模型还需要考虑入渗水分在不同层之间的分配问题。

        根据对入渗水分配过程的处理，可以再将多层模型分为两类，即简单活塞模型和包含再分配过程的模型。入渗水分进入土壤后，根据“tipping-bucket”方法从上至下对每层土壤进行填充。假设入渗水分直接将每层的“水桶”填充至稳定状态而不存在水分再分配过程时，此类水均衡模型称为简单活塞模型，如VSMB模型 (Baier and Robertson, 1996)。通常情况下取田持含水量作为稳定条件下的含水量。包含再分配过程的模型则将入渗水分的分配分为两个过程，即入渗水分的分配和水分的再分配过程。入渗水分先采用简单活塞模型的方式将各层土壤含水量填充至饱和含水量，再在重力势的作用下进行向下的排水过程，各层饱和含水量逐渐向田持逼近。采用此方案的模型有SWAT (Arnold et al., 2012)、INFIL 3.0 (FILL, 2008)、SoilWat (Holzworth et al., 2014)、DSSAT (Jones et al., 2003)、BUDGET (Raes, 2002)、Hydrobal (Bellot and Chirino, 2013)、SWB (Westenbroek et al., 2010)等。

        土壤水分所具有的势能称为土水势，一般认为非饱和带土壤水分总的土水势由重力势和基质势组成 (雷志栋等，1988)。上述的传统土壤水分运动均衡模型假设重力势远大于基质势，因而忽略基质势的作用，采用“tipping-bucket”方法描述土壤水分在重力势作用下的运动。但是该类模型并不能描述非饱和带土壤水分垂向向上的运动过程，因此只适用于非饱和带厚度较大且蒸散发作用较小的情况 (Mao et al., 2018)。此外，当土壤剖面存在土质分层时，土壤含水量会在不同土质的交界面处发生突变 (Yang et al., 2009)，传统均衡模型则避开了该问题，如APSIM (Holzworth et al., 2014)将整个剖面作为均质土壤描述。

        综上，非饱和带土壤水分运动极其复杂。Richards方程作为处理非饱和带水流运动问题的基本控制方程，具有严格的数学与物理基础。然而，Richards方程需要详尽的输入参数和较高的计算成本，限制了该方程在区域尺度实际问题的应用。基于均衡原理的简化模型由于其较高的计算效率、较好的数据适应能力、良好的模型稳定性而在区域问题中得到了广泛应用。现有的描述土壤水分运动的均衡模型主要基于“tipping-bucket”方法。然而，“tipping-bucket”方法仅能描述非饱和带土壤水分在重力势作用下的向下流动，而忽略了可能发生的土壤水分在基质势梯度下的运动，从而极大限制了模型的适用性。因此，在对区域土壤水分运动过程准确描述的基础上，对求解过程进行适当简化，发展适合于区域问题的土壤水分运动模型具有重要的意义与实用价值。

        土壤溶质运移研究的为溶质在土壤中的运移过程、规律和机理 (李保国等，2005)，土壤溶质的运移不仅与土壤水的性质有关，也与其依附的固体颗粒的物理化学生物过程有紧密的联系 (Nielsen et al., 1986)。多孔介质结构非常复杂，因此存在不同的研究方法与理论模型。根据对溶质运移过程描述的不同尺度，可以简单分为较为微观的毛管模型、平均化的对流弥散模型和宏观的随机模型。毛管模型应用相对较少，不再详述。

        随着研究的深入，许多科学家提出传统的对流弥散方程。Nielson and Biggar (1961; 1962)以及Biggar and Nielson (1962; 1963)根据实验结果提出混合置换理论，认为溶质通量是由对流、扩散和弥散的综合作用引起的，首次系统地论述了对流弥散方程的科学性和合理性 (黄领梅和沈冰，2000)。Gardner (1965)与Bresler (1973)对多孔介质与溶质间的相互作用进行了总结，并以费克定律和质量守恒原理为基础，建立了一维溶质运移方程。Bear (1972)提出一种平均化的方法来推导对流弥散方程，该方法物理概念清楚，对出现在方程中的弥散系数定量给出了与水流特性和介质特性的关系，采用了空间平均，脱离了较死板的几何体系和统计方法。对流弥散方程是溶质运移研究的核心数量模型。

        虽然ADE模型概念清晰，但模型限制因素较多，实际情况下计算误差较大；此外，ADE模型也无法解释溶质穿透曲线中的早期穿透与拖尾等物理非平衡现象 (解雪峰等，2016)。为了改进对流弥散理论而发展出了多种溶质运移模型，其示意图如图1.2所示。Coats and Smith (1964)和van Genuchten and Wierenga (1976)提出了两区模型以描述溶质运移中的物理非平衡现象，该模型将多孔介质孔隙中的水分分为不动水体与可动水体两部分，并引入两部分的比例系数与相互之间的质量交换系数2个参数。两区模型也称为可动-不可动水体模型，该模型考虑了多孔介质中不动水体的存在，能较好的描述土柱中不同位置溶质穿透曲线的提前穿透和拖尾特征，能更有效地模拟非均质介质中溶质在较大尺度上的运移过程，因此被广泛应用 (van Genuchten and Wierenga, 1976; Gao et al., 2009, 2010, 2013)。Skopp et al. (1981)在可动-不可动水体模型的基础上提出了两流区模型，两流区模型同样将多孔介质中的水分划分为两个区域，但是该模型主要依据水流速度进行不同区域的划分，且两个区域中的水流速度都不为零，如图1.2(c)所示 (王全久，2005)。两流区模型的不同区域分别采用独立的对流弥散方程描述其溶质运移过程，两区域之间可以发生溶质交换。两区模型根据孔隙尺寸划分不同区域，两流区模型根据流速划分不同区域。两流区模型与可动-不可动水体模型具有内在一致性，在特定情况下较为复杂的两流区模型可以转化为较为简单的可动-不可动水体模型 (马东豪和王全九，2004)。Gerke and van Genuchten (1993)更进一步提出了双重渗透模型。两流区模型是针对饱和的情况而言，而在双重渗透模型中的两个区域均可以为非饱和流，两区域均采用Richards方程求解其水分运动，再通过ADE方程求解其溶质运移过程，水分和溶质均可以在两区域之间相互交换。此外，还有更加复杂的三重渗透模型、流管模型等等。Köhne et al. (2009)在总结了大量已有模型的实际应用之后指出，对于田间及更小的尺度，如果仅考虑水流运动，除非有非常大的孔隙，一般而言，Richards方程的结果与分区类模型并没有较大差别，但是对于溶质而言，分区模型结果相对更好。但是，多分区模型包含有更多的参数，而在实际应用中这些参数难以通过试验或反演等方法实现，限制了分区类模型的实际应用。在流域尺度，由于数据的缺乏等原因，分区类模型尚难以实际应用。

        此外，众多国内外学者在过去的二三十年中进行了大量关于随机理论的基础研究工作 (Freeze, 1975; Jury, 1982; Schwartz et al., 1983; Gelhar and Axness, 1983; Gelhar, 1986; Dagan, 1987, 1991; 杨金忠等, 2000; Neuman and Federico, 2003; Yang et al., 2004; Dagan and Neuman, 2005; 史良胜, 2009)。随机方法提供了一种新的研究思路，但大多数研究均限于理论研究阶段。其他求解多孔介质溶质运移的方法有分数阶对流扩散方程FADE (Krepysheva et al., 2006; Meerschaert et al., 2006)、连续时间随机游走理论CTRW (Berkowitz et al., 2006; Bijeljic and Blunt, 2006)、渗透理论和其他方法 (Park et al., 2005; Raoof and Hassanizadeh, 2010; Dentz et al., 2011)。

1.2.2  饱和-非饱和水分及溶质运移耦合模型研究进展

        地表水、土壤水和地下水作为水循环系统的主要过程，长期以来却在各自相对独立的领域中发展 (王蕊等，2008; 曾献奎等，2009)。非饱和带土壤水分运动非常复杂，需要大量计算且常常缺少相应的数据，使得大尺度的地表水文模型与地下水运动模型均对非饱和带的土壤水分运动有大量的简化。地表水文模型一般过度简化非饱和带的土壤水分运动，且非常少考虑地下水的三维流动 (Shen and Phanikumar, 2010)；区域的地下水模型常常将地下水补给项等作为外部源汇项，从而避开非饱和带的计算 (Twarakavi et al., 2008)。

        随着研究的推进与实际工程的要求，需要了解非饱和带中水分和溶质分布分局部特征，越来越需要将地表水、土壤水与地下水作为整体进行研究，特别是应用于大区域地下水模型中的非饱和带土壤水分运动的计算问题。饱和带与非饱和带具有统一的控制方程，即Richards方程。因此，一些模型统一求解三维变饱和的Richards方程，以得到饱和带与非饱和带精确的土壤水分运动状态 (Pikul et al., 1974)，如FEFLOW (Diersch, 2013)、HydroGeoSphrere (Brunner and Simmons, 2012)、InHM (VanderKwaak and Loague, 2001)、MODHMS (Werner et al., 2006)等等。但是，对于大区域的饱和-非饱和带水分运动的纯三维数值模拟，一方面由于非饱和带水力特性变化剧烈，使得模型面临着节点众多、时空步长较小、计算成本过高的问题 (Van Walsum and Groenendijk, 2008; Shen and Phanikumar, 2010)；另一方面由于Richards方程非线性参数较多，实际情况中难以获得区域上的参数，限制了模型在实际中的应用。较为普遍的方式是对土壤水-地下水运动过程采用不同程度的概念简化，并将二者进行耦合。最为常用的是拟三维的计算方法 (Kuznetsov et al., 2012)。从物理过程和实际观测而言，相对于饱和带，非饱和带垂向尺度一般较小。当非饱和带不包含毛管上升区且不存在上层滞水时，非饱和带的侧向流比较微弱，一般可以将非饱和带的水流运动简化为一维垂向运动 (Markstrom et al., 2008; Kuznetsov et al., 2012; Zhu et al., 2012; Zhu et al., 2016)。从数值计算讲，饱和含水层运动缓慢，常常采用较大的时空步长进行计算，而非饱和带土壤水分变化剧烈，常常采用很小的时空步长（一般为s和cm）(朱焱, 2013;刘昭, 2016)，因此二者的时空运移尺度不匹配，造成了一系列数值求解的困难。Hatton et al. (1992)与Vertessy et al. (1993)即指出，许多澳大利亚的非饱和带水分均衡分析并不需要对整个区域进行完全三维的建模分析，而是将非饱和带近似为一维垂向运动 (Zhang et al., 1996; Ranatunga et al., 2008)。

        Pikul et al. (1974)和Niswonger et al. (2006)指出，拟三维方法是求解地下水模型，特别是大尺度问题时最有效率的方法。饱和-非饱和模型的耦合方式分为3种，分别为松散耦合、迭代耦合和全耦合 (Furman, 2008)。在一个时间步中，分别计算饱和带模型与非饱和带模型，并在二者界面交接处单向传递物理量的方法称为松散耦合。如果在二者交界面处，双向反馈传递物理量并迭代计算，直至饱和带模型与非饱和带模型在界面交界处的物理量相互匹配，则称为迭代耦合。全耦合则是将饱和带模型与非饱和带模型作为一个整体的系统统一求解，而不再对二者进行区分。Pikul et al. (1974)将垂向一维Richards方程与水平向一维的Boussiniewq方程耦合用以求解二维饱和-非饱和地下水流问题；Skaggs (1978)开发DRAINMOD模型，在非饱和带采用水均衡的方法处理水分运动，在饱和带采用Hooghoudt解和Ernst解，将两者耦合求解饱和-非饱和水流运动问题；Feddes et al. (1982)开发了SWATR模型，该模型将非饱和带的一维有限差分解与饱和带的Ernst解耦合求解饱和-非饱和水流运动问题。DRAINMOD和SWATR模型均为广泛应用的数值模型，但是它们一般适用于点尺度的求解，主要用于非饱和带与排水系统的水分运动的求解。随着计算机技术的发展，许多计算更加精细的模型被开发和应用。Havard et al. (1995)建立了LINKFLOW，该模型在非饱和带采用有限差分法计算一维垂向Richards方程，将之与饱和带的三维模型MODFLOW进行耦合。非饱和带的厚度为地表到地下水位，地下水位由饱和带模型控制，非饱和带Richards方程的求解采用predictor-corrector方法求解。

        Facchi et al. (2004)将土壤水运动模型SVAT与MODFLOW松散耦合在一起，采用GIS管理研究区所有的相关信息，用以模拟饱和-非饱和土壤水分运动。其中，SVAT模型是非饱和带的概念模型。Niswonger et al. (2006)在MODFLOW-2005的基础上开发了UZF1程序包用以模拟非饱和带的土壤水分运动及其对地下水的补给量。该模型假设非饱和带水流运动只有重力势驱动，忽略压力梯度，因此采用行波法求解简化后的Richards方程，采用Brooks-Corey模型描述非饱和水力传导度与含水量的关系，需要指定蒸发潜在深度。GSFLOW模型亦采用UZF程序包处理非饱和带水流运动 (Markstrom et al., 2008)。UZF1程序包需要指定蒸发潜在深度，不能处理非均质的问题，且仅考虑向下的水流运动，这些均极大限制了该程序包的实际应用。Seo et al. (2007)基于MODFLOW-2000开发了HYDRUS-1D程序包。该程序包以HYDRUS-1D代码为基础，移除了溶质运移与热量传输项，忽略迟滞现象，且去除了与耦合模型不相关的边界条件，将简化过的模型与MODFLOW进行松散耦合。Twarakavi et al. (2008)和Bailey et al. (2013)通过算例对比了UZF1与HYDRUS程序包的差别。Leterme et al. (2013, 2015)将该程序包应用在区域问题中，并对比了与UZF1的计算结果的差别。目前，加州大学河湾分校、比利时核研究中心与波兰格但斯克工业大学合作，采用最新版本的HYDRUS与MODFLOW-2005进行耦合，并希望加入溶质运移、优先流等更多的功能 (Szymkiewicz et al., 2018)。Beegum et al. (2018)提出了修正由于HYDRUS与MODFLOW模型在饱和-非饱和交界面水头不匹配而引起的虚假流量的方法。Van Walsum and Groenendijk (2008)开发了一个松散耦合的饱和-非饱和模型，该模型在非饱和带采用稳定解获得剖面含水量分布。Xu et al. (2012)采用SWAP模型描述非饱和带土壤水分运动，并与三维地下水模型MODFLOW-2000松散耦合。

        但是，以上所述各模型所采用的松散耦合方法存在质量误差的问题 (Kuznetsov et al., 2012)。Twarakavi et al. (2009)提出了新的改进HYDRUS程序包，该程序包在每个HYDRUS计算时间步末均更新非饱和带的厚度，以此达到质量均衡的目标。但是，这种处理导致该模型难以加入溶质运移项。为了耦合模型的质量均衡，不同的模型采用了不同的方法，最主要的是采用迭代计算的方法。Yakirevich et al. (1998)采用迭代的方式构建饱和带-非饱和带拟三维水流与溶质运移模型，将一维的Richards方程与二维的潜水方程进行耦合。Kuznetsov et al. (2012)进一步发展了该模型，其对地下水的模型采用完全三维的形式进行描述。Stoppelenburg et al. (2005)采用迭代方法将饱和带的LGM模型与非饱和带的SWAP模型进行耦合。LGM模型为拟三维 (水平向含水层二维，垂向弱透水层一维) 地下水模型，采用有限元进行离散求解。LGM模型与SWAP模型在交界面存在相互重叠的部分，因此两个模型不是在地下水位处相耦合，而是在第一层弱透水层处，SWAP模型的下边界条件选择为柯西边界，如图2.3所示。SWAP根据下边界进入的垂向流量，计算得到非饱和带对饱和带的补给流量（通过zi-1的通量与zi的通量的线性平均得到），并根据地下水位之上的单元计算给水度，将这两个量传入LGM作为其上边界条件。该方法也存在缺点，一是在两个方程水平上的迭代收敛在计算过程中非常难以实现；二是该模型迭代方法建立在特殊的地质构造上，即上层潜水，之下存在弱透水层，再之下是承压含水层，该地质条件限制了此模型的应用 (Shen and Phanikumar, 2010)。

        Shen and Phanikumar (2010)开发了基于物理基础的、可以用于中尺度 (~1000 km2) 和大尺度流域 (>5000 km2) 地表水与地下水相互耦合计算的PAWS模型。该模型对非饱和带求解垂向一维Richards方程，对饱和带求解二维水平向饱和带方程，并且有两个假设。其一是假设非饱和带中的土壤水分运动均为一维垂向运动，忽略侧向流；其二是侧向流在垂向上均匀分布在潜水含水层，如图2.4所示。该模型的计算过程为，首先从上一个时间步计算结果中提取地下水位和侧向流动通量DR；求解处于稳定情况下饱和带的解；根据地下水位的分布，求解非饱和带一维Richards方程，获得对饱和带的补给量；根据该补给量求解二维饱和带方程，重新获得地下水位与侧向流动通量DR，返回第一步重新计算。该模型避免了在方程层面上的迭代，且将三维系统拆分为低维系统，更易收敛。

        Zeng et al. (2019)提出了HYDRUS与MODFLOW的耦合迭代模型，将HYDRUS计算得到的水分通量作为二类边界条件传递给MODFLOW模型，将MODFLOW模型计算得到的水位信息作为一类边界条件传递给HYDRUS，并将二者进行迭代耦合计算，直至模型收敛。由于饱和带与非饱和带运动特征的差异，采用完全耦合方式的模型相对较少。Zhu et al. (2012)构建了完全耦合的拟三维饱和-非饱和水分及溶质运移模型Q3D (朱焱，2013)。该模型采用一维Richards方程描述水分在非饱和带的运动，而对于饱和带，将三维流速场分解为水平向二维运动和垂向一维运动，采用有限单元法求解水平向二维流速场，采用有限差分法计算垂向水流通量将各含水层水流运动进行耦合。该方法收敛性较强，结果良好，且质量可控，但编程与求解过程过于复杂，难以应用。相较于区域饱和-非饱和水流运移模型，区域饱和-非饱和溶质运移模型的研究相对较少。Morway et al. (2013)和Bailey et al. (2013)先后修改了MT3DMS和RT3D，并将之与MODFLOW及UZF程序包耦合在一起，用以计算饱和-非饱和水分及溶质运移过程。其建模思路为，将UZF与MODFLOW计算得到的水分运动信息传递给修改后的MT3DMS或RT3D，在溶质运移的计算中并不区分饱和带与非饱和带。Wu et al. (2016)则在前述模型的基础上进一步扩展了模型对溶质运移的模能力，并将新模型命名为PHT3D-UZF。Beegum et al. (2018)采用了松散耦合的方法将HYDRUS-1D与MODFLOW结合在一起，并加入了溶质运移的部分。相比于之前的HYDRUS程序包，新的耦合方法采用了更新耦合区域压力水头与溶液浓度的方法以保证边界通量的稳定和质量守恒。Zhu et al. (2013)采用完全耦合的方法开发了饱和-非饱和水分及溶质运移数值模型Q3D。此外，尚有MIKE-SHE (Refsgaard and Storm, 1995)、MODFLOW-SURFACT (Panday and Huyakorn, 2008)等商业软件可供选择。

        综上所述，现阶段虽然已有较多的各种类型的模型用以处理饱和-非饱和水分与溶质运移问题，然而尚存在一些挑战。主要有以下几个问题：首先严格控制质量守恒的饱和-非饱和水分及溶质运移模型尚需进一步发展与完善；其次是模型的稳定性问题，对于非饱和带模型而言，现有的饱和-非饱和模型多采用Richards方程描述非饱和带土壤水分运动，但是该方程需要大量的难以测量的输入参数，且在某些情况下会出现迭代不收敛、水量不均衡等问题，极大影响了模型在区域问题中的实际应用。对于饱和非饱和模型耦合方案的问题，现有大多数模型采用非迭代耦合方案，虽然不存在耦合模型的稳定性问题，但是却导致水分模型质量误差不可控，从而难以更进一步构建溶质运移模型。一些模型选择基于严格的物理基础的方程层面的迭代耦合方案，但是在实际应用过程中难以收敛 (Shen and Phanikumar, 2010)。最后，现有的模型更加关注具体的算法，较少关注模型的结构设计和前后处理模块的开发，使得模型在扩展性、实用性等方面均有所欠缺。

华北平原地下水系统变化规律研究（曹国亮）
 

5.1 区域地下水补给计算方法

        在干旱、半干旱地区，农业依赖于灌溉。因此，地下水补给量的准确估算和蒸发量的计算对实现地下水资源的可持续管理极为重要。然而，由于土壤质地和植被类型的空间变化，以及气象因素，如降水和蒸发的时空变异，地下水补给量的估算通常不容易做到。在灌溉的条件下，地下水的补给过程会更为复杂(Gee和Hillel，1988)。传统的测量地下水补给的物理方法(如蒸渗仪和渗透仪lysimeter，infiltrometer)和测量蒸散发的物理方法(如涡度相关和蒸渗仪)，以及诸多经验方法通常要求大量的试验数据，因此往往局限在小尺度应用。区域地下水补给量的估算通常采用的方法有：水均衡法、潜水位动态法(WTF)、化学和同位素方法以及数值模拟方法(Simmers，1987；Healy和Cook，2002；Scanlon等，2002；Xu和Beekman，2003)。仅应用某种方法估算地下水补给量的准确度一般很难判定，因此建议同时采用多种方法以对不同方法的结果互相验证以降低结果的不确定性(Scanlon 等，2002)。另外，由于气象因素、土地利用类型以及地表物理和生物特性的空间变异，将点状计算结果插值到区域上往往造成较大误差(Scanlon等，2002)。

        因此，对地下水补给估算方法的选择取决于研究的空间和时间尺度，以及研究区的数据情况(Scanlon等，2002；Delin等，2007)。估算区域地下水补给量的主要方法为土壤水量平衡法。其中，土壤带或根植层底界的地下水流出量作为净入渗量(net infiltration)或潜在补给量(potential recharge)。基于该种方法的数值模型很多，例如美国地质调查局(USGS)的INFIL模型(USGS，2008)，BALAN模型(Samper和García-Vera，1992)，WAVES模型(Zhang和Dawes，1998)，和美国环保局(U.S. EPA)的HELP模型(Schroeder等，1994)。虽然这些模型能够简便快速的估算地下水补给量的区域分布，前提条件是地下水潜水面要低于根植层的底界，并假设根植层底部的渗出量即为地下水补给量(潜在补给量=实际补给量)。但是认为自根植层底界至潜水面的非饱和带内的水流过程是可以忽略的是有问题的。当非饱和带或包气带厚度较大时，非饱和带将成为控制补给时间和水分传输的主要因素，忽略其中的水流过程是不可取的(Harter和Hopmans，2004；Hunt等，2008)。

        鉴于以上各种方法的优缺点，非饱和带数值模拟方法逐渐成为估计区域地下水补给时空分布的通用方法(例如，Keese等，2005；Small，2005；Liggett等，2009)。此类模型主要包括DAISY(Hansen等，1990)，UNSAT-H(Fayer，2000)，SWAP(Kroes和van Dam，2003)和HYDRUS-1D(Šimůnek等，2005)。虽然这些模型能够模拟非饱和带中的水流并估算蒸发量和地下水补给量，其应用过程中通常需要指定模型下边界条件(如潜水面位置)，没有耦合饱和带地下水流动过程。因此，其结果不能直接反映出气候变化、土地利用和土壤质地对地下水动态的影响。近年来，开发的耦合非饱和带和饱和带水流过程的数值模型，如HYDRUS-3D(Šimůnek，2006)，MODFLOW-VSF(Thom等，2006)和HydroGeoSphere(Therrien等，2005)，是对上述非饱和带数值模型的巨大改进。然而，这些模型的理论基础是刻画非饱和带水流过程的理查德(Richards’)方程，其应用主要受到计算机计算能力的限制(van Walsum和Groenendijk，2008)。而对于流域尺度的模拟计算，将非饱和带内水分运动假设为一维流动通常是合理的(Harter和Hopmans，2004)。因此，许多模型，如MIKE-SHE(Refsgaard等，1995)和MODFLOW-UZF(Niswonger等，2006)采用了非饱和带一维流动和饱和带三维流动的耦合方式以降低耦合模型的复杂程度和提高计算效率。在这些模型中，一维的理查德(Richards’)方程被进一步近似成运动波(kinematic-wave)方程并采用特征值方法求解。该求解方法的优点在于不需要对非饱和带进行网格剖分(Niswonger等，2006)。

        基于非饱和带水流过程模拟估算区域地下水补给量的主要困难在于模型参数的获取和尺度转化，如土壤的渗透系数。对于该过程，Harter和Hopman(2004)对多种地质统计方法提升土壤渗透系数的空间尺度做了细致的综述性讨论，但这些方法是基于大量的小尺度范围内的渗透系数估算结果(Vereecken等，2007)。

        作为直接测量原状土壤渗透系数的替换方法，土壤转换函数方法(PTFs)(Wösten，2001)能有效的用于区域非饱和带渗透系数的估算。利用土壤转换函数方法可以直接测量土壤特征指标(通常附带于土壤调查数据库)并估算非饱和的水力参数(Li等，2007)。许多广泛应用的土壤转换函数方法是基于统计模型的方法(Pachepsky和Rawls，2004)估算土壤水特征函数和饱和渗透系数。所利用的土壤特征指标通常为土壤颗粒分级，土壤干密度和有机成分含量(Nemes等，2003)。

        绝大多数非饱和带数值模拟软件最初的开发目的是用于农业领域研究工作，需要大量反应植被生长特征的参数。虽然土地利用和植被类型分布可以通过调查统计直接获得，遥感解译也有助于获取植被参数，但对于植被类型和土地利用不是地下水补给主控因素的地区，简化植被生长过程模拟，减少模型参数更有利于减少模型的不确定性和模型校正过程中的不唯一性。

        蒸散发(ET)，作为土壤水均衡的主要排泄项，对利用包气带数值模拟方法计算地下水补给有着重要影响(环境同位素方法估计地下水补给不需要蒸散发的计算)。上述不同的模型中对蒸散发的计算过程也各不相同。最常用的计算方式首先是基于常规气象观测数据计算潜在蒸发或参考蒸散发(ETp或ET0)(Lu等，2005)，然后利用经验公式计算裸土蒸发和植物蒸腾(如Šimůnek等，2005)。实际蒸发和植物蒸腾根据土壤含水量情况和植被信息确定，并加和得到实际蒸散发(ETa)。但是，这种方法通常要求较小的时间尺度，如日气象数据。因此，饱和带水流模型的时间步长也需要相应减小，但对于饱和带地下水流、特别是大区域地下水流模型来说，显然是不合适的。作为一种折衷方案，可以利用子时间步长，即在饱和带水流模型时间步长的基础上再进行剖分，得到的更小时间步长作为非饱和带模型的时间步长，如耦合HYDRUS和MODFLOW的程序包(Seo等，2007)。

        基于互补关系的蒸散发模型(Bouchet，1963)是另外一类用于估算区域蒸散发的常用模型。这类模型仅利用常规气象观测数据，不需要刻画陆面边界的详细信息。因此，此类模型避免了复杂的土壤-植被-大气这一复杂系统，不需要土壤含水量数据和植被参数(Hobbins等，2001)。不同的气候类型地区与区域水均衡的蒸发计算结果对比已经证实了该类方法的适用性(如，Hobbins等，2001；Xu和Singh，2005)。此外，此类模型认为陆表大气湿度包括了所有陆面过程的影响，已经考虑了灌溉对蒸发的影响(例如，Ozdogan和Salvucci，2004)。