1979—2021年鄂陵湖和扎陵湖结冰日变化及主控因素分析

黄文峰 李瑞 李志军 张庾粟 杨文焕

摘要：湖冰物候过程影响着湖泊冰封期湖泊环境的初始、发展和结束条件以及湖-气作用等。为了探究青藏高原淡水湖泊湖冰物候变化趋势及影响因素，对FLake模型进行改进，并使用MODIS温度数据进行验证，以鄂陵湖和扎陵湖为例，重建其1979—2021年结冰日数据序列，利用模型对湖泊结冰日进行敏感性分析。结果表明：改进的FLake-SELF模型对两湖的结冰日计算精度提升了35%，重建序列显示鄂陵湖和扎陵湖结冰日分别以4.5 d/（10 a）、3.8 d/（10 a）的速率推迟；气温变暖、风速降低和水汽压增加是导致两湖结冰日推迟的最主要气象因素，气温、水汽压和风速的年际差异是造成两湖结冰日年际波动的关键气象因素；相同的气候变化条件下，湖泊深度的不同是造成两湖结冰日及其变化速率差异的原因；较深的湖泊，结冰日年际波动性越大、变化速率越大，对气候变化响应更为敏感。

关键词：结冰日；气候变化；湖泊模型；鄂陵湖；扎陵湖

中图分类号：TV11

文献标志码：A

文章编号：1001-6791（2023）01-0102-13

收稿日期：2022-06-28；

网络出版日期：2022-11-07

网络出版地址：https：∥kns.cnki.net/kcms/detail∥32.1309.P.20221107.1119.002.html

基金项目：国家重点研发计划资助项目（2019YFE0197600）；冻土工程国家重点实验室开放基金资助项目（SKLFSE201813）

作者简介：黄文峰（1985—），男，安徽濉溪人，副教授，博士，主要从事冰物理与冰区生态环境研究。E-mail：huangwenfeng@chd.edu.cn

联合国政府间气候变化专门委员会（IPCC）研究表明，在过去的1个世纪里，世界各地的气温出现不同程度的上升，并预测未来一段时间内气温升高的现象还将持续^［1］。青藏高原由于自身特性，对全球气候变化的响应更为显著，升温速率是全球均值的2倍^［2-3］。高原湖泊对气候变暖响应敏感，表现在气候变化导致湖泊开始冻结及封冻历时发生改变，因此湖冰被认为是气候变化的指示因子^［4-5］。湖冰作为冰冻圈的重要组成部分^［6］，年内变化不仅影响着区域能量和物质收支平衡，还具有重要的生态和社会服务功能，体现在多个方面：湖冰的存在影响湖泊与大气直接能量交换^［7］，湖冰较高的反射率降低了进入湖泊的太阳辐射^［8］，中止了风动力驱动的湖泊混合过程，隔绝了水体与大气的直接物质交换，复氧、大气沉降被极大削弱或中断，改变湖泊微生物的生态环境使微生物群落丰富度增加^［9］。在北欧、中国黑龙江等高纬度地区，坚硬厚实的冰盖也成为跨越湖泊的天然通道等。因此，了解湖冰对气候变化的响应对于预测气候变化对湖泊及其相关生态系统的影响至关重要。

湖冰过程在气候变暖的背景下正发生变化，但不同湖泊的封冻历时变化趋势并不相同。研究发现，纳木错湖封冻历时在1978—2017年以8.0 d/（10 a）的趋势缩短^［10］；咸水湖勒斜武担湖、邦达错湖封冻历时同期却出现了延长的现象^［11］；欧洲北部湖泊封冻历时的变化趋势要高于南部地区^［12］；芬蘭Segozero湖（1950—2009年）的开始冻结日以3.7 d/（10 a）的趋势推迟，同地区其他湖泊开始冻结日60 a推迟了2～12 d不等^［13］。许多早期湖冰的研究依赖于人类对湖冰长期直接观测。在环境恶劣、人烟稀少的青藏高原，大多数湖泊没有长期的观测数据。不少学者借助卫星遥感资料对青藏高原湖冰进行研究，但能够成功用于湖冰物候过程反演的主流遥感资料时间相对较短^［14-15］，难以准确揭露湖冰物候过程的长期变化趋势。因此使用依赖于气象数据的湖泊热力学模型可以对青藏高原地区湖泊湖冰物候进行长序列的模拟，如基于相似理论发展出来的相对简单的二层模型FLake（Fresh Lake Model）^［16-17］、基于热扩散发展出来的CLIMO（Canadian Lake Ice Model）模型^［18］和WRF（Weather Research and Forecasting Model）-lake模型^［19］等。这些模型的水温模拟效果较好（误差约2～4 ℃），但因湖泊开始冻结日期对水温变化较为敏感，可引起开始冻结日计算误差高达20 d^［20］。

本文基于对影响湖泊开始冻结物理过程的分析，对FLake模型进行改进，建立FLake-SELF模型，并利用20 a的遥感监测数据验证模型；基于气象观测与再分析数据重建鄂陵湖、扎陵湖结冰日序列（1979—2021年），分析两湖开始冻结日变化趋势，辨别和评价造成湖泊开始冻结变化的关键因素，探究不同深度的湖泊对气候变化响应差异性的原因。

