蒸发条件下非饱和土水分迁移的数值模拟

熊 义 高 抗

（1.三峡大学水利与环境学院，湖北宜昌 443002；2.三峡大学三峡库区地质灾害教育部重点实验室，湖北宜昌 443002）

大气和土表的水量交换是水循环当中的重要组成部分.渗透作用使大气中的水分进入土体，蒸发蒸腾作用又使水分离开土体进入大气，由此而引起的土体含水率的变化直接影响着土体的稳定性.例如，膨胀土地基上土体遇水膨胀、失水收缩，可能引起上部建筑物的隆起、倾斜或开裂［1］.黄土的湿陷特性导致其在含水率增高时强度降低，在有荷载作用的情况下出现地基不均匀沉降的危害［2］.因此，研究自然条件下非饱和土的水分迁移规律对把握土体变化及工程应用有十分重要的意义.

长期以来，人们通常采用经典的Richards等温流方程描述土壤中液态水的运动［3］.例如，Thomas等［4］采用该方程的一维形式模拟了在气候变化下土体脱湿变干行为.吴宏伟等［5］采用该方程的二维形式分析了雨水入渗对非饱和土坡稳定性影响.然而，土壤蒸发是土壤中的水由液态水转化为气态水的过程，是一个相变的过程.土壤水分在液相和气相之间相互转换时，必然伴随着不同形式的水分（液态和气态）及热量的流动.因此，在土壤－大气连续体中分析非饱和土中水分的运动时，液态水的流动，气态水的流动和热量的流动需要同时考虑.本文以Penman－Wilson蒸发模型预测的实际蒸发量为水流量边界条件，利用较易获取的土体基本物性指标（颗粒分布、干密度、比重）预测出相应水力参数和热特性参数，同时结合基于Darcy定律和Fick定律而建立的非饱和土水－热迁移模型，分析了蒸发条件下非饱和土壤水分的运动规律，为自然环境下非饱和土水分迁移问题的探究提供参考.

1 数值模型的建立

1.1 土壤水分运动控制方程

为了建立描述土体中液态水、气态水和热量的运动方程，需引入如下假定［6］：土体的力学行为可以用其自身土颗粒、空气和水形成的连续单元体所描述，忽略由于渗透压梯度引起的液态水流动，液态水和空气的渗透系数仅仅和含水量或饱和度相关，不考虑渗透系数和吸力关系间的滞回，单元体的液态水和水蒸气局部热力学平衡始终成立，土体中任一点的温度在水的融点之上和沸点之下，水中气体的溶解不加以考虑.基于以上假设，建立了单元体中水的质量和能量守恒的方程，其一维形式表述如下：

方程（1）中有总水头和水蒸气分压两个变量，但这两个变量并不是独立的，他们之间的关系可以通过Edlefsen和Anderson［7］提出的以下关系式来描述：

式中，Pvs为纯水饱和蒸气压力（kPa）；hrair为空气相对湿度；w为水蒸汽克分子量（0.018kg／mol）；R为通用气体常数（8.314J／mol·K）；T为绝对温度（K）.

1.2 边界条件

在利用上述方程进行饱和－非饱和渗流计算时，必须给定计算模型的流量边界条件和温度边界条件.

根据数值模拟的特点，可选用土表实际蒸发量作为计算模型的流量边界条件，其值可用Penman－Wilson公式［8］计算，该公式是Wilson通过考虑土表相对湿度变化，同时引入总吸力，把土壤蒸发作为土壤湿度的连续函数而建立的：

式中，E为蒸发量（mm／day）；Γ为饱和蒸气压和温度关系曲线在温度为Ta时的斜率；Q为土表面等效净辐射量（mm／d）；v为湿度常数，一般取0.66hPa／℃，Ea＝f（u）ea（B－A）；f（u）是风函数，f（u）＝0.35×（1＋0.15u），u为风速（m／s）；ea为蒸发面上方空气中的水蒸气分压（kPa）；B为空气相对湿度的倒数，即es／ea；es为空气饱和蒸气压（kPa）；A为土表面相对湿度的倒数，即1／hr；当土表饱和时，即土表相对湿度为100%时，该公式演变为

该式即为Penman推出的计算饱和蒸发量的公式，同时可以用于计算水面的潜在蒸发量.

土体和大气间的热量交换发生在土壤表面，因此地表土温度可以作为计算模型温度的边界条件，其值可用Wilson提出的下式［6］换算：

式中，Ts为地表土的温度（K）；Ta为土表空气温度（K）.

1.3 参数确定

为了求解上述模型，必需给定材料参数及蒸发边界条件，如土水特征曲线、渗透性函数、水蒸气扩散系数、热特性参数和土表实际蒸发量.为避免传统试验测试手段周期长、费用昂贵的不足，本文根据相关气象条件（温度、相对湿度、风速）和Penman－Wilson公式预测土表实际蒸发量，利用土壤基本物性指标（颗粒分布、干密度、比重）预测上述参数，具体方法简述如下：利用颗粒分布曲线和相关物理性质指标，根据Arya－Paris［9］模型预测出不同含水量所对应的基质吸力值.在此基础上，利用Fredlund ＆Xing［10］模型对含水量－基质吸力的离散点进行拟合便得到完整土水特征曲线.

