高寒地区冻土活动层温度场数值模拟

３．２０
１．３饱和土壤冻结曲线
在等效体积热容量ｃ。的表达式中，等式右边第４项中水 的相对饱和度Ｓ。与温度ｒ之间的关系将会对热传导过程产生 重要影响。水的相对饱和度Ｓ。与温度ｒ之间所确定的曲线一 般称为土壤冻结曲线，见图１。Ｂｏｎａｃｉｎａ等证明，在进行温度场 计算过程中，水分冻结函数的实际形状并不是那么重要，但是 它一定要满足潜热的条件¨“，即
Ｎｕｍｅｒｉｃａｌ Ｓｉｍｕｌａｔｉｏｎ ｏｆ ｔｈｅ Ｔｈｅｒｍａｌ Ｒｅｇｉｍｅ ｏｆ Ｐｅｒｍａｆｒｏｓｔ Ａｃｔｉｖｅ Ｌａｙｅｒ ｉｎ Ａｌｐｉｎｅ Ｒｅｇｉｏｎｓ
ＷＡＮＧ Ｊｕｎ—ｚｈｉｌ，ＬＩＡＮＧ Ｓｉ—ｈａｉｌ，Ｗａｎ Ｌｉｌ，ＷＡＮＧ Ｈａｉ．ｚｈｏｕ２
（１．Ｓｃｈｏｏｌ ｏｆ Ｗａｔｅｒ Ｒｅｓｏｕｒｃｅｓ ａｎｄ Ｅｎｖｉｒｏｎｍｅｎｔ，Ｃｈｉｎａ Ｕｎｉｖｅｒｓｉｔｙ ｏｆ Ｇｅｏｓｃｉｅｎｃｅｓ，Ｂｅｉｊｉｎｇ １０００８３，Ｃｈｉｎａ；
生态过程…。近年来，在全球气候变暖的背景下，高寒地区冻
土发生了退化现象旧ｊ。冻土活动层温度场研究对于揭示全球 气候变暖条件下，高寒地区的水文过程和生态变化有着十分重 要的意义。 ２０世纪７０年代以前，冻土温度的研究主要是在理论分析 基础上的近似解析法和以试验、野外勘探与调查资料为基础的 统计方法∞Ｊ。２０世纪７０年代，Ｈａｒｌａｎ首先建立了水热耦合的 数学模型Ｈ１，后来，许多学者基于Ｈａｒｌａｎ的思想建立了冻土水 热耦合数值模型。国内学者很早就开展了对冻土温度的研究。 庞强强等Ｍ１根据１９９１—２０００年的地面温度观测资料并结合数 字高程模型，应用斯蒂芬公式得到了青藏高原冻土区的活动层 厚度分布图。李述训等¨ｏ应用数值方法预测了气候变暖条件 下，青藏高原多年冻土的温度状况。郝振纯等…建立了季节冻 土水热耦合迁移的数学模型，并对黄河源区站点的冻融深度进 行了模拟。 近年来，等效显热熔法越来越多地被应用于相变问题的研 究。笔者根据这一方法的原理，通过引入描述冻土活动层饱和 程度参数，对该方法进行改进，使其适用于非饱和带，并以黄河
２黄河源区季节冻土温度模拟
２．１研究区概况
黄河源区位于青藏高原东部，面积约为２０ 在４
２００ ９３０
ｋｍ２，海拔
ｍ以上，相对高差大于ｌ
０００
ｍ。源区年平均气温为
△珥＝Ｊ（妒。Ｐ。ｃ。＋妒护。ｃ．＋妒；Ｐ。ｃ。）ｄＴ—ｆＰｉ￡ｄ妒．
（３）
一３．８９℃，平均降水量为３１２．２ ｍｍ，年潜在蒸发能力为
１ ３５３．２
Ｋｅｙ ｗｏｒｄｓ：ｐｅｒｍａｆｒｏｓｔ ａｃｔｉｖｅ ｌａｙｅｒ；ｐｈａｓｅ ｃｈａｎｇｅ；ｎｕｍｅｒｉｃａｌ ｓｉｍｕｌａｔｉｏｎ；ｅｑｕｉｖａｌｅｎｔ－ａｐｐａｒｅｎｔ ｂｅａｔ ｃａｐａｃｉｔｙ ｍｅｔｈｏｄ；ａｌｐｉｎｅ
高寒地区地下水的冻结与融化直接影响区域水文过程和
源区玛多县为例，利用ＣＯＭＳＯＬ Ｍｕｈｉｐｈｙｓｉｃｓ模拟了该地区冻 土活动层的温度状况。 １
表１模拟中使用的参数 参数 ｝Ｌ隙度妒 水的比热ｃ。／（Ｊ・ｋｇ～・Ｋ。） 冰的比热ｃｉ／（Ｊ・ｋｇ．。・Ｋ－１） 固体骨架的比热ｃ。／（Ｊ・ｋｇ。１・Ｋ－１） 水的热传导率Ａ。／（Ｗ・ｍ叫・Ｋ。１） 冰的热传导率Ａ／（Ｗ・ｍ～・Ｋ－１） 固体骨架的热传导率Ａ／（Ｗ・ｍ．１・Ｋ一） 水的密度Ｐ。／（ｋｇ・ｎｌ。） 冰的密度ｐｉ／（ｋｇ・ｍ。） 固体骨架的密度ｐ。／（ｋｇ・ｍ。３） 相变潜热∥（Ｊ・ｋｇ“） 冻结温度ｒ，／℃ 相变区间△∥℃ 未冻残余饱和度Ｊｓ。。 水的相对饱和度ｏ／％ 数值
ｒｒｔ
Ｉｎ的实测地温为边界条件，以１９８０年１月的实测剖面地温
为初始条件，在假定土壤含水量不变，忽略水分对流和热源（汇） 的情况下，使用ＣＯＭＳＯＬ Ｍｕｈｉｐｈｙｓｉｃ建立数值模型。
１５
Ｌ
式中：孔为融化温度；ｔ为冻结温度；（ｐｃ）。为显热熔。
∞ ｐ 堡 胁 ∞ 船 ｍ
５
罄‘
犀
誊 墨
＊
Ｏ
［一Ｌ
妪０
赠 。
一１０
—１５
图２
１９８０－－－２０００年玛多县多年月平均地温（单位：ｍ）
图１土壤水分冻结曲线
使用ＣＯＭＳＯＬ Ｍｕｈｉｐｈｙｓｉｃｓ地球科学模块中预定义的多孔 介质热传导应用模式，通过对该模式进行一定程度地修改来描 述上述伴有相变的非饱和一维热传导模型。在实际使用过程 中，首先选取ＣＯＭＳＯＬ Ｍｕｈｉｐｈｙｓｉｃｓ中具有二阶连续导数的 ｆｌｃ２１协函数来表示冻结函数¨“，建立含有相变潜热的焓变表达 式；然后通过定义流体材料属性和固体材料属性，给出等效热 传导系数和等效体积热容量；最后将模型剖分为４５６个单元进 行求解。
２．３参数分析
在冻土温度计算中涉及的有关参数很多，这些参数的测量 较为复杂，而且测量价格昂贵，在实际科研和工程计算中每次 都进行测量就会显得十分不经济。本次模拟将冻土活动层看
．２】・
万方数据
人民黄河２０１３年第３期 成是均质各向同性的地质体，相关冻土物理参数根据冻土活动 层的土壤性质，参照文献［１５］给出（表１）。
ｃａｎ
ｓｏｕｒｃｅ
ｒｅｇｉｏｎ．Ｔｈｅ ｒｅｓｕｌｔｓ ｓｈｏｗ ｔｈａｔ ｔｈｅ ｍｏｄｉｆｉｅｄ ｅｑｕｉｖａｌｅｎｔ—ａｐｐａｒｅｎｔ ｈｅａｔ ｃａｐａｃｉｔｙ ｍｅｔｈｏｄ ｉｓ
ｅｆｆｉｃｉｅｎｔ ａｎｄ
ｂｅ ｕｓｅｄ
ｔｏ
ｓｔｕｄｙ ｔｈｅ ｐｅｒｍａｆｒｏｓｔ ａｃｔｉｖｅ ｌａｙｅｒ ｖａｒｉａｔｉｏｎ ｍｅｃｈａｎｉｓｍｓ ｏｆ ｔｈｅ ａｌｐｉｎｅ ｒｅｇｉｏｎ． ｒｅｇｉｏｎ
万方数据
人民黄河２０１３年第３期 导过程中所起的作用，采用混合理论，按照所占体积的比例构 成一个等效热参数代入控制方程中进行计算。对于饱和土壤， 等效热传导系数Ａ由水、冰、固体骨架颗粒通过体积平均共同 表示：
Ａ＝妒Ｓ。Ａ。＋妒Ｓ。Ａ．＋妒。Ａ。
热容量可以表示为 Ａ＝Ｏｔ（妒Ｓ。Ａ。＋妒ＳｉＡ，）＋妒。Ａ。
２．Ｔｈｅ Ｆｉｆｔｈ Ｔｅａｍ．Ｈｅｎａｎ Ｎｏｎｆｅｒｒｏｕｓ Ｇｅｏ—Ｅｘｐｌｏｒａｔｉｏｎ Ｂｕｒｅａｕ，Ｚｈｅｎｇｚｈｏｕ
４５００１６，Ｃｈｉｎａ）
Ａｂｓｔｒａｃｔ：Ｔｈｅｎｎａｔ ｒｅｇｉｍｅ ｏｆ ｐｅｒｍＭｒｏｓｔ ａｃｔｉｖｅ ｌａｙｅｒ ｈａｓ
ａ
ｇｒｅａｔ ｓｉｇｎｉｆｉｃａｎｃｅ
第３５卷第３期 ２０１３年３月
人民黄河
ＹＥＬＬＯＷ
ＲＩＶＥＲ
Ｖ０１．３５．Ｎｏ．３ Ｍａｒ．．２０１３
【水文・泥沙】
高寒地区冻土活动层温度场数值模拟
王俊智１，梁四海１，万
力１，王海周２
（１．中国地质大学水资源与环境学院，北京１０００８３；２．河南省有色金属地质矿产局第五地质大队，河南郑州４５００１６）
寒的气候条件决定了该区为季节冻土，并形成岛状和大片连续 多年冻土并存的分布格局Ｈ“。近年来，在全球气候变暖的背 景下，黄河源区出现了地下径流减小、植被覆盖程度降低、冻土 退化等一系列生态地质问题¨…。
２．２模型概化
黄河源区玛多县气象站１９８０－－２０００年地表以下不同深度 处的多年月平均地温曲线见图２。由图２可知，活动层温度在 ０℃上下季节性波动，且随着深度的增加，活动层温度伴有明 显的衰减和滞后现象。这是一个伴有相变的一维热传导问题。 以玛多县气象站地表以下０．０５～３．２０ ｍ为研究剖面，将 １９８０年１月至２０００年１２月设为模拟期，以地表以下０．０５ ｎｌ和
摘要：等效显热熔法是研究高寒地区冻土活动层温度场水文过程和生态变化的重要方法。通过引入描述冻土活动层 饱和程度的参数，改进了等效显热熔法，使其适用于非饱和带。在假设土壤含水量不变，忽略水分对流和热源（ｊｒ－）的前 提下，使用ＣＯＭＳＯＬ Ｍｕｌｔｉｐｈｙｓｉｃｓ模型模拟了黄河源区玛多县；东土活动层的温度状况。结果表明：改进的等效显热熔法 是合理、可信的，可以用于高寒地区冻土活动层变化规律与机制的研究。 关键词：冻土活动层；相变；数值模拟；等效显热熔法；高寒地区 中图分类号：ＴＵ４４５ 文献标志码：Ａ ｄｏｉ：１０．３９６９／ｊ．ｉｓｓｎ．１０００—１３７９．２０１３．０３．００９
ｍｍ。黄河源区特有的地理位置和地形、地貌、水文、干
式中：ｐ。＝妒５。；妒，＝筇，＝妒（１一Ｓ。）；ｃ。、ｃ，、ｃ；分别为水、冰、固
体骨架的比热容，Ｊ／（ｋｇ・Ｋ）；ｐ。、Ｐｉ、Ｐ。分别为水、冰、固体骨架 的密度；ｋｇ／ｍ３；Ｌ为水或冰发生相变转化时释放或吸收的潜 热，Ｊ／ｋｇ。 综合上述两方面，就可得到饱和土壤的等效体积热容量 Ｃ。，它是热焓相对于温度的变化率： ｃ。＝妒ｗｐ。ｃ。＋妒。ｐｉｃ。＋妒ｓｐ。ｃ。＋ｐ—ｉｊＬ斤ｄＳ－’．． （４）
（６）
ｃ。：ａ（妒。ｐ。ｃ。＋妒ｉｐ。ｃｉ）＋妒。ｐ。ｃ。＋—ｏｔｐ—ｊｉＬｒｄＳ“
是变化的，因此参数ｄ实际上是时间和空间的函数。
（７）
（２）
需要说明的是，由于不同深度、不同时刻的土壤饱和程度
式中：ｐ为孔隙度；驴。＝１一妒，表示固体颗粒占总体积的比例； Ａ。、Ａ。和Ａ。分别表示水、冰和固体骨架的热传导系数；５。和Ｓｉ 分别表示水和冰的相对饱和度，即ｓ。＋Ｓｊ＝１。 在处理相变潜热时，常将其看成是一个较大的显热熔一１。 通常的做法是在方程的右边添加一项关于相变潜热的焓 变¨…，即
收稿日期：２０１２—０９—１７ 基金项目：国享自然科学基金资助项目（４１０７２１９１）；中央高校基本稃研业务费 专项（２０ＩＯＺＹｌ２、２０１１ＹＸＬ０３９）。 作者简介：王俊智（１９８７一），男，河南郑州人，博士研究生，主要从事地下水与 环境研究工作。
Ｅ・ｍａｉｌ：ｗａｎｇｊｚ．ｃｕｇｂ＠ｇｍａｉｌ．ｔｏｍ
ｔｏ
ｃｏｎｖｅｃｔｉｏｎ
ａｎｄ ｔｈｅ
ｓｏｕｒｃｅ（ｓｉｎｋ）ｏｆ
ｈｅａｔ，ｔｈｉｓ ｐａｐｅｒ ｕｔｉｌｉｚｅｄ ＣＯＭＳＯＬ Ｍｕｈｉｐｈｙｓｉｃｓ
ｓｉｍｕｌａｔｅ ｔｈｅ ｔｈｅｒｍａｌ ｒｅｇｉｍｅ ｏｆ Ｍａｄｕｏ Ｃｏｕｎｔｙ’Ｓ