基于实测与模拟的青海湖冰厚时空变化特征*

曹 娟，姚晓军**，靳惠安,2，张调风，高永鹏，张大弘，赵全宁

(1：西北师范大学地理与环境科学学院，兰州 730070)(2：甘肃林业职业技术学院，天水 741020)(3：青海省气候中心，西宁 810000)(4：云南大学国际河流与生态安全研究院，昆明 650091)(5：青海省气象科学研究所，西宁 810000)

政府间气候变化专业委员会(IPCC)第五次评估报告(AR5)指出，过去130多年全球平均气温升高了0.85℃，主要升温时间段为1983-2012年的30年间[1]. 作为大气圈、冰冻圈和水圈的联系纽带，湖泊通过地表水的汇聚和蒸发以及地下水汇/渗流等过程参与地表水循环，并对气候变化响应敏感[2-3]. 青藏高原地区98个气象站点观测数据表明，1982-2012年平均气温升幅高达1.9℃，是全球平均升温幅度的2倍[4]. 青藏高原的自然环境和生态系统在全球占有特殊地位，该地区也被视作气候变化敏感区[5-6]，如青藏高原部分湖泊湖冰物候特征近期表现为开始冻结和完全冻结日期推迟、封冻期减少趋势[7-10]. 青海湖是维系青藏高原东北部生态安全的重要水体，在保护脆弱的高原湖泊湿地资源、阻挡西部荒漠化向东蔓延和保护鸟类及裸鲤等生物资源等方面发挥着不可替代作用. 随着西部大开发和丝绸之路经济带的建设，尤其是兰新高速铁路的开通，青海湖已成为旅游热点景区，而在冬季时有车辆掉入湖中事故发生. 因此，开展青海湖湖冰厚度研究不仅对于认识气候变暖背景下的湖冰响应规律具有重要的理论价值，而且对于制定科学合理的冬季青海湖管理措施具有现实意义.

作为湖泊冻结期的重要物理参数之一，湖冰厚度不仅反映了水-气界面能量交换强度和物质迁移过程，而且具有重要的生态价值和经济价值[11]，较其时间属性(如封冻期和消融期)难以监测[12]. 传统的冰厚测量主要采用接触式方法，如钻孔取冰、电热丝融冰以及压力传感器测厚法等，此类方法虽测量精度较高，但存在效率低、数据量少、损坏冰层等缺点[13]. 随着遥感技术快速发展和测量仪器设备的改进，一些非接触式冰厚测量方法开始出现，如基于冰和水的电导率差别原理采用电磁感应仪测量冰厚，在水下安装仰视声呐通过发射声脉冲根据回波时延差确定冰厚，该类方法仅能用于定点测量，难以反映湖面冰厚整体分布情况[14-15]；卫星数据(如CryoSat-2卫星搭载的Ku波段SIRAL)虽可用于监测大范围湖泊冰厚，但对薄冰反演效果较差且难以剔除雪对反照率的影响[16]. 从热力学角度来看，湖冰是开阔水域在大气和水的相互作用下通过内能和热能转化达到临界状态的产物，因此其厚度亦可借助热力学模型来模拟，如度日法湖冰生长模型[17]、CLIMo模型[18]、Zubov模型[19]、Mylake模型[20]等. 目前国际上对于湖冰厚度的研究主要集中在长时间序列模拟结果分析及与实测数据的对比，即模拟精度验证，或单点测量短时间序列的冰厚特征分析. 如Dibike等[20]基于NARR数据采用Mylake模型模拟了北美地区(40°～75°N)的湖泊冰厚，结果表明最大冰厚普遍分布在40 cm左右且随着湖深的增加冰厚呈现减小趋势，其中位于加拿大东西部及较高海拔地区的湖泊冰厚更大. Duguay等[18]以Barrow(北极)、Poker Flat Research Range(亚北极)、Churchill(高北森林)3个代表性地区的湖泊为研究对象，基于研究区附近的气象站点数据采用CLIMO模型进行冰厚模拟，发现随雪深变化模拟的最大冰层厚度介于159～227 cm，整个冬季平均绝对误差为2 cm. 就青海湖而言，目前系统性的冰厚观测仍很少，仅在湖区东南侧的下设水文站有冰厚观测记录，因人工观测点的单一性、湖冰形成的时空异质性，以及人工观测频次较低和观测位置不定，对于青海湖湖冰厚度时空分布及演化过程仍知之甚少. 本文基于实测钻孔冰厚数据、飞航式测冰雷达数据和较高时空分辨率的Terra MODIS遥感影像，采用度日法湖冰生长模型模拟青海湖冰厚变化，开展冰厚精度评价，并分析青海湖2000-2019年冰厚时空变化规律，从而为深入研究青海湖湖冰生消机理提供参考.

