陸洋 王曉春

(南京信息工程大學(xué)海洋科學(xué)學(xué)院, 南京 210044)

提要 融池是影響北極海冰變化的重要因子。但在目前廣泛使用的海冰模式CICE6.0中, 融池模擬與觀測(cè)存在較大差異。CICE6.0中的CESM融池參數(shù)化方案可以計(jì)算融池覆蓋率與深度, 但其中的部分重要參數(shù)存在一定的經(jīng)驗(yàn)性, 影響了融池模擬。本研究為CICE6.0海冰模式的CESM融池參數(shù)化方案開發(fā)了伴隨模式, 利用CICE6.0海冰模式、融池伴隨模式和L-BFGS極小化算法, 使用一年冰及多年冰區(qū)域的MODIS衛(wèi)星融池覆蓋率觀測(cè)數(shù)據(jù), 對(duì)CESM參數(shù)化方案中的融池縱橫比參數(shù)進(jìn)行了分時(shí)段分區(qū)的參數(shù)估計(jì)。結(jié)果表明: 使用伴隨模式可以有效地對(duì)融池參數(shù)進(jìn)行估計(jì); 得到的融池參數(shù)隨時(shí)空變化, 符合融池過程時(shí)空變化強(qiáng)烈的特征; 估計(jì)的參數(shù)值用于模擬, 減小了融池覆蓋率的模擬誤差, 與MODIS觀測(cè)更為一致。

0 引言

北極夏季, 海冰和冰上積雪開始融化形成融水, 融水在自身重力作用下沿著冰表面流動(dòng), 匯聚于冰表面低洼處形成大小融池。由于較高的反照率, 海冰可以反射大部分的入射短波輻射。海冰融化將會(huì)使得更多的熱量被極區(qū)吸收, 從而促使更多的海冰融化, 稱為“海冰-反照率”正反饋[1]。融池的反照率低于海冰, 融池的形成將降低海冰反照率, 增強(qiáng)正反饋, 起到加速海冰融化的作用。因此, 融池在極區(qū)氣候系統(tǒng)中起到重要作用, 使得其成為近年來氣候科學(xué)家關(guān)注的熱點(diǎn)。隨著北極海冰融化加劇[2], 一年冰比重顯著增加[3-4], 而一年冰平整的表面有利融池形成后在水平方向擴(kuò)展。盡管缺乏長期觀測(cè)數(shù)據(jù), 但一般認(rèn)為, 由于一年冰的增加, 北極的融池覆蓋率正在增加[5]?？梢灶A(yù)料, 融池將對(duì)未來北極地區(qū)的氣候產(chǎn)生更明顯的影響。因此, 融池及其在海冰模式中的參數(shù)化研究具有十分重要的意義。

由于北極地區(qū)自然環(huán)境的限制, 融池的現(xiàn)場(chǎng)觀測(cè)一直比較缺乏。融池覆蓋率(Melt Pond Fraction, MPF)、反照率是觀測(cè)的重點(diǎn)。北極海表面熱平衡收支觀測(cè)計(jì)劃(Surface Heat Budget of the Arctic Ocean, SHEBA)執(zhí)行期間曾對(duì)海冰反照率進(jìn)行了長時(shí)間連續(xù)觀測(cè)[6], 測(cè)量得到的融池反照率隨時(shí)間迅速變化, 變化范圍達(dá)到0.1～0.65, 這與融池自身的發(fā)展演化密切相關(guān)。Eicken等[7]基于觀測(cè)指出季節(jié)尺度融池發(fā)展存在4個(gè)大致的階段: (1)融池形成階段; (2)融水流失階段; (3)融池穩(wěn)定擴(kuò)展階段; (4)融池重新凍結(jié)階段。不同階段的融池行為和控制機(jī)制不同。衛(wèi)星遙感為融池觀測(cè)提供了新的數(shù)據(jù)。2012年, R?sel等[8]利用MODIS(Moderate Resolution Imaging Spectroradiometer)光學(xué)傳感器反射率數(shù)據(jù)反演了融池覆蓋率, 提供了2000—2011年每年5—9月中8天一次的北極融池覆蓋率產(chǎn)品。近年來, 開始利用微波遙感對(duì)融池進(jìn)行觀測(cè)[9-11], 但尚未發(fā)布產(chǎn)品?，F(xiàn)場(chǎng)觀測(cè)和衛(wèi)星觀測(cè)均表明, 融池過程具有明顯的時(shí)空變化特征, 海冰模式中的融池參數(shù)化需要考慮這種特征。

融池在海冰模式中的參數(shù)化是非常重要的研究課題。由于融池在空間及時(shí)間上的明顯變化,根據(jù)冰表面特征直接指定固定的冰面反照率的做法, 無法很好考慮融池, 可能造成較大的模擬偏差[6]。為在模式中對(duì)反照率進(jìn)行符合實(shí)際的計(jì)算,需要發(fā)展融池模式, 將其納入海冰模式中。融池觀測(cè)為融池模式的發(fā)展奠定了基礎(chǔ), 目前已有幾種復(fù)雜程度不同的融池模式用在海冰模式中。早期的融池模式[12-14]為一維模式或者二維模式, 往往與現(xiàn)有氣候系統(tǒng)模式中的海冰模式不匹配, 不能很好地納入海冰模式中。Holland等[15]提出了一個(gè)簡(jiǎn)單的融池參數(shù)化方案, 該方案可以計(jì)算融水體積, 并根據(jù)給定的融池縱橫比(融池深度與覆蓋率之比)分配融水體積。CICE海冰模式(Los Alamos Sea Ice Model)采用了該方案并稱為CESM(Community Earth System Model)方案。Flocco等[16-17]基于海冰模式中的厚度分布函數(shù)(Ice Thickness Distribution, ITD)建立了更為復(fù)雜的融池方案, 考慮了冰面地形對(duì)融水分配的影響,利用海冰厚度來確定冰面地形。在該方案中, 融水產(chǎn)生后首先覆蓋最低的冰面, 然后覆蓋次低的冰面, 逐層累加。此外, 該方案還考慮了由于海冰微孔隙結(jié)構(gòu)所導(dǎo)致的融池水滲流等過程。CICE引入該方案稱為地形方案(Topography Scheme,TOPO)。Hunke等[18]綜合了CESM方案和TOPO方案, 提出平整冰方案(Level Ice Scheme, LVL),利用CICE中已有的平整冰變量, 假定融池只能存在于平整冰上, 并且假定融池體積的變化存在縱橫比, 從而分配融水體積。王傳印等[19]結(jié)合觀測(cè)數(shù)據(jù)對(duì)三種方案在CICE5.0中的表現(xiàn)進(jìn)行了較為詳細(xì)的對(duì)比研究, 認(rèn)為CESM方案的融池凍結(jié)條件更為合理, 但存在多年冰上積雪不能完全融化、融池被低估的問題。