1 研究区概况

鄂陵湖、扎陵湖是青藏高原大型淡水湖泊的典型代表（图1），位于青藏高原北部、黄河源头地区（34°46′03″N —35°05′21″N，97°01′27″E—97°54′21″E）。

鄂陵湖海拔为4 269 m，常年平均水深为18.9 m，最大水深为31 m，南北长约32.3 km，东西宽约31.6 km，形似葫芦状，南宽北窄，湖水面积约679 km²，是青藏高原第一大淡水湖泊，且水量正以0.08 Gt/a的速率上升^［21］，年均气温为-3.6 ℃，10月至次年4月的平均气温都在0 ℃以下，湖泊每年11月末12月初开始冻结，次年4月下旬湖冰开始消融。最大冰厚可达0.6 m以上，封冻历时长达半年以上^［22-23］。

扎陵湖位于鄂陵湖西侧、上游，海拔为4 294 m，平均水深为8.9 m，最大水深为13 m，位于湖心偏东北侧；该湖南北宽约21.6 km，东西长约35 km，形似大贝壳，面积约526 km²，冰期为每年11月中下旬到次年4月中旬，封冻历时与鄂陵湖相当。两湖为姊妹湖，最近距离仅为11 km，湖区气候条件及气候变化特征基本一致，但开始冻结、完全消融日期和最大冰厚却存在明显差异（图1）。

2 研究资料与方法

2.1 数据来源与处理

2.1.1 气象数据

研究需要的气温、风速、水汽压和短波辐射均来自中国区域高时空分辨率地面气象要素驱动数据集^［24］，该数据集以Princeton 再分析资料、GLDAS 资料、GEWEX-SRB 辐射资料以及TRMM降水资料为背景场，融合了中国气象局常规气象观测数据制作而成。数据覆盖范围为1979年1月1日至2015年12月31日，时间分辨率为3 h，空间分辨率为0.1°。该数据可从国家青藏高原科学数据中心（https：∥www.tpdc.ac.cn）获得。已有研究证明该数据在中国陆面过程和水文过程研究中表现良好^{［17，25］}。云量数据采用欧洲中期天气预报中心（ECMWF）ERA-Interim数据集（https：∥apps.ecmwf.int/datasets/）进行补充^［26］，时间跨度为1979—2015年，时间分辨率为6 h，空间分辨率为0.125°。

2016—2021年数据用ERA-5数据集补充，时间分辨率为1 h，空间分辨率为0.1°（气温、风速、短波辐射）和0.25°（水汽压、云量）。2015年和2016年气温变幅较大，衔接误差达到2.7 ℃，因此2016—2021年气温数据用距离研究区最近的玛多县气象站和中国区域高时空分辨率地面气象要素驱动数据集的气温序列进行粗略校准，即将2016—2019年ERA-5提供的鄂陵湖气温增加2.0 ℃、扎陵湖增加1.5 ℃，使数据的衔接较为平整。其他数据衔接误差较小，其影响可以忽略。

2.1.2 遥感数据

本文基于多种遥感产品提取鄂陵湖、扎陵湖湖冰物候信息，主要包括地表温度产品MOD11A1以及雪覆盖产品MOD10A1、MYD10A1经过处理之后获得的青藏高原湖冰范围数据集^［27］。该数据集利用MODIS归一化差雪指数（NDSI）和传统的雪图算法获取晴天条件下的湖泊日冰量和覆盖范围；利用湖泊表面条件的时空连续性，通过一系列算法对云覆盖条件下的湖冰范围进行重新确定，以降低云污染的影响^［28］。本文选用2002—2015年的每日湖冰数据，数据从国家青藏高原科学数据中心下载。

湖面温度 （Lake Surface Water Temperature，LSWT）选用MODIS地表温度数据产品MOD11A1，空间分辨率为1 km，时间分辨率12 h（北京时间10：30和22：30每日2次）。有研究表明MOD11A1产品有较高的精度，与试验地区地表温度相差1.49 K^［29］。选取2002—2021年时间段，共14 586幅，数据量较大，本文采用MTR（Modis Reprojection Tool）工具进行批量投影转换和拼接处理，用水温提取程序^［30］，对湖泊水面上所有的数据点进行平均计算，将白天均值和夜晚均值再次求平均，取用作每日湖面温度。

2.2 湖冰物候过程特征值

湖冰物候指湖泊冻结、消融的时间，常用以下关键物候指标来表征：开始冻结（Freeze-up start，FUS，湖冰首次出现）、完全冻结（Freeze-up end，FUE，湖冰完全覆盖湖面）、开始消融（Break-up start，BUS，湖冰开始融化、破裂）、完全消融（Break-up end，BUE，湖冰完全消失）。湖泊冰情的演变情况可以由上述4个物候节点反映。例如，封冻历时（Freeze duration，FD）是湖泊开始冻结至完全消融的一段时间。

