土壤水动力学SWD土壤水运动的数值解法和软件应用

土壤水动力学(SoilWaterDynamics)

第9章土壤水动力学的数值求解与软件应用

(Numericalsimulation&applicationofSoilWaterdynamics)毛晓敏中国农业大学水利与土木工程学院Telmail:maoxiaomin@第9章土壤水动力学

的数值求解与软件应用土壤水分运动的数值计算方法常用土壤水模拟软件应用举例9.1土壤水分运动的数值计算方法定解问题及常用数值计算方法Richards方程的差分离散边界条件处理土壤水分运动参数取值差分方程的求解9.1.1定解问题及常用数值计算方法定解问题：Richards方程边界条件初始条件求解方法：解析法、半解析法：数值方法：解析法：求出解的函数表达式:精度高，但有很大局限性:简化、特殊情况。数值法：求用数值表示的有限个离散点和离散时段上的近似解:具有广泛适用性；容易程序化，数值模拟软件多；随着计算机技术的提高和数值算法的改进得以迅速发展。常用数值计算方法：全离散化方法：时间和空间进行离散，利用微分方程推导出物理量在节点值的代数方程组，求解方程组得到微分方程在节点上的近似解有限差分法(FDM)：用差商替代导数控制容积法(CVM)：积分方程离散，守恒有限元(FEM)：借助基函数，直接求函数的近似解半离散方法：时间或空间离散有限解析法（FAM）线法(MOL)有限差分法－利用Taylor级数用差商代替导数：前向差分：后向差分：中心差分：一阶导数的近似表达式：

前差后差有限差分法的离散：i＝12。。。。。。。。Ij＝12

。

。

。

。

。

。

。Jfi,jΔx中心差若采用中心差的形式，二阶导数如何表示？有限元法－用简单的插值函数代替每个单元上未知函数分布，形成一系列代数方程组的方法

有限元法的离散：有限元法变化形式较多，如：1）从权函数的选择来说，有配置法、矩量法、最小二乘法和伽辽金法；2）从计算单元网格的形状来划分，有三角形网格、四边形网格和多边形网格；3）从插值函数的精度来划分，又分为线性插值函数和高次插值函数等。

fjGalerkin法引入基函数ψi

并采用ψi作权函数，使得：

举例—Galerkin法：对于：可用矩阵表示为：式中H、dH/dt分别为由h1，h2,…,hn和dh1/dt,dh2/dt,…dhn/dt组成的列矢量；D为导水矩阵，P为贮水矩阵。若采用隐式方法：通过求解大型稀疏矩阵方程，得到Hk+1数值计算的收敛性与稳定性收敛性：Δt、Δx→0时，差分方程解→微分方程的解稳定性：计算误差有界数值弥散问题：对流占优情况下容易出现数值弥散，利用特征线法等来克服9.1.2Richards方程的差分离散定解问题：初始：地表瞬时湿润时间差分：显式差分格式：α=0——不稳定中心差分(Crank-Nicholson)格式：α=0.5—可能不稳定隐式差分格式：α=1——稳定参数取值：时段初(β=0)时段中(β=0.5)时段末(β=1)——效果较好，需要迭代9.1.3边界条件处理一类边界：θ0=θs二类边界（通量边界）：直接差分半节点水量平衡虚拟节点9.1.4土壤水分运动参数取值时间：显式：用时段初含水率剖面计算参数预报校正：以时段初值求解方程，用新的含水率校正参数，再次求解迭代：半节点土壤水分运动参数：算术平均几何平均调和平均…9.1.5差分方程的求解对于非显式差分，得到方程：代入边界条件得到三对角方程组如下，用追赶法（Thomas算法）求解求解步骤求解土壤水分运动方程时，首先要给出有关土壤水分运动参数的表达式、初始条件与边界条件，然后按一定的时间、空间步长进行离散。在每一离散时段内，进行以下迭代计算：(1)以时段初的土壤含水率作为时段末含水率的预报值。(2)计算相关的土壤水分运动参数，利用追赶法求解土壤水分运动的差分方程组（三对角方程组），得到时段末含水率。(3)如果时段末含水率的预报值和计算值相差不大，可以认为迭代收敛，完成本时段计算，进入下一时段继续求解；否则，以时段末含水率的计算值作为新的预报值，返回步骤(2)继续进行迭代。在求解过程中，也可以根据各时段迭代情况对时间步长进行调整。一般情况下，如果迭代次数较少，表明时段内含水率变化不大，可以增大时间步长；否则，如果迭代次数较多，需要减小时间步长。这种自动调节步长的方法既可以保证模拟计算的精度，也可以有效地提高计算速度。

