渤海海冰演化的三维数值模型

张庆河，张 娜,

(1. 天津大学水利工程仿真与安全国家重点实验室，天津 300072；

2. 天津城市建设学院土木工程学院，天津 300384)

渤海海冰演化的三维数值模型

张庆河1，张 娜1,2

(1. 天津大学水利工程仿真与安全国家重点实验室，天津 300072；

2. 天津城市建设学院土木工程学院，天津 300384)

基于非结构化网格、有限体积三维海洋模型FVCOM建立了海洋-海冰耦合的渤海海冰数值模拟系统．模型中考虑了海冰动力过程和海冰热力过程，讨论了适合渤海的冰厚度分类规则、临界冰厚、短波辐射在冰内的透射率以及冰的热传输系数的合理取值．利用美国国家环境预报中心(NCEP)一小时一次的再分析气象数据作为驱动对2003—2004年冬季渤海海冰过程进行了模拟．模拟的渤海温盐计算值与实测值吻合较好，模拟的渤海冰盖分布与同期卫星遥感结果接近，模拟的辽东湾初冰日、终冰日的海冰面积、海冰厚度及其最大值和最大值对应的时刻与实测值基本一致．模拟与实测结果对比表明，建立的三维海冰数值模型能比较合理地反映渤海动力过程、气象条件影响下的温盐场过程以及渤海海冰的生消演化过程，具有较高的计算精度和计算效率，适合用于渤海海冰中长期数值模拟.

渤海海冰；有限体积海洋模型(FVCOM)；非结构化网格；数值模拟

海冰是较高纬度地区海域出现的自然现象，我国的渤海和黄海北部海域，冬季常受西伯利亚寒潮的影响，气温降低，每年都会形成海冰，且容易造成灾害[1]．1969年2—3月间渤海曾出现历史上罕见的冰情；2010年冬，渤海及黄海北部又发生了近30年来最严重冰情，海水结冰导致航道阻塞，船只、海上设施和海岸工程遭到破坏，港口码头封冻，水产养殖受损[2]．鉴于海冰灾害给人类在海上的生命财产安全以及生产活动造成的严重影响和危害，深入研究海冰生消演变规律、建立和完善海冰预警预报机制，对于海洋防灾减灾具有十分重要的意义．

海冰数值模拟是研究海冰生消变化规律的一种重要方法，在与气象预报结合的基础上可以进一步进行海冰预报．我国早期的海冰数学模型工作主要集中于海冰动力学的研究[3-4]，只适宜做非常短期的海冰预报，而海冰的生长与消融是同时受动力条件和热力条件影响的，特别是在进行较长期预报时必须综合考虑热力学影响．近年来随着冰-海耦合技术的不断提升，海冰数值模拟的时间尺度也从短期发展到中长期．唐茂宁等[5]以MM5数值产品作为气象输入条件，利用结构化网格的POM模式建立海冰-海洋耦合数值模式同时模拟动力学和热力学过程，对2003—2004年冬季渤海海冰过程进行了模拟；刘煜等[6]基于冰-海耦合模型，模拟了从1997—1998年冬季至2008—2009年冬季渤海海冰的季节性变化．在国外，Herman等[7]基于结构化网格的POM模型，同时考虑海冰热力学和动力学，模拟了波罗的海海冰的季节性变化；Gao等[8]基于有限体积法的FVCOM三维海洋模型开发了非结构化网格海冰-海洋耦合数值模型，该模型同时考虑海冰动力学和热力学过程，成功应用于北极海域海冰的数值模拟．为了深入研究具有复杂几何形状的渤海的海冰演化规律，发展渤海海冰预报的数值模型，笔者将基于FVCOM模型建立渤海海冰演化的海冰-海洋耦合三维数值模式，并收集资料验证了模型的合理性．

1 数值模型