1 研究区概况

青海湖地处青藏高原东北部(36°32′～37°15′N，99°36′～100°46′E)，跨青海省刚察、海晏和共和3县，东西长约109 km，南北宽约40 km(图1). 根据青海省第一次水利普查公布的2013年测量数据，当水位在3193.5 m时，青海湖水域面积为4294 km2，平均水深18.3 m，最大水深26.6 m，蓄水量785.2×108m3，湖水呈弱碱性，pH值为9.23，含盐量为14.13 g/L[21]，是中国最大的内陆湖和咸水湖. 青海湖是新构造运动导致周围山地上升隆起，外泄通道堵塞形成的构造断陷湖，湖中现有海心山和三块石两个湖心岛，湖区东面自北向南依次分布着尕海、沙岛湖、海晏湾和耳海4个子湖. 青海湖入湖径流较大的河流主要有布哈河、沙柳河以及哈尔盖河等，其中布哈河径流量最大，约占入湖径流总量的60%[22]. 湖区处于我国东部季风区、西北干旱区和西南高寒区的交汇地带，属典型的高原半干旱高寒气候[23]，年平均气温介于-1.0～1.5℃，并由东南向西北逐渐降低，年降水量为350～400 mm，其中夏季降水达60%以上，年蒸发量可达1300～2000 mm. 因湖水中含有无机盐类，青海湖冻结温度比0℃稍低，通常于每年12月中旬至1月初开始冻结，12月底或次年1月初完全被冰层覆盖，3月下旬冰层逐渐消融，至4月初或中旬完全消融. 1958-1983年沙陀寺水文站观测结果表明，青海湖冰厚一般为0.5 m左右，最厚达0.7 m，且由岸边向湖心方向逐渐变薄[24]. 基于遥感数据的监测结果表明2000-2016年青海湖封冻期和完全封冻期平均为(88±20)和(77±20)d,湖冰存在期和湖冰消融期持续时间平均为(108±18)和(10±6)d[10].

图1 青海湖及湖冰厚度测区(点)分布

2 数据与方法

2.1 湖冰钻孔数据与飞航式雷达测厚数据

2019年1月13日-3月24日期间，在青海湖开展湖冰钻孔测厚和飞航式测冰雷达测厚试验，测厚样点分布及测飞区域如图1所示. 出于工作人员安全考虑和受无人机电池电量限制，测厚样点及测飞区域绕湖区周边布设. 为保障钻孔数据与飞航式测冰雷达测厚数据时空一致性，两项试验同时同地开展，其中钻孔取冰厚度采用L型数字式量冰尺(精度为0.01 mm)测量(图2a)；飞航式测冰雷达测厚采用IGPR-30雷达系统，以无人机为载具平台，并配备激光测距、GPS差分定位和高清摄像头等模块，设备雷达中心频率为400 MHz，采样时间间隔为2 ps，探测精度可达mm级，冰层探测厚度>6 m(图2b). 该设备的冰厚测量原理是雷达向下发射高频电磁波，当到达空气-冰界面和冰-水界面时发生反射，反射的电磁波信号被天线接收，假设雷达波发射至各界面和返回被天线接收所有时间相同，则冰厚可根据双程走时、冰层介电常数及电磁波在空气中的传播速度计算得到. 2019年1月13日-14日，在湖区西北侧及鸟岛附近各布设一个200 m×200 m观测场，并按50 m等间距在冰面布设标记点，首先采用飞航式测冰雷达对各标记点冰厚进行测量，然后钻孔取冰测量冰厚并用于校正飞航式测冰雷达模拟冰厚的介电常数. 2019年1月30日-31日，在仙女湾、哈尔盖河入湖口、尕海以及湖区东部进行钻孔取冰测厚(5个测点)和雷达测厚(6个测区)试验. 2019年3月1日-3日，在泉吉、鸟岛、布哈河口以及湖区西南部泉湾、峡湾、黑马河、三江源祭海台、二郎剑等地共测点10个，雷达测厚区共8个. 2019年3月22日-24日，在二郎剑、黑马河、布哈河口、鸟岛等地共测飞7个区域. 为更好地了解不同方向湖冰厚度分布情况，飞航式测冰雷达测量时随机采用多种航线(图2c～f)，并在每个测区各测飞2次冰厚. 以上湖冰厚度实测数据用于拟合度日法湖冰生长模型参数及模拟结果误差分析.

