陆洋 王晓春

(南京信息工程大学海洋科学学院, 南京 210044)

提要 融池是影响北极海冰变化的重要因子。但在目前广泛使用的海冰模式CICE6.0中, 融池模拟与观测存在较大差异。CICE6.0中的CESM融池参数化方案可以计算融池覆盖率与深度, 但其中的部分重要参数存在一定的经验性, 影响了融池模拟。本研究为CICE6.0海冰模式的CESM融池参数化方案开发了伴随模式, 利用CICE6.0海冰模式、融池伴随模式和L-BFGS极小化算法, 使用一年冰及多年冰区域的MODIS卫星融池覆盖率观测数据, 对CESM参数化方案中的融池纵横比参数进行了分时段分区的参数估计。结果表明: 使用伴随模式可以有效地对融池参数进行估计; 得到的融池参数随时空变化, 符合融池过程时空变化强烈的特征; 估计的参数值用于模拟, 减小了融池覆盖率的模拟误差, 与MODIS观测更为一致。

0 引言

北极夏季, 海冰和冰上积雪开始融化形成融水, 融水在自身重力作用下沿着冰表面流动, 汇聚于冰表面低洼处形成大小融池。由于较高的反照率, 海冰可以反射大部分的入射短波辐射。海冰融化将会使得更多的热量被极区吸收, 从而促使更多的海冰融化, 称为“海冰-反照率”正反馈[1]。融池的反照率低于海冰, 融池的形成将降低海冰反照率, 增强正反馈, 起到加速海冰融化的作用。因此, 融池在极区气候系统中起到重要作用, 使得其成为近年来气候科学家关注的热点。随着北极海冰融化加剧[2], 一年冰比重显著增加[3-4], 而一年冰平整的表面有利融池形成后在水平方向扩展。尽管缺乏长期观测数据, 但一般认为, 由于一年冰的增加, 北极的融池覆盖率正在增加[5]。可以预料, 融池将对未来北极地区的气候产生更明显的影响。因此, 融池及其在海冰模式中的参数化研究具有十分重要的意义。

由于北极地区自然环境的限制, 融池的现场观测一直比较缺乏。融池覆盖率(Melt Pond Fraction, MPF)、反照率是观测的重点。北极海表面热平衡收支观测计划(Surface Heat Budget of the Arctic Ocean, SHEBA)执行期间曾对海冰反照率进行了长时间连续观测[6], 测量得到的融池反照率随时间迅速变化, 变化范围达到0.1～0.65, 这与融池自身的发展演化密切相关。Eicken等[7]基于观测指出季节尺度融池发展存在4个大致的阶段: (1)融池形成阶段; (2)融水流失阶段; (3)融池稳定扩展阶段; (4)融池重新冻结阶段。不同阶段的融池行为和控制机制不同。卫星遥感为融池观测提供了新的数据。2012年, Rösel等[8]利用MODIS(Moderate Resolution Imaging Spectroradiometer)光学传感器反射率数据反演了融池覆盖率, 提供了2000—2011年每年5—9月中8天一次的北极融池覆盖率产品。近年来, 开始利用微波遥感对融池进行观测[9-11], 但尚未发布产品。现场观测和卫星观测均表明, 融池过程具有明显的时空变化特征, 海冰模式中的融池参数化需要考虑这种特征。

融池在海冰模式中的参数化是非常重要的研究课题。由于融池在空间及时间上的明显变化,根据冰表面特征直接指定固定的冰面反照率的做法, 无法很好考虑融池, 可能造成较大的模拟偏差[6]。为在模式中对反照率进行符合实际的计算,需要发展融池模式, 将其纳入海冰模式中。融池观测为融池模式的发展奠定了基础, 目前已有几种复杂程度不同的融池模式用在海冰模式中。早期的融池模式[12-14]为一维模式或者二维模式, 往往与现有气候系统模式中的海冰模式不匹配, 不能很好地纳入海冰模式中。Holland等[15]提出了一个简单的融池参数化方案, 该方案可以计算融水体积, 并根据给定的融池纵横比(融池深度与覆盖率之比)分配融水体积。CICE海冰模式(Los Alamos Sea Ice Model)采用了该方案并称为CESM(Community Earth System Model)方案。Flocco等[16-17]基于海冰模式中的厚度分布函数(Ice Thickness Distribution, ITD)建立了更为复杂的融池方案, 考虑了冰面地形对融水分配的影响,利用海冰厚度来确定冰面地形。在该方案中, 融水产生后首先覆盖最低的冰面, 然后覆盖次低的冰面, 逐层累加。此外, 该方案还考虑了由于海冰微孔隙结构所导致的融池水渗流等过程。CICE引入该方案称为地形方案(Topography Scheme,TOPO)。Hunke等[18]综合了CESM方案和TOPO方案, 提出平整冰方案(Level Ice Scheme, LVL),利用CICE中已有的平整冰变量, 假定融池只能存在于平整冰上, 并且假定融池体积的变化存在纵横比, 从而分配融水体积。王传印等[19]结合观测数据对三种方案在CICE5.0中的表现进行了较为详细的对比研究, 认为CESM方案的融池冻结条件更为合理, 但存在多年冰上积雪不能完全融化、融池被低估的问题。

