毛乌素沙地冻融期气态水迁移机理及影响因素

2022-05-10高万德陈云飞卢玉东刘秀花

水科学进展 2022年2期

郑策，高万德，陈云飞，卢玉东，刘秀花

(1. 长安大学旱区地下水文与生态效应教育部重点实验室，陕西西安 710064；2. 长安大学水利与环境学院，陕西西安 710064)

毛乌素沙地是中国四大沙地之一，面积达4.22万km2，区内水资源短缺，属于黄河流域生态脆弱地带[1- 2]。同时，由于地处季节性冻土带，受冻融循环过程影响，冬季时土壤水、热传输过程变得复杂，包气带液态水与气态水运移对区内许多工程建设以及农业生产活动造成较大影响，但目前研究中对于气态水运移过程的理解较为薄弱。厘清气态水运移机理及其影响因素，有助于更好地理解冻融期沙地水文过程机理[3]。

长期以来，学者们对于包气带水分运移过程开展大量研究，关于液态水的运移机理较为明确。随着研究深入，以Philip和de Vries所提出的理论为基础，许多试验与模拟研究成果均表明气态水同样是包气带水分迁移过程的重要组成部分[4- 6]。然而，相关研究主要围绕土壤未冻结的情况取得了丰硕的成果，针对冻融过程中的研究较少[7- 9]。一方面是由于传统的非饱和水分运移研究中往往忽略了负温条件下气态水的迁移，只考虑液态水运移过程。对于湿润地区或许可以忽略气态水的影响，但对于地下水位埋深较大的干旱半干旱地区而言是不合适的[10- 12]。另一方面，由于野外研究的复杂性，部分水汽热耦合研究采用室内土柱试验与数值模拟相结合的方法来开展，因而在研究中存在一些简化现象，如试验与模型中土柱上部设置为定水头或者定通量边界、土柱剖面设置为均质等。通过该研究虽然可以了解水分运移过程，但实际上与野外情况存在明显区别，不能很好地阐明水分运移机理，并且所建立的数值模型也不一定适用于野外复杂条件，容易出现数值计算问题[13- 15]。

对于诸如毛乌素沙地之类的干旱半干旱地区，在土壤未冻结时期，水汽通量已被证明是总水分通量的重要组成部分，对水分运移影响不应被忽略[16- 17]。通过厘清冻融期内温度梯度驱动下气态水运移过程，可为揭示旱区沙地包气带水文循环机理提供重要理论依据，并进一步为毛乌素沙地生态系统脆弱区的植被恢复以及生态系统管理提供科学依据。为此，本文以毛乌素沙地为研究对象，结合野外原位监测与数值模拟手段，对冻融期包气带水汽热耦合运移过程展开研究，分析气态水的分布特征、运移规律及其在总水分通量中占比，探讨水汽运移对土壤水分分布的影响。

1 研究区概况与试验研究方法

1.1 研究区概况

本研究的原位监测试验场选择地下水位埋深为9 m的毛乌素沙地榆林风沙滩区，位于国家林业局榆林荒漠生态站中(38°23′19″ N，109°42′29″ E)，如图1所示。由多年气象资料可知，研究区内主要气候特征为冷热剧变及日照强烈，年平均气温为6.4 ℃，且年温差与日温差变幅较大；年平均降雨量与蒸发量分别为360.2 mm及2 343.4 mm，其中，降雨量在年内分布不均，超过60%的降雨集中在7—9月，年际变化较大[18]。试验场内土地平坦，只在夏季雨水较多时有零星沙柳生长。通常在每年11月至次年3月间，土壤会经历冻融循环过程，最大冻结深度会超过1 m。

图1 试验点位置与场地布置Fig.1 Location of the experimental site and field instrumentation layout

表1 研究区不同层位土壤物理性质