在遥感数据出现或利用以前，湖泊开始冻结通常是指观测点视线范围内首次出现湖冰的日期，完全消融是指观测点视线范围内无冰的最后日期。对于小湖而言该类定义基本与实际结冰过程相符，但对于大、中型湖泊而言，该类定义与湖泊的冰物候存在较大差异。即便是遥感技术应用于湖冰过程监测后，因其重访期、空间分辨率和精度问题仍不能精确地解译湖泊开始冻结、完全消融的确切时间^［31］。

本文选择湖面温度降至低于冰点（淡水湖泊的冰点约为0 ℃）作为湖泊开始冻结日；相似地，定义湖面温度开始升至冰点的日期作为湖冰完全消融日。尽管该定义与实际开始冻结日和完全消融日有一定出入，但可直观刻画整个湖泊结冰与冰层消失的平均状态。通过与湖冰发展遥感数据的比对（图2，该数据集在湖泊封冻历时前后内首次识别出湖冰的日子为FUS、最后一次识别湖冰的日子为BUE，s_bia为偏差），发现鄂陵湖、扎陵湖湖面开始冻结和完全消融日湖面冰层覆盖率基本分别维持在5%±6%（均值±标准差，最小值不低于0）和17%±25%，这与遥感数据通常定义的开始冻结日、完全消融日（或开湖日）十分接近，說明了该定义用于指示湖冰物候过程的合理性。

2.3 模型的说明与验证

2.3.1 FLake模型及改进

通过选用湖面平均温度作为湖冰物候指标，可简化计算，将湖泊温度变化作为垂向一维问题。文章采用一维淡水湖泊模型FLake，采用相对简单的分层方案，在浅湖适用性较好^［16］。FLake将湖泊分为2层，上层湖水充分混合、温度趋于等温，称其为混合层（图3），与大气相接；下层为温跃层，采用自相似理论对其温度结构进行参数化，与湖底相接。自相似理论本质是温跃层中的量纲一温度剖面可以通过“量纲一深度的通用函数”以合理的精度进行参数化^［32］，即：

式中：t为时间；z为深度；θ（z，t）为深度z处的水温；θ_s（t）为水面温度（此时等于混合层水温）；θ_b（t）为湖底温度；h（t）为混合层深度；d为湖水深度；ψ_θ（ζ）为深度比值ζ=［z-h（t）］/［d-h（t）］的函数，满足边界条件ψ_θ（0）=0和ψ_θ（1） =1。该模型可计算雪层、湖冰层、水层和湖底沉积层的温度结构，并将自相似理论应用于所有介质层。

湖面温度受季节影响出现周期变化，春季、夏季湖面升温，秋季、冬季湖面温度下降，湖泊冻结就是水面温度降至0 ℃的过程。对于淡水湖泊来说，水在约4 ℃时密度最大。进入秋季，湖泊表面降温过程最为迅速，其降温过程可以以4 ℃为分界线分为3个阶段：高于4 ℃时湖水在风力和重力的作用下进行局部混合，4 ℃时湖水在自身重力和风力的作用下达到整体混合，低于4 ℃时水体仅在风动力强迫下局部混合或者无垂向混合。试算发现，FLake模型计算的结冰日一般晚于观测值14～20 d，存在较大误差。通过分析FLake模型4～0 ℃降温过程中的水温垂向曲线发现，可能是由于过于简单的分层方案使水温垂向曲线刻画比较粗略，湖泊散失一定热量时湖表降温幅度偏小，湖表降温速率偏低，使结冰日计算值晚于观测值。为了减小模拟误差，需要对驱动水温垂向结构和表温的物理过程模拟更加精细化，因此尝试将刻画4～0 ℃阶段更为详细的SELF模型融入FLake模型中^［33］，即：降温过程中，当FLake计算湖表温度将至4 ℃后，4～0 ℃阶段使用SELF模型计算，形成FLake-SELF联合模型，以期改善结冰日的模拟情况。

SELF模型将每日的降温过程分为混合阶段和降温阶段2个阶段（图4，虚线为上一时刻水温）。

（1） 混合阶段。湖水在风的作用下进行混合，使混合层内水温趋于一致，即风的动能输入用于搅动水体混合而改变水体势能，可表示为

式中：η为风能传递系数；P_w（t）为风能;E_P为水体势能;ΔE_P为水体混合之后势能的改变量；

θ_mix（t）为混合层水温。式（2）是关于h（t）单一变量的方程，由此计算混合层深度。

（2） 降温阶段。降温阶段仅作用于混合层，水温计算如下：

式中：Δθ_s（t）为水面与前次混合层水温的差值；E_c（t）为单位时间内水面能量交换；c_p为淡水比热容；ρ₀（t）为淡水密度。

2.3.2 模型设置

FLake、SELF模型驱动数据要求相似，均含有短波辐射、风速、气温、水汽压和云量这5个气象要素，2个模型均使用中国区域高时空分辨率地面气象要素驱动数据集和ERA再分析数据计算湖水表面温度。对FLake-SELF模型进行初始化，水深设置为平均水深（鄂陵湖18.9 m，扎陵湖8.9 m），初始日期设置为夏季1979年6月1日。由于缺少初始混合层温度和厚度，初始水温和混合层厚度利用重复迭代的方法获得，在模拟过程中用1979年的驱动数据反复迭代模拟10 a，使湖泊温度不再发生变化，以消除模型初始值的影响^［32］。鄂陵湖、扎陵湖水的消光系数（λ）根据经验公式估算^［34］，模型其他参数依赖内部的经验公式和数值数据估算，湖泊依赖性较小，不需要重新调整^［35］。计算消光系数的经验公式为

