基于能量平衡的融雪期积雪深度模拟

黄东浩，王 鹏，周丽丽，孙宏超

（沈阳农业大学水利学院/辽宁省水土流失防控与生态修复重点实验室，沈阳 110161）

东北黑土区地处我国中高纬度地区，雪期长度约8个月[1]，降雪量占年降水量7%～25%，由此产生的融雪径流量和输沙模数分别占全年的13.3%～24.9%和5.8%～27.7%[2]，对该区域水土可持续利用具有重要影响。随着全球气候变暖，该区融雪期提前与土壤冻融循环次数增多，从而加剧土壤侵蚀[3-4]；
同时消融早期的融雪径流明显增加，导致后期积雪减少，改变了流域融雪径流年内分配[5]，对流域内水资源的管理、规划及可持续利用造成困难[6]。因此，准确预报积雪消融过程对该区水资源管理及融雪侵蚀研究具有重要的指导意义。

积雪消融预报主要依靠融雪模型，由于模型输出表达方式的不同，模型的复杂程度也有所不同[7]，如CROCUS[8]模型可以模拟雪内部结构变化结构特征。就目前而言，国内外积雪消融预报模型主要分为两类：一类是基于度日因子的概念模型；
另一类是基于物理过程的能量模型。度日因子模型由于输入参数少、计算简单而被广泛应用[9]，但度日模型存在的缺点[10]为：(1)随着时间的推移，模拟精度明显下降，不适用于长时间序列的融雪模拟；
(2)不能定量描述积雪内部消融状况。针对以上不足，一些学者进行了模型优化，AKBARI 等[11]通过优化积雪覆盖面积和引入流量沉降系数，来提高度日模型的模拟精度；
穆艾塔尔·赛地等[12]基于地形因子将流域分带后引入度日模型来模拟积雪消融变化。但就本质而言度日因子模型是一个数学统计模型，由于积雪消融过程中雪表面时空差异较大，该模型无法描述积雪内部的消融过程，若有降雨输入时还会降低精度[13]。而基于能量平衡的融雪模型可定量描述融雪过程，并渐成研究趋势[14]，模型输入参数要求也较高[15]，目前研究方向主要有：(1)将能量平衡融雪模型与其他模型集成，孟现勇等[16]将基于能量平衡的融雪模型耦合到SWAT模型融雪模块，提升了雪深、径流量的模拟精度；
CHEN 等[17]利用SNTHERM-MEMLS 耦合模型来模拟雪深、雪温，取得了较好的模拟结果。(2)提高模型参数的准确性，闻昕等[18]利用机器学习的方法，从实测数据自动提取融雪模型参数来提高模型精度。(3)建立分布式融雪模型[19]，ABBOTT等[20]采用分布式的能量平衡融雪模型（SHE）来反应融雪在空间上的差异性；
BROWN等[21]开发了分布式融雪模型UEBGRID在喜马拉雅地区进行模拟试验，发现模型对雪深有较好的表现。(4)结合地形与覆盖物的参数，开发一种基于时间和空间的融雪模型，这种模型在融化过程中考虑了降雨、辐射大气湍流等因素[22]。

近年来，关于东北黑土区积雪消融过程的研究中取得了一定进展，为融雪径流预测提供了科学依据。杨鑫等[23]利用黑龙江鹤山农场10年气象资料采用线性回归方法分析积雪深度和融雪径流深的关系。田琳[24]基于遗传算法对松花江白山流域构建积雪消融变化经验方程，并预测径流量。整体而言，东北黑土区也多以遥感的方法和统计学方法对积雪深度进行反演模拟[25]，太阳辐射对积雪消融过程的影响以及积雪融化的物理过程等方面考虑较少。