海冰的数值模拟涉及到三维水动力学、海温、盐度以及海冰热力学、海冰动力学控制方程的求解及其耦合，图1给出了海洋-海冰耦合数值模型的组成情况和各子模型间的变量传递．下面对涉及的模型进行简要介绍．

图1 海冰数值模型的组成Fig.1 Composition of the sea ice numerical model

1.1 FVCOM三维海洋模型

FVCOM是描述海洋运动的三维模型，其基于直角坐标系的原始控制方程如下．

动量方程

连续性方程

海温方程

盐度方程

密度方程

式中：x、y和z分别代表坐标系的东向、北向和垂向坐标轴；u、v和w分别为x、y和z轴的速度分量；T为海温；S为盐度；ρ为密度；p为压强；f为科氏力系数；g为重力加速度；Km为垂向涡黏系数；Kh为垂向涡黏热扩散系数；Fu(Fv)、FT和FS分别代表水平动量项、水平热扩散项和水平盐扩散项．控制方程(1)～(7)可实现垂向上的σ坐标变换，即

式中σ在[-1,0]之间变化．在海表面z=b，σ=0，在海底zH=-，1σ=-；D为总水深，DHb=+，H和b分别为相对于0z=的底部水深和自由表面高度.

1.2 海冰模型的基本方程

海冰模型的基本方程主要包括海冰热力学方程、海冰动力学方程、海冰传输方程和海冰本构方程．

1.2.1 海冰热力学的控制方程

1) 来自大气到达海冰上表面的净热通量

式中：NETF为净热通量；SHF为感热通量；LHF为潜热通量；LWF↑为冰面向大气的长波辐射；LWF↓为大气逆辐射；SWF↓为来自大气到达海冰表面的总太阳辐射；α为冰面反照率；0I为短波辐射在冰内的透射率．

2) 海水结冰的判断条件

当海水温度低至临界冻冰温度fTSμ=-时开始形成海冰，μ为海水的冻结温度与盐度比例，μ= 0.054．用Ffrzmlt判断冻结或融化，其表达式为

式中：wT为海表面温度；wρ为海水的密度；pc为海水的比定压热容；mixh为海洋混合层深度；TD为海冰热力学的计算时间步长．当frzmlt0F＞时海水开始结冰，此时，冰底面(指海冰下表面)温度为临界冻冰温度fT.

3)海冰融化的判断条件

冰底面的增长和消融是由冰底面的传导热通量与海洋热通量之间的净热通量控制的，用BOTF表示，其表达式为

式中：hc为热传输系数；*u为海冰和海水之间的相对速度．

1.2.2 海冰动力学的控制方程

直角坐标系下，描述海冰运动的动量方程为

式中：m为单位面积的冰和雪的混合质量；us和vs为冰速在x和y方向上的分量；σij为内部应力张量，下标1、2和i、j分别代表x和y方向上的分量；(τax,τay)和(τwx,τwy)分别代表x和y方向上的海表面风应力和海表面水应力．表面风应力和表面水应力的矢量形式可写为

式中：c为值从0～1的海冰密集度；su为表面速度矢量；u、C、ρ分别代表速度矢量、拖曳力系数和密度；下标a和w分别代表空气和水．

海冰本构关系采用Hunke的弹黏塑性(EVP)模型[9]，其方程形式为

式中：E为由弹性模量定义的弹性参数；下标kk表示Einstein求和；η为剪切黏度；ζ为体积黏度；p为冰强度；δij为Kronecker函数，i=j，δij=1，i≠j，δij=0；ε˙ij为冰应变率．

1.2.3 海冰传输方程

海冰热力学和动力学模型求解的基本问题就是计算冰厚度随空间和时间的演化分布，这2个过程可通过海冰传输方程耦合起来，其标准对流方程形式为

1.3 模型的数值离散