不同差分处理方法的模拟结果比较(a)C－N格式(b)隐式格式(c)特征差分格式降水入渗过程：降水强度0.2427mm/min，历时180min(a)中心差分格式(b)全隐式及指数加权隐式格式降雨入渗试验模拟结果（降雨强度0.05374cm/min，历时140min）9.2常用土壤水模拟软件软件概述：土壤水模拟软件：一般基于饱和-非饱和土壤水分运动方程，采用有限差或有限元进行数值求解。常用软件：（一般不仅仅局限于土壤水分运动模拟，有的也包括地下水运动、溶质运移、热迁移、甚至作物生长等。)一维软件：SWAP、HYDRUS系列、SHAW等二维/三维软件：SWAP由瓦赫宁根大学集成的SWAP(Soil—Water—Atmosphere—Plant)模型是由详尽完善的土壤水运移、溶质迁移、热量传输、土壤蒸发、植物蒸腾和作物生长子模块组成。HYDRUS1D、HYDRUS2D/3D系列软件是美国农业部盐土实验室开发的模拟饱和-非饱和介质中水分、热、溶质运移的模型。其中HYDRUS1D在该系列模型中主要进行垂向一维条件下的模拟。模型采用水头形式的Richards方程，因此，可以对饱和-非饱和区进行统一模拟，可以处理诸如有压入渗、地下水位线波动情况下的区域水流运动问题。SHAW是由美国农业部Flerchinger等人建立的，该模式是一个最初发展模拟土壤冻结和融化的一维模式，通过对土壤冻结、融雪、土壤温度、水量、蒸腾、水平衡影响的处理，

可以准确预测土壤的冻结／融化过程．二维/三维软件

1)专注于地表产汇流过程的模型，如SWAT；2）专注于干旱区区域水均衡的模型，如清华大学开发的干旱区散耗型水文模型；3）专注于区域地下水运动的模型，如FEFLOW，HYDRUS2D/3D，MODFLOW，SUTRA等。地表产汇流模型SWAT模型是近年来发展起来的分布式水文模型，可用于模拟预测各种管理措施及气候变化对水资源供给的影响，评价流域非点源污染现象。模型陆面模拟组成包括水文、气候、侵蚀、土壤温度、植物生长、营养物质、杀虫剂、以及土地利用和水资源管理。干旱区散耗型水文模型

由于描述流域水文循环的分布式水文模型多考虑了产汇流过程，并把产汇流模拟结果作为检验模拟准确性的重要标准。这与干旱区绿洲的情形是很不相同的。绿洲水文循环的特点是降水稀少、蒸发强烈，水资源多依靠上游河渠输水并经过各级渠道分配到田间，最终消耗于田间蒸散发。对于这种情况，清华大学开发了相应的干旱区散耗型水文模型。区域地下水运动模型

MODFLOW是由美国地质调查局

(USGS)开发的孔隙介质中三维地下水流数值模拟的模型，其源代码、执行文件和使用说明书都可以在其官方网站上免费下载。MODFLOW自从开发以来，已经经过数次更新，从起初的MODFLOW88更新为MODFLOW96、MODFLOW2000、MODFLOW2005等。以MODFLO为基础与地表水、溶质运移软件等耦合产生一些新的模拟软件，MODFLOW也被嵌入一些用户界面中（VP，GMS，PMwinetc.)，使其具有方便的前后处理界面。FEFLOW是由德国WASY水资源规划系统研究所研制开发的地下水模型软件包，可以进行复杂二维和三维饱和、非饱和带水流（包括密度流）、热量运移和溶质对流弥散的模拟，同时可以模拟溶质迁移过程中发生的线性非线性吸附反应和一阶动力学化学反应等。SUTRA由USGS开发，采用有限元法进行二维/三维、饱和/非饱和、变密度情况下地下水流动和溶质/能量运移的模拟]。我国已有报道介绍该软件，但应用尚不多见。国外多采用该模型进行海水入侵方面的模拟研究。土壤水模拟软件HYDRUS

HYDRUS1D、HYDRUS2D/3D系列软件

模型采用水头形式的Richards方程，可以对饱和-非饱和区进行统一模拟，可以处理诸如有压入渗、地下水位线波动情况下的区域水流运动问题。在给定充分供水条件下的地面蒸发和作物蒸腾量条件下，该模型可以模拟土壤水盐胁迫下的实际根系吸水（作物蒸腾）和地面实际蒸发过程。可以处理多种灵活的边界条件（大气边界、自由排水、定水位、动水位等），具有良好的用户界面和灵活的输入输出功能。模型采用有限元方法进行数值离散求解。9.3HYDRUS应用举例（1）

应用软件模拟前要考虑：计算区域和计算时段？土壤水分运动参数？如何给定初始条件？设起初为零通量稳定剖面，下边界含水率