模型的气象驱动数据来自于中国区域高时空分辨率地面气象要素驱动数据集^［24］和ERA数据集产品^［26］，模拟结果用MODIS遥感数据进行验证，以均方根误差（E_RMS）、s_bia作为评估指标。

式中：m为模拟结果；

o为遥感结果；m為模拟结果平均值；o为遥感结果平均值；下标i为第i年（1979年为第1年），i=1，2，…，n，n为模拟数量。

3 结果与讨论

3.1 鄂陵湖和扎陵湖湖冰物候特征分析

鄂陵湖和扎陵湖2002—2015年不同遥感监测方案的开始冻结日和完全融化日序列如图5所示。扎陵湖开始冻结多始于11月中旬，次年4—5月完全融化，多年平均开始冻结日为319 d，完全消融日为470 d（366 d为次年1月1日），封冻历时为151 d，标准差为18 d；鄂陵湖多于11月底12月初开始冻结，次年5—6月完全消融，多年平均开始冻结日为334 d，完全消融日为486 d，封冻历时151 d，标准差为12 d。鄂陵湖开始冻结晚于扎陵湖，但是两者封冻历时没有明显差别。2002—2014年13个冰封期内两湖湖冰物候变化趋势不明显。不同的方法比较可以看出，由湖面温度获得湖泊冰封期变化浮动小于由湖冰范围得到的冰封期，这表明由温度数据得到的湖泊冰封期更加稳定。

3.2 FLake-SELF模型的验证

已有的遥感数据序列可用于检验FLake-SELF的模拟效果，图6给出了鄂陵湖、扎陵湖非冰期湖面温度模拟结果。模型计算的湖面温度与MODIS监测数据高度一致（相关系数R²=0.84），鄂陵湖s_bia=3.3 ℃、E_RMS=4.2 ℃，扎陵湖s_bia=3.6 ℃、E_RMS=4.3 ℃，两湖模拟效果类似。Crosman等^［36］研究发现MODIS的观测值要低于实际观测值，并且云存在导致MODIS观测温度有较大的负偏差。Wen等^［37］认为，FLake模型在湖泊封冻历时内忽略了太阳辐射透过冰层对湖水的加热效应，低估了冰封期内湖泊的热容量及水温，因此模拟的完全消融日期晚于观测值，且相差较大，鄂陵湖E_RMS=17.4 d，扎陵湖E_RMS=23.5 d;但是开始冻结日模拟结果与观测值相差较小，鄂陵湖E_RMS=13.4 d，扎陵湖E_RMS=9.2 d。

改进的FLake-SELF模型开始冻结日模拟效果较FLake模拟结果更接近遥感观测值（图7）。鄂陵湖FLake模拟开始冻结日最大误差出现在2004年，为36 d；新模型同年误差为23 d，减少了13 d。扎陵湖Flake模拟最大误差出现在2019年，为24 d；新模型为19 d，减少了5 d。扎陵湖的FLake开始冻结日E_RMS=14.2 d，鄂陵湖E_RMS=20.5 d；新模型扎陵湖开始冻结日模拟精度提升了35.2%，鄂陵湖提升了34.6%。可见，用新模型来重建鄂陵湖、扎陵湖开始冻结日期序列，可以更好地探究开始冻结日期变化规律。

3.3 鄂陵湖和扎陵湖开始冻结日变化趋势

利用改进的FLake-SELF模型重构两湖1979—2021年开始冻结日序列，如图8所示。多年平均而言，鄂陵湖开始冻结出现在340 d，扎陵湖为322 d，鄂陵湖开始冻结日晚于扎陵湖约18 d，与遥感观测结果一致（图5）。

基于Mann-Kendall非参数检验法的开始冻结日变化趋势分析表明，1979—2021年鄂陵湖、扎陵湖开始冻结日推迟变化趋势均达到99%的显著性水平。43 a中，鄂陵湖开始冻结日期已从11月下旬变化到12月中旬，扎陵湖从11月中旬变化到12月上旬，两湖开始冻结日期均明显推迟。鄂陵湖开始冻结日平均推迟4.5 d/（10 a），扎陵湖平均推迟3.8 d/（10 a），鄂陵湖开始冻结日延后速率略大于扎陵湖。两湖开始冻结日年际变化有明显的一致性（R²=0.82），最早开始冻结日均出现在1985年，可能是两地相似的气候条件造成的，鄂陵湖、扎陵湖在1985年年均气温分别为-5.77、-6.01 ℃，这是43 a内观测到的最低年均气温。

3.4 影响湖泊开始结冰的主要气象因素