但是, 以上融池參數(shù)化方案均建立在一定的假設(shè)上, 由于對(duì)某些過程缺乏認(rèn)識(shí), 包含了一些不確定性較大的參數(shù), 降低了海冰模式的準(zhǔn)確性和可靠性[20]。對(duì)于這部分不確定參數(shù), 一般可以通過敏感性試驗(yàn)的方法來研究和減小其不確定性。但這種方法有很強(qiáng)的局限性, 一方面氣候系統(tǒng)是復(fù)雜系統(tǒng), 存在大量反饋過程, 另一方面在待研究參數(shù)較多的情況下, 這種方法所需要的計(jì)算代價(jià)較為高昂[21]。伴隨模式(Adjoint Model,ADM)為融池參數(shù)化方案中參數(shù)的估計(jì)提供了另外一種可能性[22-23]。伴隨模式是原始非線性模式的切線性模式(Tangent Linear Model, TLM)的轉(zhuǎn)置, 利用伴隨模式進(jìn)行一次積分可以得到模式結(jié)果對(duì)所有模式變量(包括邊界條件、初始條件和模式參數(shù))的導(dǎo)數(shù), 進(jìn)而可以定量地對(duì)參數(shù)進(jìn)行調(diào)整。伴隨模式已經(jīng)廣泛用于模式參數(shù)估計(jì)和初始場(chǎng)估計(jì)中。比如, Liu等[24]利用MITgcm(Massachusetts Institute of Technology General Circulation Model)海冰-海洋耦合模式及其伴隨模式進(jìn)行了海冰模式中海冰密集度初始場(chǎng)優(yōu)化的研究,得到的海冰密集度初始場(chǎng)和實(shí)際更接近, 說明伴隨模式可以用于對(duì)海冰模式初始場(chǎng)進(jìn)行調(diào)整。Fenty等[22]利用MITgcm及其伴隨模式進(jìn)行了2004年海冰-海洋耦合狀態(tài)估計(jì), 同化了海洋水文資料和海冰密集度資料, 相比未同化資料的試驗(yàn), 對(duì)海冰和海洋狀態(tài)的估計(jì)都更加準(zhǔn)確。Kim等[23]利用CICE模式的伴隨模式對(duì)模式中21個(gè)熱力和動(dòng)力學(xué)參數(shù)同時(shí)進(jìn)行了敏感性分析, 并利用計(jì)算得到的梯度信息對(duì)參數(shù)進(jìn)行了試驗(yàn)性調(diào)整,驗(yàn)證了伴隨方法在海冰參數(shù)估計(jì)和優(yōu)化中的可行性。以上研究表明, 伴隨模式已經(jīng)在海冰模式參數(shù)估計(jì)中得到了廣泛應(yīng)用, 但是這些研究中很少考慮融池參數(shù)。由于融池過程對(duì)于海冰模式的重要性以及融池模式中存在許多不確定參數(shù), 有必要利用伴隨模式考慮重要的融池參數(shù), 并通過優(yōu)化算法對(duì)其進(jìn)行估計(jì), 以減小海冰模式在融池模擬中的誤差。

CESM融池參數(shù)化方案的融池縱橫比參數(shù)是基于SHEBA一年冰觀測(cè)數(shù)據(jù)的經(jīng)驗(yàn)參數(shù)[25]。由于多年冰上融池傾向于在深度方向上增長, 和一年冰融池存在顯著不同, 這個(gè)固定值參數(shù)對(duì)整個(gè)北極海冰的適用性較低。因此, 本文的研究目標(biāo)為通過伴隨模式方法, 估計(jì)海冰模式中CESM融池參數(shù)化方案的融池縱橫比參數(shù), 以改善該融池方案對(duì)融池覆蓋率的模擬結(jié)果。第1節(jié)介紹本文所用的數(shù)據(jù)、CICE6.0海冰模式和CESM融池參數(shù)化方案。第2節(jié)介紹本文所構(gòu)建的融池伴隨模式、L-BFGS極小化算法、伴隨模式參數(shù)估計(jì)算法以及利用控制試驗(yàn)的模擬數(shù)據(jù)對(duì)參數(shù)估計(jì)算法的驗(yàn)證結(jié)果。第3節(jié)針對(duì)多年冰及一年冰海域, 利用MODIS融池觀測(cè)數(shù)據(jù)作為觀測(cè)約束, 對(duì)融池參數(shù)進(jìn)行最優(yōu)估計(jì), 并將其用于模擬, 檢驗(yàn)其對(duì)融池模擬的影響。第4節(jié)為本文結(jié)論與討論。

1 數(shù)據(jù)、CICE6.0海冰模式及融池參數(shù)化方案

1.1 融池覆蓋率與海冰密集度觀測(cè)數(shù)據(jù)