本研究在吉林省梅河口市吉兴小流域的水土保持监测站完成数据的监测与采集，吉兴小流域面积约为15 km2，地理坐标为125°28′13″～125°31′34″E、42°10′26″～42°14′32″N，海拔高度范围为393～970 m，该区地形起伏较大，总体呈南高北低的趋势，为东北低山丘陵区。属于典型的温带大陆性季风气候，多年平均降水量约为708.8 mm，冬季积雪深度约为300 mm，降雪期集中在12 月到次年2 月份，降雪日数在6～32 d，为典型的季节性冻融区。年平均气温4.6 ℃，全年≥10 ℃积温达到2 732 ℃，多年平均日照时数2 556 h。

在研究区吉兴小流域站点架设PC400 型雪与土壤多参数系统采集土壤温度、雪深、土水势数据；
气象生态环境监测仪采集太阳辐射、风速风向、土壤湿度、空气温度、相对湿度数据；
使用Snow Fork雪特性分析仪对雪液态含水量、雪密度进行数据采集，观测时间间隔为60 min，试验仪器基本情况见表1。

表1 试验仪器情况Table 1 Experimental equipment situation

2.1 模型计算原理

模型基于能量平衡计算（图1为雪面能量平衡示意图），通过输入短波辐射、空气温度、风速（2 m）和相对湿度等状态变量来驱动模型运行。模型的能量过程用于模拟融雪、再冻结以及积雪的热量变化，模拟计算后可以得到雪面净辐射通量、显热通量、潜热通量、总能量和雪深。

图1 雪面能量平衡示意图[26]Figure 1 Schematic diagram of the energy balance of the snow surface

2.1.1 基本能量平衡方程 能量平衡模型：

式中：Qnet为模拟时段Δt内积雪层单位面积上总能量的净输入量(W·m-2)；
Qr为雪面净辐射通量(W·m-2)；
Qh为显热通量(W·m-2)；
Qle为潜热通量(W·m-2)；
Qp为融雪期降水所带来的热通量(W·m-2)；
Qg为来自土壤的地热通量(W·m-2)。

雪面净辐射通量（Qr）的计算：

式中：Qs为入射的短波辐射(W·m-2)；
α为雪面反射率；
εs为积雪比辐射率；
Ta为空气温度(℃)；
Ts为雪面温度(℃)；
εαc为有云覆盖下的大气比射率；
σ为STEFANBOHZMAN 常数。雪面反射率是影响雪表面能量平衡的重要参数，采用Noah陆面模式的积雪反照率方案[27]：

式中：αmax为新雪反照率；
t为距离降雪的天数；
A和B为常数，分别取0.94 和0.82。本研究中采取对于雪面温度的计算，Abdelghani发现当温度在零下5 ℃时，雪面温度和空气温度存在线性关系[28]，因此雪面温度采用高黎明等[29]建立的线性关系式：

显热通量和潜热通量采用在雪面能量平衡的运用较好的空气动力学梯度方法[30]，该模型计算简单且计算表现较好。

显热通量（Qh）的计算：

式中：Qh＞0时，能量从大气传向雪层，雪面获得能量；
Qh＜0时，能量从雪层传向大气，雪面损失能量。ρ为空气密度；
Cp为空气的比热容；
ras为积雪层表层与附近的参考高度间的空气动力学阻力[31]：

式中：χ为参考高度；
μ为参考高度处的风速；
Zm为动力阻抗系数；
Zh为热量和水汽阻抗系数；
k为VON KARMAN常数。

潜热通量（Qle）的计算：

式中：Qle＞0时，能量从大气传向雪层；
Qle＜0时，能量从雪层传向大气；
λv为升华潜热；
Rd为干空气常数；
e(Ta)为空气水汽压；
e(Ts)为雪面水汽压；

e(Ta)和e(Ts)的计算：

esvp(Ta)为饱和水汽压，采用泰登公式[32]。

esvp(Ta)的计算：

降雨携带输入能量(Qp)：

式中：ρw为水的密度；
Cw为水的比热容；
Pr为降雨量；
Qp＜0表示从雪层吸热。在模拟计算时段内基本都为固态降水，因此本研究不考虑降雨带来的热量。