以扎陵湖为例，1979—2021年每10 a年均气温增长0.62 ℃、水汽压增长17 Pa、风速减少0.2 m/s、短波辐射减少0.2 W/m²、云量增加0.01，各气象要素均发生不同程度变化。通过设置敏感性试验来识别引起结冰日变化的关键因素。将气象驱动数据进行组合开展5种虚拟试验，即假定1979—2021年只保留其中1个气象要素取实测序列，其他气象要素43 a间均用1980年数据替换而保持年过程不变。如短波辐射（R）组，其气温、水汽压、风速、云量1979—2021年均维持1980年的年变化过程不变，仅有短波辐射取43 a的实际观测数据。以此类推，定义气温（T_a）组、水汽压（e）组、风速（W）组和云量（C）组。驱动数据的背景值选择1980年，是因为该年各气象因素年均值与1979—1982年间的平均值相差最小，仅为7%，可以用于表征扎陵湖模拟初期的气候状态。

模拟结果表明（图9），气温、水汽压、风速、云量和短波辐射单独变化时可造成扎陵湖开始冻结日每10 a分别推迟2.1 d、推迟1.4 d、推迟0.9 d、推迟0.6 d和提前0.5 d。横向比较可以发现，气温变化对扎陵湖开始冻结日影响最大，其引发的开始冻结日变化与实际变化相关性最高（相关系数R²=0.81），其次为水汽压（R²=0.79）和风速（R²=0.52）。另一方面，重建序列的方差（D）和变差系数（C_V）分别为6.3 d和1.9%，气温组试验结果的方差和变异系数分别为4.3 d和1.3%，与重建序列最为接近，表明其对开始冻结日年际差异性贡献最大;其次为风速（D=2.8 d，C_V=0.9%）和水汽压（D=2.3 d，C_V=0.7%）。可見，气温、风速、水汽压是扎陵湖开始冻结日年际差异性和推迟的关键气象因素。

3.5 不同深度湖泊开始冻结日对气候变化响应的差异

对比发现，扎陵湖平均气温、水汽压、风速、短波辐射和云量分别比鄂陵湖低0.5 ℃、低5 Pa、低0.026 m/s、高0.83 W/m²和低0.001，相差极小，两湖气象条件基本一致。扎陵湖气象要素平均每10 a变化6.51%，鄂陵湖平均每10 a变化6.56%，两湖气象要素变化趋势基本相同。为了消除气候差异对两湖开始冻结日的影响，将两湖气候条件互换模拟，结果表明互换气候条件之后扎陵湖开始冻结日变化仅为1.28 d，因此，可以推测气象要素不是造成鄂陵湖、扎陵湖结冰日期表现不同的主要原因。为了探究两湖结冰日不同的原因，考虑湖泊形态特点，设计了5组不同湖深（5、10、20、30、40 m）的模型模拟，以分析湖泊深度对开始冻结日在相同气候变化下的影响。

结果表明，相同气候变化情势下，随湖泊深度的增加开始冻结日逐渐推迟（图10）。平均湖深为5、10、20、30、40 m时近40 a平均开始冻结日期分别出现在311、324、339、347、355 d。湖深20 m开始冻结日比10 m时晚出现14.8 d，这与卫星观测的鄂陵湖（18.9 m）比扎陵湖（8.9 m）平均结冰日（2002—2021年）晚14 d的结果相近。因此可推断主要是两湖深度的不同造成了两湖开始冻结日的差异。同时还发现，湖深若为40 m时，2017年则可能出现不结冰的现象。

不同深度湖泊过去43 a推迟日数也存在明显差异。湖深5、10、20、30、40 m时近40 a开始冻结日期推迟速率分别为每10 a推迟3.5、3.7、5.2、6.0、7.0 d。随着湖深的增加结冰日推迟天数不断增大。43 a湖深10 m与20 m推迟天数相差1.5 d，这也与鄂陵湖、扎陵湖开始冻结日重建结果（相差0.7 d）相近。

不同深度的湖泊开始冻结日年际波动幅度不同，5、10、20、30、40 m深的湖泊在过去43 a的時间序列里开始冻结日期标准差分别为6.3、6.6、7.9、9.0、10.2 d，随着深度的增加开始冻结日年际波动幅度也在增加。

3.6 湖泊开始冻结日敏感性分析

IPCC第6次报告指出，当前地球气候变暖是前所未有的，预计到2100年根据不同温室气体排放气温将升温1.3～5 ℃^［1］。根据中国1960—1990年30 a的短波辐射量研究，发现全国大部分地区太阳辐射量降低了15%左右^［38］。研究指出西南地区1969—2009年风速以每10 a减小0.24 m/s的速率显著降低^［39］。因此，对气温、短波辐射、风速3个气象因素进行敏感性分析。

以扎陵湖气候要素为背景值，各气象要素分3次进行敏感性分析（气温分别增加1.5、3.0、4.5 ℃；短波辐射分别减少5%、10%、15%；风速分别减少0.4、0.8、1.2 m/s）（图11）。结果表明：气温组开始冻结日变化最为明显，增温4.5 ℃时平均开始冻结日从初始条件的322 d推迟至332 d，推迟了10 d；其次为短波辐射组，辐射量减少15%时平均开始冻结日期较初始条件提前了5 d；风速组开始冻结日变化不明显，风速减少1.2 m/s时平均开始冻结日期较初始条件仅推迟了2.3 d。气温的增加和风速的减小均削弱了湖泊与大气的热量交换过程，造成了开始冻结日推迟的现象；短波辐射是湖水的主要能量来源，因此短波辐射的减小造成了湖泊蓄能降低，湖面更容易达到冰点，开始冻结日提前。