MODIS融池覆蓋率數(shù)據(jù)為德國漢堡大學(xué)研發(fā)的產(chǎn)品[8], 下載自漢堡大學(xué)綜合氣候數(shù)據(jù)中心網(wǎng)站: http://icdc.zmaw.de。該數(shù)據(jù)時(shí)間跨度為2000—2011年, 每年5月9日—9月6日, 時(shí)間分辨率為8天一次?？臻g范圍為180°W～180°E、60°N～90°N, 空間分辨率為12.5km×12.5km。衛(wèi)星數(shù)據(jù)和現(xiàn)場(chǎng)觀測(cè)驗(yàn)證數(shù)據(jù)的均方根誤差(Root Mean Squared Error, RMSE)的范圍在3.8%～11.2%之間。MODIS是光學(xué)傳感器, 因此反演的融池覆蓋率中包含云層所引起的數(shù)據(jù)缺測(cè)。這些缺測(cè)主要發(fā)生在2000年、2001年、2002年和2007年的80°N以北區(qū)域。較小的缺測(cè)區(qū)域已經(jīng)在數(shù)據(jù)發(fā)布前進(jìn)行了插值, 但忽略了面積大于12.5 km2的區(qū)域[26]。

本文使用的海冰密集度數(shù)據(jù)為英國氣象局哈德萊中心(Met Office Hadley Centre)的海冰和海表面溫度數(shù)據(jù)集HadISST1[27]。HadISST1整合了來自歷史觀測(cè)和衛(wèi)星觀測(cè)的海冰數(shù)據(jù), 提供了自1871年開始, 空間分辨率為1°×1°的全球逐月海冰密集度。

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

本文使用的CICE6.0海冰模式由美國Los Alamos國家實(shí)驗(yàn)室研發(fā)[28], 建立在海冰厚度分布函數(shù)[29]的基礎(chǔ)上, 是一個(gè)大尺度熱力-動(dòng)力學(xué)模式,考慮了成脊等次網(wǎng)格尺度機(jī)械過程的影響[30]。模式中使用位移極點(diǎn)網(wǎng)格或者三極網(wǎng)格來處理北極點(diǎn)的奇點(diǎn)問題[31]。ICEPACK是CICE6.0的一個(gè)子模塊[32], 刻畫了海冰熱力學(xué)和生物地球化學(xué)等網(wǎng)格元內(nèi)的垂向物理過程(Column Physics)。ICEPACK通過計(jì)算大氣-海冰邊界層和海洋-海冰邊界層的能量收支來更新冰的溫度, 并以此確定冰的增長或融化。本文使用CICE6.0為ICEPACK提供的單獨(dú)的驅(qū)動(dòng)程序, 獨(dú)立運(yùn)行ICEPACK來模擬單點(diǎn)(垂直柱)的海冰狀態(tài)變化。

海冰模式的大氣強(qiáng)迫數(shù)據(jù)使用基于日本氣象廳(Japan Meteorological Agency, JMA)55年再分析資料JRA-55的海洋-海冰模式大氣強(qiáng)迫數(shù)據(jù)集JRA55-do[33]。本文研究的時(shí)間段為2005年, 所用變量包括10米處的氣溫、10米處的比濕、10米處的風(fēng)矢量、向下短波輻射、向下長波輻射和降水。JRA55-do數(shù)據(jù)的時(shí)間間隔為3 h, 本文中將其線性插值到1 h間隔, 以匹配ICEPACK的1 h的積分步長。

1.3 CESM融池參數(shù)化方案

CESM融池參數(shù)化方案是CICE6.0海冰模式中3個(gè)融池參數(shù)化方案之一, 和文獻(xiàn)[15]中的方案有一些微小的不同之處。它可以計(jì)算融池的覆蓋率與深度, 描述了融池從生成到凍結(jié)的過程。該方案的核心是建立了融池深度和覆蓋率之間的線性關(guān)系, 認(rèn)為深度隨覆蓋率線性增長, 即:

式中,hp表示融池深度,δp為融池縱橫比,ap表示融池覆蓋率。在CESM方案的默認(rèn)設(shè)置中δp為0.8。融池水來自于海冰融水、積雪融水和降雨:

其中為模式中第 個(gè)時(shí)間步的單位面積融池體積(m),ρi和ρs分別為冰和雪的密度(kg·m-3),ρw為融水的密度(kg·m-3),frain為降水率(kg·m-2·s-1),Δt為時(shí)間步長(s), Δhi和Δhs分別為在Δt內(nèi)海冰頂層的融化量(m)和積雪的融化量(m)。部分融水會(huì)從海冰邊緣、冰上裂縫和海豹呼吸洞等流失入海,r為融水進(jìn)入融池的比例:

其中ai為海冰密集度。CICE6.0中rmin和rmax的默認(rèn)值分別為0.15和0.7。而在文獻(xiàn)[15]中, 該比例為:r=0.15+0.7ai。氣溫降低時(shí), 融池逐漸重新凍結(jié), 體積減小:

其中r2為系數(shù)0.01,Tp為參考溫度-2℃,Tsfc為海冰表面溫度(℃)。

2 融池伴隨模式及融池縱橫比參數(shù)估計(jì)算法

2.1 融池方案的伴隨模式及其驗(yàn)證

