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

第３５卷第３期２０１３年３月人民黄河ＹＥＬＬＯＷＲＩＶＥＲＶ０１．３５．Ｎｏ．３Ｍａｒ．．２０１３【水文・泥沙】高寒地区冻土活动层温度场数值模拟王俊智１，梁四海１，万力１，王海周２（１．中国地质大学水资源与环境学院，北京１０００８３；２．河南省有色金属地质矿产局第五地质大队，河南郑州４５００１６）摘要：等效显热熔法是研究高寒地区冻土活动层温度场水文过程和生态变化的重要方法。通过引入描述冻土活动层饱和程度的参数，改进了等效显热熔法，使其适用于非饱和带。在假设土壤含水量不变，忽略水分对流和热源（ｊｒ－）的前提下，使用ＣＯＭＳＯＬＭｕｌｔｉｐｈｙｓｉｃｓ模型模拟了黄河源区玛多县；东土活动层的温度状况。结果表明：改进的等效显热熔法是合理、可信的，可以用于高寒地区冻土活动层变化规律与机制的研究。关键词：冻土活动层；相变；数值模拟；等效显热熔法；高寒地区中图分类号：ＴＵ４４５文献标志码：Ａｄｏｉ：１０．３９６９／ｊ．ｉｓｓｎ．１０００—１３７９．２０１３．０３．００９ＮｕｍｅｒｉｃａｌＳｉｍｕｌａｔｉｏｎｏｆｔｈｅＴｈｅｒｍａｌＲｅｇｉｍｅｏｆＰｅｒｍａｆｒｏｓｔＡｃｔｉｖｅＬａｙｅｒｉｎＡｌｐｉｎｅＲｅｇｉｏｎｓＷＡＮＧＪｕｎ—ｚｈｉｌ，ＬＩＡＮＧＳｉ—ｈａｉｌ，ＷａｎＬｉｌ，ＷＡＮＧＨａｉ．ｚｈｏｕ２（１．ＳｃｈｏｏｌｏｆＷａｔｅｒＲｅｓｏｕｒｃｅｓａｎｄＥｎｖｉｒｏｎｍｅｎｔ，ＣｈｉｎａＵｎｉｖｅｒｓｉｔｙｏｆＧｅｏｓｃｉｅｎｃｅｓ，Ｂｅｉｊｉｎｇ１０００８３，Ｃｈｉｎａ；２．ＴｈｅＦｉｆｔｈＴｅａｍ．ＨｅｎａｎＮｏｎｆｅｒｒｏｕｓＧｅｏ—ＥｘｐｌｏｒａｔｉｏｎＢｕｒｅａｕ，Ｚｈｅｎｇｚｈｏｕ４５００１６，Ｃｈｉｎａ）Ａｂｓｔｒａｃｔ：ＴｈｅｎｎａｔｒｅｇｉｍｅｏｆｐｅｒｍＭｒｏｓｔａｃｔｉｖｅｌａｙｅｒｈａｓａｇｒｅａｔｓｉｇｎｉｆｉｃａｎｃｅｔＯｏｎｒｅｖｅａｌｉｎｇｈｙｄｒｏｌｏｇｉｃａｌｐｒｏｃｅｓｓｅｓａｎｄｅｃｏｌｏｇｉｃａｌｃｈａｎｇｅｓｉｎＡｌｐｉｎｅａｒｅ・ｇｉｏｎｓ．Ｔｈｅｅｑｕｉｖａｌｅｎｔ—ａｐｐａｒｅｎｔｈｅａｔｃａｐａｃｉｔｙｗａｓａｎｉｍｐｏｒｔａｎｔｍｅｔｈｏｄｗａｓｓｔｕｄｙｔｈｉｓｉｓｓｕｅ．Ｂｙｉｎｔｒｏｄｕｃｉｎｇｔｏｐａｒａｍｅｔｅｒｄｅｓｃｒｉｂｉｎｇｔｈｅｓａｔｕｒａｔｉｏｎｏｆｔｈｅａｃｔｉｖｅｌａｙｅｒ，ｔｈｅｅｑｕｉｖａｌｅｎｔ—ａｐｐａｒｅｎｔｈｅａｔｃａｐａｃｉｔｙｍｅｔｈｏｄｔｕｒｅｍｏｄｉｆｉｅｄｂｅｃａｐａｂｌｅｏｆｕｎｓａｔｕｒａｔｅｄｚｏｎｅ．Ｏｎｔｈｅｐｒｅｍｉｓｅｏｆｉｇｎｏｒｉｎｇｓｏｉｌｍｏｉｓ—ｔｏｃｏｎｖｅｃｔｉｏｎａｎｄｔｈｅｓｏｕｒｃｅ（ｓｉｎｋ）ｏｆｈｅａｔ，ｔｈｉｓｐａｐｅｒｕｔｉｌｉｚｅｄＣＯＭＳＯＬＭｕｈｉｐｈｙｓｉｃｓｓｉｍｕｌａｔｅｔｈｅｔｈｅｒｍａｌｒｅｇｉｍｅｏｆＭａｄｕｏＣｏｕｎｔｙ’Ｓｐｅｒｍａ—ｆｒｏｓｔａｃｔｉｖｅｌａｙｅｒｉｎｔｈｅＹｅｌｌｏｗＲｉｖｅｒｃｒｅｄｉｔａｂｌｅ．Ｉｔｃａｎｓｏｕｒｃｅｒｅｇｉｏｎ．Ｔｈｅｒｅｓｕｌｔｓｓｈｏｗｔｈａｔｔｈｅｍｏｄｉｆｉｅｄｅｑｕｉｖａｌｅｎｔ—ａｐｐａｒｅｎｔｈｅａｔｃａｐａｃｉｔｙｍｅｔｈｏｄｉｓｅｆｆｉｃｉｅｎｔａｎｄｂｅｕｓｅｄｔｏｓｔｕｄｙｔｈｅｐｅｒｍａｆｒｏｓｔａｃｔｉｖｅｌａｙｅｒｖａｒｉａｔｉｏｎｍｅｃｈａｎｉｓｍｓｏｆｔｈｅａｌｐｉｎｅｒｅｇｉｏｎ．ｒｅｇｉｏｎＫｅｙｗｏｒｄｓ：ｐｅｒｍａｆｒｏｓｔａｃｔｉｖｅｌａｙｅｒ；ｐｈａｓｅｃｈａｎｇｅ；ｎｕｍｅｒｉｃａｌｓｉｍｕｌａｔｉｏｎ；ｅｑｕｉｖａｌｅｎｔ－ａｐｐａｒｅｎｔｂｅａｔｃａｐａｃｉｔｙｍｅｔｈｏｄ；ａｌｐｉｎｅ高寒地区地下水的冻结与融化直接影响区域水文过程和源区玛多县为例，利用ＣＯＭＳＯＬＭｕｈｉｐｈｙｓｉｃｓ模拟了该地区冻土活动层的温度状况。１生态过程…。近年来，在全球气候变暖的背景下，高寒地区冻土发生了退化现象旧ｊ。冻土活动层温度场研究对于揭示全球气候变暖条件下，高寒地区的水文过程和生态变化有着十分重要的意义。２０世纪７０年代以前，冻土温度的研究主要是在理论分析基础上的近似解析法和以试验、野外勘探与调查资料为基础的统计方法∞Ｊ。２０世纪７０年代，Ｈａｒｌａｎ首先建立了水热耦合的数学模型Ｈ１，后来，许多学者基于Ｈａｒｌａｎ的思想建立了冻土水热耦合数值模型。国内学者很早就开展了对冻土温度的研究。庞强强等Ｍ１根据１９９１—２０００年的地面温度观测资料并结合数字高程模型，应用斯蒂芬公式得到了青藏高原冻土区的活动层厚度分布图。李述训等¨ｏ应用数值方法预测了气候变暖条件下，青藏高原多年冻土的温度状况。郝振纯等…建立了季节冻土水热耦合迁移的数学模型，并对黄河源区站点的冻融深度进行了模拟。近年来，等效显热熔法越来越多地被应用于相变问题的研究。笔者根据这一方法的原理，通过引入描述冻土活动层饱和程度参数，对该方法进行改进，使其适用于非饱和带，并以黄河・２０・伴有相变的非饱和一维热传导数学模型１．１热传导控制方程忽略水分的对流，地温丁随时间ｔ的变化可以用经典热传导方程来描述。８ｏ：ｖ（ＡＶＴ）＋ｈ＝Ｃ。警Ｕ‘（１）式中：Ａ为热传导系数，Ｗ／（ｎｌ・Ｋ）；ｈ为热源（热汇）；Ｃ。为体积热容量，Ｊ／（Ｋ・ｉｎ３）。１．２等效显热熔法等效显热熔法就是在处理热参数时考虑各个组分在热传收稿日期：２０１２—０９—１７基金项目：国享自然科学基金资助项目（４１０７２１９１）；中央高校基本稃研业务费专项（２０ＩＯＺＹｌ２、２０１１ＹＸＬ０３９）。作者简介：王俊智（１９８７一），男，河南郑州人，博士研究生，主要从事地下水与环境研究工作。Ｅ・ｍａｉｌ：ｗａｎｇｊｚ．ｃｕｇｂ＠ｇｍａｉｌ．ｔｏｍ万方数据人民黄河２０１３年第３期导过程中所起的作用，采用混合理论，按照所占体积的比例构成一个等效热参数代入控制方程中进行计算。对于饱和土壤，等效热传导系数Ａ由水、冰、固体骨架颗粒通过体积平均共同表示：Ａ＝妒Ｓ。Ａ。＋妒Ｓ。Ａ．＋妒。Ａ。（２）式中：ｐ为孔隙度；驴。＝１一妒，表示固体颗粒占总体积的比例；Ａ。、Ａ。和Ａ。分别表示水、冰和固体骨架的热传导系数；５。和Ｓｉ分别表示水和冰的相对饱和度，即ｓ。＋Ｓｊ＝１。在处理相变潜热时，常将其看成是一个较大的显热熔一１。通常的做法是在方程的右边添加一项关于相变潜热的焓变¨…，即△珥＝Ｊ（妒。Ｐ。ｃ。＋妒护。ｃ．＋妒；Ｐ。ｃ。）ｄＴ—ｆＰｉ￡ｄ妒．（３）式中：ｐ。＝妒５。；妒，＝筇，＝妒（１一Ｓ。）；ｃ。、ｃ，、ｃ；分别为水、冰、固体骨架的比热容，Ｊ／（ｋｇ・Ｋ）；ｐ。、Ｐｉ、Ｐ。分别为水、冰、固体骨架的密度；ｋｇ／ｍ３；Ｌ为水或冰发生相变转化时释放或吸收的潜热，Ｊ／ｋｇ。综合上述两方面，就可得到饱和土壤的等效体积热容量Ｃ。，它是热焓相对于温度的变化率：ｃ。＝妒ｗｐ。ｃ。＋妒。ｐｉｃ。＋妒ｓｐ。ｃ。＋ｐ—ｉｊＬ斤ｄＳ－’．．（４）１．３饱和土壤冻结曲线在等效体积热容量ｃ。的表达式中，等式右边第４项中水的相对饱和度Ｓ。与温度ｒ之间的关系将会对热传导过程产生重要影响。水的相对饱和度Ｓ。与温度ｒ之间所确定的曲线一般称为土壤冻结曲线，见图１。Ｂｏｎａｃｉｎａ等证明，在进行温度场计算过程中，水分冻结函数的实际形状并不是那么重要，但是它一定要满足潜热的条件¨“，即ｒｒｔＬ＝Ｊ（ｐｃ）。ｄＴ（５）ｏｔ式中：孔为融化温度；ｔ为冻结温度；（ｐｃ）。为显热熔。∞堡胁罄‘犀∞誊墨船＊ｍＯ［一Ｌ图１土壤水分冻结曲线此外，在等效体积热容量Ｃ。的表达式中，由于水的相对饱和度ｓ。与温度７１之间的关系是导数关系，因此，所选取的饱和土壤冻结函数要更加光滑。１．４非饱和等效显热熔法对于冻土活动层为非饱和带的情况，还需要对上述理论进行改进。由于气体的热传导系数和密度相对于固体或液体一般小２～３个数量级，在实际热传导模拟过程中，气体所产生的热效应往往可以忽略，只需要考虑水、冰和固体骨架对实际热传导的贡献。因此，对于非饱和情形，引进参数Ｏｌ描述冻土活动层的饱和程度，即水和冰体积之和所占孔隙的比例，来剔除气体的影响。因此，非饱和情形下的等效热传导率和等效体积万方数据热容量可以表示为Ａ＝Ｏｔ（妒Ｓ。Ａ。＋妒ＳｉＡ，）＋妒。Ａ。（６）ｃ。：ａ（妒。ｐ。ｃ。＋妒ｉｐ。ｃｉ）＋妒。ｐ。ｃ。＋—ｏｔｐ—ｊｉＬｒｄＳ“（７）需要说明的是，由于不同深度、不同时刻的土壤饱和程度是变化的，因此参数ｄ实际上是时间和空间的函数。２黄河源区季节冻土温度模拟２．１研究区概况黄河源区位于青藏高原东部，面积约为２０９３０ｋｍ２，海拔在４２００ｍ以上，相对高差大于ｌ０００ｍ。源区年平均气温为一３．８９℃，平均降水量为３１２．２ｍｍ，年潜在蒸发能力为１３５３．２ｍｍ。黄河源区特有的地理位置和地形、地貌、水文、干寒的气候条件决定了该区为季节冻土，并形成岛状和大片连续多年冻土并存的分布格局Ｈ“。近年来，在全球气候变暖的背景下，黄河源区出现了地下径流减小、植被覆盖程度降低、冻土退化等一系列生态地质问题¨…。２．２模型概化黄河源区玛多县气象站１９８０－－２０００年地表以下不同深度处的多年月平均地温曲线见图２。由图２可知，活动层温度在０℃上下季节性波动，且随着深度的增加，活动层温度伴有明显的衰减和滞后现象。这是一个伴有相变的一维热传导问题。以玛多县气象站地表以下０．０５～３．２０ｍ为研究剖面，将１９８０年１月至２０００年１２月设为模拟期，以地表以下０．０５ｎｌ和３．２０Ｉｎ的实测地温为边界条件，以１９８０年１月的实测剖面地温为初始条件，在假定土壤含水量不变，忽略水分对流和热源（汇）的情况下，使用ＣＯＭＳＯＬＭｕｈｉｐｈｙｓｉｃ建立数值模型。１５１０５ｐ妪０赠。一１０—１５图２１９８０－－－２０００年玛多县多年月平均地温（单位：ｍ）使用ＣＯＭＳＯＬＭｕｈｉｐｈｙｓｉｃｓ地球科学模块中预定义的多孔介质热传导应用模式，通过对该模式进行一定程度地修改来描述上述伴有相变的非饱和一维热传导模型。在实际使用过程中，首先选取ＣＯＭＳＯＬＭｕｈｉｐｈｙｓｉｃｓ中具有二阶连续导数的ｆｌｃ２１协函数来表示冻结函数¨“，建立含有相变潜热的焓变表达式；然后通过定义流体材料属性和固体材料属性，给出等效热传导系数和等效体积热容量；最后将模型剖分为４５６个单元进行求解。２．３参数分析在冻土温度计算中涉及的有关参数很多，这些参数的测量较为复杂，而且测量价格昂贵，在实际科研和工程计算中每次都进行测量就会显得十分不经济。本次模拟将冻土活动层看．２】・人民黄河２０１３年第３期成是均质各向同性的地质体，相关冻土物理参数根据冻土活动层的土壤性质，参照文献［１５］给出（表１）。表１模拟中使用的参数参数｝Ｌ隙度妒水的比热ｃ。／（Ｊ・ｋｇ～・Ｋ。）冰的比热ｃｉ／（Ｊ・ｋｇ．。・Ｋ－１）固体骨架的比热ｃ。／（Ｊ・ｋｇ。１・Ｋ－１）水的热传导率Ａ。／（Ｗ・ｍ叫・Ｋ。１）冰的热传导率Ａ／（Ｗ・ｍ～・Ｋ－１）固体骨架的热传导率Ａ／（Ｗ・ｍ．１・Ｋ一）水的密度Ｐ。／（ｋｇ・ｎｌ。）冰的密度ｐｉ／（ｋｇ・ｍ。）固体骨架的密度ｐ。／（ｋｇ・ｍ。３）相变潜热∥（Ｊ・ｋｇ“）冻结温度ｒ，／℃相变区间△∥℃未冻残余饱和度Ｊｓ。。水的相对饱和度ｏ／％数值活动层对地表气温的变化十分敏感。结合１９８０－－２０００年该地区的气温年平均变化曲线（图５）可知，１９８０－－１９８５年，气温呈现下降趋势，活动层温度下降，冻土有向深部发展的趋势；１９８５－－１９９０年，气温呈现回升趋势，活动层温度升高，冻土退化；１９９０－－１９９５年，气温呈现下降趋势，活动层温度下降，冻土得到一定程度的恢复；由于２０００年之前的２ａ气温较高，因此２０００年地温较１９９５年有小幅升高，在滞后作用的影响下，冻土轻度退化。１９８０－－２０００年气温呈现上升趋势，活动层温度升高，冻土发生了退化。嘣一撇啪删哪呲一圳。¨卯３模拟结果分析３．１地温曲线采用改进的等效显热熔法，利用ＣＯＭＳＯＬＭｕｌｔｉｐｈｙｓｉｃｓ模拟得到了不同深度、不同时刻的冻土活动层温度。为了检验模拟效果，分别选取地表以下０．４０、０．８０、１．６０ｎｌ的模拟结果与实测数据进行对比，见图３。总体来说，模拟结果还是不错的，与实测数据相比仅在局部存在微小差异。ｉ◆薹薹◆０旺５１０ｇ燕Ｉｓ鞋＝２０Ｚ５３．０２３４５ｉ６７８９１０月份ｃ：１ｑ９０年ｎ０５（媳。薹飘＾Ｍ＾小＾．＾八八＾＾从＾』＼ＭＭ，菲礴：［０毒”卿２０２ｓ３．０藜≈一：☆蔼丽而杰荫而蔚高再打南面菇葫而２（。髋。——宴测值一一一模拟值图３地表以下０．４０、０．８０、１．６０ｉｎ处地温实测值与模拟值２３４５６７８９１０月份０ｆｄｊ：９９５年３．２地温场等值线不同时刻、不同深度的地温场等值线图是分析地温状况的有效工具。实测地温数据采样点有限，用其绘制出的地温场等值线图分辨率较低，采用模拟结果则不存在这一缺憾。选取１９８０年、１９８５年、１９９０年、１９９５年、２０００年为特征年，利用模拟数据绘制该地区地温场等值线，见图４，并对１９８０－－２０００年的地温状况进行分析。地温场等值线表明：黄河源区玛多县地区一般在８月、９月活动层地温开始降低，１０月底至翌年４月自浅表向深部逐渐冻结，６月左右冻结深度达到最大；一般在４月初活动层自地表开始向深部融化，１０月初融化深度达到最大。・２２・图４Ｏ５１．０郭２。！ｊ３０１３４５６７８９】０１ｌ１２月份ｆｅ）２０００年１９８０、１９８５、１９９０、１９９５年和２０００年地温场等值线（单位：℃）从地温场等值线图上，还可以得到活动层每年的最大冻结深度。万方数据人民黄河２０１３年第３期年份图５１９８０－－２０００年玛多县气温年平均变化曲线由于毛细水压力和溶解在孔隙水中盐分的存在，因此冻土活动层的冻结温度往往低于纯水的冰点¨“。考虑到上述因素，以一０．５℃为活动层土壤冻结临界温度，结合图４可知，黄河源区玛多县地区冻土活动层的最大冻结深度一般在６月出现，深度约为２～３ｍ。１９８０－－２０００年，由于活动层温度升高，因此冻土发生了退化。将模拟得到的最大冻结深度与实测数据相比较可知（图６），两者趋势一致，相关系数为０．４４７。巨挂毯诺蚂－蛞Ｋ《年份图６１９８０－－２０００年玛多县冻土活动层年最大冻结深度４讨论冻土活动层的冻结融化过程受自身物理化学性质和地表气温等因素的共同影响。活动层是一个非常复杂的地质体，其含水量、孔隙度、土壤热参数等随时间和空间不断变化，而模型中对这些参数的处理却较为简单，导致模拟结果与实测数据相比存在一定偏差。这是该模型改进时需要考虑的地方。总体上看，该模型虽然存在一定不足，但是对于冻土活动层温度场的研究来说是十分有效的，可以用于高寒地区冻土活动层变化规律与机制的研究，如黄河源区活动层的季节冻融与年际变化过程；冻土活动层年最大冻结深度的变化过程；全球气候变暖条件下黄河源区的冻土退化等。５结语冻土活动层温度场的研究对于揭示全球气候变暖条件下，高寒地区的水文过程和生态变化有着十分重要的意义。等效万方数据显热熔法是研究高寒地区冻土活动层温度场水文过程和生态变化的重要方法。由于冻土活动层多为非饱和带，需要引入描述活动层饱和程度的参数来改进等效显热熔法。基于改进的等效显热熔法，使用ＣＯＭＳＯＬＭｕｈｉｐｈｙｓｉｃｓ建立了伴有相变的一维热传导数值模型模拟了黄河源区玛多县冻土活动层的温度状况，取得了较好的效果。由于冻土活动层的含水量、孔隙度、土壤热参数等是随时间和空间不断变化的，而模型中对这些参数的处理却较为简单，因此模拟结果与实测数据相比存在一定偏差。这些是将来模型改进时需要考虑的。总体看来，该模型是合理、可信的，可以用于高寒地区冻土活动层变化规律与机制的研究。参考文献：［Ｉ］曹文炳，万力，周训，等．黄河源区冻结层上水地质环境影响研究［Ｊ］．水文地质工程地质，２００３，３０（６）：６—１０．［２］周训，金晓梅，梁四海，等地下水科学专论［Ｍ］．北京：地质出版社，２０１０．［３］费里德曼，冻土温度状况计算方法［Ｍ］．徐学祖，译．北京：科学出版社，１９８２．１４１Ｈａｒｌａｎ．Ａｎａ／ｙｓｉｓｏｆＣｏｕｐｌｅｄＨｅａｔ—ＦＭｄＴｒａｎｓｐｏｒｔｉａＰａｒｔｉａＵｙＦｒｏｚｅｎＳｏｉｌ［Ｊ］．ＷａｔｅｒＲｅｓｏｕｒｃｅＲｅｓｅａｒｃｈ，１９７３，９（５）：１３１４—１３２３．［５ｊ庞强强，李述训，吴通华，等．青藏高原冻土区活动层厚度分布模拟［Ｊ］．冰川冻土，２００６，２８（３）：３９０—３９５．［６］李述训。程国栋．气候变暖条件下青藏高原高温冻土热状况变化趋势数值模拟［Ｊ］．冰川冻土，１９９６（增刊１）：１９０—１９６［７］郝振纯，张晓鹏，张磊磊，等．气候变暖下黄河源区冻土变化的数值模拟［Ｊ］．黑龙江水专学报，２００９，３６（３）：ｔ００—１０４［８］ＣａｒｓｌａｗＨＳ，ＪａｅｇｅｒＪＣ．ＣｏｎｄｕｃｔｉｏｎｏｆＨｅａｔｉｎＳｏｌｉｄｓ［Ｍ］．Ｏｘｆｏｒｄ：ＣｌａｒｅｎｄｏｎＰｒｅｓｓ，１９４７，ｆ９］ＧｏｏｄｒｉｃｈＬＥ．ＡｎＩｎｔｒｏｄｕｃｔｏｒｙＲｅｖｉｅｗｏｆＮｕｍｅｒｉｃａｌＭｅｔｈｏｄｓｆｏｒＧｒｏｕｎｄＴｈｅｒ．ｍａｌＲｅ，ｍｅＣ，，ｄｃｕｌａｔｉｏｎｓ［Ｒ］．Ｏｔｔａｗａ：ＮａｔｉｏｎａｌＲｅｓｅａｒｃｈＣｏｕｎｃｉｌｏｆＣａｎａｄａＤｉｖｉｓｉｏｎｏｆＢｕｉｌｄｉｎｇＲｅｓｅａｒｃｈ，１９８２．［１０］ＭｏｔｔａｇｈｙＤ，ＲａｔｈＶＬａｔｅｎｔＨｅａｔＥｆｆｅｃｔｓｉｎＳｕｂｓｕｒｆａｃｅＨｅａｔＴｒａｎｓｐｏｒｔＭｏｄｅｌ・ｉｎｇａｎｄＴｈｅｉｒＩｍｐａｃｔｏｎＰａｌａｅｏ—ＴｅｍｐｅｒａｔｕｒｅＲｅｃｏｎｓｔｒｕｃｔｉｏｎｓ［Ｊ］．Ｇｅｏｐｈｙｓｉ—ｃａｌＪｏｕｍａｌＩｎｔｅｒｎａｔｉｏｎａｌ，２００６，１６４（１）：２３６—２４５．［ｔ１ｊＢｏｎａｅｉｎａＧ，Ｃｏｍｉｎｉ．ＯｎｔｈｅＳｏｌｕｔｉｏｎｏｆｔｈｅＮｏｎｌｉｎｅａｒＨｅａｔＣｏｎｄｕｃｔｉｏｎＥｑｕａ—ｔｉｏｎｓｂｙＮｕｍｅｒｉｃａｌＭｅｔｈｏｄｓ［Ｊ］．ＩｎｔｅｒｎａｔｉｏｎａｌＪｏｕｒｎａｌｏｆＨｅａｔａｎｄＭａｓｓＴｒａｎｓ－ｆｅｒ，１９７３，１６：５８１—５８９［１２］金会军，王绍令，吕兰芝，等．黄河源区冻土特征及退化趋势［Ｊ］．冰川冻土，２０１０，３２（１）：１０一１７［１３］赵云云，赵其华黄河源头多年冻土退化原因及变化趋势［Ｊ］．人民黄河，２００９，３１（６）：１０—１２．［１４］中仿科技．专业数值分析系统ＣＯＭＳＯＬＭｕｌｔｉｐｈｙｓｉｃｓ［Ｊ］．ＣＡＤ／ＣＡＭ与制造业信息化，２００８（９）：４０—４４．［】５］徐学祖，王家澄，张立新．冻土物理学［Ｍ］．北京：科学出版社，２００１．【责任编辑吕艳梅】・２３・