FVCOM水动力、温盐方程以及海冰模型控制方程均采用有限体积法进行数值离散，潮流与温盐模型通过内外模分离的数值离散方法，二维外模数值格式为基于非结构化三角形网格的有限体积法，将连续性方程、动量方程在控制体区域积分后，通过改进的四阶龙格库塔方法求解，计算垂向平均的潮位和潮流．三维内模的动量方程和温盐方程采用简单的显式与隐式相结合的差分格式，对流项采用二阶精度的龙格库塔时间推进格式，垂向扩散使用隐式格式求解，计算三维潮流和三维温盐，海冰模型的数值离散采用控制方程的积分形式．关于详细的离散格式参见文献[10]．

2 渤海模型参数讨论与选取

2.1 水动力、温盐模型的参数

渤海水动力、温盐模型的主要计算参数及取值见表1．

表1 模型主要计算参数取值Tab.1 Values of main parameters in the model

2.2 海冰模型的参数

2.2.1 冰厚分类规则

海冰传输方程通过离散每个网格内的冰厚种类求解，不同冰厚种类的冰增长率是不一样的，因此冰厚的分类规则对海冰模拟有重要影响．20世纪70年代，Thorndike等[11]将极地海冰视为包含无冰海水和厚冰脊的混合物，并首次提出冰厚度的分布函数．21世纪初，Lipscomb[12]给出了适合于北冰洋等极寒海域的冰厚分类规则，通常包括无冰海水、一年冰、多年厚冰以及厚冰脊．渤海海冰的特点是冰期短(仅冬季有)、冰层薄(常冰年冰厚一般为0.20～0.40,m，最大冰厚0.60,m左右)，灾情最重的年份其最大堆积冰高度不超过3,m[13]．如果将北极海冰冰厚的分类规则应用于渤海，则会使薄冰与厚冰合并，不能准确模拟．敏感性测试结果表明，冰厚度的分类比较细时，模拟结果较合理，因此本文采用适合渤海的冰厚分类规则为：0、0.10,m、0.30,m、0.60,m、1.00,m、1.50,m、2.10,m和2.80,m．敏感性计算表明，对于比该规则更细的划分，模拟结果并没有改善．

2.2.2 临界冰厚

临界冰厚的设置是为了保证热力学计算的稳定性，是影响冰模型计算结果的一个重要参数．在热力学模式中，Hibler[14]提出在一个计算网格内存在2类冰厚分布，将薄冰和无冰海水归为一类，厚冰归为一类，针对北极海冰模拟给出了区分2类冰的临界冰厚为0.5,m．本文根据渤海海冰总体偏薄的特点，选取临界冰厚为0.01,m．

2.2.3 热传输系数

热传输系数hc和垂向涡黏热扩散系数hK对于模式的计算结果影响较大，尤其是热传输系数直接影响着海冰融化的速度．McPhee[15]在北极地区测得的适合该海区的海冰热传输系数为0.006，Josberger[16]在白令海和弗拉姆海峡测量得到的为2×10-4～8×10-4；Haapala等[17]在对波罗的海海冰的研究中取ch为2.8×10-4；Shirasawa等[18]在加拿大Resolute Passage海的实验结果ch为2.3×10-3，本文选取ch为2.2×10-5[19]．敏感性测试研究结果表明，越小则海冰融化速度越慢．垂向涡黏热扩散系数hK的选取关系到模式的计算稳定性，多次数值试验结果表明，选取1.0×10-4较合理．

2.2.4 短波辐射透射率

进入冰内的太阳短波辐射能缩短融冰期、阻止海冰厚度增加，是影响冰模型计算结果的一个重要参数．海冰内部可能含有大量气泡、卤水泡和泥沙等杂质，因此太阳辐射在海冰内的透射率并不是一个恒定值．根据Beer定律[20]，本文的透射系数为

式中：I(z)为冰表面以下深度为z处的太阳短波辐射；Fswabs为到达冰面除去反射部分的太阳辐射，Fswabs=(1-α)FSW↓；i0对于近红外辐射取值为0，对于可见光辐射取值为0.7；κi表示冰内太阳辐射的衰减系数，对于可见光波长设定为1.4,m-1．