本研究的野外工作开始于2017年4月，通过选取典型剖面，从地表至地下200 cm处每20 cm用环刀取3个重复土样，在实验室内利用Bettersize 2000激光粒度分析仪(丹东百特仪器有限公司)测定土壤颗粒组成信息，通过烘干法测定土壤重度，测定结果见表1。随后，分别于地表以下10 cm、20 cm、50 cm、100 cm、130 cm、200 cm处安装Hydra Probe Ⅱ(Stevens公司)土壤含水量与温度监测探头，地表以下20 cm、100 cm、200 cm处安装HFP01(Hukseflux公司)热通量监测探头以及地表以下900 cm处安装SDP(Stevens公司)地下水位监测探头。所有传感器通过包气带水热监测井壁孔洞插入土壤相应深度中，并连接至CR1000数据传感器(Campbell公司)，监测频率为1 min，每10 min记录1次数据。需要注意的是，Hydra Probe II传感器测量含水量的原理是基于频域分解法，通过分解土壤中复介电常数的实部与虚部来确定土壤体积含水量，其中液态水的介电常数为78，远高于土壤基质、冰、空气等其他组分，因此在土壤冻结后该传感器所测量数值为液态水含水量，即未冻水含量[19]。同时，在场地内安装有Decagon公司的气象要素监测设备，可监测气温与相对湿度(VP- 3传感器)、降雨(ECRN- 100传感器)、风速(Davis Cup传感器)等信息。

1.2 水汽热耦合运移方程

(1) 水分运移方程。考虑温度、气态水以及冰的共同影响时，对于垂向一维水流运动，其控制方程可用修改后的Richards方程来表示[20]：

(1)

式中：θl、θv和θi分别为土壤中液态水含水量、气态水含水量以及含冰量，cm3/cm3；ρl、ρv和ρi分别为液态水、气态水及冰密度，g/cm3；t为时间，h；z为垂向坐标轴(向上为正)，cm；h为基质势，cm；T为温度，K；qlh、qlT、qvh和qvT分别表示等温液态水、非等温液态水、等温气态水以及非等温气态水通量，cm/h；Klh和Kvh分别表示等温液态水与气态水水力传导度，cm/h；KlT和KvT分别表示非等温液态水与气态水水力传导度，cm2/(K·h)。可以看出，气态水运移受基质势梯度与温度梯度驱动，数值分别与相应的水力传导度及驱动力梯度有关，如下所示：

(2)

(2) 热传导方程。对于一维垂向热传导过程，考虑气态水与冰的影响后，控制方程可表示为如下形式[20]：

(3)

式中：Cp为土壤的体积热容量，J/(cm3·K)，可根据土壤中固、液、气各相乘以各自体积分数求和得到，其中Cw和Cv分别为土壤中液态水及气态水的体积热容量；L0与Lf分别为蒸发及冻结潜热，J/kg；λ为热传导系数，J/(cm·h3·K)；ql和qv分别为液态水及气态水通量，cm/h；S为源汇项，h-1。

(3) 土壤水力参数。由基质势梯度及温度梯度引起的液态水与气态水水力传导度可分别由下式计算得到：

(4)

式中：Ks为饱和导水率，cm/h；Se为有效饱和度；l和m为经验参数；Gwt为描述土壤水分特征曲线受温度影响的参数；γ和γ0分别为土壤水的表面张力和25 ℃时的表面张力；D为土壤中水汽扩散系数，cm2/h；ρsv为饱和水汽密度，g/cm3，由温度数据计算得到；M为水的摩尔质量；g为重力加速度；R为气体常数；η为水汽扩散增强因子；Hr为相对湿度，计算公式如下：

(5)

可以看出，相对湿度主要与温度及基质势有关，而与饱和水汽密度乘积即为水汽密度(ρv=Hrρsv)。公式(4)中其余所涉及参数具体计算过程可参考Hansson等[20]。需要注意的是，在不考虑冻融过程时，水、热传输方程中存在有3个未知量，分别为含水量、基质势以及温度，其中含水量与基质势之间的关系可以用土壤水分特征曲线来表示。当考虑冻融作用影响时，方程中引入了新的未知量，即含冰量。此时，可以通过Clapeyron方程建立冻土中基质势与温度之间的联系：

(6)

式中：Tf为液态水冻结温度，K。将上式代入土水特征曲线方程中(如van Genuchten方程)，可以得到θl与T之间的关系，即土壤水分冻结曲线。

1.3 数值模型建立及验证

1.3.1 模型建立