既然气温是造成开始冻结日变化最重要的因素，则根据IPCC气温预测结果，可对两湖未来开始冻结日的变化进行大致预测。以2015年气温为基础，将5 ℃升温幅度平均地分配在2016—2100年，构建出未来气温变化。其他气象要素保持不变，均采用与近5 a均值最接近的年份数据，即水汽压采用2017年数据、风速和云量采用2020年数据、短波辐射采用2021年数据。使用该数据计算的2016—2021年开始冻结日与上面构建的开始冻结日序列相比，误差仅为4.3 d，说明这种处理的气象数据具有一定的合理性。预测结果表明，两湖的开始冻结日仍维持推迟现象，鄂陵湖推迟速率为2.9 d/（10 a）、扎陵湖为1.6 d/（10 a）；到2100年，鄂陵湖开始冻结日为373 d，扎陵湖开始冻结日为349 d，相比于近5 a，开始冻结日将分别推迟23 d和12 d。

4 结论

本文将计算预冻结期的SELF模型融入FLake模型，构建了计算湖泊开始冻结日的FLake-SELF模型，对研究区开始冻结日进行了模拟和重建，研究了青藏高原鄂陵湖、扎陵湖1979—2021年开始结冰日的演变特征，并通过敏感性试验识别、分析了湖泊开始结冰的主控因素，主要结论如下：

（1） 将SELF模型融入FLake模型来模拟淡水湖结冰日，并以MODIS湖面温度数据作为验证。新模型对鄂陵湖、扎陵湖开始冻结日模拟精度均有提升。利用该模型重建了两湖1979年以来的开始冻结日序列，发现两湖开始冻结日均呈现出不同程度的推迟（鄂陵湖推迟速率为4.5 d/（10 a）、扎陵湖为3.8 d/（10 a））。

（2） 模型试验揭示，1979年以来，气温、风速、水汽压的年际波动是鄂陵湖、扎陵湖开始冻结日年际波动的关键因素，气温升高、水汽压升高、风速降低是造成湖泊开始冻结日推迟的关键原因。

（3） 两湖气象、气候条件及变化趋势基本相同，其开始冻结日及其推迟速率的差异主要是由两湖深度不同造成的。相同的气候条件下深（淡水）湖开始冻结日对气候变化的响应更为敏感。

本文发展的淡水湖泊开始冻结日计算模型为监测困难、条件恶劣地区的湖冰物候变化的研究提供了一种模拟方法，可用于长序列开始冻结日重建、湖冰物候变化机制分析和未来变化预测，但模型未能模拟湖冰完全消融日，未来仍需进一步完善模型的结构与参数化方案。

参考文献：

［1］IPCC. Climate change 2022：summary for policymakers ［M/OL］. ［2022-05-24］. https：∥www.ipcc.ch/report/ar6/wg3/downloads/report/IPCC_AR6_WGIII_SPM.pdf

［2］YANG K，WU H，QIN J，et al. Recent climate changes over the Tibetan Plateau and their impacts on energy and water cycle：a review［J］. Global and Planetary Change，2014，112：79-91.

［3］季芳，范林峰，匡星星，等. 青藏高原多年凍土退化对蒸散发的影响［J］. 水科学进展，2022，33（3）：390-400. （JI F，FAN L F，KUANG X X，et al. Role of permafrost degradation on evapotranspiration on the Qinghai-Tibet plateau［J］. Advances in Water Science，2022，33（3）：390-400.（in Chinese））

［4］GRANT L，VANDERKELEN I，GUDMUNDSSON L，et al. Attribution of global lake systems change to anthropogenic forcing［J］. Nature Geoscience，2021，14（11）：849-854.

［5］刘吉平，金晶，李晟铭，等. 1995—2015年全球不同气候带典型冰川补给型湖泊动态变化［J］. 水科学进展，2021，32（1）：1-9. （LIU J P，JIN J，LI S M，et al. Dynamic changes of typical glacier-fed lakes in different climatic zones around the world from 1995 to 2015［J］. Advances in Water Science，2021，32（1）：1-9.（in Chinese））

［6］丁永建，张世强，吴锦奎，等. 中国冰冻圈水文过程变化研究新进展［J］. 水科学进展，2020，31（5）：690-702. （DING Y J，ZHANG S Q，WU J K，et al. Recent progress on studies on cryospheric hydrological processes changes in China［J］. Advances in Water Science，2020，31（5）：690-702.（in Chinese））

［7］HUANG W F，ZHAO W，ZHANG C，et al. Sunlight penetration dominates the thermal regime and energetics of a shallow ice-covered lake in arid climate［J］. The Cryosphere，2022，16（5）：1793-1806.