但是, 以上融池参数化方案均建立在一定的假设上, 由于对某些过程缺乏认识, 包含了一些不确定性较大的参数, 降低了海冰模式的准确性和可靠性[20]。对于这部分不确定参数, 一般可以通过敏感性试验的方法来研究和减小其不确定性。但这种方法有很强的局限性, 一方面气候系统是复杂系统, 存在大量反馈过程, 另一方面在待研究参数较多的情况下, 这种方法所需要的计算代价较为高昂[21]。伴随模式(Adjoint Model,ADM)为融池参数化方案中参数的估计提供了另外一种可能性[22-23]。伴随模式是原始非线性模式的切线性模式(Tangent Linear Model, TLM)的转置, 利用伴随模式进行一次积分可以得到模式结果对所有模式变量(包括边界条件、初始条件和模式参数)的导数, 进而可以定量地对参数进行调整。伴随模式已经广泛用于模式参数估计和初始场估计中。比如, Liu等[24]利用MITgcm(Massachusetts Institute of Technology General Circulation Model)海冰-海洋耦合模式及其伴随模式进行了海冰模式中海冰密集度初始场优化的研究,得到的海冰密集度初始场和实际更接近, 说明伴随模式可以用于对海冰模式初始场进行调整。Fenty等[22]利用MITgcm及其伴随模式进行了2004年海冰-海洋耦合状态估计, 同化了海洋水文资料和海冰密集度资料, 相比未同化资料的试验, 对海冰和海洋状态的估计都更加准确。Kim等[23]利用CICE模式的伴随模式对模式中21个热力和动力学参数同时进行了敏感性分析, 并利用计算得到的梯度信息对参数进行了试验性调整,验证了伴随方法在海冰参数估计和优化中的可行性。以上研究表明, 伴随模式已经在海冰模式参数估计中得到了广泛应用, 但是这些研究中很少考虑融池参数。由于融池过程对于海冰模式的重要性以及融池模式中存在许多不确定参数, 有必要利用伴随模式考虑重要的融池参数, 并通过优化算法对其进行估计, 以减小海冰模式在融池模拟中的误差。

CESM融池参数化方案的融池纵横比参数是基于SHEBA一年冰观测数据的经验参数[25]。由于多年冰上融池倾向于在深度方向上增长, 和一年冰融池存在显著不同, 这个固定值参数对整个北极海冰的适用性较低。因此, 本文的研究目标为通过伴随模式方法, 估计海冰模式中CESM融池参数化方案的融池纵横比参数, 以改善该融池方案对融池覆盖率的模拟结果。第1节介绍本文所用的数据、CICE6.0海冰模式和CESM融池参数化方案。第2节介绍本文所构建的融池伴随模式、L-BFGS极小化算法、伴随模式参数估计算法以及利用控制试验的模拟数据对参数估计算法的验证结果。第3节针对多年冰及一年冰海域, 利用MODIS融池观测数据作为观测约束, 对融池参数进行最优估计, 并将其用于模拟, 检验其对融池模拟的影响。第4节为本文结论与讨论。

1 数据、CICE6.0海冰模式及融池参数化方案

1.1 融池覆盖率与海冰密集度观测数据