在融雪期间，雪面温度一般不超过0 ℃，且此时土层为冻土层，雪层底部的热传导度比较低，因此地热通量对积雪消融的影响基本可以忽略不计。

融雪深度h计算：

式中：h为融雪深度(cm)；
S为雪的溶解热；
ρs为积雪密度。

能量闭合率[33]：

式中：CR为能量闭合率。CR值越接近1，能量闭合率越好。

模型评价：

纳什系数（Nash）：

平均绝对误差（MAE）：

式中：n为观测次数；
H0为模拟值；
Hp为观测值；
Have为观测平均值；
纳什系数（Nash）反映模拟值与观测值的比值与1∶1 曲线的靠近程度，纳什系数（Nash）越接近1，相对误差的绝对值越小，模拟精度就越高；
平均绝对误差（MAE）越接近0，模拟的效果就越好。

2.2 数据来源及模型参数取值

本研究选取2017年2月23日-3月10日和2018年3月6日-3月23日的积雪消融过程为分析时段。模型建立的基本数据主要为野外实测数据、雪面反照率初始值采用中国全球陆面再分析40年产品（CRA/Land），模型中参数取值见表2。

表2 参数取值Table 2 Parameter values

3.1 能量辐射模拟结果分析

积雪消融是热传导过程，而净辐射通量是反应地表能量、动量、水分输送过程的重要参数[34]，常用于计算积雪消融过程。通过模拟2017-2018 年融雪期能量(图2)发现，净辐射昼夜变化很大，日间为正值，最高可达175.6 W·m-2；
夜间为负值，最低可达-79.4 W·m-2。经计算，2017 年和2018 年净辐射通量日均值占日总能量均值的变化趋势存在差异，即2017呈逐渐减小趋势；
2018年变化趋势为先减小后增大。2017年净辐射通量日平均值为117.8 W·m-2，占总能量的74.7%，最高达80%；
2018年日均值为122.5 W·m-2，占总能量的67.4%，最高达到93%。可见净辐射通量是主导能量平衡的主要部分。总能量的输入呈周期性的变化，日间净输入为正，11点达到最大值，两年消融期的净能量峰值均在150 W·m-2以上，2018年的峰值达到200 W·m-2；
夜间基本上为负值，最低可达-76.65 W·m-2。

图2 2017年和2018年净辐射通量、总能量模拟值Figure 2 Simulated values of net radiative flux and total energy from 2017 to 2018

然而受湍流交换差异的影响，2018 年净辐射通量对积雪消融能量的贡献较2017 年有所下降，降幅达7.2%。湍流交换是雪表面能量交换过程的主要形式，并受到显热通量和潜热通量影响。其中，显热通量受温度控制，体现近地面空气与雪层表面的热交换。潜热通量则取决于水气压梯度、空气气压梯度，表现为雪的升华。显热通量和潜热通量是积雪消融的主要影响因子，受天气类型的影响较大。风速平均值均在2 m·s-1以下，2018 年风速最大值达到6.9 m·s-1；
温度平均值基本保持在0 ℃以下，2018 年温度最高值达到16.1 ℃。可见，2018年风速、温度的平均值和最大值以及变化幅度皆均比2017年高（表3和图3）。

图3 201年和2018年温度、风速实测值；
显热通量、潜热通量模拟值Figure 3 Measured values of temperature,wind speed,simulated values of sensible heat flux and latent heat flux from 2017 to 2018

表3 2017年和2018年温度、风速的平均值与最大值Table 3 Average and maximum values of temperature and wind speed from 2017 to 2018