本文為對(duì)融池模式中的參數(shù)進(jìn)行估計(jì), 研發(fā)了融池模式的伴隨模式。構(gòu)造一個(gè)數(shù)值模式的伴隨模式有兩種路徑: 第一種路徑是基于解析的伴隨模式方程組離散得到伴隨模式, 第二種路徑是從原模式的代碼轉(zhuǎn)換得到。前者只適用于比較簡(jiǎn)單的模式, 后者則可用于復(fù)雜模式。在第二種路徑中, 可以借助自動(dòng)微分(Automatic Differential,AD)軟件工具來自動(dòng)生成伴隨模式代碼[34], 這對(duì)于在設(shè)計(jì)之初考慮了數(shù)據(jù)同化問題并為之預(yù)先優(yōu)化了代碼結(jié)構(gòu)的模式(如MITgcm)是方便的。另外一些模式, 其伴隨模式需要直接轉(zhuǎn)換(人工逐行編寫伴隨模式代碼, 其遵循的編碼基本規(guī)范和自動(dòng)微分軟件是相同的)得到, 如日本氣象廳氣象研究所(Meteorological Research Institute, MRI)的海冰模式伴隨模式[35]。本文所使用的CICE6.0模式并沒有為開發(fā)伴隨模式進(jìn)行特別的代碼優(yōu)化,無法使用自動(dòng)微分工具生成伴隨模式[36]。因此,本文采用了直接轉(zhuǎn)換的方式, 這樣可以在一定程度上保證伴隨模式的正確性。

如1.3節(jié)所述, CESM融池方案在物理和數(shù)值計(jì)算上由4個(gè)部分組成: (1)首先計(jì)算由于冰雪的融化和降水所帶來的融水體積增量; (2)其次, 計(jì)算融水流失比例, 按比例扣除融水得到最終進(jìn)入融池的融水體積增量; (3)根據(jù)氣溫和冰表面溫度,計(jì)算縮減后的融池水體積; (4)根據(jù)融池水體積,通過指定的融池縱橫比參數(shù), 同時(shí)考慮海冰厚度對(duì)融池深度的限制, 計(jì)算融池覆蓋率和深度。根據(jù)這樣的特點(diǎn), 本文開發(fā)伴隨模式的簡(jiǎn)要步驟如下: 首先, 確定融池方案的輸入變量和輸出變量,選擇融池縱橫比δp作為輸入變量, 輸出變量為融池方案中依賴于δp參數(shù)的變量; 其次, 逐行線性化融池方案得到切線性模式; 最后, 將切線性模式算子轉(zhuǎn)置得到伴隨模式算子。本文伴隨模式中的基本態(tài)是重新計(jì)算的。對(duì)于模式中的分支結(jié)構(gòu)語句(比如IF…ELSE…END)的處理方法是: 對(duì)兩個(gè)分支的語句分別進(jìn)行線性化, 分支結(jié)構(gòu)的判斷式若涉及輸入變量, 則用其基本態(tài)代替。對(duì)于融池方案的4個(gè)部分, 分別寫出伴隨代碼, 最終將計(jì)算順序顛倒組成完整的融池伴隨模式。伴隨模式編寫的具體規(guī)范和細(xì)節(jié)可以參考MM5伴隨模式系統(tǒng)的技術(shù)報(bào)告[37]。

得到的切線性模式和伴隨模式必須進(jìn)行正確性檢驗(yàn)才能使用。當(dāng)擾動(dòng)足夠小時(shí), 由正向模式控制積分及其擾動(dòng)積分計(jì)算得到的擾動(dòng)大小應(yīng)該和切線性模式計(jì)算結(jié)果非常接近。因此, 可以利用正向模式, 按下式來檢驗(yàn)切線性模式的正確性:

其中,Q是正向模式,z是正向模式輸入變量向量,α為一尺度系數(shù),h是輸入變量的擾動(dòng),代表從正向模式控制積分及其擾動(dòng)積分的輸出所計(jì)算的擾動(dòng)大小。P是切線性模式,Ph代表將h作為輸入運(yùn)行切線性模式得到的擾動(dòng)大小。在實(shí)際計(jì)算中,h應(yīng)該足夠小, 使得可以忽略其高階項(xiàng),同時(shí)應(yīng)該大于機(jī)器舍入誤差。當(dāng)尺度系數(shù)變小時(shí),正向模式得到的擾動(dòng)值與切線性模式得到的擾動(dòng)值之比R趨向于1。

可以利用切線性模式, 按照下式檢驗(yàn)伴隨模式的正確性:

其中,P是切線性模式,dz是切線性模式的輸入變量向量,PT是伴隨模式。將dz作為輸入, 運(yùn)行切線性模式得到Pdz。將Pdz作為輸入, 運(yùn)行伴隨模式PT得到PTPdz。＜＞表示向量內(nèi)積。若上式在機(jī)器誤差精度內(nèi)成立, 則伴隨模式代碼通過驗(yàn)證。

按照公式( 5 )對(duì)切線性模式進(jìn)行檢驗(yàn)。其中尺度系數(shù)α從1.0減小到0.1, 擾動(dòng)量h保持為0.2。結(jié)果如圖 1所示, 隨著α向0逼近,R逼近1, 說明切線性模式的代碼是正確的。按公式(6)對(duì)伴隨模式進(jìn)行了檢驗(yàn), 其中dz=αh, 尺度系數(shù)α和擾動(dòng)向量 的設(shè)置同切線性模式檢驗(yàn)。伴隨模式檢驗(yàn)結(jié)果如圖 2所示, 隨著α逼近0, 公式(6)兩邊之差的絕對(duì)值逐漸逼近0, 說明伴隨模式的代碼是正確的。

圖1 CESM融池切線性模式檢驗(yàn)結(jié)果.α為擾動(dòng)尺度系數(shù), R(公式5)為正向模式積分得到的擾動(dòng)大小與切線性模式積分得到的擾動(dòng)大小的比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融池伴隨模式檢驗(yàn)結(jié)果.α為擾動(dòng)尺度系數(shù),Δ(公式6兩邊之差的絕對(duì)值) 為切線性模式計(jì)算的擾動(dòng)大小與伴隨模式計(jì)算的擾動(dòng)大小之差的絕對(duì)值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極小化算法及目標(biāo)函數(shù)