［8］HUANG W F，ZHANG Z，LI Z J，et al. Under-ice dissolved oxygen and metabolism dynamics in a shallow lake：the critical role of ice and snow［J］. Water Resources Research，2021，57（5）：e2020WR027990.

［9］GARNAUD C，MACKAY M，FORTIN V. A one-dimensional lake model in ECCC′s land surface prediction system［J］. Journal of Advances in Modeling Earth Systems，2022，14（2）：e2021MS002861.

［10］GUO L N，ZHENG H X，WU Y H，et al. Responses of lake ice phenology to climate change at Tibetan Plateau［J］. IEEE Journal of Selected Topics in Applied Earth Observations and Remote Sensing，2020，13：3856-3861.

［11］XIONG C，LEI Y H，QIU Y B. Contrasting Lake ice phenology changes in the Qinghai-Tibet plateau revealed by remote sensing［J］. IEEE Geoscience and Remote Sensing Letters，2021，18（12）：2132-2136.

［12］NEWTON A M W，MULLAN D J. Climate change and Northern Hemisphere Lake and river ice phenology from 1931—2005［J］. The Cryosphere，2021，15（5）：2211-2234.

［13］EFREMOVA T，PALSHIN N，ZDOROVENNOV R. Long-term characteristics of ice phenology in Karelian lakes［J］. Estonian Journal of Earth Sciences，2013，62（1）：33-41.

［14］TUTTLE S E，ROOF S R，RETELLE M J，et al. Evaluation of satellite-derived estimates of lake ice cover timing on linnévatnet，kapp linné，svalbard using in situ data［J］. Remote Sensing，2022，14（6）：1311.

［15］LANG J H，MA Y M，LI Z G，et al. The impact of climate warming on lake surface heat exchange and ice phenology of different types of lakes on the Tibetan Plateau［J］. Water，2021，13（5）：634.

［16］MIRONOV D V. Parameterization of lakes in numerical weather prediction. description of a lake model［R］. Offenbach：Deutscher Wetterdienst，2008.

［17］宋興宇，文莉娟，李茂善，等. 不同湖泊模式对青藏高原典型湖泊适用性对比研究［J］. 高原气象，2020，39（2）：213-225. （SONG X Y，WEN L J，LI M S，et al. Comparative study on applicability of different lake models to typical lakes in Qinghai-Tibetan Plateau［J］. Plateau Meteorology，2020，39（2）：213-225.（in Chinese））

［18］QIU B，HUANG A N，SHI X L，et al. Implementation and evaluation of an improved lake scheme in Beijing climate center atmosphere-vegetation interaction model［J］. Journal of Geophysical Research：Atmospheres，2020，125（9）：e2019JD031272.

［19］WU Y，HUANG A N，LAZHU，et al. Improvements of the coupled WRF-Lake model over Lake Nam Co，Central Tibetan Plateau［J］. Climate Dynamics，2020，55（9/10）：2703-2724.

［20］ROBINSON A L，ARIANO S S，BROWN L C. The influence of snow and ice albedo towards improved lake ice simulations［J］. Hydrology，2021，8（1）：11.

［21］龙笛，李雪莹，李兴东，等. 遥感反演2000—2020年青藏高原水储量变化及其驱动机制［J］. 水科学进展，2022，33（3）：375-389. （LONG D，LI X Y，LI X D，et al. Remote sensing retrieval of water storage changes and underlying climatic mechanisms over the Tibetan Plateau during 2000—2020［J］. Advances in Water Science，2022，33（3）：375-389.（in Chinese））

［22］WEN L J，WANG C，LI Z G，et al. Thermal responses of the largest freshwater lake in the Tibetan Plateau and its nearby saline lake to climate change［J］. Remote Sensing，2022，14（8）：1774.

［23］WANG M X，WEN L J，LI Z G，et al. Mechanisms and effects of under-ice warming water in Ngoring Lake of Qinghai-Tibet Plateau［J］. The Cryosphere，2022，16（9）：3635-3648.

［24］YANG K，HE J. China meteorological forcing dataset （1979—2015）［DS］. Lhasa：National Tibetan Plateau Data Center，2016. ［DOI：10.3972/westdc.002.2014.db. CSTR：18406.11.westdc.002.2014.db］

［25］CHEN Y Y，YANG K，HE J，et al. Improving land surface temperature modeling for dry land of China［J］. Journal of Geophysical Research Atmospheres，2011，116（D20）：D20104.

［26］DEE D P，UPPALA S M，SIMMONS A J，et al. The ERA-Interim reanalysis：configuration and performance of the data assimilation system［J］. Quarterly Journal of the Royal Meteorological Society，2011，137（656）：553-597.

［27］QIU Y B. The Tibet Plateau river lake ice range / coverage data set V1.0（2002—2018）［DS］. Lhasa：National Tibetan Plateau Data Center，2019. ［DOI：10.11922/sciencedb.744］

［28］QIU Y B，XIE P F，LEPPRANTA M，et al. MODIS-based daily lake ice extent and coverage dataset for Tibetan Plateau［J］. Big Earth Data，2019，3（2）：170-185.