MODIS融池覆盖率数据为德国汉堡大学研发的产品[8], 下载自汉堡大学综合气候数据中心网站: http://icdc.zmaw.de。该数据时间跨度为2000—2011年, 每年5月9日—9月6日, 时间分辨率为8天一次。空间范围为180°W～180°E、60°N～90°N, 空间分辨率为12.5km×12.5km。卫星数据和现场观测验证数据的均方根误差(Root Mean Squared Error, RMSE)的范围在3.8%～11.2%之间。MODIS是光学传感器, 因此反演的融池覆盖率中包含云层所引起的数据缺测。这些缺测主要发生在2000年、2001年、2002年和2007年的80°N以北区域。较小的缺测区域已经在数据发布前进行了插值, 但忽略了面积大于12.5 km2的区域[26]。

本文使用的海冰密集度数据为英国气象局哈德莱中心(Met Office Hadley Centre)的海冰和海表面温度数据集HadISST1[27]。HadISST1整合了来自历史观测和卫星观测的海冰数据, 提供了自1871年开始, 空间分辨率为1°×1°的全球逐月海冰密集度。

1.2 CICE6.0海冰模式及ICEPACK模块

本文使用的CICE6.0海冰模式由美国Los Alamos国家实验室研发[28], 建立在海冰厚度分布函数[29]的基础上, 是一个大尺度热力-动力学模式,考虑了成脊等次网格尺度机械过程的影响[30]。模式中使用位移极点网格或者三极网格来处理北极点的奇点问题[31]。ICEPACK是CICE6.0的一个子模块[32], 刻画了海冰热力学和生物地球化学等网格元内的垂向物理过程(Column Physics)。ICEPACK通过计算大气-海冰边界层和海洋-海冰边界层的能量收支来更新冰的温度, 并以此确定冰的增长或融化。本文使用CICE6.0为ICEPACK提供的单独的驱动程序, 独立运行ICEPACK来模拟单点(垂直柱)的海冰状态变化。

海冰模式的大气强迫数据使用基于日本气象厅(Japan Meteorological Agency, JMA)55年再分析资料JRA-55的海洋-海冰模式大气强迫数据集JRA55-do[33]。本文研究的时间段为2005年, 所用变量包括10米处的气温、10米处的比湿、10米处的风矢量、向下短波辐射、向下长波辐射和降水。JRA55-do数据的时间间隔为3 h, 本文中将其线性插值到1 h间隔, 以匹配ICEPACK的1 h的积分步长。

1.3 CESM融池参数化方案

CESM融池参数化方案是CICE6.0海冰模式中3个融池参数化方案之一, 和文献[15]中的方案有一些微小的不同之处。它可以计算融池的覆盖率与深度, 描述了融池从生成到冻结的过程。该方案的核心是建立了融池深度和覆盖率之间的线性关系, 认为深度随覆盖率线性增长, 即:

式中,hp表示融池深度,δp为融池纵横比,ap表示融池覆盖率。在CESM方案的默认设置中δp为0.8。融池水来自于海冰融水、积雪融水和降雨:

其中为模式中第 个时间步的单位面积融池体积(m),ρi和ρs分别为冰和雪的密度(kg·m-3),ρw为融水的密度(kg·m-3),frain为降水率(kg·m-2·s-1),Δt为时间步长(s), Δhi和Δhs分别为在Δt内海冰顶层的融化量(m)和积雪的融化量(m)。部分融水会从海冰边缘、冰上裂缝和海豹呼吸洞等流失入海,r为融水进入融池的比例:

其中ai为海冰密集度。CICE6.0中rmin和rmax的默认值分别为0.15和0.7。而在文献[15]中, 该比例为:r=0.15+0.7ai。气温降低时, 融池逐渐重新冻结, 体积减小:

其中r2为系数0.01,Tp为参考温度-2℃,Tsfc为海冰表面温度(℃)。

2 融池伴随模式及融池纵横比参数估计算法

2.1 融池方案的伴随模式及其验证