2017 年的风速和温度变幅较2018 年小，因此显热通量和潜热通量的也比较稳定。2017 年和2018 年显热通量平均值为正值25.6 W·m-2和35.6 W·m-2，潜热通量为负值为-8.9 W·m-2和-9.0 W·m-2，计算结果较为接近。但整体表现为2018年要比2017年更大。显热通量的值为正时，空气向雪面输送热量，这是因为空气温度高于雪面温度，会形成稳定的逆温层[35]。潜热通量为负值时，说明此时雪表面正在以升华和蒸发的形式消融，向空气输送热量。2018年3月15日至23日出现了潜热通量为正值的时段，这是因为此时饱和水气压大于空气水汽压，温差和风速也较小，湍流交换量小，不利于雪的升华和蒸发，此时水汽发生凝结，表现为空气向雪面输送热量[36]。

在冬季，地表与大气之间存在着复杂的物质和能量交换过程。其中一个重要的指标是能量闭合率，它反映了地表所接收到的净辐射与地表释放出去的感热、潜热和土壤热之间是否达到平衡。理想情况下，这两者应该相等，但实际观测中往往存在一定程度的不平衡。图4为积雪消融过程中步长为1 h数据计算的能量闭合率日平均值。由图4 可知，2017-2018 年能量闭合率随积雪深度的减小而减小，这是因为雪的反照率比裸露地表高，当雪深减少时，裸露面积增加，反照率降低，会导致能量闭合率变低，能量闭合率最大可达0.25～0.67。2018 年能量闭合率要优于2017 年，2017 年能量闭合率低的原因可能是消融周期时间长，积雪覆盖的天数较多，导致积雪表面吸收的太阳辐射量减少，因为积雪会反射大量的入射太阳辐射[37]，这将导致能量闭合率值降低。此外，较长的消融周期还会导致积雪表面形成冰层，增加雪表面反射率也会致使能量闭合率较低。

图4 2017年和2018年雪深与能量闭合率关系Figure 4 Relationship between snow depth and energy closure rate from 2017 to 2018

3.2 积雪深度模拟

由能量辐射模拟结果可知，净辐射能量是积雪消融的主要驱动力，净辐射通量在日间为正值，因此积雪的消融也主要发生在日间。图5 描述了2017 年3 月7日的能量变化和积雪的消融状况，该天能量与积雪变化明显，能够较好反映积雪消融过程。由图5 可知，净辐射通量是先增大后减小的趋势，在12 点达到最大值152 W·m-2，夜间为负值。显热通量和潜热通量11 点达到了极值，分别为25.34 W·m-2和-24.44 W·m-2，但显热通量和潜热通量出现了时间上明显的滞后性，滞后约1 h，这是因为能量主要来自太阳辐射，热传输速度不同步。在出现极值后至23∶00，显热通量先减小后增大，潜热通量先增大后减小。积雪从8∶00 开始消融，在12∶00 雪深急剧减小，融雪速率达到峰值，17∶00 后积雪消融开始变缓，直至18∶00完全停止。整体看来，净辐射通量、显热通量、潜热通量集中在9∶00-15∶00变化，共计7 h，与雪深变化的时间相吻合。

图5 净辐射通量、显热通量、潜热通量、雪深的日变化值Figure 5 Diurnal variation of net radiation flux,sensible heat flux,latent heat flux and snow depth

2017-2018 年融雪深度实测值与模拟值对比，以及1∶1 线图见图6。由图6 可知，2017-2018 年融雪期积雪深度的模拟值和实测值总体趋势较为一致。2017 年实测值与模拟值整体均匀分布在1∶1 线两侧，模拟结果较好（Nash=0.97；
MAE=0.36 cm）；
2018 年在雪深为小于15 cm 时，模拟值较高于实测值（Nash=0.83；
MAE=1.46 cm），尽管2018 年模拟结果可满足精度要求，但模拟值出现偏高或者偏低的情况，其主要原因归结为：（1）本研究采用的是单层融雪模型，将融雪层作为一维的整体考虑，当雪层能量为正且输入达到临界状态时，积雪才开始消融，而实际上整个雪层能量输入还未达到临界值时，表层可能会先消融，此时表层液态水受重力作用向下迁移，使模拟值比实际值略微偏高[38]；
（2）另外，暴露在空气中的雪表面产生较强的湍流交换也是影响雪深变化的一个重要原因。