由于操作方便且功能完善，Hydrus- 1D软件被广泛应用于包气带水汽热耦合运移研究中，但在标准版的Hydrus- 1D软件中并未考虑冰- 水相变过程，只能被用于土壤未冻结情况。虽然Hansson等[20]提出了适用于冻融过程Hydrus- 1D程序，但该程序只适用于较为简单的情景，当所建立的模型复杂时容易遇到数值计算问题，并且在该程序中并未涉及到水汽运移及其相变过程，因而总体上应用较少。在该程序基础上，以冻融过程包气带水汽热耦合运移理论为指导，利用Fortran语言编写相应程序，修改后的程序充分考虑了冻融条件下土壤液态水、气态水运移以及热传导过程。同时，对程序中可能引起数值计算不收敛问题进行改进，主要包括基于有效能量理论修正了相变时由表观热容增大引起的温度计算误差，以及将含冰量计算模块加入每个时间节点计算出新的土壤温度[21]。由于水汽热耦合过程高度非线性以及边界条件复杂性，因此采用数值法求解，其中，水分运移方程通过有限差分法进行时空离散，热传导方程分别采用有限差分法和有限元法对时间及空间剖面进行离散，求解时依次计算水分运移与热传导部分。以水分运移方程为例，求解时步骤如下：首先根据前一个时间步长土壤水、热计算数据更新传导度、容水度等水力参数；随后离散方程并构建矩阵方程组求解最新时间步长含水量与基质势值，采用迭代法反复求解直到精度达到预先设定要求；最后计算液态水与气态水通量并将计算结果传入热传导模块。

利用该程序建立包气带水- 汽- 冰- 热耦合模型，选取2017年11月至2018年3月为研究期，收集原位监测土壤数据与气象数据。模拟时以小时为单位，模拟期共计3 624 h，模拟剖面为地表至地下200 cm。模型初始水、热条件采用实测的含水量与温度数据，并通过线性插值法获得每个节点数值。水分运移上边界受气象因素控制，由降雨与蒸发条件决定，下边界由于距离地下水位(约为9 m)较远，因此忽略地下水的影响。热传导上、下边界则分别采用计算与实测得到的土壤温度值作为边界条件。采用van Genuchten- Mualem模型计算非饱和液态水水力传导度，所需参数通过野外入渗试验以及ROSETTA模型联合确定，最终采用的土壤水力参数见表2。

表2 研究区不同层位水力参数

1.3.2 模型验证

研究期内实测与模拟的土壤未冻水含量与温度对比变化情况如图2所示，图中也列出了定量化评价指标决定系数(R2)与均方根误差(ERMS)的计算结果。

由图2可以看出，10 cm、50 cm及100 cm处土壤水分在冻结期内明显下降，并且变化过程存在明显滞后现象，表明土壤冻结现象发生，所建立的模型很好地反映出了这一过程。由于冬季降雨稀少(降雨与降雪总和为20.8 mm，且主要以降雨形式发生在土壤融化之后)且蒸发较弱，浅层土壤水分变化受外界因素影响较小，因此模拟误差较小，如10 cm处R2与ERMS值分别为0.96及0.004 cm3/cm3。相对来说，100 cm误差值较大，ERMS值为0.009 cm3/cm3，造成这一现象主要是由于80 cm处上侧与下侧土壤性质发生变化，引起了水力参数的改变，造成模型误差增大，同样的现象也出现在Yu等[22]研究中。通过在所研发的程序中嵌入参数反演模块，优化所使用的参数，可在后续研究中进一步提高模型精度。另一方面，在研究期内，土壤温度变化与气温变化保持一致，呈现出先降后升的趋势，温度变幅随土壤深度增大而减小。所有深度温度模拟结果均较好，平均R2与ERMS值分别为0.97及0.4 ℃。总体来说，模型运行期内并未出现数值计算问题，且土壤水、热变化均反映出较好的拟合精度，表明该模型适用于评价研究区冻融期土壤水分运移过程。

图2 实测与模拟含水量及温度对比Fig.2 Comparisons of simulated and measured unfrozen water content and soil temperature

1.3.3 典型模拟时段选取