本文为对融池模式中的参数进行估计, 研发了融池模式的伴随模式。构造一个数值模式的伴随模式有两种路径: 第一种路径是基于解析的伴随模式方程组离散得到伴随模式, 第二种路径是从原模式的代码转换得到。前者只适用于比较简单的模式, 后者则可用于复杂模式。在第二种路径中, 可以借助自动微分(Automatic Differential,AD)软件工具来自动生成伴随模式代码[34], 这对于在设计之初考虑了数据同化问题并为之预先优化了代码结构的模式(如MITgcm)是方便的。另外一些模式, 其伴随模式需要直接转换(人工逐行编写伴随模式代码, 其遵循的编码基本规范和自动微分软件是相同的)得到, 如日本气象厅气象研究所(Meteorological Research Institute, MRI)的海冰模式伴随模式[35]。本文所使用的CICE6.0模式并没有为开发伴随模式进行特别的代码优化,无法使用自动微分工具生成伴随模式[36]。因此,本文采用了直接转换的方式, 这样可以在一定程度上保证伴随模式的正确性。

如1.3节所述, CESM融池方案在物理和数值计算上由4个部分组成: (1)首先计算由于冰雪的融化和降水所带来的融水体积增量; (2)其次, 计算融水流失比例, 按比例扣除融水得到最终进入融池的融水体积增量; (3)根据气温和冰表面温度,计算缩减后的融池水体积; (4)根据融池水体积,通过指定的融池纵横比参数, 同时考虑海冰厚度对融池深度的限制, 计算融池覆盖率和深度。根据这样的特点, 本文开发伴随模式的简要步骤如下: 首先, 确定融池方案的输入变量和输出变量,选择融池纵横比δp作为输入变量, 输出变量为融池方案中依赖于δp参数的变量; 其次, 逐行线性化融池方案得到切线性模式; 最后, 将切线性模式算子转置得到伴随模式算子。本文伴随模式中的基本态是重新计算的。对于模式中的分支结构语句(比如IF…ELSE…END)的处理方法是: 对两个分支的语句分别进行线性化, 分支结构的判断式若涉及输入变量, 则用其基本态代替。对于融池方案的4个部分, 分别写出伴随代码, 最终将计算顺序颠倒组成完整的融池伴随模式。伴随模式编写的具体规范和细节可以参考MM5伴随模式系统的技术报告[37]。

得到的切线性模式和伴随模式必须进行正确性检验才能使用。当扰动足够小时, 由正向模式控制积分及其扰动积分计算得到的扰动大小应该和切线性模式计算结果非常接近。因此, 可以利用正向模式, 按下式来检验切线性模式的正确性:

其中,Q是正向模式,z是正向模式输入变量向量,α为一尺度系数,h是输入变量的扰动,代表从正向模式控制积分及其扰动积分的输出所计算的扰动大小。P是切线性模式,Ph代表将h作为输入运行切线性模式得到的扰动大小。在实际计算中,h应该足够小, 使得可以忽略其高阶项,同时应该大于机器舍入误差。当尺度系数变小时,正向模式得到的扰动值与切线性模式得到的扰动值之比R趋向于1。

可以利用切线性模式, 按照下式检验伴随模式的正确性:

其中,P是切线性模式,dz是切线性模式的输入变量向量,PT是伴随模式。将dz作为输入, 运行切线性模式得到Pdz。将Pdz作为输入, 运行伴随模式PT得到PTPdz。＜＞表示向量内积。若上式在机器误差精度内成立, 则伴随模式代码通过验证。

按照公式( 5 )对切线性模式进行检验。其中尺度系数α从1.0减小到0.1, 扰动量h保持为0.2。结果如图 1所示, 随着α向0逼近,R逼近1, 说明切线性模式的代码是正确的。按公式(6)对伴随模式进行了检验, 其中dz=αh, 尺度系数α和扰动向量 的设置同切线性模式检验。伴随模式检验结果如图 2所示, 随着α逼近0, 公式(6)两边之差的绝对值逐渐逼近0, 说明伴随模式的代码是正确的。

图1 CESM融池切线性模式检验结果.α为扰动尺度系数, R(公式5)为正向模式积分得到的扰动大小与切线性模式积分得到的扰动大小的比Fig.1.Correctness check of CESM melt pond scheme tangent linear model.The α is a scaling factor.R (defined in Equation 5) is the ratio between the norm of perturbation from forward model integration and that from tangent linear model integration

图2 CESM融池伴随模式检验结果.α为扰动尺度系数,Δ(公式6两边之差的绝对值) 为切线性模式计算的扰动大小与伴随模式计算的扰动大小之差的绝对值Fig.2.Correctness check of CESM melt pond scheme adjoint model.The α is a scaling factor.Δ (the norm of difference between the left-hand side and the right-hand side of Equation 6) is the difference between the perturbation from the tangent linear model and that from the adjoint model