3 渤海海冰模拟结果验证与分析

3.1 模型设置

3.1.1 空间网格与时间步长

采用上面建立的模型对渤海海冰进行了从2003年11月3日0：00至2004年2月29日23：00长达4个月的数值模拟．海冰热力学模拟的时间步长与内模式的时间步长一致，为20,s．为了保证EVP本构方程求解弹性波过程的稳定性，海冰动力学模拟的时间步长设为热力学时间步长的1/2，为10,s；外模式的时间步长为2,s．水平方向使用非结构化三角形网格进行划分，并在渤海湾、辽东湾和莱州湾近岸将网格局部加密，水平方向的最小网格步长为0.01°，最大网格步长为0.055°．渤海局部(辽东湾)及渤海测站布置的网格划分如图2所示．垂向采用σ坐标共划分了10层，并在水面附近加密．

3.1.2 近岸动边界处理

根据实际的沿岸冰分布形态，在FVCOM三维干湿动边界处理技术的基础上，设定潮间带海冰不随涨落潮移动，只要海温低于临界冻冰温度，则网格点干出时，依然有冰覆盖，冰密集度与冰厚为网格点干出前一刻的数值，干出点网格的冰速为0．上述做法可以保证海冰求解的稳定性，也接近实际情况．

3.1.3 模型初、边值及外部驱动条件

模型初值、边界条件及外部驱动条件如下．

(1) 模型计算的初始三维海温场由美国国家环境预报中心(NCEP)发布的气候预测系统(CFSR)一小时再分析海温数据插值而成，水平方向空间分辨率为0.3°[21]．初始盐度场参考1998年10月实测渤海近岸秋季表层、底层盐度分布[22]，结合NCEP一小时盐度数据插值而成．

(2) 开边界潮位由潮汐预报软件Chinatide提供[23]．

(3) 模型输入条件是指风、气温、相对湿度、云量、短波辐射、长波辐射和海表面气压等全场逐时驱动参数，其中气象参数用来计算感热通量、潜热通量和净热通量等各种热通量．外部驱动场来自NCEP一小时再分析数据，该数据的水平方向空间分辨率为0.3°，时间间隔为1,h．

图2 渤海局部(辽东湾)网格划分及渤海测站布置Fig.2 Computational meshes for part of the Bohai Sea(the Liaodong Bay)and the observation stations in the Bohai Sea

3.2 模拟与实测结果的比较

为了验证模型的合理性，在渤海海域收集了计算时间内的潮汐预报表潮位资料和各种实测资料，其中包括鲅鱼圈和塘沽潮汐预报表潮位资料，葫芦岛和龙口实测海温资料，大连盐度实测资料，MODIS卫星遥感资料以及JZ20-2平台冰厚资料，图2中标出了资料收集点的位置．以下给出模拟结果与所收集资料的比较情况．

3.2.1 潮位模拟与潮汐表结果比较

水动力模型计算了从2003年11月1日0：00至2004年2月29日23：00共4个月的渤海潮波，鲅鱼圈和塘沽2个站点的潮位计算值与天文潮潮汐表的对比结果见图3和图4，图中时间从2003年11月20日0：00至2003年11月30日0：00．由图3和图4可知，2个验潮站的计算潮位与潮汐表潮位吻合较好，水动力结果的合理性为模拟渤海温盐及海冰奠定了基础．

图3 鲅鱼圈计算水位与潮汐表对比Fig.3 Comparison between the calculated values and the water elevation from the tide table in Bayuquan

3.2.2 温盐模拟与实测结果比较和分析

图5所示为葫芦岛、龙口海温模拟值与实测值的对比情况，图6所示为大连盐度模拟值与实测值的对比情况．图5和图6中开始时间为2003年11月9日，结束时间为2004年2月28日；受海面结冰影响，实测数据只到2003年11月30日．由图5和图6可知，验证点温盐模拟结果与实测结果有很好的一致性.