根据模拟结果中0 ℃所处深度，绘制了冻融期内土壤冻结深度曲线(见图3)，可将冻融过程划分为4个时期，分别为① 未冻结期(11月1—18日及3月13—31日)：土壤未经历冻融过程；② 初始冻结期(11月19日至12月4日)：浅层土壤出现明显夜间冻结白天融化现象，该时期内冻结深度最大不超过30 cm；③ 向下冻结期(12月5日至翌年2月17日)：浅层土壤完全冻结不再经历冻融循环过程，冻结深度不断增大；④ 融化期(2月18日至3月12日)：冻结锋运移到最深处后，随后由于气温逐渐回升，冻结深度逐渐减小。

在上述4个时期中，分别选取典型时段，具体分析小时尺度上土壤水分、水汽密度以及水汽通量的分布特征及变化规律，典型时段选取具体信息见表3。同时，表3中列出了所选取时段小时尺度土壤未冻水含量及温度的模拟值与实测值误差计算结果，该结果为时段内4个典型深度(10 cm、50 cm、100 cm及200 cm)计算平均值。可以看出，未冻水含量与温度拟合结果均较好，再次表明所建立的模型适用于毛乌素沙地包气带水汽热耦合运移研究。

2 结果分析与讨论

2.1 不同阶段土壤液态和固态水分布特征

在所选取的4个时段内(见表3)，土壤未冻水含量(即液态水含水量)时空分布特征如图4所示，冻融过程中含水量剖面分布存在有如下2个显著特点。

图4 不同时段未冻水含量时空分布Fig.4 Spatial- temporal distribution of unfrozen water content in different periods

首先，P2时段浅层、P3时段上部以及P4时段中部冻层内未冻水含量均明显低于其余未冻结部分，且温度越低时未冻水含量越小，最小值仅为0.01 cm3/cm3。其中P2时段浅层10 cm范围内土壤受气温日变化影响，存在白天融化夜间冻结现象。由P4时段等值线变化情况可以看出融化锋运移规律，向下与向上的融化锋分别从42 cm移动至54 cm以及从98 cm移动至93 cm。其次，剖面80～160 cm土层内含水量要高于其上下层位，这除了受土壤颗粒组成影响外，还与气态水迁移过程有关(具体分析见下文)。此外，对比P1与P4时段水分剖面可以看出，冻融过程会改变剖面含水量分布情况，尤其是造成深层含水量明显降低，如在这2个时段200 cm处平均未冻水含量分别为0.091 cm3/cm3及0.077 cm3/cm3。相应地，在土壤冻结后随着未冻水含量减小，含冰量逐渐增大，冻层内不同深度含冰量变化情况如图5所示(其中并未考虑P1时段)。对于浅层2 cm与10 cm而言，受气温日变化影响，在P2时段含冰量日变幅较大，在P3时段由于土壤完全冻结导致其变化较小，含冰量呈现缓慢增大趋势，最大值达到了0.062 cm3/cm3；对于50 cm而言，由于在P4时段融化锋穿过该深度，因此含冰量显著减小，由0.047 cm3/cm3减小至0。

图5 不同深度含冰量随时间变化情况Fig.5 Variations of ice content with time at different depths

表4 不同深度平均相对湿度 %

2.2 不同阶段水汽密度时空分布特征

由于现有手段无法直接监测土壤水汽通量，其数值通常可由水汽密度计算得到，而水汽密度的计算与土壤相对湿度有关，表4中列出了不同深度相对湿度变化情况。可以看出，在土壤未冻结时，孔隙中相对湿度几乎处于饱和状态；当土壤冻结后，受温度与基质势减小影响，相对湿度也会减小。根据土壤相对湿度计算结果，进一步绘制出P1—P4时段内水汽密度时空分布特征，如图6所示。

图6 不同时段水汽密度时空分布Fig.6 Spatial- temporal distribution of vapor density in different periods

由图6可以看出，水汽密度时空分布特征与温度变化联系紧密。对于浅层40 cm范围，受温度日变化影响，水汽密度存在明显的日变化规律，随深度增大变幅减弱。如在P1时段，表层2 cm处水汽密度最小值与最大值分别为5.4×10-6g/cm3与1.1×10-5g/cm3，而20 cm处变幅仅为6.9×10-6～9.1×10-6g/cm3。此外，在P2与P3时段，受土壤冻结作用影响，表层水汽密度在白天时增幅要弱于P1与P4时段。对于40 cm以下范围，水汽密度受气温变化影响较小，并不存在日变化现象，其中，在P1—P3时段数值呈缓慢减小趋势，在P4时段略有增大。相对来说，下部土层的水汽密度要高于上部，其中最大值出现在P1时段200 cm处，为1.2×10-5g/cm3，这与该点处温度(14.1 ℃)为剖面中最高有关。此外，受气温季节性变化影响，土壤水汽密度同样呈现出季节性变化规律。如50 cm处水汽密度呈现出先减小后增大的趋势(4个时段内分别为8.7×10-6g/cm3、5.6×10-6g/cm3、4.6×10-6g/cm3及4.9×10-6g/cm3)，而在更深层则由于P4时段融化过程未完全结束仍处于减小状态。