正向模式、伴隨模式和極小化算法的結(jié)合可以實(shí)現(xiàn)對(duì)模式參數(shù)的最優(yōu)估計(jì)。本文中所使用的極小化算法為 L-BFGS算法(Limited-memory Broyden-Fletcher-Goldfarb-Shanno Algorithm), 它是一種擬牛頓法, 可以應(yīng)用于大規(guī)模優(yōu)化計(jì)算[38]。算法按照下述公式對(duì)被估計(jì)參數(shù)進(jìn)行迭代:

其中,Xk為被估計(jì)參數(shù)的第k次迭代值,為第k次迭代的被估計(jì)參數(shù)目標(biāo)函數(shù)的梯度, 可以通過伴隨模式給出,ρk為第k次迭代的步長, 由極小化算法自動(dòng)調(diào)整。迭代終止條件為: (1)迭代次數(shù)大于2000次或者(2)達(dá)到指定精度, 滿足其中,為目標(biāo)函數(shù)梯度的范數(shù),ε的值在本文中設(shè)為10-5,為被估計(jì)參數(shù)的范數(shù)。

本文的研究目標(biāo)為估計(jì)融池縱橫比參數(shù), 減小融池覆蓋率模擬與觀測(cè)之間的誤差, 所以本文定義目標(biāo)函數(shù)為:

其中,J為目標(biāo)函數(shù),為第i個(gè)觀測(cè)時(shí)間點(diǎn)MODIS觀測(cè)的融池覆蓋率,為對(duì)應(yīng)第i個(gè)觀測(cè)時(shí)間點(diǎn)模式模擬的融池覆蓋率,N為觀測(cè)次數(shù)。目標(biāo)函數(shù)J代表了模式和觀測(cè)之間的誤差, 當(dāng)目標(biāo)函數(shù)取得極小值時(shí), 就得到了參數(shù)的估計(jì)。

2.3 海冰融池參數(shù)估計(jì)算法及其驗(yàn)證

為了對(duì)融池參數(shù)進(jìn)行參數(shù)估計(jì), 需要構(gòu)建一套有效的參數(shù)估計(jì)算法。該算法框架設(shè)計(jì)如下:第一步, 設(shè)定參數(shù)的初始估計(jì)值; 第二步, 進(jìn)行正向模式積分, 得到融池的基本態(tài), 將其保存,同時(shí)保存融池中間強(qiáng)迫數(shù)據(jù)以驅(qū)動(dòng)融池伴隨模式,所謂融池中間強(qiáng)迫數(shù)據(jù)指作為融池模式輸入的5個(gè)海冰狀態(tài)參量(海冰頂層融化速率、積雪融化速率、海冰密集度、海冰體積和海冰表面溫度)和降水變量在模式積分中每一時(shí)間步的值; 第三步,將正向積分得到的融池覆蓋率數(shù)據(jù)與觀測(cè)數(shù)據(jù)根據(jù)公式(8)計(jì)算目標(biāo)函數(shù); 第四步, 利用融池中間強(qiáng)迫數(shù)據(jù), 進(jìn)行伴隨模式的反向積分, 得到伴隨變量的值; 第五步, 計(jì)算目標(biāo)函數(shù)的梯度; 第六步, 將計(jì)算得到的目標(biāo)函數(shù)值和梯度輸入LBFGS極小化算法, 極小化算法據(jù)此來自動(dòng)調(diào)整最優(yōu)步長及更新參數(shù)的估計(jì)值, 并判斷是否達(dá)到指定的迭代終止條件, 若達(dá)到條件, 則終止程序并返回最終的參數(shù)估計(jì)值, 反之則從第二步開始利用最新的參數(shù)估計(jì)值進(jìn)行新的迭代循環(huán)。算法流程如圖 3所示。

圖3 融池參數(shù)估計(jì)算法流程圖Fig.3.Flow chart of the melt pond parameter estimation algorithm

驗(yàn)證本文參數(shù)估計(jì)算法的方法和前人工作[39-41]相似, 采取了理想試驗(yàn)的方式, 具體步驟為: 首先選取CICE6.0中的默認(rèn)值0.8作為預(yù)先給定的δp參數(shù)值, 進(jìn)行正向模式控制積分, 得到模擬的融池覆蓋率逐小時(shí)數(shù)據(jù), 將此數(shù)據(jù)作為融池的“觀測(cè)”數(shù)據(jù)保存以供使用。然后設(shè)定參數(shù)的初始估計(jì)值為一隨機(jī)值0.7, 進(jìn)行參數(shù)估計(jì)。檢驗(yàn)伴隨模式及相應(yīng)的極小化算法能否得到參數(shù)默認(rèn)值0.8。

參數(shù)估計(jì)程序迭代了5次后結(jié)束, 得到的參數(shù)δp估計(jì)值為0.8。圖4a為參數(shù)δp在迭代過程中的變化。可以看到, 在迭代了2次之后參數(shù)估計(jì)值基本上接近預(yù)先設(shè)置的值0.8。圖4b顯示了目標(biāo)函數(shù)J的值在迭代過程中的變化, 可以看出,目標(biāo)函數(shù)在迭代前2步下降較快, 第2步已經(jīng)接近零。圖4c為目標(biāo)函數(shù)梯度范數(shù)在迭代過程中的變化, 這里目標(biāo)函數(shù)梯度范數(shù)在兩次迭代后就變得非常接近于零, 表明目標(biāo)函數(shù)將達(dá)到其極值。

圖4 參數(shù)估計(jì)中被估計(jì)δp參數(shù) (a), 目標(biāo)函數(shù) (b) 和目標(biāo)函數(shù)梯度范數(shù) (c) 隨迭代次數(shù)的變化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