图6 2017年和2018年实测雪深与模拟雪深对比Figure 6 Comparison of measured and simulated snow depth from 2017 to 2018

在积雪日变化方面，张娟等[39]对三江源地区研究表明，积雪消融时间发生在10∶00-19∶00，历时9 h；
桑婧等[40]在额尔古纳的研究则发现开始融雪时间在2∶00，消融时间近17 h。然而，本研究区积雪消融时间集中在净辐射通量、显热通量、潜热通量日内集中变化时间段，融化时间在9∶00-15∶00，历时7 h。在积雪消融天数方面，周杨等[41]在青藏高原玛多地区研究结果表明集中积雪消融天数为12～18 d，而本研究区融雪时间在9～15 d。综上所述，可以看出各地区之间积雪消融时段差别明显。

在模拟积雪深度变化方面，东北黑土区以往的研究多以空气温度等指标研究积雪变化规律，此方法属于统计学的方法，不能描述雪面能量收支特征。而本研究通过能量平衡原理建立的积雪消融模型能够更加准确描述雪面能量收支与积雪深度变化的特征。该模型可为季节性冻融区建立融雪模型的建立提供理论依据。本模型的不足之处在于假设积雪的密度与温度为均匀分布，忽略了雪层内部温度和密度的变化，因此在特定的环境下可能会产生误差，在后续的研究中可以将雪层分层模拟，考虑积雪内部的非均匀性，进而更准确地模拟积雪融化过程。

积雪消融主要是与大气间物质能量交换的过程，但由于试验地观测资料局限和过程的复杂性，净辐射通量、显热通量、潜热通量不是直接观测得到，而是采用了前人研究经验对所涉及能量进行估算，这在今后工作中应当增加必要的能量观测，进一步分析积雪与天气系统之间的变化状况。

本研究基于能量平衡原理建立单层融雪模型对融雪期雪深进行分析与模拟，得出以下结论：净辐射通量、显热通量、潜热通量是研究区融雪能量来源，融雪期间呈周期性的变化。其中净辐射通量是主要的能量来源，2017 年和2018 年积雪消融观测期，净辐射通量平均值分别为117.8 W·m-2和122.5 W·m-2，占到总能量的74.7%～67.4%，而湍流交换分别占到了18.8%～25.8%。能量闭合率随雪深的增大而变大，最大值为0.25～0.67，积雪覆盖会使能量闭合率偏低。使用单层融雪模型模拟的积雪深度变化，2017年和2018年纳什系数（Nash）分别达到0.97和0.83；
平均绝对误差（MAE）分别为0.36 cm和1.46 cm。该模型能较好地模拟融雪深度的变化。

猜你喜欢 雪深融雪潜热 一元复始万象更新优雅(2020年2期)2020-04-30自融雪沥青混合料路用性能及融雪能力试验研究工程与建设(2019年1期)2019-09-03Effect of moxibustion combined with acupoint application on enteral nutrition tolerance in patients with severe acute pancreatitisNursing Communications(2019年3期)2019-08-30道岔融雪设备的管理与维护铁道通信信号(2018年8期)2018-11-10一种基于频率与相位算法的雷达雪深探测方法成都信息工程大学学报(2018年3期)2018-08-29高原冬季雪深与重庆夏季降水的年际关系研究成都信息工程大学学报(2018年2期)2018-06-27工业革命时期蒸汽动力的应用与热力学理论的关系科学文化评论(2018年2期)2018-05-14铁路防灾雪深图像采集的设计和实现铁道通信信号(2016年5期)2016-06-01融雪毯儿童故事画报·发现号趣味百科(2016年1期)2016-02-25青藏高原东部夏季降水凝结潜热变化特征分析成都信息工程大学学报(2015年6期)2015-12-02