不同初始含水率砂土水-汽-热耦合迁移数值模拟

郭 毅，晋 华，桂金鹏

(太原理工大学水利科学与工程学院，太原 030024)

干旱半干旱地区，昼夜温差较大，温度波动引起的液态水流动和水汽输送会对土壤水分的时空分布产生显著影响[1,2]。通常情况下，干旱区近地表土壤表层的含水量很低，水汽通量成为总水分通量的重要组成部分[3-5]。而且，水汽迁移同时携带了大量能量，进而影响了土壤内的温度分布。例如，Cahill等研究发现在裸土顶层2 cm处，40%到60%的热通量传递是由水蒸气迁移产生的[6]。土壤的水热状态会对发生在土壤内物理，生物和化学过程产生重要影响[7]。因此，有必要对干旱区土壤水-汽-热耦合迁移规律展开研究。

Penman[8]首次提出了基于Fick定律的土壤中水蒸气迁移理论。随后，Philip和de Vries[9]建立了温度梯度和含水率梯度双重作用下的水-汽-热耦合迁移数学模型，即PDV模型，并引入了水汽扩散增强因子解释温度梯度驱动的水蒸气通量。Milly[10]改进了PDV模型，在模型使用基质吸力代替了体积含水率，考虑了土壤的非均质性。Nassar等[11]对PDV模型进行了扩展，提出了土壤水分、热量和溶质耦合作用的传输方程。Bittelli等[12]进行了裸土条件下热量，液态水和水汽传输的实验研究和数值模拟。国内学者也在该领域做了大量的研究。例如，蔡树英[13]通过室内土柱蒸发试验，对PDV模型进行了验证，并与等温模型进行了对比。Liu Wei等[14]进行了自然环境下热量、水分和气体传递的数值模拟，分析了环境温度、相对湿度、太阳辐射和风速等参数对水热迁移的影响。王华军等[15]对砂土和壤土进行了考虑水蒸气迁移的数值模拟研究，发现两种土质在低饱和度时，水蒸汽迁移效应明显。

综上，国内外学者已经进行了很多水-汽-热迁移的相关研究，但不同初始含水率下水蒸气对土壤水热迁移影响的定量研究尚未见到，因此，有必要开展该方面的研究。本文以一端有热源影响的一维土柱为研究对象，利用Hydrus-1D软件将非饱和土壤水分运动方程和热量传递方程进行了耦合迭代计算，对不同初始含水率的砂土进行了有无水蒸气影响对比模拟，分析水蒸气迁移对整体水分迁移和温度传递的影响，以期为进一步研究土壤水-汽-热耦合传输规律提供参考。

1 理论模型

1.1 物理模型

本次物理模型以100 cm长土柱为研究对象，对土柱左端施加恒温热源，在距热源3 cm处设置第1个观测点，然后分别以5、10、10、15、15和20 cm间隔的布置其他6个观测点，研究土柱含水率和温度的变化规律，物理模型如图1所示。为了简化模型，做出如下假设：①土样为均质和各向同性的多孔介质；②不考虑土壤的形变；③只考虑水平方向的水汽热传递；④不考虑热辐射，只考虑热传导和热对流。

图1 一维土柱水-汽-热耦合迁移模型示意图Fig.1 One-dimensional model of coupled water-vapor-heat migration in soil

1.2 数值模型

1.2.1 水分运动方程

非饱和土壤水分传输有两种方式，一种是液态水的传输，一种是气态水即水蒸气的传输，驱动力有基质势梯度和温度梯度，水分迁移方程可用如下方程表述[16]:

(1)

式中：θ是总体积含水率，θ=θl+θv；θl是液态水体积含水率；θv是水蒸气体积含水率；KLh是基质势梯度驱动的液态水等温水力传导率；KLT是温度梯度驱动的液态水非等温水力传导率；Kvh是基质势梯度驱动的水蒸气等温水力传导率；KvT是温度梯度驱动的水蒸气非等温水力传导率；T是温度；t是时间；h是负压水头。x表示距热源的距离，原点位于土柱最左端，向右为正。

土壤水分特征曲线是描述负压水头与含水率的关系曲线，对分析和评价土壤水分的保持和运动有重要的意义。该特征曲线和液态水等温水力传导率KLh可由VG模型表述[17]。

(2)

m=1-1/n,n>1

(3)

(4)

KLT，Kvh和KvT由下列经验公式求得[19-21]：

(6)

(7)

(8)

式中：GWT是表征土壤水分特征曲线的温度依赖性的参数，取值7[22]，γ是土壤水的表面张力，γ=75.6-0.1425T-2.38×10-4T2；γ0是25 ℃时候土壤水表面张力；D是水蒸气扩散度；ρw是液态水密度；ρvs是饱和水汽密度；Hr是相对湿度；M是水分子的摩尔质量；g是重力加速度；R是气体常数；η是水蒸气扩散增强因子[23]。