土体在非饱和情况下，其渗透系数随着其自身含水量的变化而变化.在已知土体颗粒分布、比重及干密度的基础上，利用Kozeny－Carman［11］模型预测出土样饱和渗透系数，同时结合Fredlund ＆Xing模型预测所得的土水特征曲线相关参数，便可得到完整的渗透性函数曲线.

根据Philip和deVries［12］模型的描述，在土体孔隙率（可由比重和干密度求得）已知的前提下，水蒸气扩散系数可由其与土体含水量和环境温度的函数关系求得.

土体的热特性参数主要指土壤的体积比热和热传导率，其值取决于土体中固相物质（土颗粒）、液相物质（孔隙水）和气相物质（孔隙气）所占的体积比，以及他们各自的热容和热传导率.水蒸气的潜化热是指单位质量的液体转变为相同温度的水蒸气时吸收的热量，其值可以利用一个和温度呈线性的方程描述［6］：Lv＝4.186×103×（607－0.7T）.

2 模型验证

2.1 几何模型求解及材料参数预测

为了验证上述方程及参数预测方法的可靠性，将该理论用于模拟Wilson土柱蒸发试验实测情况.采用有限单元法来求解偏微分方程中含水量及温度变化率，利用Vadose／w建模并实现数值模拟.

试验用土为均质砂土，其颗粒分布曲线如图1所示，比重为2.67，土柱模型干容重为1.64g／cm3，高度为30cm，直径为10cm.总共划分300个四边形有限元单元，沿高度划分为30个等份，沿直径划分10等份.

图1 颗粒分布曲线

由土壤基本物性指标（颗粒分布、干密度、比重），结合上述相关参数的预测方法可获取模型计算所需参数.其完整的土水特征曲线如图2所示，渗透性函数相关参数见表1，土样热传导率和体积比热容与含水率关系曲线如图3所示.

图2 土水特征曲线

图3 热传导率函数曲线和体积比热容函数曲线

表1 Fredlund ＆Xing模型预测土水特征曲线和渗透性函数相关参数

影响蒸发的主要外部气象因素主要有风速、温度、湿度.因此，为了求解瞬时水流量边界和温度边界，需要给定实时地气象参数.Wilson土柱蒸发试验条件为：平均环境温度为38℃，相对湿度接近10%，风速为0m／s.在此环境条件下土柱模型从饱和含水量开始持续蒸发40d，这一过程当中，土柱没有任何水分的补充，土柱水分的蒸发量即为土柱质量的减少量，其值采用电子天平即时获取.

2.2 计算结果分析

土柱蒸发试验蒸发量实测值与计算值对比结果如图4所示.从图4中可以看出，在前4d，土柱蒸发量实测值与计算值均稳定在7mm／d的潜在蒸发量附近.但从第4天末开始，土柱蒸发量迅速降低，在不到两周的时间里，其值从7mm／d降低到1mm／d.随后，土柱的蒸发量缓慢降低并趋于平稳，在超过两周的时间里，其值从1mm／d降低到0.5mm／d.通过对比不难发现，在为期40d的蒸发试验中，蒸发量实测值和计算值的变化趋势基本一致.同时，试验及数值模拟结果也验证了非饱和土蒸发呈现3个阶段：稳定阶段、减速阶段和残余阶段.

图4 蒸发量实测值与计算值比较

图5为蒸发试验的不同时刻土柱纵剖面体积含水量分布情况.从图中可以看出，在蒸发试验进行到第4天末，土柱不同深度处体积含水量均低于初始饱和含水量，但此时不同深度的含水量降低幅度较一致，因此在纵剖面上并没有明显的体积含水量梯度出现.在这一时刻，剖面含水量计算值和实测值吻合较好.但随着蒸发试验的持续进行，剖面含水量继续降低，在第10天末，土柱表层体积含水量接近于0，而表层以下10cm处的体积含水量为19%，这标志着在土柱表面有含水量极低的干燥土层的形成.此后，干燥面不断的向下推移，在第40天末，干燥面推移至土表以下10cm处，此时的干燥土层厚度达到10cm.通过对比可以发现，在蒸发进行到第10天末及第40天末，土柱纵剖面含水量的计算值和实测值出现了明显的差异，特别是在干燥土层以下相对湿润的置，出现这一现象的原因有可能是随着含水量的降低，预测的土体渗透系数减小得过快，因而在某一含水量时，输水能力的预测值小于其真实值，在模拟计算时，同一单元节点处液态水流量小于实际情况，随着深度的增加，这一现象更为突出，这也从侧面说明了土体渗透系数函数对干燥土层以下土壤水分的运动起决定性的作用.而在干燥土层内部，土体渗透系数因含水量的急剧降低而趋近于0.