2.2 L-BFGS极小化算法及目标函数

正向模式、伴随模式和极小化算法的结合可以实现对模式参数的最优估计。本文中所使用的极小化算法为 L-BFGS算法(Limited-memory Broyden-Fletcher-Goldfarb-Shanno Algorithm), 它是一种拟牛顿法, 可以应用于大规模优化计算[38]。算法按照下述公式对被估计参数进行迭代:

其中,Xk为被估计参数的第k次迭代值,为第k次迭代的被估计参数目标函数的梯度, 可以通过伴随模式给出,ρk为第k次迭代的步长, 由极小化算法自动调整。迭代终止条件为: (1)迭代次数大于2000次或者(2)达到指定精度, 满足其中,为目标函数梯度的范数,ε的值在本文中设为10-5,为被估计参数的范数。

本文的研究目标为估计融池纵横比参数, 减小融池覆盖率模拟与观测之间的误差, 所以本文定义目标函数为:

其中,J为目标函数,为第i个观测时间点MODIS观测的融池覆盖率,为对应第i个观测时间点模式模拟的融池覆盖率,N为观测次数。目标函数J代表了模式和观测之间的误差, 当目标函数取得极小值时, 就得到了参数的估计。

2.3 海冰融池参数估计算法及其验证

为了对融池参数进行参数估计, 需要构建一套有效的参数估计算法。该算法框架设计如下:第一步, 设定参数的初始估计值; 第二步, 进行正向模式积分, 得到融池的基本态, 将其保存,同时保存融池中间强迫数据以驱动融池伴随模式,所谓融池中间强迫数据指作为融池模式输入的5个海冰状态参量(海冰顶层融化速率、积雪融化速率、海冰密集度、海冰体积和海冰表面温度)和降水变量在模式积分中每一时间步的值; 第三步,将正向积分得到的融池覆盖率数据与观测数据根据公式(8)计算目标函数; 第四步, 利用融池中间强迫数据, 进行伴随模式的反向积分, 得到伴随变量的值; 第五步, 计算目标函数的梯度; 第六步, 将计算得到的目标函数值和梯度输入LBFGS极小化算法, 极小化算法据此来自动调整最优步长及更新参数的估计值, 并判断是否达到指定的迭代终止条件, 若达到条件, 则终止程序并返回最终的参数估计值, 反之则从第二步开始利用最新的参数估计值进行新的迭代循环。算法流程如图 3所示。

图3 融池参数估计算法流程图Fig.3.Flow chart of the melt pond parameter estimation algorithm

验证本文参数估计算法的方法和前人工作[39-41]相似, 采取了理想试验的方式, 具体步骤为: 首先选取CICE6.0中的默认值0.8作为预先给定的δp参数值, 进行正向模式控制积分, 得到模拟的融池覆盖率逐小时数据, 将此数据作为融池的“观测”数据保存以供使用。然后设定参数的初始估计值为一随机值0.7, 进行参数估计。检验伴随模式及相应的极小化算法能否得到参数默认值0.8。

参数估计程序迭代了5次后结束, 得到的参数δp估计值为0.8。图4a为参数δp在迭代过程中的变化。可以看到, 在迭代了2次之后参数估计值基本上接近预先设置的值0.8。图4b显示了目标函数J的值在迭代过程中的变化, 可以看出,目标函数在迭代前2步下降较快, 第2步已经接近零。图4c为目标函数梯度范数在迭代过程中的变化, 这里目标函数梯度范数在两次迭代后就变得非常接近于零, 表明目标函数将达到其极值。

图4 参数估计中被估计δp参数 (a), 目标函数 (b) 和目标函数梯度范数 (c) 随迭代次数的变化Fig.4.Change of the δp parameter (a), cost function (b) and norm of cost function gradient (c) with iteration number in parameter estimation process

为了检验初始估计值对参数估计结果的影响,进行了若干次试验, 每次的初始估计值不同, 结果如表1所示。表1中第3列为参数的初始估计值, 从0.1变化到1.0。第4列为参数的最终估计值。从表中可以看出, 在不同初始估计值下, 参数估计算法都能得到“正确的参数估计值”0.8。同时, 迭代次数均在10次以下, 具有较高的计算效率。以上计算验证了伴随模式和参数估计的算法。结果表明, 参数估计算法是有效的。