水蒸气扩散度可以定义为:

D=τθaDa

(9)

式中：τ是土体孔隙曲率系数；Da是自由大气扩散系数；θa是被空气所填充的孔隙率。τ和Da可用下列给出的公式求得[24]。

(10)

(11)

ρvs、Hr、η分别由下列经验公式求得[24,10,23]：

ρvs=

(12)

(13)

(14)

式中：fc是土壤中黏土的质量分数。

1.2.2 热量传递方程

非饱和土壤中热量传输分为固体骨架热传导、液态水对流换热、水蒸气对流换热和液态水的蒸发潜热四部分，基于能量守恒定律，热量传输可用如下公式表述[16]：

(15)

土壤表观比热容和表观热导率可以下列经验公式求得[25,26]：

CP(θ)=Cnθn+Cwθw+Cvθv

(16)

λ(θ)=λ0(θ)+βtcw|q|

(17)

式中：θ是体积分数；下标n,w,v分别代表固液气三相；βt为热弥散系数；λ0(θ)为多孔介质热导率，可由公式λ0(θ)=b1+b2θ+b3θ0.5[27]求得。

1.2.3 初始条件与边界条件

初始条件：在初始状态下，土柱内各处温度、含水率相等，即：

T(x,t)=T0t=0

(18)

θ(x,t)=θ0t=0

(19)

边界条件：在x=0处

T(x,t)=TH

(20)

(21)

式中：TH为热源温度。

在x=L处

T(x,t)=Tc

(22)

(23)

式中：Tc为环境温度。

式(1)、式(15)与上述初始条件和边界条件构成非饱和土壤一维水平水-汽-热耦合迁移数学模型。

1.3 数值计算

非饱和土壤水-汽-热迁移方程组具有非线性和耦合性，无法求得解析值，故利用Hydrus-1D软件对水热传递方程进行离散化，进而求得数值解。模型求解区域为100 cm长土柱，空间计算步长为1 cm，模拟时间为1 440 min，初始计算时间步长为1 min，最小计算时间步长为0.01 min，最大计算时间步长为5 min。

1.4 数值模拟方案

基于非饱和土壤水-汽-热耦合迁移数值模型，设置以下6种模拟情景，模拟研究不同初始含水率下水蒸气对土壤整体水分迁移和温度变化的影响，模拟方案如表1所示。为了保证模拟结果适用的广泛性，所选土质为Hydrus软件内置的砂土，砂土VG模型参数如表2所示。

表1 模拟方案表Tab.1 Simulation scheme table

表2 砂土VG模型参数Tab.2 Parameters of VG model for sand

2 模拟结果与分析

考虑到砂土含水率和温度在距热源较近测点变化剧烈，在离热源较远的测点变化较小，因此本文选择在3和8 cm测点的温湿度动态变化进行研究。同时，由于含水率的绝对变化量相对较小，为方便对比水蒸汽在非饱和土壤水汽热迁移中的贡献，对含水率和温度变化均进行了无量纲处理，具体公式如下[16]：

无量纲土壤含水率：

(26)

无量纲土壤温度：

(27)

2.1 含水率变化模拟结果与分析

图2是初始含水率分别为0.08、0.11、0.13、0.15、0.20和0.30 cm3/cm3时砂土含水率随时间变化曲线。可以看到，无论是否考虑水蒸气迁移，在距热源3和8 cm处均会形成一个短暂的含水率峰值，该峰值现象的出现与陈红兵、王华军等[28,29]人的研究结果一致，这是因为在温度梯度和含水率梯度双重作用下，水分产生了迁移，在该点处形成了短暂的水分滞留；在考虑水蒸气迁移后，含水率变化更加剧烈，表现为含水率峰值与模拟结束时刻土壤最终含水率之间的差值变大。下面以模拟结束时刻距热源3 cm处的最终含水率为例，分析不同初始含水率下水蒸气迁移对整体水分迁移的影响。当砂土初始含水率为0.08 cm3/cm3时，仅考虑液态水迁移，含水率减少0.039，考虑水蒸气迁移，含水率减少0.154，含水率的减少量增幅为294.87%；当砂土初始含水率为0.11 cm3/cm3，仅考虑液态水迁移，含水率减少0.091，考虑水蒸气迁移后，含水率减少0.130，含水率减少量增幅为42.86%；当砂土初始含水率为0.13 cm3/cm3时，仅考虑液态水迁移，含水率减少0.133，考虑水蒸气迁移，含水率减少0.148，含水率减少量增幅为11.28%；当砂土初始含水率为0.15 cm3/cm3，仅考虑液态水迁移，测点含水率减少0.159，考虑水蒸气迁移，含水率减少0.169，含水率减少量增幅为6.29%；当砂土初始含水率为0.20 cm3/cm3，仅考虑液态水迁移，含水率减少0.208，考虑水蒸气迁移后，含水率减少0.209，含水率减少量增幅为0.48%；当砂土初始含水率为0.30 cm3/cm3，仅考虑液态水迁移，含水率减少0.222，考虑水蒸气迁移，含水率减少0.222，含水率减少量增幅为0%。进一步地，图3是不同初始含水率下，考虑水蒸气迁移后3 cm处含水率减少量的增幅图。可以看到：初始含水率越低，水蒸气对土壤最终含水率的影响越大，初始含水率为0.08 cm3/cm3时，考虑水蒸气后，砂土含水率减少量增幅极大，这表明水蒸气迁移在整体水分迁移中占据了主导作用，随初始含水率的增加到0.11 cm3/cm3，砂土含水率减少量增幅迅速减小，当初始含水率增加到0.20 cm3/cm3时，水蒸气迁移对整体水分迁移的影响可以忽略。