图4 塘沽计算水位与潮汐表对比Fig.4 Comparison between the calculated values and the water elevation from the tide table in Tanggu

图5 葫芦岛和龙口海表面温度计算值与实测值Fig.5 Comparison between the calculated and measured sea surface temperatures for Huludao and Longkou

图6 大连海表面盐度计算值与实测值Fig.6 Comparison between the calculated and measured sea surface salinities for Dalian

为了直观观察渤海在海冰发生期的海温平面分布以及沿水深的垂向分布，图7给出了2004年1月25日模拟海表面温度场．图8以JZ20-2平台所在位置为例给出不同日期垂向海温分布模拟结果．从渤海温、盐的计算结果可知，由于冬季降水很少，径流较弱，风速较大引起海表面蒸发加速，致使盐度略有上升，但总体变化不大．冬季海温主要受冷空气的影响呈现温度不断降低的趋势，但海温的日变化较气温平稳，计算结果为冬季表层海温南高北低(如图5所示)，沿岸低外海高(如图7所示)．由于冬季大风使得表层与底层混合均匀，底层海温分布与表层一致(如图8所示)，这与以往的研究成果相吻合[24]．可以认为，温、盐的计算值与实测值吻合较好，为海冰模拟奠定了基础．

图7 2004年1月25日渤海海表面温度场分布Fig.7 Sea surface temperature distribution of the Bohai Sea on Jan. 25th，2004

图8 JZ20-2平台不同日期垂向海温分布Fig.8 Vertical distribution of sea water temperature on different dates at JZ20-2 platform position

图9 2003年12月7日渤海海冰数值模拟冰盖分布与MODIS遥感影像冰盖分布对比Fig.9 Comparison between the calculated ice cover and remote sensing image of MODIS in the Bohai Sea on Dec. 7th，2003

3.2.3 海冰模拟与卫星遥感结果比较

渤海海冰冰盖模拟的结果与MODIS卫星遥感的同期冰盖分布对比如图9～图12所示．根据MODIS遥感影像监测，2003—2004年冬季为常冰偏轻年，初冰日为2003年12月6日，终冰日2004年2月24日，冰期76天．模式计算的初冰日与遥感监测的日期吻合，模式计算的终冰日为2004年2月23日，只比遥感监测的日期提前1,d．

图13以辽东湾为例给出了海冰面积变化过程的模拟结果与卫星遥感结果对比．从MODIS卫星遥感结果看，辽东湾最大海冰面积发生于2004年1月25日，约为13,814,km2，整个渤海最大海冰面积发生在2004年1月25日，约为20,402,km2．模式计算的海冰面积最大时刻与遥感监测的时间一致，但最大面积较监测结果略偏小，辽东湾最大海冰面积约为12,432,km2，渤海最大海冰面积约为18,361,km2．

图10 2004年1月25日渤海海冰数值模拟冰盖分布与MODIS遥感影像冰盖分布对比Fig.10 Comparison between the calculated ice cover and remote sensing image of MODIS in the Bohai Sea on Jan. 25th，2004

图11 2004年2月2日渤海海冰数值模拟冰盖分布与MODIS遥感影像冰盖分布对比Fig.11 Comparison between the calculated ice cover and remote sensing image of MODIS in the Bohai Sea on Feb. 2nd，2004

图12 2004年2月14日渤海海冰数值模拟冰盖分布与MODIS遥感影像冰盖分布对比Fig.12 Comparison between the calculated ice cover and remote sensing image of MODIS in the Bohai Sea on Feb. 14th，2004

图13 辽东湾海冰数值模拟面积与MODIS 遥感影像面积对比Fig.13 Comparison between the calculated sea ice area and the remote sensing data of MODIS in the Liaodong Bay