表1 不同初始估计值下海冰模式融池方案的参数估计结果Table 1.Results of melt pond parameter estimation under different first-guess values

3 北极多年冰、一年冰融池纵横比参数估计

大量融池观测数据统计表明[20,42], 融池覆盖率与融池所在海冰类型有关, 多年冰表面较为粗糙, 受到冰脊等地形限制, 倾向于形成覆盖面积小、深度大的融池。而一年冰表面平整, 易于形成覆盖面积大、深度浅的融池。因此, 多年冰上融池覆盖率较低, 一年冰上融池覆盖率高。本文所优化的参数为融池纵横比参数, 它是融池深度和融池覆盖率的比值, 与海冰类型具有密切关系。为此, 在验证了参数估计算法的可靠性之后,利用MODIS观测数据, 分别对北极多年冰海域和一年冰海域, 选择代表性位置, 进行参数估计。

3.1 研究位置选取

多年冰为历经至少一个夏天仍然存在的海冰。据此, 本研究中利用海冰密集度资料HadISST1, 确定多年冰范围为在3—9月一直被海冰覆盖的海域, 一年冰范围则为3月被海冰覆盖、9月无海冰覆盖的海域。图 5为2005年的3月和9月月平均海冰密集度的二维分布, 从图中可以看出9月边缘海的海冰损失最为严重, 俄罗斯以北的近海海域几乎已无海冰存在。因此, 本文选择的多年冰研究位置位于加拿大北极群岛以北的120°W、80°N, 在图5b中为白色实心正方形所在位置, 一年冰研究位置位于喀拉海的90°E、78°N, 在图5b中为白色空心正方形所在位置。

图5 2005年3月(a)和9月(b)的北极海冰密集度.(b)中白色实心及空心正方形分别为融池纵横比参数估计的多年冰区域(120°W, 80°N)及一年冰区域(90°E, 78°N)Fig.5.Sea ice concentration of March 2005 (a) and September 2005 (b).The white filled square in (b) is the multi-year ice region (120°W, 80°N) and the white hollow square is the first-year ice region (90°E, 78°N) for melt pond aspect ratio parameter estimation

3.2 融池参数估计试验设计

本文主要研究伴随模式在海冰融池参数估计中的应用, 为此设计了以下两组试验。(1)多年冰参数估计试验。使用多年冰位置处2005年的MODIS观测数据约束模式, 先将MODIS数据插值到1 h, 然后将目标函数定义为逐小时的模拟与观测误差平方和。融池纵横比参数的初始估计值设定为0.7。考虑到在融池的发育演化过程中融池覆盖率和融池深度的变化不完全同步[7], 对于融池发展的不同阶段, 模式中应该使用不同的参数。为此, 根据MODIS融池覆盖率资料特点, 将插值后的MODIS数据从5月9日—9月6日, 每隔8天分为一段, 共分为15个时间段, 每个时间段估计一个参数。以多年冰位置处的JRA55-do数据作为强迫场驱动ICEPACK对单点海冰状态进行模拟, 模拟时间为2005年1月1日—5月9日, 积分步长为3600 s, 初始海冰密集度为100%,初始海冰平均厚度为2.79 m, 初始融池覆盖率和深度分别为0%和0 m, 热力学方案选择糊状层方案[43], 短波辐射方案选择Delta-Eddington方案[44],融池参数化方案选择CESM方案, 使用的ITD有5个厚度类别(见表2), 冰层垂向上分为等间距的7层, 雪层为1层, 其他为模式的默认设置。从5月9日开始分段积分, 每个时间段内, 参数估计的流程与2.3中相似, 并直接使用了MODIS观测数据, ICEPACK正向积分8天, 融池伴随模式反向积分8天, 将该段积分得到的海冰-融池状态作为下一时间段正向积分的初始条件。分段积分完成后, ICEPACK接着从2005年9月6日积分到2006年1月1日。(2)一年冰参数估计试验。一年冰试验使用一年冰位置处2005年的MODIS观测数据约束模式, 使用一年冰位置处2005年的JRA55-do数据作为大气强迫, 而初始海冰平均厚度则设定为1.80 m。其余设置与多年冰试验相同。

表2 本文所使用海冰厚度类别Table 2.Ice thickness distribution used in our experiment