為了檢驗(yàn)初始估計(jì)值對(duì)參數(shù)估計(jì)結(jié)果的影響,進(jìn)行了若干次試驗(yàn), 每次的初始估計(jì)值不同, 結(jié)果如表1所示。表1中第3列為參數(shù)的初始估計(jì)值, 從0.1變化到1.0。第4列為參數(shù)的最終估計(jì)值。從表中可以看出, 在不同初始估計(jì)值下, 參數(shù)估計(jì)算法都能得到“正確的參數(shù)估計(jì)值”0.8。同時(shí), 迭代次數(shù)均在10次以下, 具有較高的計(jì)算效率。以上計(jì)算驗(yàn)證了伴隨模式和參數(shù)估計(jì)的算法。結(jié)果表明, 參數(shù)估計(jì)算法是有效的。

表1 不同初始估計(jì)值下海冰模式融池方案的參數(shù)估計(jì)結(jié)果Table 1.Results of melt pond parameter estimation under different first-guess values

3 北極多年冰、一年冰融池縱橫比參數(shù)估計(jì)

大量融池觀測(cè)數(shù)據(jù)統(tǒng)計(jì)表明[20,42], 融池覆蓋率與融池所在海冰類型有關(guān), 多年冰表面較為粗糙, 受到冰脊等地形限制, 傾向于形成覆蓋面積小、深度大的融池。而一年冰表面平整, 易于形成覆蓋面積大、深度淺的融池。因此, 多年冰上融池覆蓋率較低, 一年冰上融池覆蓋率高。本文所優(yōu)化的參數(shù)為融池縱橫比參數(shù), 它是融池深度和融池覆蓋率的比值, 與海冰類型具有密切關(guān)系。為此, 在驗(yàn)證了參數(shù)估計(jì)算法的可靠性之后,利用MODIS觀測(cè)數(shù)據(jù), 分別對(duì)北極多年冰海域和一年冰海域, 選擇代表性位置, 進(jìn)行參數(shù)估計(jì)。

3.1 研究位置選取

多年冰為歷經(jīng)至少一個(gè)夏天仍然存在的海冰。據(jù)此, 本研究中利用海冰密集度資料HadISST1, 確定多年冰范圍為在3—9月一直被海冰覆蓋的海域, 一年冰范圍則為3月被海冰覆蓋、9月無海冰覆蓋的海域。圖 5為2005年的3月和9月月平均海冰密集度的二維分布, 從圖中可以看出9月邊緣海的海冰損失最為嚴(yán)重, 俄羅斯以北的近海海域幾乎已無海冰存在。因此, 本文選擇的多年冰研究位置位于加拿大北極群島以北的120°W、80°N, 在圖5b中為白色實(shí)心正方形所在位置, 一年冰研究位置位于喀拉海的90°E、78°N, 在圖5b中為白色空心正方形所在位置。

圖5 2005年3月(a)和9月(b)的北極海冰密集度.(b)中白色實(shí)心及空心正方形分別為融池縱橫比參數(shù)估計(jì)的多年冰區(qū)域(120°W, 80°N)及一年冰區(qū)域(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 融池參數(shù)估計(jì)試驗(yàn)設(shè)計(jì)

本文主要研究伴隨模式在海冰融池參數(shù)估計(jì)中的應(yīng)用, 為此設(shè)計(jì)了以下兩組試驗(yàn)。(1)多年冰參數(shù)估計(jì)試驗(yàn)。使用多年冰位置處2005年的MODIS觀測(cè)數(shù)據(jù)約束模式, 先將MODIS數(shù)據(jù)插值到1 h, 然后將目標(biāo)函數(shù)定義為逐小時(shí)的模擬與觀測(cè)誤差平方和。融池縱橫比參數(shù)的初始估計(jì)值設(shè)定為0.7?？紤]到在融池的發(fā)育演化過程中融池覆蓋率和融池深度的變化不完全同步[7], 對(duì)于融池發(fā)展的不同階段, 模式中應(yīng)該使用不同的參數(shù)。為此, 根據(jù)MODIS融池覆蓋率資料特點(diǎn), 將插值后的MODIS數(shù)據(jù)從5月9日—9月6日, 每隔8天分為一段, 共分為15個(gè)時(shí)間段, 每個(gè)時(shí)間段估計(jì)一個(gè)參數(shù)。以多年冰位置處的JRA55-do數(shù)據(jù)作為強(qiáng)迫場(chǎng)驅(qū)動(dòng)ICEPACK對(duì)單點(diǎn)海冰狀態(tài)進(jìn)行模擬, 模擬時(shí)間為2005年1月1日—5月9日, 積分步長為3600 s, 初始海冰密集度為100%,初始海冰平均厚度為2.79 m, 初始融池覆蓋率和深度分別為0%和0 m, 熱力學(xué)方案選擇糊狀層方案[43], 短波輻射方案選擇Delta-Eddington方案[44],融池參數(shù)化方案選擇CESM方案, 使用的ITD有5個(gè)厚度類別(見表2), 冰層垂向上分為等間距的7層, 雪層為1層, 其他為模式的默認(rèn)設(shè)置。從5月9日開始分段積分, 每個(gè)時(shí)間段內(nèi), 參數(shù)估計(jì)的流程與2.3中相似, 并直接使用了MODIS觀測(cè)數(shù)據(jù), ICEPACK正向積分8天, 融池伴隨模式反向積分8天, 將該段積分得到的海冰-融池狀態(tài)作為下一時(shí)間段正向積分的初始條件。分段積分完成后, ICEPACK接著從2005年9月6日積分到2006年1月1日。(2)一年冰參數(shù)估計(jì)試驗(yàn)。一年冰試驗(yàn)使用一年冰位置處2005年的MODIS觀測(cè)數(shù)據(jù)約束模式, 使用一年冰位置處2005年的JRA55-do數(shù)據(jù)作為大氣強(qiáng)迫, 而初始海冰平均厚度則設(shè)定為1.80 m。其余設(shè)置與多年冰試驗(yàn)相同。

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