由公式(5)可知，温度与基质势是影响相对湿度与水汽密度的最主要因素，进一步计算温度与基质势变化对水汽密度影响。结果表明，当基质势不变时，温度升高1 ℃，水汽密度增幅约为6%；当温度不变时，基质势变化10倍与100倍时，水汽密度变化量分别仅为0.1%及1.3%。整体上来看，不论对于浅层还是深层土壤，水汽密度的变化与温度联系紧密，而受基质势影响较小。

2.3 不同阶段水汽通量变化规律

气态水在包气带中运移主要受基质势梯度与温度梯度影响，其中，由基质势梯度驱动的等温气态水通量仅在表层处较大，随深度增大数值显著减小，一般情况下小于由温度梯度驱动的非等温气态水通量2～5个数量级，因此本节通过选取剖面浅部(2 cm)、上部(10 cm)、中部(50 cm)以及下部(150 cm)典型深度，分析qvT随时间变化情况(图7)。

图7 不同深度非等温气态水通量日变化Fig.7 Diurnal variations of the thermal vapor flux at different depths

由图7可以看出，受温度变化影响，浅层2 cm和10 cm处通量呈现出明显的日变化规律。以2 cm变化情况为例，在每日8:00—18:00由于表层温度较高，因此通量方向向下，在其他时刻则在温度梯度驱动下向表层运移，不论是向上或向下的通量，其数值均随深度增大而减小。与qvh不同，冻结后qvT值则要小于未冻结情况，如10 cm处向下的通量最大值在P1时段达到-2.1×10-4cm/h，在P2与P3时段分别仅为-3.4×10-5cm/h及-5.3×10-6cm/h。对于50 cm以下深度，由于研究期内地表温度较低，随深度增大温度升高(与图6水汽密度剖面分布相似)，温度梯度方向向上，因此通量方向向上，水汽主要向上部运移。正是由于深层土壤气态水(以及液态水)在基质势梯度与温度梯度驱动下，不断向冻结锋处运移，引起了未冻水含量的持续降低。

作为气态水运移的驱动力，温度梯度(Tg，℃/cm，为剖面相邻两点温度差)不仅影响水汽通量大小，同时也控制着水汽运移方向。由图7可知，2 cm处水汽通量变化幅度较大，以该深度为例，Tg与qvT关系如图8所示。可以看出，不同时段计算得到两者相关系数均较高(0.88≤R≤0.98)，表明不论是向上或向下的通量均与温度梯度有很好的相关性。当Tg相等时，不同时段水汽通量值并不相同，这主要是由于除了梯度以外，传导度同样也会影响通量大小(见公式(2))。在土壤冻结后，水汽密度降低等因素导致气态水水力传导度减小，因此在同一温度梯度下，水汽通量值也小于土壤未冻结状态。

图8 不同时段2 cm处温度梯度与非等温气态水通量散点分布Fig.8 Scatter plots of the temperature gradient and thermal vapor flux in different periods at the 2 cm depth

2.4 水汽通量占比及驱动力分析

在研究期内，qvh在未冻结层位总水分通量中的占比趋近于零，当土壤冻结后其占比有所增大，在所选择的时段内最大值可达10%。相比较而言，qvT在总水分通量中占比较大且变化明显，因此本节重点分析qvT占比变化及其驱动力，结果如图9所示。