3.3 多年冰上融池纵横比参数的估计

多年冰参数估计试验结果如图6所示。从图中可以看出, 该点观测的融池覆盖率存在两个相邻峰值。自5月9日开始, 融池覆盖率逐渐增大,到6月26日, 开始快速增长, 7月20日达到第一个峰值, 约为30%。随后开始下降, 但是在8月初很快又上升到第二个峰值, 约为35%, 之后开始持续的下降过程。使用默认参数0.8的模拟与MODIS观测相差较大, 模拟的融池覆盖率偏小。尤其是在6月底至7月20日之间, 观测值大于20%, 然而模拟的融池覆盖率接近于零, 未能模拟出观测中融池覆盖率的上升, 模拟与8天一次的MODIS观测的均方根误差为16.69 %。使用估计参数模拟时, 对于每一个分段, 融池纵横比参数使用该段的估计值, 初始条件使用上一段的积分终值。使用估计的参数值后, 模拟的融池覆盖率增大, 与MODIS观测更为接近, 均方根误差为7.19%, 与使用默认参数值的结果相比减小了56.93%。由此可见, 使用估计的融池参数, 可以有效地降低融池覆盖率的模拟误差。但是, 并非每个时间段的模拟结果都能得到优化, 比如8月13日之后的时间。这可能是由于这个时间段其他参数对融池的影响更大。比如, 过低的气温导致融水过早减少, 此时仅仅调整融池纵横比这一个参数是不够的。

图6 多年冰区域使用估计的融池参数模拟的融池覆盖率(红色线), 使用默认参数模拟的融池覆盖率(蓝色线)以及样条插值后的MODIS观测(黑色线, 方块标注的是MODIS原始观测)之间的对比Fig.6.Comparison between simulated melt pond fraction using estimated pond parameter (red), simulated melt pond fraction using the default parameter (blue) and the MODIS observation (black, The MODIS pond fraction observation is interpolated using spline interpolation, and the squares mark the original MODIS observations) in MYI region.

多年冰参数估计值有明显的时间变化。融冰期初期, 5月9日—6月2日, 参数值较大, 达到3左右。6月2日—6月26日, 参数值下降到大约1左右。6月26日—8月29日, 参数值进一步下降到0.05以下。融冰期末尾, 8月29日—9月6日的时间段内, 参数值又上升到0.8的默认值。整体来看, 多年冰上估计的融池参数从融冰初期至融冰末期呈下降趋势。这与使用默认参数值模拟的融池覆盖率与观测的差异有关, 参数估计算法不断调整参数值以使得模拟和观测的误差减小。在5月和6月, 使用默认参数值模拟的融池覆盖率略高于观测值, 因而估计的纵横比参数值大于默认值。这是因为融池纵横比增大, 按照公式( 1 )将增加融池深度, 在融水体积不变的情况下, 在进一步积分中将减少融池面积, 起到修正融池覆盖率模拟的效果。在7月和8月, 使用默认参数值模拟的融池覆盖率远低于观测值, 估计的参数值很小。这是因为融池纵横比减小, 按照公式( 1 )将减少融池深度, 在融水体积不变的情况下, 在进一步积分中将增加融池面积, 起到修正融池覆盖率模拟的效果。当然这样的过程只能在其他因子不变的情况下发生。

3.4 一年冰上融池纵横比参数的估计

一年冰参数估计试验结果如图7所示。从图中可以看出, 该点的MODIS融池覆盖率序列存在一个峰值, 在6月上旬增加到15%左右, 随后下降为零。使用默认参数的模拟与MODIS观测相差较大, 模拟的融池覆盖率偏大, 在观测融池覆盖率降低为零后, 模拟融池覆盖率仍然为15%以上, 持续时间过长, 直到9月初仍有融池存在,模拟与观测的均方根误差为6.31%。使用估计的参数值后, 模拟的融池覆盖率明显下降, 与MODIS观测吻合得更好, 均方根误差为4.60%,减小了27.10%。由此可见, 使用估计的融池参数,可以有效地降低融池覆盖率的模拟误差。注意到6月18日之后, 观测融池覆盖率降低为零, 模拟的融池覆盖率未能抓住这种特征。这可能是由于实际中一年冰完全融化, 而模式中海冰仍然存在。使用估计的融池参数, 也不能完全消除这种误差。