3.3 多年冰上融池縱橫比參數(shù)的估計(jì)

多年冰參數(shù)估計(jì)試驗(yàn)結(jié)果如圖6所示。從圖中可以看出, 該點(diǎn)觀測(cè)的融池覆蓋率存在兩個(gè)相鄰峰值。自5月9日開始, 融池覆蓋率逐漸增大,到6月26日, 開始快速增長, 7月20日達(dá)到第一個(gè)峰值, 約為30%。隨后開始下降, 但是在8月初很快又上升到第二個(gè)峰值, 約為35%, 之后開始持續(xù)的下降過程。使用默認(rèn)參數(shù)0.8的模擬與MODIS觀測(cè)相差較大, 模擬的融池覆蓋率偏小。尤其是在6月底至7月20日之間, 觀測(cè)值大于20%, 然而模擬的融池覆蓋率接近于零, 未能模擬出觀測(cè)中融池覆蓋率的上升, 模擬與8天一次的MODIS觀測(cè)的均方根誤差為16.69 %。使用估計(jì)參數(shù)模擬時(shí), 對(duì)于每一個(gè)分段, 融池縱橫比參數(shù)使用該段的估計(jì)值, 初始條件使用上一段的積分終值。使用估計(jì)的參數(shù)值后, 模擬的融池覆蓋率增大, 與MODIS觀測(cè)更為接近, 均方根誤差為7.19%, 與使用默認(rèn)參數(shù)值的結(jié)果相比減小了56.93%。由此可見, 使用估計(jì)的融池參數(shù), 可以有效地降低融池覆蓋率的模擬誤差。但是, 并非每個(gè)時(shí)間段的模擬結(jié)果都能得到優(yōu)化, 比如8月13日之后的時(shí)間。這可能是由于這個(gè)時(shí)間段其他參數(shù)對(duì)融池的影響更大。比如, 過低的氣溫導(dǎo)致融水過早減少, 此時(shí)僅僅調(diào)整融池縱橫比這一個(gè)參數(shù)是不夠的。

圖6 多年冰區(qū)域使用估計(jì)的融池參數(shù)模擬的融池覆蓋率(紅色線), 使用默認(rèn)參數(shù)模擬的融池覆蓋率(藍(lán)色線)以及樣條插值后的MODIS觀測(cè)(黑色線, 方塊標(biāo)注的是MODIS原始觀測(cè))之間的對(duì)比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.

多年冰參數(shù)估計(jì)值有明顯的時(shí)間變化。融冰期初期, 5月9日—6月2日, 參數(shù)值較大, 達(dá)到3左右。6月2日—6月26日, 參數(shù)值下降到大約1左右。6月26日—8月29日, 參數(shù)值進(jìn)一步下降到0.05以下。融冰期末尾, 8月29日—9月6日的時(shí)間段內(nèi), 參數(shù)值又上升到0.8的默認(rèn)值。整體來看, 多年冰上估計(jì)的融池參數(shù)從融冰初期至融冰末期呈下降趨勢(shì)。這與使用默認(rèn)參數(shù)值模擬的融池覆蓋率與觀測(cè)的差異有關(guān), 參數(shù)估計(jì)算法不斷調(diào)整參數(shù)值以使得模擬和觀測(cè)的誤差減小。在5月和6月, 使用默認(rèn)參數(shù)值模擬的融池覆蓋率略高于觀測(cè)值, 因而估計(jì)的縱橫比參數(shù)值大于默認(rèn)值。這是因?yàn)槿诔乜v橫比增大, 按照公式( 1 )將增加融池深度, 在融水體積不變的情況下, 在進(jìn)一步積分中將減少融池面積, 起到修正融池覆蓋率模擬的效果。在7月和8月, 使用默認(rèn)參數(shù)值模擬的融池覆蓋率遠(yuǎn)低于觀測(cè)值, 估計(jì)的參數(shù)值很小。這是因?yàn)槿诔乜v橫比減小, 按照公式( 1 )將減少融池深度, 在融水體積不變的情況下, 在進(jìn)一步積分中將增加融池面積, 起到修正融池覆蓋率模擬的效果。當(dāng)然這樣的過程只能在其他因子不變的情況下發(fā)生。

3.4 一年冰上融池縱橫比參數(shù)的估計(jì)

一年冰參數(shù)估計(jì)試驗(yàn)結(jié)果如圖7所示。從圖中可以看出, 該點(diǎn)的MODIS融池覆蓋率序列存在一個(gè)峰值, 在6月上旬增加到15%左右, 隨后下降為零。使用默認(rèn)參數(shù)的模擬與MODIS觀測(cè)相差較大, 模擬的融池覆蓋率偏大, 在觀測(cè)融池覆蓋率降低為零后, 模擬融池覆蓋率仍然為15%以上, 持續(xù)時(shí)間過長, 直到9月初仍有融池存在,模擬與觀測(cè)的均方根誤差為6.31%。使用估計(jì)的參數(shù)值后, 模擬的融池覆蓋率明顯下降, 與MODIS觀測(cè)吻合得更好, 均方根誤差為4.60%,減小了27.10%。由此可見, 使用估計(jì)的融池參數(shù),可以有效地降低融池覆蓋率的模擬誤差。注意到6月18日之后, 觀測(cè)融池覆蓋率降低為零, 模擬的融池覆蓋率未能抓住這種特征。這可能是由于實(shí)際中一年冰完全融化, 而模式中海冰仍然存在。使用估計(jì)的融池參數(shù), 也不能完全消除這種誤差。

圖7 一年冰區(qū)域使用估計(jì)的融池參數(shù)模擬的融池覆蓋率(紅色線), 使用默認(rèn)參數(shù)模擬的融池覆蓋率(藍(lán)色線)以及樣條插值后的MODIS觀測(cè)(黑色線, 方塊標(biāo)注的是MODIS原始觀測(cè))之間的對(duì)比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.