图5 不同深度处体积含水量实测值与计算值比较

图6～7分别为1cm深度处单元节点和8cm深度处节点水流量情况.从图中可以看出，在蒸发前期各点处水分的迁移以液态形式流动为主.但随着含水量的降低，单元点处的渗透系数不断减小，从而导致该点处液态水通量迅速降低.从某一时刻起，水蒸气扩散成为该点处水分迁移的主要形式.以1cm深度处节点为例，在前3天，液态水流量维持在7×10－8m／s附近，水蒸汽流量接近0（此时土体饱和度高，水蒸气扩散系数小）.从第4天起，液态水流量开始快速下降，而水蒸气流量随着水蒸气扩散系数的增大而增加.在第6天，液态水流量降低至5.5×10－9m／s，而此时水蒸气流量为2.0×10－8m／s，水蒸气的流动占据主导地位.

图6 1cm深度处节点液态水和水蒸气随时间变化情况

图7 8cm深度处节点液态水和水蒸气通量随时间变化情况

图8为不同时刻土体纵剖面温度分布情况.从图8可以看出在蒸发初期，土体温度迅速降低，在第1天末，距离土表1cm深度处土层温度为30.6℃，距离土表25cm深度处土层的温度为35.7℃，可见靠近土表层位置处温度降低的幅度大于底层的.但随着时间的推移，土体不同深度处的温度均开始回升，在第42天，土体纵剖面温度重新回到了初始值附近.出现上述现象的原因是液态水转变成水蒸气需要消耗热量.在试验初期，土表保持7mm／d潜在蒸发强度，土表蒸发需要消耗大量的热量，此时汽化潜热成为单元体热量变化的主要形式.而在第4天后，随着蒸发量的迅速降低，所需消耗的热量减小，热传导在单元体热量变化中的作用逐渐凸显出来.因此，在该作用下土体温度逐渐变得均匀并再次接近于环境温度.同时还不难发现，土柱纵剖面上最低点的温度和干土层的向下推移的位置较一致，这也从另一个角度说明了蒸发面的位置随着蒸发的进行而不断向土体内部移动.

图8 不同深度处温度实测值与计算值比较

3 结 语

利用土体基本物性指标预测了非饱和土水－热迁移模型计算所需的水力参数、热特性参数，同时引入气象环境参数预测了土表实际蒸发量和温度值，模拟了蒸发条件下非饱和土中水分及热量的迁移规律.

通过该理论预测得到的非稳定蒸发量、纵剖面含水量分布、纵剖面温度分布与Wilson室内蒸发试验实测吻合较好.这说明采用上述模型及相关参数预测方法来模拟蒸发条件下土壤水分的迁移是可行的.

模拟结果表明渗透系数的急剧降低导致干土层的出现，从而导致土壤水分的运动形式发生了根本的转变，即在干土层内部水蒸气扩散成为土壤水分向上运动的主要形式.而此时在水力梯度下的渗透作用几乎可以忽略，土壤的实际蒸发面也不再是土壤表面.

［1］ 柯尊敬，黄绍铿.膨胀土工程性质的研究总报告［R］.南宁：广西大学土木系，1995.

［2］ 王铁行，李 宁，谢定义.土体水热力耦合问题研究意义、现状及建议［J］.岩土力学，2005，26（3）：487－492.

［3］ 陈建斌，孔令伟，赵艳林，等.非饱和土的蒸发效应与影响因素分析［J］.岩土力学，2007，28（1）：36－40.

［4］ Thomas H R，Rees S W.The Numerical Simulation of Seasonal Soil Drying in an Unsaturated Clay Soil［J］.International Journal for Numerical＆Analytical Methods Geomechanics，1993，17（1）：119－132.

［5］ 吴宏伟，陈守义，庞宇威.雨水入渗对非饱和土坡稳定性影响的参数研究［J］.岩土力学，1999，20（1）：1－14.

［6］ Wilson G W.Soil Evaporation Fluxes for Geotechnical Engineering Problems［D］.Saskatoon：University of Saskatoon，1990.

［7］ Edlefsen N E，Anderson A B C.Thermodynamics of Soil Moisture［J］.Hilgardia，1943，15：31－39.

［8］ Wilson G W，Fredlund D G，Barbour S L.Coupled Soilatmosphere Modeling for Soil Evaporation［J］.Canadian Geotechnical Journal，1994，31：151－161.

［9］ Arya L M，Paris J FA.Physicoempirical Model to Predict the Soil Moisture Characteristic from Particle－Size Distribution and Bulk density Data［J］.Soil Science Society of America Journal，1981，45：1023－1030.

［10］Fredlund D G，Xing A.Equations for the Soil－Water Characteristic Curve［J］.Canadian Geotechnical Journal，1994，31：521－532.

［11］Carman P C.Permeability of Saturated Sands，Soils and Clays［J］.Journal of Agricultural Science，1939，29：26.

［12］De Vries D A.Thermal Properties of Soils［M］.Amsterdam：North－Holland Publishing Company，1963.