就渤海湾和莱州湾而言，遥感监测的渤海湾和莱州湾于2004年1月22日开始有冰盖分布，2004年2月2日渤海湾和莱州湾内海冰全部融化；模式计算的渤海湾于2004年1月20日开始结冰，1月30日冰期结束；莱州湾于2004年1月21日开始结冰，1月31日冰期结束，模式计算的渤海湾与莱州湾冰期仅有短短的10,d，与遥感监测的日期吻合较好．

图14显示了JZ20-2平台附近海域的冰厚模拟值与实测值[3]对比情况，图中初始时刻为2003年12月31日0：00．由图14可知，JZ20-2平台的模拟值与实测值的初冰日、终冰日吻合较好，模拟的冰厚最大值以及最大值对应的时刻、冰厚变化趋势与实测较为一致．

图14 JZ20-2平台冰厚模拟值与实测值对比Fig.14 Comparison between the calculated and the measured sea ice thicknesses at JZ20-2 platform

4 结 论

基于FVCOM三维海洋模型建立了海洋-海冰耦合的渤海海冰数值模拟系统，讨论了适合渤海的冰厚分类规则、临界冰厚和短波辐射在冰内的透射率以及冰的热传输系数的合理取值，利用NCEP一小时再分析气象数据作为驱动对2003—2004年冬季渤海海冰过程进行了模拟，得出如下结论：

(1) 葫芦岛与龙口2个近岸实测站点的海表面温度、大连的海表面盐度的验证结果表明计算值与实测值吻合较好．由2004年1月25日渤海海表面温度分布及JZ20-2平台位置不同时刻垂向海温分布结果可知，冬季表层海温呈南高北低、沿岸低外海高特征．由于冬季大风使得表层与底层混合均匀，海温沿垂向分布非常均匀，底层海温与表层几乎相同．

(2) 海冰模拟结果与MODIS卫星遥感的同期冰盖分布、海冰面积及JZ20-2平台附近海域的冰厚实测值对比，初冰日、终冰日吻合较好，模拟的冰面积和冰厚最大值及最大值对应的时刻，以及冰面积、冰厚变化趋势均与实测较为一致．

(3) 模拟和实测结果比较表明，本文建立的渤海三维海冰数值模型能比较合理地反映渤海动力过程、气象条件影响下的温盐场过程，再现了渤海海冰的发生、发展和消融过程．基于FVCOM三维海洋模型，综合考虑海冰热力学和动力学过程的海洋-海冰耦合模式，具有与气象预报相结合进行海冰预测的潜力，非结构化网格可实现对复杂岸线模拟及对关心区域局部加密，未来为精细化模拟特别是港口海冰发展模拟打下了基础，值得深入研究和推广．

［1］ 王相玉，袁本坤，商 杰，等. 渤黄海海冰灾害与防御对策[J]. 海岸工程，2011，30(4)：46-55.

Wang Xiangyu，Yuan Benkun，Shang Jie，et al. Disaster of sea ice and defense countermeasures in the Bohai Sea and the Yellow Sea[J]. Coastal Engineering，2011，30(4)：46-55(in Chinese).

［2］ 孙 劭，苏 洁，史培军. 2010年渤海海冰灾害特征分析[J]. 自然灾害学报，2011，20(6)：87-93.

Sun Shao，Su Jie，Shi Peijun. Features of sea ice disaster in the Bohai Sea in 2010[J]. Journal of Natural Disasters，2011，20(6)：87-93(in Chinese).

［3］ 季顺迎，岳前进，赵 凯. 渤海海冰动力学的质点网格法数值模拟[J]. 水动力学研究与进展，2003，18(6)：748-760.

Ji Shunying，Yue Qianjin，Zhao Kai. Numerical simulation of sea ice dynamics with particle-in-cell approach in the Bohai Sea [J]. Journal of Hydrodynamics，2003，18(6)：748-760(in Chinese).