一年冰上估計(jì)的參數(shù)值較大, 達(dá)到10以上,和默認(rèn)參數(shù)值0.8存在較大的差異。一年冰上使用默認(rèn)參數(shù)值模擬的融池覆蓋率高于觀測(cè)值, 因而估計(jì)的參數(shù)值較大。對(duì)于海洋生態(tài)模式, 已有研究[45]指出, 相比固定參數(shù)以及僅隨空間或時(shí)間變化的參數(shù), 使用伴隨同化方法得到的隨空間和時(shí)間變化的生態(tài)參數(shù)更為優(yōu)越, 符合生態(tài)機(jī)制。對(duì)于本文研究的融池參數(shù), 由于融池過程在時(shí)間及空間上的明顯變化, 隨時(shí)間和空間變化的參數(shù)可能更為合理。

4 結(jié)論與討論

海冰模式中的融池參數(shù)化方案存在某些經(jīng)驗(yàn)性的參數(shù), 使用伴隨模式方法可以對(duì)這些參數(shù)進(jìn)行估計(jì)以減小不確定性。本文針對(duì)CICE6.0海冰模式中的CESM融池參數(shù)化方案進(jìn)行了參數(shù)估計(jì)。首先, 根據(jù)伴隨模式代碼生成的規(guī)則構(gòu)建了融池方案的伴隨模式。接著, 給出了使用伴隨模式估計(jì)融池參數(shù)的算法流程。這一算法包括融池方案、融池方案的伴隨模式及L-BFGS極小化算法。利用控制試驗(yàn)數(shù)據(jù)作為觀測(cè)約束, 進(jìn)行了該參數(shù)估計(jì)算法的驗(yàn)證。最后, 利用MODIS觀測(cè)資料作為約束條件, 對(duì)CESM方案的參數(shù)進(jìn)行了分時(shí)間段、分區(qū)域的估計(jì)。主要結(jié)論如下。

1.本文所采用的代碼直接轉(zhuǎn)換方法可以正確地得到融池參數(shù)化方案的切線性模式和伴隨模式。構(gòu)造的伴隨模式可以正確地計(jì)算目標(biāo)函數(shù)的梯度。使用正向模式、伴隨模式和極小化算法組成的參數(shù)估計(jì)算法可以用于對(duì)CICE6.0的CESM融池參數(shù)化方案中的參數(shù)進(jìn)行估計(jì)。

2.使用伴隨模式方法得到的參數(shù)估計(jì)值是隨時(shí)間和空間變化的, 區(qū)別于CICE6.0中的默認(rèn)固定參數(shù)值。在北極夏季不同時(shí)間段, 多年冰海域和一年冰海域的融池參數(shù)估計(jì)值具有較大差異。這與融池過程的時(shí)空變化有著密切關(guān)系。

3.使用得到的最優(yōu)參數(shù)估計(jì)值進(jìn)行模擬, 對(duì)選取樣本區(qū)域融池的模擬誤差減小。在兩個(gè)單點(diǎn)的模擬中, 相比使用默認(rèn)參數(shù)值的試驗(yàn), 多年冰上模擬融池覆蓋率和MODIS觀測(cè)數(shù)據(jù)的均方根誤差減小了56.93%, 一年冰上均方根誤差減小了27.10%。

融池的時(shí)空變化受到輻射、氣溫、降水、海冰類型等多種變量的影響, 具有復(fù)雜的機(jī)制。本文主要以融池縱橫比參數(shù)作為估計(jì)和優(yōu)化的目標(biāo),研發(fā)了融池伴隨模式。一年冰上估計(jì)的參數(shù)大于多年冰上估計(jì)的參數(shù)。這一點(diǎn)不同于融池在一年冰上更淺的特征。這可能是由于融池伴隨模式只考慮了一個(gè)參數(shù)融池縱橫比, 并且在一年冰及多年冰區(qū)域都只選擇了一個(gè)點(diǎn)進(jìn)行參數(shù)估計(jì)試驗(yàn),為了修正一年冰上模擬融池覆蓋率偏大的特征,必須增大融池縱橫比參數(shù), 以增加融池深度, 進(jìn)而減小融池覆蓋率。這在數(shù)學(xué)上是合理的, 在物理上卻不一定反映了這一點(diǎn)的真實(shí)情況。但本研究的目標(biāo)是以融池為例探討伴隨模式在海冰模式參數(shù)估計(jì)中的應(yīng)用, 本文研究結(jié)果表明伴隨模式方法可以應(yīng)用于海冰融池參數(shù)的估計(jì), 并得到了優(yōu)化的融池縱橫比參數(shù), 為進(jìn)一步使用伴隨模式估計(jì)融池參數(shù)化方案中的其他參數(shù), 從整體上優(yōu)化融池和海冰模式的模擬提供了基礎(chǔ)。

融池縱橫比決定了融水在融池面積和融池深度兩者之間的分配, 但只是CESM方案中重要參數(shù)之一, 方案中的融水流失比參數(shù)則可以直接影響融池體積。由于各物理過程之間的復(fù)雜相互作用, 下一步應(yīng)改進(jìn)伴隨模式, 比如發(fā)展二維的海冰伴隨模式, 對(duì)多個(gè)重要參數(shù)同時(shí)進(jìn)行估計(jì)和調(diào)整以達(dá)到融池模擬最優(yōu)的目的。此外, CESM方案的物理過程比較簡(jiǎn)單, 針對(duì)物理過程更完備的方案研發(fā)其伴隨模式, 并進(jìn)行參數(shù)估計(jì), 將會(huì)改善海冰模式中融池過程的描述, 提高海冰模式的準(zhǔn)確性。