［29］HULLEY G C，HOOK S J，BALDRIDGE A M. Investigating the effects of soil moisture on thermal infrared land surface temperature and emissivity using satellite retrievals and laboratory measurements［J］. Remote Sensing of Environment，2010，114（7）：1480-1493.

［30］王敏，梁雨华，马建行，等. 基于MODIS数据的2015年夏季全球大型湖泊表面温度格局初析［J］. 湿地科学，2017，15（4）：571-581. （WANG M，LIANG Y H，MA J H，et al. Preliminary study on pattern of surface temperature of large lakes in the world in summer in 2015 based on MODIS LST product［J］. Wetland Science，2017，15（4）：571-581.（in Chinese））

［31］REED B，BUDDE M，SPENCER P，et al. Integration of MODIS-derived metrics to assess interannual variability in snowpack，lake ice，and NDVI in Southwest Alaska［J］. Remote Sensing of Environment，2009，113（7）：1443-1452.

［32］MIRONOV D，HEISE E，KOURZENEVA E，et al. Implementation of the lake parameterisation scheme FLake into the numerical weather prediction model COSMO［J］. Boreal Environment Research，2009，15（2）：77-95.

［33］TOFFOLON M，CORTESE L，BOUFFARD D. SELF v1.0：a minimal physical model for predicting time of freeze-up in lakes［J］. Geoscientific Model Development，2021，14（12）：7527-7543.

［34］HKANSON L. Models to predict Secchi depth in small glacial lakes［J］. Aquatic Sciences，1995，57（1）：31-53.

［35］BALSAMO G，DUTRA E，STEPANENKO V，et al. Deriving an effective lake depth from satellite lake surface temperature data：a feasibility study with MODIS data［J］. Boreal Environment Research，2010，15：178-190.

［36］CROSMAN E T，HOREL J D. MODIS-derived surface temperature of the Great Salt Lake［J］. Remote Sensing of Environment，2009，113（1）：73-81.

［37］WEN L J，LYU S H，KIRILLIN G，et al. Air-lake boundary layer and performance of a simple lake parameterization scheme over the Tibetan Highlands［J］. Tellus A：Dynamic Meteorology and Oceanography，2016，68（1）：31091.

［38］虞海燕，劉树华，赵娜，等. 我国近59年日照时数变化特征及其与温度、风速、降水的关系［J］. 气候与环境研究，2011，16（3）：389-398. （YU H Y，LIU S H，ZHAO N，et al. Variation characteristics of the sunshine duration and its relationships with temperature，wind speed，and precipitation over recent 59 years in China［J］. Climatic and Environmental Research，2011，16（3）：389-398.（in Chinese））

［39］张志斌，杨莹，张小平，等. 我国西南地区风速变化及其影响因素［J］. 生态学报，2014，34（2）：471-481. （ZHANG Z B，YANG Y，ZHANG X P，et al. Wind speed changes and its influencing factors in Southwestern China［J］. Acta Ecologica Sinica，2014，34（2）：471-481.（in Chinese））

Analysis on the variation of ice-on date of Lakes Ngoring and Gyaring from 1979—2021 and its influencing factors

The study is financially supported by the National Key R&D Program of China （No.2019YFE0197600） and the Open Fund of State Key Laboratory of Frozen Soil Engineering（No.SKLFSE 201813）.

HUANG Wenfeng^1，2，LI Rui^1，2，LI Zhijun³，ZHANG Yusu^1，2，YANG Wenhuan⁴

（1. Key Laboratory of Subsurface Hydrology and Ecological Effect in Arid Region （Ministry of Education），Chang′an University，Xi′an 710054，China;2. School of water and Environment，Chang′an University，Xi′an 710054，China;3. State Key Laboratory of Coastal and Offshore Engineering，Dalian University of Technology，Dalian 116024，China;4. School of Energy and Environment，Inner Mongolia University of Science & Technology，Baotou 014010，China）

Abstract：Lake ice phenology affects the lake environment and lake-air interaction through the ice-covered period. In order to investigate the trends and influencing factors of lake ice phenology in freshwater lakes on the Tibetan Plateau，the FLake model was employed after improvement. The model validation used MODIS surface temperature data based reconstruction of the ice-on date time series in 1979—2021 for Lakes Ngoring and Gyaring. Sensitivity experiments were conducted using the improved model. The results showed that the present FLake-SELF model improved the accuracy of the ice-on date calculation by 35%，and the reconstructed sequences showed that the ice-on dates of Lakes Ngoring and Gyaring have been delayed by 4.5 d/（10 a） and 3.8 d/（10 a），respectively. Increasing air temperature，decreasing wind speed，and increasing vapor pressure are the most important meteorological factors leading to the delay of ice-on date，and interannual variability in these factors were the most important reasons for the interannual fluctuations of ice-on dates in the two lakes. In the same climate conditions，the different lake depth was the reason for the differences in ice-on date and their delay rate between the two lakes. The deeper the lake was，the greater was the interannual fluctuation in ice-on date，the greater was the rate of change of ice-on date，and the more sensitive was the ice-on date to climate change.

Key words：ice-on date;climate change;lake model;Lake Ngoring;Lake Gyaring