［4］ 王瑞学，季顺迎，岳前进. 海冰动力学数值模拟中改进的PIC方法[J]. 计算力学学报，2006，23(4)：440-446.

Wang Ruixue，Ji Shunying，Yue Qianjin. Modified PIC method for numerical simulation of sea ice dynamics[J]. Chinese Journal of Computational Mechanics，2006，23(4)：440-446(in Chinese).

［5］ 唐茂宁，刘钦政，刘 煜，等. 渤海海冰季节演变的数值模拟[J]. 海洋预报，2010，27(2)：48-52.

Tang Maoning，Liu Qinzheng，Liu Yu，et al. The numerical simulation of sea ice in the Bohai Sea[J]. Marine Forecasts，2010，27(2)：48-52(in Chinese).

［6］ Liu Yu，Liu Qinzheng，Su Jie，et al. Seasonal simulations of a coupled ice-ocean model in the Bohai Sea and North Yellow Sea since the winter of 1997/1998[C]//Proceedings of the Twenty-First(2011) International Offshore and Polar Engineering Conference. Maui，Hawaii，USA，2011：19-24.

［7］ Herman A，Jedrasik J，Kowalewski M. Numerical modelling of thermodynamics and dynamics of sea ice in the Baltic Sea[J]. Ocean Science，2011，7：257-276.

［8］ Gao Guoping，Chen Changsheng，Qi Jianhua，et al. An unstructured-grid，finite-volume sea ice model：Development，validation，and application[J]. J Geophys Res，2011，116(9)：2010JC006688.

［9］ Hunke E C. Viscous-plastic sea ice dynamics with the EVP model：Linearization issues[J]. J Comput Phys，2001，170：18-38.

［10］ Chen C S，Cowles G，Beardsley R C. An Unstructured Grid，Finite-Volume Coastal Ocean Model：FVCOM User Manual[R]. USA：School of Marine Science and Technology，University of Massachusetts-Dartmouth，2006.

［11］ Thorndike A S，Rothrock D A，Maykut G A，et al. The thickness distribution of sea ice[J]. J Geophys Res，1975，80(33)：4501-4513.

［12］ Lipscomb W H. Remapping the thickness distribution in sea ice models[J]. J Geophys Res，2001，106(C7)：13989-14000.

［13］ 白 珊，刘钦政，李 海，等. 渤海的海冰[J]. 海洋预报，1999，16(3)：1-9.

Bai Shan，Liu Qinzheng，Li Hai，et al. Sea ice in the Bohai Sea of China[J]. Marine Forcasts，1999，16(3)：1-9(in Chinese).

［14］ Hibler W D. A dynamic thermodynamic sea icemodel[J]. J Phys Oceanogr，1979，9：815-846.

［15］ McPhee M G. Turbulent heat flux in the upper ocean under sea ice [J]. J Geophys Res，1992，97(C4)：5365-5379.

［16］ Josberger E G. Bottom ablation and heat transfer coefficients from the 1983 marginal ice zone experiment [J]. J Geophys Res，1987，92(C7)：7012-7016.

［17］ Haapala J，Lepparanta M. Simulating the Baltic Sea ice season with a coupled ice：Ocean model[J]. Tellus A，1996，48：622-643.

［18］ Shirasawa K，Ingram G R. Currents and turbulent fluxes under the first-year sea ice in resolute passage，Northwest Territories，Canada[J]. J of Marine Syst，1997，11：21-32.

［19］ 苏 洁，吴辉碇，刘钦政，等. 渤海冰-海洋耦合模式[J]. 海洋学报，2005，27(1)：19-26.

Su Jie，Wu Huiding，Liu Qinzheng，et al. A coupled ice-ocean model for the Bohai Sea[J]. Acta Oceanologica Sinica，2005，27(1)：19-26(in Chinese).