图2 不同初始含水率条件下砂土含水率随时间变化曲线Fig.2 The moisture variation curve of sand with time under different initial moisture content

图3 砂土不同初始含水率下考虑水蒸气迁移后含水率减少量增幅图Fig.3 The increase of water content reduction of sand after considering water vapor migration under different initial water content

2.2 温度变化模拟结果与分析

图4是热源温度为50 ℃，初始含水率分别为0.08、0.11、0.13、0.15、0.20和0.30 cm3/cm3时温度随时间变化曲线。对比有无水蒸气对温度模拟结果的影响，这里以模拟120 min后距热源3 cm处的温度变化为例进行分析：当初始含水率为0.08 cm3/cm3时，与仅有液态水迁移相比，考虑水蒸气迁移后，模拟点处温度从0.800增加到0.828，增幅为3.5%；当初始含水率为0.11 cm3/cm3，与仅有液态水迁移相比，考虑水蒸气迁移后，模拟点处温度从0.801增加到0.823，增幅为2.75%；初始含水率为0.13 cm3/cm3，与仅有液态水迁移相比，考虑水蒸气迁移后，模拟点处温度从0.799增加到0.819，增幅为2.50%；初始含水率为0.15 cm3/cm3，与仅有液态水迁移相比，考虑水蒸气迁移后，温度从0.798上升到0.816，增幅为2.26%；含水率为0.20 cm3/cm3，与仅有液态水迁移相比，考虑水蒸气迁移后，温度从0.796上升到0.808，增幅为1.51%；含水率为0.30 cm3/cm3，与仅有液态水迁移相比，考虑水蒸气迁移后，模拟点处温度从0.792上升到0.796，增幅为0.51%。进一步地，图5是不同初始含水率下，考虑水蒸气迁移后3 cm处砂土在120 min时刻温度增幅图。可以看到：水蒸气的迁移加快了砂土温度上升的速度，且随着初始含水率增加，水蒸气对温度提升幅度不断减小。砂土中水分迁移和热量传递是个耦合的过程，砂土温度变化的趋势与含水率变化是密切相关，当初始含水率较少时，土体孔隙的气相成分所占比例较大，水蒸气的迁移较其在高含水率中的迁移阻力较小，水蒸气迁移对土壤含水率影响很大，故砂土温度上升过程差异较大。随着含水率增加，水蒸气迁移对整体水分迁移贡献越来越小，考虑和不考虑水蒸水的模拟差别越来越小。当含水率增加到0.20 cm3/cm3时，模拟点处砂土温度的变化几乎一致，可以忽略水蒸气迁移对温度变化的影响。

图4 不同初始含水率条件下砂土温度随时间变化曲线Fig.4 The temperature variation curve of sand with time under different initial moisture content

图5 砂土不同初始含水率考虑水蒸气迁移后温度增幅图Fig.5 The temperature increase diagram of sand with different initial moisture content after considering water vapor migration

3 结 论

本文基于Hydrus-1D软件构建了热源作用下一维非饱和土壤水-汽-热耦合迁移模型，进行了多种初始含水率下有无水蒸气的对比模拟，可以得到以下结论。

(1)与仅考虑液态水迁移相比，考虑水蒸气迁移后会导致含水率迁移加剧，当初始含水率为0.08 cm3/cm3时，水蒸气迁移在整体水分迁移中占主导作用，随着初始含水率增加0.11 cm3/cm3，水蒸气迁移对整体水分迁移影响大幅减小，当初始含水率继续增加到0.20 cm3/cm3，水蒸气对整体水分迁移的效果可忽略。

(2)与仅考虑液态水迁移相比，考虑水蒸气迁移后会导致土壤温度升高速度提升，且初始含水率越低，水蒸气对温度的提升幅度越大，随着初始含水率不断增加，水蒸气对土壤温度的影响不断减小直至消失。

(3)在对低含水率砂土进行水热迁移模拟时，有必要将水蒸气因素考虑进去提高模拟精度，而对于含水率较高的砂土，可仅考虑液态水迁移以简化计算。