图7 一年冰区域使用估计的融池参数模拟的融池覆盖率(红色线), 使用默认参数模拟的融池覆盖率(蓝色线)以及样条插值后的MODIS观测(黑色线, 方块标注的是MODIS原始观测)之间的对比Fig.7.Comparison between simulated melt pond fraction using estimated pond parameter (red), simulated melt pond fraction using the default parameter (blue) and the MODIS observation (black, The MODIS pond fraction observation is interpolated using spline interpolation, and the squares mark the original MODIS observations) in FYI region.

一年冰上估计的参数值较大, 达到10以上,和默认参数值0.8存在较大的差异。一年冰上使用默认参数值模拟的融池覆盖率高于观测值, 因而估计的参数值较大。对于海洋生态模式, 已有研究[45]指出, 相比固定参数以及仅随空间或时间变化的参数, 使用伴随同化方法得到的随空间和时间变化的生态参数更为优越, 符合生态机制。对于本文研究的融池参数, 由于融池过程在时间及空间上的明显变化, 随时间和空间变化的参数可能更为合理。

4 结论与讨论

海冰模式中的融池参数化方案存在某些经验性的参数, 使用伴随模式方法可以对这些参数进行估计以减小不确定性。本文针对CICE6.0海冰模式中的CESM融池参数化方案进行了参数估计。首先, 根据伴随模式代码生成的规则构建了融池方案的伴随模式。接着, 给出了使用伴随模式估计融池参数的算法流程。这一算法包括融池方案、融池方案的伴随模式及L-BFGS极小化算法。利用控制试验数据作为观测约束, 进行了该参数估计算法的验证。最后, 利用MODIS观测资料作为约束条件, 对CESM方案的参数进行了分时间段、分区域的估计。主要结论如下。

1.本文所采用的代码直接转换方法可以正确地得到融池参数化方案的切线性模式和伴随模式。构造的伴随模式可以正确地计算目标函数的梯度。使用正向模式、伴随模式和极小化算法组成的参数估计算法可以用于对CICE6.0的CESM融池参数化方案中的参数进行估计。

2.使用伴随模式方法得到的参数估计值是随时间和空间变化的, 区别于CICE6.0中的默认固定参数值。在北极夏季不同时间段, 多年冰海域和一年冰海域的融池参数估计值具有较大差异。这与融池过程的时空变化有着密切关系。

3.使用得到的最优参数估计值进行模拟, 对选取样本区域融池的模拟误差减小。在两个单点的模拟中, 相比使用默认参数值的试验, 多年冰上模拟融池覆盖率和MODIS观测数据的均方根误差减小了56.93%, 一年冰上均方根误差减小了27.10%。

融池的时空变化受到辐射、气温、降水、海冰类型等多种变量的影响, 具有复杂的机制。本文主要以融池纵横比参数作为估计和优化的目标,研发了融池伴随模式。一年冰上估计的参数大于多年冰上估计的参数。这一点不同于融池在一年冰上更浅的特征。这可能是由于融池伴随模式只考虑了一个参数融池纵横比, 并且在一年冰及多年冰区域都只选择了一个点进行参数估计试验,为了修正一年冰上模拟融池覆盖率偏大的特征,必须增大融池纵横比参数, 以增加融池深度, 进而减小融池覆盖率。这在数学上是合理的, 在物理上却不一定反映了这一点的真实情况。但本研究的目标是以融池为例探讨伴随模式在海冰模式参数估计中的应用, 本文研究结果表明伴随模式方法可以应用于海冰融池参数的估计, 并得到了优化的融池纵横比参数, 为进一步使用伴随模式估计融池参数化方案中的其他参数, 从整体上优化融池和海冰模式的模拟提供了基础。

融池纵横比决定了融水在融池面积和融池深度两者之间的分配, 但只是CESM方案中重要参数之一, 方案中的融水流失比参数则可以直接影响融池体积。由于各物理过程之间的复杂相互作用, 下一步应改进伴随模式, 比如发展二维的海冰伴随模式, 对多个重要参数同时进行估计和调整以达到融池模拟最优的目的。此外, CESM方案的物理过程比较简单, 针对物理过程更完备的方案研发其伴随模式, 并进行参数估计, 将会改善海冰模式中融池过程的描述, 提高海冰模式的准确性。