在土壤未冻结时(图9(a))，2 cm处qvT的占比为10%～20%，随深度增大占比逐渐减小，在50 cm以下范围占比仅为5%。由于水力传导度在剖面内变化不大，因此，qvT减小主要是受温度梯度(在表层变化剧烈，而随深度增大显著减小)变化影响。当土壤冻结后(图9(b)、图9(c))，qvT占比显著增大，如在P2时段的2 cm、10 cm以及P3时段的2 cm、10 cm、50 cm处，其数值均超过90%。在这些深度，占比的变化规律与含冰量的变化(见图5)保持一致，即含冰量增大时水汽通量占比也增大，造成这一现象主要是由于土壤冻结后冰的存在使液态水通量显著减小。此时虽然受温度梯度与非等温气态水水力传导度减小影响，qvT通量值小于未冻结时期，但占比却显著增大，表明气态水主导冻层内水分运移过程。同样的，在P4时段50 cm处，水汽占比不断减小是由于融化锋穿过该深度引起液态水通量增大造成的。通过以上分析看出，在所选择的典型时段内水汽通量在土壤冻结后占比较大。对于整个冻融期，虽然冻层内气态水在多数时刻均占据主导地位，但事实上仍然是液态水通量对总水分通量影响较大。以10 cm处为例，该深度在11月下旬至2月下旬期间土壤冻结，期间绝大多数时刻水汽通量占比均超过90%，但在整个研究期内(11月至3月)，气态水累积通量对总水分通量的贡献率只有约8%，这主要是因为该深度在融化时液态水通量显著增大，会超过气态水通量1～3个数量级。尽管如此，气态水在冻融过程中持续向浅层运移，对冻层内总含水量增加以及剖面水分分布均有较明显的影响，在研究中不可被忽略。

2.5 讨论

上述结果表明，在冻融时期水汽运移主要受温度梯度驱动，水汽通量在剖面内随深度增大而显著减小，其中浅层通量值较大主要与以下2个因素有关[5]：① 受蒸发作用影响，表层土壤含水量低，水汽扩散率大；② 地表与大气间水汽密度梯度较大。在土壤冻结后，虽然水汽通量数值要小于非冻结期，但由于冰的出现造成液态水通量的显著减小(2～5个数量级)，因此水汽通量在总水分通量中占比较高，影响反而更大。由于包气带中水汽密度多数时刻均处于接近饱和状态，因此只要温度发生变化，水汽便会相应蒸发或者凝结。对于干旱半干旱地区而言，包气带含水量较低且冬季降雨稀少，在季节性冻融过程中气态水运移对于剖面土壤水分分布影响较大。

图9 不同时段非等温气态水通量所占比例变化Fig.9 Variations in the ratio of isothermal vapor flux during different periods

随着气温降低，在11月初表层土壤温度成为剖面中较低点，包气带中水汽密度梯度方向向上(图6(a))，驱动气态水从深层向表层运移并最终逸散到空气中；随着气温不断降低，浅层土壤逐渐冻结，向上运移的气态水在冻层内凝结，造成冻层内总含水量不断增大，并最终冻结成冰，使季节冻土成为富冰冻土，而深层土壤中液态水含量则不断减小[23]；当土壤融化过程开始后，包气带内部冻层内水汽密度明显小于其上部与下部(图6(d))，此时气态水同时从两侧源源不断向该层内运移并凝结；随着气温不断回升，上部土壤温度梯度逐渐增大，由表层向下的水汽通量也随之增大，这也是造成80～100 cm范围内含水量要高于上、下层位的主要原因(图4)。万力等[24]在额济纳旗绿洲地区同样发现该现象，认为在季节冻土形成过程中气态水的迁移促进了高含水量带的形成。在冻融过程中水汽凝结与聚集现象对于维系荒漠植被生态系统有重要意义，可为诸如沙柳、樟子松等多年生典型沙区植被度过漫长的冬季干旱期提供重要水资源。

相比于Hansson等[20]提出的Hydrus- 1D冻融程序，本研究中所使用的修改后的程序不仅优化了以往数值计算中存在的问题，同时也很好地阐明了气态水在冻融过程中的运移机制。通过在程序中加入参数反演模块优化所使用的水热参数，有望进一步提升模型拟合精度[25]。在厘清裸地条件下气态水运移机理及其影响后，以数值模拟手段为基础，并结合水、汽稳定同位素技术，明确不同深度气态水的来源及其对区内典型植被的支撑作用，可为建立合理的种植模式及生态系统脆弱区的植被恢复提供依据。