1、HYDRUS-1D软件应用软件应用马文翠马文翠软件简介软件简介vHYDRUS-1D是美国农业部盐渍土实验室开发的模拟非饱和土壤中的水、热、溶质运移的软件，它在模拟土壤中水分运动、盐分、污染物(如农药)和养分(如土壤氮素)运移方面得到广泛应用。软件分为1维、2维、3维三种，分别命名维Hydrus-1D、Hydrus-2D、Hydrus-3D。由于非饱和土壤水主要是1维垂向运动的形式，因此Hydrus-1D的应用非常广泛。问题描述问题描述根系层1.细砂壤土2.中砂1 m1.5 m1.5 mv以天津一块农田为背景，对2012年冬小麦收割期农田尺度下土壤水分运动进行数值分析v土壤剖面从地面算起的深度：
2、300cmv土壤分层：分两层v准备模拟哪个时间段水分变化：030dv根系分布：1mv初始地下水位：地面以下6m软件界面软件界面前处理前处理后处理后处理模拟计算模拟计算主要过程主要过程水分运移方程水分运移方程v将坐标原点选在地面，取z轴向下为正，则一维饱和非饱和带水分运移基本方程为Richards方程：v式中：为含水率；k（）为水力传导度：h为负压水头；S 为植物根系吸水量，对裸露区为0。(1)hKStzz几何信息几何信息模型有2个土层长度单位是cm土壤剖面的深度是300 cm土壤分区长度 L分区水量 L流入或流出量 LT-1平均压力水头L实际表面通量 LT-1(入渗/蒸发:-/+)剖面底实际通
3、量LT-1(流入/出:+/-)整个区域的水质量平衡绝对误差L整个区域的水质量平衡相对误差%质量平衡信息质量平衡信息输出信息输出信息输出5组结果，每6天一组选择输出时间点迭代参数迭代参数最多迭代次数含水量迭代精度压力水头迭代精度增大步长迭代次数信号缩小步长的迭代次数信号步长增大比例步长缩小比例最小吸力间距最大吸力间距用于生成水分特征曲线的散点值上、下边界条件上、下边界条件定压力水头定通量大气边界，可积水大气边界，产流变水头变水头/通量定压力水头定通量变水头变通量自由排水深度排水渗透面水平排水初始条件使用压力水头使用含水量输入Etp和LAI消光系数最大积水厚度消光系数消光系数利用比尔定律计算消光系
4、数利用比尔定律计算消光系数:式中：TPAR是冠层底测得的瞬时光合有效辐射(molm-2s-1)，PAR是冠层顶测得的瞬时光合有效辐射，LAI是观测日的叶面积指数。利用冠层上、下测得的光合有效辐射资料，根据其变化规律拟合方程求积分，得到入射PAR的日总量(molm-2d-1)和透射的PAR日总量，再根据当日的叶面积指数来计算日平均消光系数。定义：定义：植物群体的日同化量大都表示为冠层截获光合有效辐射(PAR)量的函数，而群体对PAR的截获量可通过到达冠层顶部的PAR在通过作物群体时的衰减进行计算，衰减程度用消光系数(K)表示。它取决于太阳角度、植被类型及叶片空间分布特征。v 植被蒸腾量和土壤蒸发
5、量分开处理，HYDRUS-1D推荐使用一个经验公式来把潜在蒸散量分割为蒸腾潜力和土壤蒸发潜力：v 其中ETp为潜在蒸散量（cm/d），Tp为潜在蒸腾量(cm/d)，Ep为土壤潜在蒸发量(cm/d)，LAI是叶面积指数，k为消光系数。潜在蒸散量处理潜在蒸散量处理上边界通量上边界通量随时间变化边界条件随时间变化边界条件v 即地面土壤达到最干燥状态时的压力水头。理论上，土壤十分干燥时，吸力很大，而液态孔隙水的压强很小，与空气湿度保持平衡关系，有v hA为最小压力水头，Hr为空气绝对湿度，RT/Mg为空气的摩尔气体常数。公式需要近地面空气湿度，一般可取饱和水汽湿度，因2 cm以下土壤空气湿度往往饱和，
6、随温度发生变化，可据近地面气温变化推算地表土壤空气湿度，再换算成压力水头。HYDRUS中需输入hA的绝对值，缺省值为 v hCritA=|hA|=106 cmv 此数值只对土壤蒸发起作用。建议：hCritA 对应的土壤含水量应至少比残余含水量大0.005，在模拟根系吸水时，hA 还应低于P3。否则，当根系吸水的临界压力水头和地面蒸发的最小压力水头满足时,会导致回流现象，不合理。最小压力水头值最小压力水头值(hCritA)时间信息时间信息时间单位时间步长信息模拟30 d内的变化上边界是随时间变化的，每天一组数据，共30组数据自动处理蒸腾量在每天24小时内的变化其中Tp(t)是瞬时潜在蒸腾量，t为
7、时间。模型假设早上6点之前以及晚上18点-24点的蒸腾量总和只占全天蒸腾量的1%。HYDRUS-1D 可以使用一个经验公式来处理每天24小时潜在蒸腾量的变化，设某天的潜在蒸腾量为 （cm/d），则pT蒸腾量处理蒸腾量处理土壤水力模型土壤水力模型单孔介质模型双重介质模型双重介质模型双重孔隙度 双重渗透性吸湿和疏干滞后过程模型主要在溶质运移问题中使用van Genuchten-Mualem公式处理水分特征曲线 s 饱和含水量-r 残余含水量-a,m,n，l-经验参数 1/L,-,-,-Se 有效含水量-Ks 渗透系数 L/TK(h)压力水头为 h时的非饱和导水率 L/T水力参数水力参数编号编号 r
8、 s (cm 1)n Ks(cm/d)l土壤数据库Sandy loamsandr：土壤残余含水量s：土壤饱和含水量Alpha、N：土壤持水函数中的参数L-1Ks：饱和导水率，即渗透系数 LT-1l：导水函数中的曲度参数-神经网络预测模块神经网络预测模块99-101%0.5-2.00.0-1.0根吸水根吸水水分胁迫模型盐分胁迫模型根系吸盐模型补偿吸水域值（临界湿润度）其中，S(x)是吸水强度函数(cm/(cmd)，(h)是水分胁迫函数，b(x)是根系吸水分配（密度）函数，描述根系吸水的空间变异，Tp是潜在蒸腾量(cm/d)，Ta是实际蒸腾量(cm/d)，h为土壤压力水头(cm)，LR为根系层的深
9、度。根吸水模型根吸水模型根吸水参数根吸水参数 P3 P2L P2H Popt P0压力水头(pressure head,cm)=1=0r2Hr2L数据库吸水强度(cm/d)胁迫函数两个临界吸水强度v确定根系随深度的分布，b(x)满足：v其中，LR为根系层厚度,x是深度。本算例中根系层厚度1 m，假设根系分布函数为线性，并有v ,x1v容易证明上式满足积分为1的条件。选择Conditions/Root Distribution，使用Edit condition，划定根系层范围(100 cm)，顶部数值设为0，底部设为2，让程序自动进行线性差值。xxb2)(根系吸水分配函数根系吸水分配函数b(x)
10、根系分布根系分布根系分布函数使用工具条使用工具条ConditionsRoot Distribution初始条件初始条件v选择Conditions/Initial Condition。本算例中初始地下水位位于地下600 cm，假设土壤剖面初始状态是静力平衡态，则模型底部压力水头为-300 cm，地面压力水头为-600 cm。使用Edit Condition，设置顶部和底部压力水头，并让程序自动插值形成初始条件。加入观察点加入观察点v选择Conditions/Observation Points，使用Insert工具条添加若干观察点。土壤剖面数据汇总土壤剖面数据汇总运行模型，查看结果运行模型，查看
11、结果-20000-16000-12000-8000-40000051015202530Time daysN1N2N3N4N5Observation Nodes:Pressure Heads0.040.060.080.100.120.140.160.180.20051015202530Time daysN1N2N3N4N5Observation Nodes:Water Content-300-250-200-150-100-5000.000.050.100.150.20Theta-T0T1T2T3T4T5Profile Information:Water Content-300-250-200-150-100-500-20000-16000-12000-8000-40000h cmT0T1T2T3T4T5Profile Information:Pressure Head00.10.20.30.40.508001600240032004000|h|cmM1M2Hydraulic Properties:Theta vs.h