图2 青海湖冰厚测量照片与无人机航飞路线

图3 ERA5气候数据与观测值的温度对比

2.2 ERA5气候再分析资料

采用2000-2019年ERA5 Climate Reanalysis数据集中青海湖地区距地面2 m处气温格点数据作为研究区气温背景资料，参考青海湖历年开始冻结和湖冰完全消融日期[10]，截取该数据集中的每年10月1日至翌年4月30日数据. ERA5 Climate Reanalysis数据集是欧洲中期天气预报中心(ECMWF)于2017年发布的全球气候再分析资料第五代产品，提供了对全球近期气候的综合数字描述，时间分辨率为12 h，空间分辨率为0.125°×0.125°，该数据从欧洲中期天气预报中心网站(https://www.ecmwf.int)下载获得. 因ERA5 Climate Reanalysis数据集为NC格式且文件较多，利用Python语言编写脚本程序，将文件格式批量转换为GRID格式并计算逐格点日均温数据.

距青海湖最近的国家基准站仅有刚察气象站(图1)，通过将2000-2019年每年10月1日至翌年4月30日站点实测气温数据与该站点所在ERA5气候再分析数据格点对应分析，发现二者趋势较为一致(图3). 分析结果表明实测数据与ERA5气候数据相关系数为0.88(P<0.05)，均方根误差(RMSE)为5.17℃，两组数据间具有较好的相关性且相对误差较小. 相较于用该气象站点数据来表征青海湖整个湖面气温数据，ERA5气候再分析数据分辨率较高，能够更客观地反映青海湖湖面气温变化.

2.3 光学卫星遥感影像

为了解青海湖湖面冻结状况及湖冰空间分布情况，从美国国家航空航天局陆面分布式数据中心(NASALPDA)网站(https://data.giss.nasa.gov)共下载2888景MODIS MOD09GQ数据产品，该数据时间分辨率为1 d，空间分辨率为250 m. 此外，从美国地质勘探局(USGS)网站(http://glovis.usgs.gov)共下载19景2001-2019年青海湖封冻期且云量覆盖较小的Landsat TM/ETM+/OLI遥感影像，用于人工目视解译提取青海湖湖冰范围，以评估基于MODIS MOD09GQ数据提取的湖冰面积占比精度. 鉴于MODIS MOD09GQ数据量较大且需要进行影像拼接，首先利用MRT工具对其进行批量拼接和几何校正，然后通过编写的Python脚本程序提取湖冰并计算其面积占比.

2.4 度日法湖冰生长模型

湖冰可被视作水-气界面的绝缘膜，湖冰生消改变大气与湖水之间的热量、动量和物质交换[25]，如湖冰在生长过程中由于热能损失可改变湖面周围空气的热量. 早在19世纪，Stefan[26]就认为冰的生长过程是一种特殊的热传导过程，冰底结冰释放的热量在线性温度梯度情况下通过冰传导出去. Leppäranta[25]在其研究中指出，Stefan提出的理想状态下度日法湖冰生长模型的应用前提是不考虑热惯性和内部热源，冰层底部水中没有热通量且冰层顶部的温度已知，在模拟过程中只考虑温度变化对冰厚的影响，则冰体每日的生长速率为：

ρi·L·dH/dt=ki(Tf-T0)/H

(1)

式中，H为冰体厚度；ρi为冰体密度；L为冰体融化潜热；ki为冰体导热系数；T0为冰体上表面温度；Tf为冰体下表面温度. 若设t=0时H=H0为初始条件，则式(1)的解析式为：

(2)

(3)

(4)

若时间间隔t为1 d，S则被定义为水体处于冻结期的负积温(℃·d)，冰厚H计算可简化为：

(5)

3 结果与分析

3.1 基于实测数据的青海湖冰厚时空变化特征

表1为2019年1-3月青海湖实测冰厚数据，测区编号如图1所示. 由表1可知，在观测时段内，各采样点的冰厚总体呈增长趋势，但各月份冰厚增长速率不同. 1月30日测量数据显示湖区平均冰厚为40.2 cm，3月2日平均冰厚为43.5 cm，3月22日为51.4 cm，整个3月平均增长速率为0.30 cm/d，较2月份(0.12 cm/d)增长迅速. 但结果受该时间段内所测区域冰厚影响较大(可能所测地区均为冰厚较大或较小区域). 空间尺度上，实测冰厚大体呈现北厚南薄、东厚西薄特点. 就各测区平均增长速率而言，测区1(2019-01-30)、3(2019-01-31)、7(2019-03-01)仅有同一时间段内采样记录(未计算平均增长速率)，平均冰厚分别为42.7、46.4、49.2 cm. 有两段时间采样记录的测区有9个，可分为3个主要变化区：(1)冰厚快速增长区包括测区6、13、9，冰厚平均增长速率均大于0.30 cm/d，分别为0.46、0.34、0.32 cm/d. (2)冰厚平缓增长区平均增长速率均分布在0.25 cm/d左右，包括测区11(0.28 cm/d)、10(0.25 cm/d)、8(0.21 cm/d). (3)冰厚缓慢增长区包括测区12、14，冰厚平均增长速率仅分别为0.15、0.11 cm/d. 测区2、4、5、15将3个采样时间均包括在内，其对应平均冰厚增长速率依次减小，分别为0.30、0.24、0.12、0.08 cm/d.

表1 2019年1-3月青海湖实测冰厚数据

“航线长度”表示飞航式测冰雷达测飞长度，“-”表示钻孔测点无长度.

3.2 2019年青海湖冰厚模拟结果及精度评估

图4为2018-2019年冬季青海湖逐日气温与冰厚模拟值变化. 由该图可知，2018年11月8日青海湖开始冻结，当日湖面平均气温为-6.44℃，2019年3月29日冰厚达到最大值时气温为-1.21℃. 连续负值气温持续天数为141 d，累计负积温-1537.64℃，最大冰厚为51.05 cm. 将湖区相同位置的冰厚实测值与模拟值对应分析，两者之间的相关系数R2为0.92(P<0.05). 在模拟期内，青海湖冰厚平均增长速率为0.34 cm/d，负积温绝对值每平均累加10℃冰厚增加0.25 cm，其中2018年11月8日-11月21日及2018年12月1日-12月11日的冰厚日平均增长速率分别为0.71、0.59 cm/d，均大于0.50 cm/d；2019年2月21日-3月29日的平均增长速率为0.14 cm/d；其余时间段内湖冰厚度日平均增长速率介于0.20～0.50 cm/d之间.

图4 2018-2019年冬季青海湖气温变化与冰厚模拟值

为检验度日法湖冰生长模型模拟的青海湖冰厚精度，对模拟结果与校正后的雷达测冰厚度数据进行对比(图5)，结果表明平均误差为±2 cm，但湖区不同部位冰厚模拟精度有所差异，其中误差较大区域主要是布哈河和哈尔盖河入湖口处，其次是湖区南侧和西南侧. 如1月13日鸟岛附近冰厚误差可达10.45 cm，这可能与湖冰形成初期厚度较薄以及布哈河淡水补给有关，布哈河水流汇入湖中不仅通过热量交换使湖水温度有所上升，而且可通过动力过程(如漩涡、湍流)使冰层底部不断被侵蚀，从而导致模拟数值大于湖冰实际厚度. 随着湖冰厚度增加，模拟精度有所改善，如鸟岛附近同一测区在1月30日和3月1日误差分别为0.42 cm和1.35 cm；湖区其他部位冰厚误差大多介于1～3 cm之间. 图5的另一个显著特征是，度日法湖冰生长模型模拟的冰厚值在3月中旬前普遍高估，而之后则存在低估现象，如哈尔盖河入湖口处冰厚误差在1月31日为3.28 cm，3月24日则为-2.23 cm，这表明湖冰生消不仅与气温相关，而且受湖面风速、太阳辐射、水深、湖流等要素影响[10]. 整体而言，在缺乏多要素观测条件下，采用度日法湖冰生长模型模拟的湖冰厚度虽与实测冰厚有所差异，但其精度亦在可接受范围内，这为基于气温观测数据预测青海湖湖冰厚度提供了一种捷径.

图5 2019年青海湖冰厚模拟及误差

3.3 2000-2019年青海湖冰厚时空变化特征

3.3.1 青海湖湖冰厚度时间变化特征 图6为2000-2019年青海湖年际(当年11月至翌年3月)及各月平均冰厚变化图. 在年际尺度上，近20年青海湖平均冰厚呈波动变化趋势，平均冰厚介于32～37 cm. 根据平均冰厚大小，可将此时段划分为3个阶段：(1)2000-2007年冰厚稳定期，此时段青海湖平均冰厚为34.65 cm，最大变幅为2.26 cm；(2)2008-2016年冰厚剧烈波动期，平均冰厚由2008年的32.78 cm快速增长至2010年的36.43 cm；之后湖冰厚度稳定处于高值，2010-2014年青海湖平均厚度为36.18 cm，最大变幅仅1.38 cm，其中2013年湖冰平均厚度(37.05 cm)达到近20年最大值；2015-2016年湖冰平均厚度出现骤降，两年平均冰厚仅32.87 cm，其中2016年为研究时段平均冰厚最小(32.08 cm)的年份；(3)2017-2018年冰厚恢复期，此时段青海湖湖冰平均厚度为34.79 cm，略高于第一阶段.

由图6亦可知，青海湖在11月份平均冰厚较小(10.61 cm)，12月份和1月份平均冰厚快速增加，可分别 达24.45和37.33 cm，增厚速率分别为0.45和0.41 cm/d；之后，2月份(0.29 cm/d)和3月份(0.14 cm/d)冰厚变化速率较小，这两个月的平均冰厚分别为45.55和50.03 cm. 这表明青海湖在冻结前期(12月-翌年1月)冰厚增长迅速，而中期(2-3月)冰厚增长变得缓慢. 湖冰厚度在各月的差异直接导致其平均冰厚在年际上的不同，如2010-2014年平均冰厚最大主要与这一时段1-3月冰厚普遍较大有关.

图6 2000-2019年青海湖封冻期年际及各月平均冰厚变化

3.3.2 青海湖湖冰厚度空间变化特征 基于2000-2019年ERA5 Climate Reanalysis气温数据，采用度日法湖冰生长模型计算同时段青海湖湖冰厚度可知，近20年间青海湖多年平均冰厚在空间上分布不均，总体上介于23～41 cm之间. 受ERA5 Climate Reanalysis气温数据空间分辨率影响，模拟的湖冰平均厚度与实际冰厚必然存在差异，如后者在空间上应为连续变化，但其仍能反映出青海湖湖冰厚度在空间上的显著差异性，即整体呈现北厚南薄、东厚西薄的特点(图7). 青海湖冰厚空间分布格局与其冻结消融空间模式密切相关. 已有研究表明，青海湖首先在邻近海晏湾的东侧边缘开始冻结，之后湖区东北部、西北部开始封冻，湖冰由湖岸逐渐扩张至湖心[10]. 相应地，冻结越早和冻结期越长的区域其湖冰越厚. 此外，湖冰厚度还与湖泊水深、入湖河流径流及湖下热流分布有关，如湖区西侧虽冻结较早，但布哈河河水汇入使该区域冰厚较小；而湖区西南侧(三江源祭海台周边)湖冰最薄的原因主要是该区域分布有温泉，高温水流涌入使湖水温度高于其他湖区，并不断融化湖冰，使湖冰厚度难以增大.

图7 2000-2019年青海湖平均冰厚空间分布及变化趋势

青海湖各区域2000-2019年平均冰厚的变化幅度不同. 以湖泊质心为中心将青海湖划分为4个区域，从左上方开始按顺时针方向分别称为Ⅰ、Ⅱ、Ⅲ、Ⅳ区(图7). 就多年平均变化速率而言，Ⅱ区的平均变化速率最快(-0.046 cm/a)，Ⅳ区变化速率最慢(-0.026 cm/a)，Ⅰ区与Ⅲ区平均变化速率分别为-0.031、-0.041 cm/a. 年际变化速率最大的区域位于湖区东南角(-0.051 cm/a)，最小的区域则位于西南角(-0.019 cm/a)，各区的变化趋势均为负值，表明2000-2019年青海湖各区域平均冰厚呈减小趋势，整体而言湖区西部变化幅度小于东部. 这是由于青海湖湖区东侧没有受到淡水补给以及盛行西北风的扰动，其冰厚的变化更能体现出随气候变化的较大变化幅度，而湖区西侧冰层受到此类干扰因素较多，在一定程度上减缓了气候变化所导致的冰厚变化，从而表现出更为稳定的变化幅度.

3.4 青海湖冰厚变化与湖冰物候特征的关系

3.4.1 2000-2019年青海湖湖冰物候特征 图8为2000-2019年每年11月至翌年4月青海湖逐日湖冰面积占比及日均温变化曲线，显然湖泊冻结过程较为缓慢，而湖冰消融过程则较为迅速. 具体而言，每年11月份青海湖湖冰面积占比介于1%～5%且波动剧烈，这主要与湖冰初生过程中冰层较薄，新生湖冰或碎冰极易受气温骤升或大风天气影响有关，导致湖冰在形成初期出现反复冻结消融现象[11]. 至12月初，当湖冰面积占比达到5%时，湖冰趋于稳定并随着温度降低湖冰范围逐渐扩大直至完全冻结(12月下旬或1月上旬). 3月中下旬，随着气温上升湖冰开裂消融，至4月上旬或下旬湖冰完全消融.

图8 2000-2019年青海湖湖冰面积占比及日均温变化

前人在研究青海湖湖冰物候特征时以湖冰面积占比达到10%作为湖泊开始冻结和完全消融日期判别依据[10]，此阈值最早由Reed等[28]在研究北美阿拉斯加西南地区湖冰物候特征时所提出，其认为当湖冰面积占比小于10%时湖冰会出现反复消融冻结现象. 由图8可知，当青海湖的湖冰面积占比接近5%时，湖冰冻结或消融发展趋势将不会改变，因此将湖冰面积占比连续3 d达到5%定义为湖泊开始冻结日期和湖冰完全消融日期，由开始冻结日期至湖冰面积占比首次达到95%的时间间隔定义为冻结时长，由湖冰面积占比最后达到95%日期至完全消融日期定义为消融时长，湖冰面积占比介于95%～100%的时长则为完全封冻时长. 统计表明，近20年青海湖平均开始冻结和完全冻结日期分别为12月9日和1月6日，平均开始消融和完全消融日期分别为3月23日和4月10日. 但是，在个别年份青海湖冻结和消融日期存在较大差异，如2016-2017年开始冻结和完全冻结均最晚，分别为2016年12月25日和2017年1月23日，这应与该年11-12月气温偏高有关，根据距青海湖最近的刚察气象站观测资料，2016年11-12月平均气温为-7.77℃，较2000-2019年同期气温偏高2.63℃. 开始冻结和完全冻结最早日期分别为2012年11月30日和2005年12月25日，开始消融和完全消融最早日期均出现在2016年，分别为3月15日和3月22日，最晚日期分别为2012年3月31日和2011年4月16日. 尽管2013-2017年青海湖开始冻结时间较晚(多为12月中旬)、开始消融和完全消融时间较早(分别出现在3月中旬和3月末)，但整体来看并未发生显著变化，4个时间节点的变幅均介于4～6 d.

3.4.2 冰厚变化与湖冰物候特征的关系 从冻结时长、完全冻结时长和消融时长变化(图9)来看，青海湖冻结时长明显大于消融时长，其平均时长分别为29和11 d，完全冻结时长平均为75 d，最大和最小时长分别为94 d(2005-2006年)和50 d(2016-2017年). 除个别年份外，青海湖冻结时长和消融时长虽有所波动，但整体变化不大；完全冻结时长在2014-2017年有明显的减少趋势，且2004-2006年也较小，与之相对应的刚察气象站1-3月平均气温也较高(普遍高于-7.80℃)，较2000-2019年1-3月平均气温高0.66℃. 与采用相同数据和方法而阈值为10%的研究结果[10]对比，两种结果表现出一定差异性，具体表现为平均开始冻结日期有所提前(相差7 d)，完全冻结日期和开始消融日期基本一致(相差1 d)，完全消融日期推后(相差6 d). 平均冻结时长较10%阈值的分析结果长14 d，完全冻结时长短2 d，消融时长仅长1 d.

湖冰物候(封冻期长短)是影响湖泊冰厚变化的重要因素之一. 根据分析青海湖各年冻结时长、完全冻结时长、消融时长、封冻期(湖泊开始冻结至完全消融[29])和湖冰年平均厚度可知(图9)，湖冰厚度与4个时期相关系数R2分别为-0.05、0.13、0.06、0.23(P值均小于0.05)，其中冰厚与完全冻结时长及封冻期的相关性较大，这是因为湖冰存在的时期越长，意味着湖泊水面被冰层覆盖的时间越长，除冰-气间的热量交换外，湖冰表面的高反射率使冰体吸收的太阳辐射减少，并进一步加剧近地面气温降低和冰-气热量交换增强，最终导致湖冰冻结深度更大. 表明冰厚垂直方向的增厚主要集中于完全封冻期，而冻结期与消融期则表现为湖冰在湖面水平方向的冻结扩张或消融破裂，故而与厚度变化相关性较低.

图9 2000-2019年青海湖湖冰物候与平均冰厚的关系

4 讨论

当湖泊表面被冰覆盖时，湖冰厚度增长主要取决于热力条件和水力条件，其中下冰面-湖水界面的热交换主要是湖冰底部与水流之间的湍流热交换，冰体内部通过导热完成热量传递，冰体表面则与大气发生热交换，并受气温、风速、短波和长波辐射的吸收反射及云量等条件影响. 研究表明，冰层垂向热传导通量是驱动冰层生长的主要能量源，随着冰厚的增加冰内垂向热传导通量逐渐下降. 当冰厚为12～30 cm时，湖冰上冰层垂向热传导通量随气温的变化存在明显高频波动特性，而冰厚大于30 cm的中下冰层波动特性明显衰减[30]. 青海湖每年11月至翌年1月气温温度梯度较大，冰层较薄，大气与湖水之间的热量传输剧烈，冰内垂向热传导通量波动明显，这导致湖冰的厚度增长速率更快；而2-3月温度梯度变化较小，当冰厚大于35 cm时，冰内垂向热传导通量波动明显衰减导致冰层生长速率减缓.

湖泊水域范围通常具有波动性，如在丰水期水位上升、面积增大，而在枯水期则呈相反趋势. 受降水和冰川融水增加以及蒸发减少等影响，青藏高原大多数地区湖泊面积近期呈增加趋势[9, 31]. 已有研究表明，1973-2018年青海湖水位和面积总体呈现先下降(减少)后上升(增加)趋势，并大致以2004年为界，2004-2018年湖泊水位上升2.07 m，面积增加19.59 km2，2018年的水域面积基本与20世纪70年代末期的面积持平[32]. 根据2019年夏季考察，青海湖面积仍在持续增加，鸟岛附近道路、停车场及周边草场均被淹没. 通常，湖泊水深越大，其相应的热容也越大，并直接影响湖冰的生消周期[9]. 通过对比本研究时段内各年青海湖冰厚与水位变化趋势发现，两者之间并不存在相关性，如2000-2004年青海湖水位持续下降，2004年之后水位快速上升，但其冰厚在对应时段内并没有发生类似的变化趋势，这可能与青海湖盐度较小及平均水深较浅有关，并导致冰下与湖水热量交换强度微弱.

随着青海湖面积变化，近年来湖岸线亦发生动态变化. 已有研究表明，近45年来岸线发生变化最大的区域主要是青海湖东岸的沙岛、西岸的鸟岛和铁布卡湾及北岸沙柳河入口区域，变化时间点与湖泊面积退缩及扩张时间一致[32]. 以2004年为界，前期湖岸线向湖内推进，尤其导致沙岛湖与青海湖主体脱离，2004年之后湖岸线持续向外后退，其中2017-2018年是相邻年间岸线变化最大的年份[31]. 青海湖冻结初期较薄的湖冰会被吹至湖岸边，沿着湖岸形成岸冰，而这主要与湖岸线形状和局地盛行风向有关[28]. 湖岸线较为复杂的区域(如湖区东部海晏湾、尕海、沙岛及湖区西部的铁卜卡湾)最先出现岸冰，且最早形成连续的固定冰，上述区域亦为青海湖空间平均冰厚较厚的区域. 虽然研究时段内岸线变化与平均冰厚变化并没有相同的变化趋势，但岸线变化最大的年份(2017-2018年)平均冰厚较前一年增加2.46 cm，亦是近20年来冰厚变化较大的年份(除2010-2011年较前一年增加2.93 cm).

湖冰厚度变化不仅影响温室气体排放强度[33]，而且通过控制湖面反射率和湖-气间物质与能量交换影响湖区热量和能量收支平衡[3,34-35]，并通过改变湖水温度分层控制湖泊生物地球化学的水生生物食物链进程[36]. 同时对于交通运输，渔业冬捕、冬季旅游和地方文化认同具有重要意义[37-38]. 近10余年青海湖冬季旅游持续增热，如地方部门举办冰面徒步穿越青海湖活动，个别游客也会私自下湖游玩. 尽管青海湖国家级自然保护区管理局已出台相关政策规范冬季冰上旅游安全，但每年仍有车辆陷入湖中事故发生，这一方面与青海湖范围广阔难以有效监管有关，另一方面也与游客对青海湖冰厚认识不足有关. 上节分析表明青海湖通常在每年1月中下旬平均冰厚可达35 cm左右，其承载能力也会得到极大增强，因此从1月下旬至3月上旬开展冰面活动较为安全. 但需要注意的是，因青海湖冰厚在空间上的异质性，在一些湖区应加强游客管理，尤其是在湖区西南侧.

5 结论

本文基于ERA5 Climate Reanalysis气温数据集、MODIS MOD09GQ数据产品和2019年湖冰钻孔测厚数据及雷达测厚数据，应用度日法湖冰生长模型计算青海湖2019年湖冰厚度，并对2000-2019年青海湖冰厚时空变化进行分析，得到以下结论：

1)分析2019年1-3月的实测数据可知，1月30日湖区平均冰厚为40.2 cm，3月2日平均冰厚测量值为43.5 cm，3月22日为51.4 cm，3月平均增长速率为0.30 cm/d，较2月份(0.12 cm/d)增长迅速. 实测冰厚空间整体呈现北厚南薄，东厚西薄的特点. 就冰厚平均增长速率而言，测区1、3、7仅有同一时间段内观测记录，有两段时间内观测记录的9个测区，主要有冰厚快速增长区(测区6、13、9)、冰厚平缓增长区(测区11、10、8)、冰厚较慢增长区(测区12、14). 测区2、4、5、15包含3段采样时间，其对应平均冰厚增长速率依次减小.

2)2018年11月-2019年3月湖冰厚度平均增长速率为0.34 cm/d，与实测数据相比模拟冰厚平均误差为±2 cm，但在湖区不同区域模拟精度有所差异，其中误差较大区域为布哈河口、哈尔盖河入湖处，以及湖区西南部以及湖区南部. 同时，冰厚模拟数值存在3月中旬前高估而之后低估现象.