［20］ Hunke E C，Lipscomb W H. CICE：the Los Alamos Sea Ice Model Documentation and Software User’s Manual Version 4. 0 LA-CC-06-012[R]. USA：T-3 Fluid Dynamics Group，Los Alamos National Laboratory，2008.

［21］ 张 寅，罗亚丽，管兆勇. NCEP全球预报系统在ARM SGP站点预报大气温度、湿度和云量的检验[J]. 大气科学，2012，36(1)：170-184.

Zhang Yin，Luo Yali，Guan Zhaoyong. Temperature，relative humidity，and cloud fraction predicted by the NCEP global forecast system at the ARM SGP site during 2001—2008：Comparison with ARM observations[J]. Chinese Journal of Atmospheric Sciences，2012，36(1)：170-184(in Chinese).

［22］ 程济生. 黄渤海近岸水域生态环境与生物群落[M].青岛：中国海洋大学出版社，2004.

Cheng Jisheng. The Ecological Environment and Biological Communities in the Coastal Waters of Bohai Sea and Yellow Sea[J]. Qingdao：Ocean University of China Press，2004(in Chinese).

［23］ 李孟国，郑敬云. 中国海域潮汐预报软件Chinatide的应用[J]. 水道港口，2007，28(1)：65-68.

Li Mengguo，Zheng Jingyun. Introduction to Chinatide software for tide prediction in China seas[J]. Journal of Waterway and Harbor，2007，28(1)：65-68(in Chinese).

［24］ 宋文鹏. 渤海冬、夏季温盐场结构及其海流特征分析[D]. 青岛：中国海洋大学海洋环境学院，2009.

Song Wenpeng. The Analysis of the Structure of T-S and the Current Characteristics in Baohai Sea During Winter and Summer[D]. Qingdao：College of Physical and Environmental Oceanography，Ocean University of China，2009(in Chinese).

A Three-Dimensional Numerical Model of Sea Ice Evolution in the Bohai Sea

Zhang Qinghe1，Zhang Na1,2
(1. State Key Laboratory of Hydraulic Engineering Simulation and Safety，Tianjin University，Tianjin 300072，China；2. School of Civil Engineering，Tianjin Institute of Urban Construction，Tianjin 300384，China)

An ice-ocean coupled sea ice numerical simulation system of the Bohai Sea，including the dynamic processes and thermodynamic processes of sea ice，was developed based on the unstructured grid three-dimensional finitevolume coastal ocean model(FVCOM). Reasonable values of ice thickness，critical thickness for thinnest ice，ratio of shortwave radiation penetrating into the ice and heat transfer coefficient were discussed. The sea ice evolution of the Bohai Sea in 2003—2004 winter was simulated with the hourly time-series atmospheric data of NCEP(National Centers for Environmental Prediction)as a driving force. The simulated results of sea water temperature and salinity of the Bohai Sea are in good agreement with the measured data. The simulated ice coverage of the Bohai Sea is close to the image obtained by satellite remote sensing in the same period. The simulated sea ice area and sea ice thickness of Liaodong Bay are basically the same as the measured data for the start date，the end date，the maximum value and the corresponding time. It is shown that the presented three-dimensional numerical model of sea ice can reasonably simulate the process of hydrodynamic，sea water temperature and salinity，and the evolution of sea ice forced by meteorological conditions. The model can be applied to medium- and long-term sea ice numerical simulation with high accuracy and computational efficiency.

sea ice in the Bohai Sea；finite-volume coastal ocean model；unstructured grid；numerical simulation

O242；P731.15

A

0493-2137(2013)04-0333-09

DOI 10.11784/tdxb20130408

2012-12-04；

2012-12-13.

国家自然科学基金创新研究群体科学基金资助项目(51021004)；天津市自然科学基金重点资助项目(12JCZDJC30200).

张庆河（1966— ），男，博士，教授，qhzhang@tju.edu.cn．

张 娜，coastlab@163.com.