3)2000-2019年青海湖年平均冰厚为32～37 cm，年际变化呈波动趋势，其中2000-2007年为冰厚稳定期，2008-2016年为冰厚剧烈波动期，2017-2018年为冰厚恢复期. 青海湖冻结前期冰厚增长迅速，12月与1月湖冰增长速率分别高达0.45和0.41 cm/d；2月后冰厚增长速率放缓，2月和3月的增长速率分别为0.29和0.14 cm/d，这与冰层垂向热传导通量在不同深度冰层内随气温变化的波动频率密切相关. 青海湖湖冰平均厚度在空间上存在显著差异，整体表现为北厚南薄、东厚西薄的空间格局，且湖区西部冰厚变化幅度小于东部.

4)近20年青海湖平均开始冻结和完全冻结日期分别为12月9日和1月6日，平均开始消融和完全消融日期分别为3月23日和4月10日，平均冻结时长和消融时长分别为29和11 d，完全封冻时长平均为75 d. 湖冰平均厚度变化与完全封冻时长以及封冻期相关性较高，而与冻结时长及消融时长相关性较小，表明冰层增厚主要集中于完全封冻期，冻结期与消融期则为湖冰在湖面的冻结扩张或消融破裂.

湖冰演化是一个极为复杂的过程，尽管本文模拟了2000-2019年青海湖冰厚空间分布特征及变化趋势，但因度日法湖冰生长模型仅考虑气温对冰厚的影响，并未考虑冰体内部导热过程及冰下的热力与水力作用，不能模拟消融期的冰厚变化及未能考虑积雪存在对冰厚模拟的影响，这使得模拟结果与真实湖冰厚度存在一定差异，在未来模拟青海湖湖冰厚度时应加以考虑. 同时，受试验观测区域限制，对青海湖中间大片区域的冰厚模拟结果仍有待验证.

致谢：在开展湖冰测厚试验期间，青海湖国家级自然保护区管理局和刚察县气象局工作人员、大连中睿科技发展有限公司陈洁博士及云南大学祁苗苗博士给予了大量帮助，在此致以谢忱.