基于GEOS-Chem V12.6.3的全球CO2濃度同化系統的構建

霍 霄, 王茂華, 張欽偉, 魏 崇, 黃永健, 顧倩榮

(1.中國科學院上海高等研究院 碳數據與碳評估研究中心, 上海 201210; 2.中國科學院大學, 北京 100049; 3.中國科學院上海高等研究院 低碳轉化科學與工程重點實驗室, 上海 201210;4.中國科學院腦科學與智能技術卓越創新中心, 上海 200031)

自1750年工業革命開始,大量的CO2因人類對化石能源的大規模使用而被排入大氣,使得大氣中CO2濃度從工業革命前的約278 mL/m3上升到2019年的約410.5 mL/m3[1]。大氣中CO2濃度的上升導致輻射強迫的變化,是引起全球變暖的主要原因[2]。因此,估算出大氣中CO2濃度及其變化情況,對更好地理解全球碳循環機理,研究碳源、碳匯的分布,從而為減緩全球變暖提供科學支持具有重要的意義。

國際社會目前已發射了多顆衛星監測大氣CO2濃度變化,包括SCIAMACHY[3]、GOSAT[4]、GOSAT-2[5]、OCO-2[6]、OCO-3[7]和TanSat[8]。除此之外,如MicroCarb[9]、GeoCARB[10]等其他衛星也將于2022年左右進行發射。衛星觀測覆蓋范圍廣、分布均勻,為大氣中CO2濃度和CO2通量估算研究的開展,提供了扎實的數據基礎。

基于不同的大氣化學傳輸模型和數據同化算法,已開發出多個全球碳同化系統,包括CarbonTracke[11]、PCTM-Adjoint[12]、LMDZ-4DVar[13]、TM5-4DVar[14]和GEOS-Chem Adjoint[15]。同化衛星觀測能較好地提升CO2模擬濃度的準確度和反演通量的不確定度[14,16-21]。Liu等[22]通過對Aqua衛星搭載的AIRS(atmospheric infrared sounder)載荷的XCO2數據進行同化,生成全球6小時間隔的三維大氣CO2濃度場。Basu等[14]采用TM5傳輸模型和4D-Var相結合的方法,在4°×6°的模式分辨率下,通過同化GOSAT觀測對海洋和生態系統碳通量進行估計。Wang等[23]采用PCTM傳輸模型和Bayes綜合反演相結合的方法,在2°×2.5°的模式分辨率下,分別對原位觀測、GOSAT觀測以及融合的原位觀測和GOSAT觀測數據進行了同化,并分析了對應的濃度模擬結果和通量反演結果。模型誤差能顯著影響同化效果,全球碳通量的估算需綜合考慮不同模型的實驗結果[24-26]。Crowell等[27]采用不同的大氣化學模型和同化算法,分別同化了OCO-2衛星觀測和原位觀測,同化結果顯示由這兩種觀測估計的全球總碳匯相近。Basu等[26]采用觀測系統模擬實驗的方法,評估了ACTM、LMDZ等五種全球大氣化學模型對于同化OCO-2偽觀測的影響。

目前更多的研究集中于同化GOSAT衛星數據,并且多同化系統結果相互比較對于合理估算碳通量具有重要的作用。因此本研究以GEOS-Chem V12.6.3版本為前向模式,開發了一個能同化OCO-2衛星XCO2數據的全球大氣CO2濃度四維變分同化系統。首先,借鑒在大氣化學傳輸模式GEOS-Chem[28]V9.0.2基礎上開發的伴隨模型GEOS-Chem Adjoint V35[15]版本的思路和方法開發了觀測算子伴隨模塊、積云對流伴隨模塊、行星邊界層伴隨模塊和平流伴隨模塊,并利用有限差分的方法對這4個伴隨模塊進行了驗證,然后分別進行了針對2018年的模擬實驗和同化OCO-2衛星XCO2數據的實驗,實驗結果與獨立的觀測比較顯示,同化衛星XCO2數據,能顯著提高對全球大氣CO2濃度估計的準確性。

1 同化系統的構建

1.1 GEOS-Chem前向模式

GEOS-Chem[28]是一個全球三維大氣化學傳輸模式。本研究以內部模塊結構如圖1所示的GEOS-Chem V12.6.3版本作為同化系統的前向模式。在同化過程中,系統CO2濃度狀態量的初始值和先驗CO2通量由HEMCO模塊[29]讀入,在由FlexGrid模塊讀入的氣象場的驅動下,先經過平流模塊,然后經過通量添加模塊疊加上由先驗CO2通量引起的CO2濃度變化量之后,再經過行星邊界層和積云對流模塊的向前積分,得到前向模式對全球大氣CO2濃度分布的模擬結果。

圖1 GEOS-Chem V12.6.3版本的內部模塊結構圖

1.2 同化計算過程

在GEOS-Chem前向模式的基礎上,構建了一個全球大氣CO2濃度同化系統,見圖2,采用四維變分的方法對OCO-2衛星的XCO2數據進行同化,圖2的虛線框內是同化系統的伴隨計算過程。同化系統中各模塊的說明見表1。

表1 全球大氣CO2濃度同化系統中主要模塊說明

四維變分的同化過程,是通過求解式(1)表示的代價函數的代價及其梯度,然后將求得的代價函數的代價和梯度輸入進梯度優化算法中去調整對前向模式的系統初始狀態的估計,再通過迭代計算的過程,得到使代價函數取得最小值時的系統初始狀態值,即為同化優化后的結果[31]。四維變分的詳細計算過程說明如下。

J(x0)=Jo(x0)+Jb(x0)

(1)

式中:J(x0)表示系統CO2濃度狀態量的初始值為x0時代價函數的代價,可以分解為由觀測代價Jo(x0)和背景代價Jb(x0)構成。Jo(x0)表示由觀測值yi和系統狀態xi之間的偏差引起的代價,Jb(x0)表示由系統初始狀態x0和背景值xb之間的偏差引起的代價。Jo(x0)和Jb(x0)的表達式為:

(2)

(3)

式中:N為當前同化窗口中同化時刻的總數;H為觀測算子;R為觀測誤差協方差;B為背景誤差協方差。

為提高計算代價函數梯度的效率,采用了伴隨的方法[32],先求出式(2)觀測代價函數的梯度?Jo(x0),然后求出式(3)中背景代價函數的梯度?Jb(x0),最后加和后得到總的代價函數的梯度?J(x0)。

1.2.1觀測代價函數的梯度?Jo(x0)

同化系統中的CO2濃度狀態量xi,表示的是大氣層中每個三維格點處的CO2點濃度。OCO-2衛星的XCO2數據,表示的是從地表到大氣層頂端各大氣分層中處的CO2濃度,經加權平均后得到的柱濃度[33]。兩者的物理意義不同,因此在計算式(2)表示的觀測代價函數的梯度?Jo(x0)時,首先要構建一個觀測算子H,用于將系統CO2濃度狀態量xi,轉換成與衛星觀測yi有相同物理意義的大氣CO2柱濃度。觀測算子H的表達式[19,34]為:

(4)

式(2)中的R是OCO-2衛星的觀測誤差協方差矩陣。為了減少OCO-2觀測數據之間的相關性,將R簡化為對角陣,我們參照Crowell等[36]的方法,先用式(5)～(6)將OCO-2衛星在時間間隔1 s內的XCO2數據和相應的不確定度σ進行平均,得到1 s均值,然后再用式(5)～(6)將時間間隔10 s內的1 s均值再進行平均,得到10 s均值。

(5)

(6)

經式(5)～(6)處理之后,相鄰2個OCO-2衛星XCO2的10 s均值數據之間,在空間上相距約為80 km,可以認為它們之間沒有相關性[36]。因此對于XCO2的10 s均值數據,觀測誤差協方差矩陣R可以簡化成對角陣,按照式(7)～(8)表示的殘差計算方法[37]即可得到R。

(7)

εj=H(xj)-yj-b

(8)

使用伴隨方法求觀測代價函數的梯度?Jo(x0)的過程,可以分解為先求出同化窗口中當前同化時刻tn對應的伴隨項λ(tn),然后求出時刻tn的前一時刻tn-1對應的伴隨項λ(tn-1),再用相同的方法向后遞歸,直到求出初始時刻t0對應的伴隨項λ(t0),λ(t0)即為所求的梯度?Jo(x0)[38-39]。觀測代價函數在當前同化時刻tn的伴隨項λ(tn)的計算公式為

λ(tn)=HTR-1(H(xn)-yn)

(9)

在GEOS-Chem前向模式中,系統狀態量xi從時刻ti遷移到時刻ti+1的過程,可以表示成式(10)[40-41],因此當0≤i

xi+1=Μconv(Μpbl(Μadvc(xi)+Edt))

(10)

(11)

式中:Μconv、Μpbl、Μadvc分別為GEOS-Chem前向模式中的積云對流、行星邊界層和平流模塊;dt為時刻ti和時刻ti+1之間的時間隔,即前向模式的積分時間步長;Edt為在一個積分時間步長內,輸入模式的先驗CO2通量。

分別為同化系統中,與前向模式中的Μconv、Μpbl、Μadvc模塊相對應的伴隨模塊。

式(11)等號右邊的第2項,表示通過前向模式的伴隨過程,將時刻ti+1的伴隨項λ(ti+1)向后積分到時刻ti。利用式(11),一直到向后遞歸到同化周期初始時刻t0時得到的λ(t0),即為觀測代價函數的梯度?Jo(x0)[38]。

1.2.2背景代價函數的梯度?Jb(x0)

式(3)中的B是CO2背景濃度的誤差協方差矩陣,因前向模式采用2°×2.5°較大尺度的網格空間,為簡化背景誤差協方差矩陣B的計算,參照Park等[42]的方法,采用對角陣來近似B,并參照Bannister等[43]的方法,將CO2濃度的方差作為B中對角線處的值,具體的方差值采用Villalobos等[20]的研究結果為4 mL/m3。在得到B之后,將背景代價函數Jb(x0)對初始狀態x0求導,即可得到背景代價的梯度?Jb(x0),具體公式為:

?Jb(x0)=B-1(x0-xb)

(12)

1.2.3CO2濃度狀態更新

在根據式(1)求得代價函數J(x0)的值,和將式(11)～(12)的結果加和得到代價函數的梯度?J(x0)之后,調用擬Newton優化算法LBFGS-B[44-46]程序,對當前同化過程的CO2濃度狀態的初始值進行更新。LBFGS-B算法對初始值進行更新的具體過程為[44]:

x0,new=x0-ηH-1?J(x0)

(13)

式中:x0,new為更新后的CO2濃度狀態初始值;H為代價函數J(x0)對于CO2濃度狀態初始值x0的Hesse矩陣;η為LBFGS-B算法內部的計算步長。

2 實驗與數據

2.1 實驗設計

為檢驗同化系統中伴隨計算過程的準確性,及同化OCO-2衛星數據對提高大氣CO2濃度估計準確性的作用,分別設計了伴隨驗證、模擬和同化3個實驗。

2.1.1伴隨模塊驗證實驗

表2 伴隨驗證實驗的計算方法

2.1.2模擬實驗

模擬實驗僅用來自MERRA-2[47]的氣象數據驅動GEOS-Chem前向模式,模擬全球大氣CO2濃度的變化情況,不對OCO-2衛星的XCO2數據進行同化。模擬實驗的時間范圍是2017年10月至2018年12月,水平分辨率為2°×2.5°,垂直分辨率為47層,積分時間步長為10 min。GEOS-Chem模擬時的初始CO2濃度來自CT2019[48](CarbonTracker 2019)數據集,先驗CO2通量中的人為通量部分來自ODIAC 2019[49],陸地生態通量部分來自CASA[50]模型的GFED4.1s,野火通量和海洋通量來自CT2019中的野火和海洋后驗通量。

2.1.3同化實驗

同化實驗中GEOS-Chem模式的參數設置與模擬實驗完全相同。在同化實驗中,利用構建的同化系統,同化OCO-2衛星天底模式的XCO2數據,得到優化后的全球大氣CO2濃度。同化周期為2017年10月至2018年12月,其中2017年10月至2017年12月為模式的啟動期,同化窗口為6 h,同化積分時間步長為10 min。在將同化結果與地面觀測數據對比做誤差分析時,不包括啟動期時間內的部分。

2.2 實驗數據

2.2.1用于同化的OCO-2衛星的XCO2數據

OCO-2衛星是NASA在2014年發射的一顆大氣CO2濃度觀測衛星,運行在約705 km高度的太陽同步軌道上,觀測的空間分辨率為1.29 km×2.25 km[51-52]。在同化實驗中使用的是OCO-2第9版的XCO2Lite數據產品,先從Lite數據產品中篩選出質量標記為“好”(xco2_quality_flag=0)的天底觀測模式(operation_mode=0)的XCO2數據[35],然后再按照式(5)～(6)進行10 s均值處理,2018年全年數據的處理結果見圖3。

圖3 經質量篩選和10 s均值處理后的2018年OCO-2衛星天底模式的XCO2數據

2.2.2用于驗證的地面觀測數據

用于驗證實驗結果的大氣CO2柱濃度觀測數據來自TCCON的GGG2014[53]數據集。研究中使用了27個在2018年內有觀測數據的TCCON站點的數據。為了分析實驗結果在內陸、海洋和海岸3種不同下墊面下的誤差,將27個TCCON站點分為內陸、海洋和海岸3種類型。這3種類型的分類規則為:1)在大陸上,并且離最近海岸線距離大于100 km為內陸型;2)在島上,并且離最近海岸線距離小于100 km為海洋型;3)在大陸上,并且離最近海岸線距離小于100 km為海岸型。用于驗證實驗結果的地面和航飛觀測數據來自obspack_co2_1_CARBON-TRACKER_CT2019_2020-01-16[54]。研究中使用了來自Obspack的77個地面站點和12個航飛站點的觀測數據。77個地面站點中有53個氣瓶采樣(flask)方式的站點,16個可編程氣瓶采樣(programmable flask packages, PFP)方式的站點和8個原位觀測(in situ)方式的站點。12個航飛站點都是PFP觀測方式。

2.2.3用于參照對比的CT2019B數據

CT2019B[48]是CarbonTracker同化系統[11,55]發布的2019版本全球大氣CO2同化相關的數據集,研究中使用CT2019B數據集中全球3°×2°垂直層數為25層的大氣CO2濃度數據產品,與實驗結果進行參照和對比。

2.3 實驗結果評價指標

采用平均誤差ME(mean error)、平均絕對誤差MAE(mean absolute error)、均方根誤差RMSE(root mean square error)和相關系數CORR(correlation)指標來評價模擬、同化實驗結果和用于驗證的地面觀測數據之間的偏差。這4個評價指標的表達式分別為:

(14)

(15)

(16)

3 實驗結果分析與討論

3.1 伴隨模塊驗證實驗

圖4是用有限差分法,驗證同化系統中觀測算子、積云對流、行星邊界層和平流伴隨模塊計算正確性的實驗結果。圖4(a)～(d)分別對應4個伴隨模塊的驗證結果,橫軸表示用伴隨模塊計算出的梯度,豎軸表示用對應的有限差分法計算出的梯度。圖4中,除觀測算子伴隨模塊驗證結果擬合直線的斜率為0.97,平流伴隨模塊驗證結果擬合直線的R2為0.98之外,其它伴隨模塊驗證結果擬合直線的斜率和R2都為1.00,證明了4個伴隨模塊計算過程的正確性。

圖4 同化系統伴隨模塊驗證實驗結果

3.2 模擬和同化實驗

3.2.1與TCCON的對比

采用Wunch等[52]的方法,將實驗結果、CT2019B與全球27個TCCON站點的觀測數據進行對比。先以TCCON站點為中心,對位于南緯25°以北的TCCON站點,劃出一個經度×緯度為10°×5°的區域,對位于南緯25°以南的TCCON站點,劃出一個經度×緯度為120°×20°的區域,對落在該區域中的模擬、同化實驗的結果和CT2019B的數據分別進行平均,然后再將模擬、同化實驗的平均結果和CT2019B的均值數據轉換成柱濃度數據,與TCCON觀測數據進行對比。表3中列出了模擬、同化實驗的結果與TCCON觀測數據進行對比分析的結果,同時作為參照,也列出了CT2019B的對比結果,圖5(a)和圖6(a)分別展示了全球對比結果和三種TCCON類型對比結果的箱形圖。

表3 模擬、同化實驗、CT2019B與TCCON觀測之間誤差分析結果

箱型內的實線和虛線分別表示中位數和均值,箱型的上下兩條橫邊分別表示上四分位數(Q3)和下四分位數(Q1),箱型上下兩條橫邊之間的距離表示四分位間距(IQR),箱型外的上下兩條短橫線分別表示Q3+1.5×IQR和Q1-1.5×IQR值的位置。

圖6 模擬、同化實驗結果和CT2019B與TCCON和地面觀測對比結果的箱形圖

從表3和圖5(a)可以看出,同化實驗結果和模擬實驗相比,與TCCON觀測的ME和RMSE,分別從0.62 mL/m3和1.42 mL/m3,下降到0.37 mL/m3和1.25 mL/m3,改善了40.32%和11.97%,這與Wang等[56]同化OCO-2衛星數據得到的與TCCON觀測0.80 mL/m3的ME大致相符,也與O’Dell等[48]得出的OCO-2衛星XCO2數據與TCCON之間的ME和RMSE分別為0.30 mL/m3和1.08 mL/m3的結論基本相符,說明同化結果與TCCON之間的偏差,基本達到了OCO-2衛星XCO2數據的程度。

從表3和圖6(a)可以看出,對于海洋和內陸這2種類型的TCCON站點,同化實驗結果比模擬實驗結果都有明顯的改善,ME分別從0.67 mL/m3和1.06 mL/m3降到了0.51 mL/m3和0.74 mL/m3,改善了23.88%和30.19%,RMSE分別從1.33 mL/m3和1.50 mL/m3降到了1.17 mL/m3和1.25 mL/m3,改善了12.03%和16.67%。但是對于海岸類型的TCCON站點,模擬實驗的ME和RMSE分別為-0.74 mL/m3和1.27 mL/m3,同化實驗的ME和RMSE分別為-0.86 mL/m3和1.31 mL/m3,同化實驗的結果并沒有改善。作為參照的CT2019B與同化結果的情況相同,與海岸TCCON站點之間的ME和RMSE分別為-0.94 mL/m3和1.36 mL/m3,明顯大于與海洋和內陸TCCON站點之間的ME和RMSE。

這種現象的原因在于,雖然OCO-2衛星天底模式的XCO2數據,并不包括對海洋上空的觀測,但因為海洋上空的大氣狀態較為一致,所以同化系統利用OCO-2衛星的數據,在對陸地上同化結果進行約束的同時,通過前向模式仍能較好地約束海洋上空的結果[57],然而對于海陸交界的海岸區域,在受復雜多變海陸風影響的同時,OCO-2衛星天底模式的XCO2數據在海岸區域又非常稀疏,因此難以通過同化對海岸區域的結果進行有效的約束。對于這3種類型的TCCON站點,模擬實驗結果和同化實驗結果與TCCON觀測都有一個高的相關性,然而在Reunion、Anmyeondo和Eureka站點相關性是較差的,與之對應的是CT2019B和這3個TCCON站點的觀測也有一個較差的相關性。

3.2.2與地面站和航飛觀測的對比

在與Obspack[54]中的77個地面站和12個航飛觀測數據進行對比時,在時間上選取與觀測數據的觀測時間最接近的實驗結果和CT2019B數據,在空間上將觀測數據插值到最近鄰的GEOS-Chem模式和CT2019B的格點上,然后將該格點處的實驗結果和CT2019B的數據,線性插值到觀測數據對應的氣壓分層后,再進行對比分析。實驗結果、CT2019B數據與地面站觀測數據對比分析的結果見表4,對應的箱形圖見圖5(b)和圖6(b),與航飛觀測數據對比分析的結果見表5,對應的垂直廓線的對比結果見圖5(c)和圖7。

表5 模擬、同化實驗、CT2019B與航飛觀測之間的誤差分析結果

圖7 模擬、同化實驗、CT2019B與12個航飛觀測站點之間垂直廓線的對比

同化實驗結果和模擬實驗相比,總體上有較明顯的改善,其中與地面站和航飛觀測之間的平均誤差ME,分別從0.71 mL/m3和0.93 mL/m3,下降到0.41 mL/m3和0.51 mL/m3,改善了42.25%和45.15%。

從表4和圖6(b)可以看出,對于氣瓶采樣和原位觀測的地面觀測數據,同化實驗結果與模擬實驗相比都有改善,ME分別從0.70 mL/m3和0.74 mL/m3下降到0.40 mL/m3和0.44 mL/m3,改善幅度為42.86%和40.54%,RMSE分別從2.63 mL/m3和4.56 mL/m3下降到2.52 mL/m3和4.44 mL/m3,改善幅度為4.18%和2.63%。對于PFP采樣的觀測數據,與模型實驗相比,同化的效果不明顯,與此相對應的是,CT2019B的數據與PFP采樣的觀測數據之間的ME、MAE和RMSE,也明顯高于與氣瓶采樣和原位觀測之間的對應值。

從表5可以看出,除Brig站點之外,同化實驗結果與其它11個航飛觀測站點觀測數據之間的ME,都要比模擬實驗有明顯的改善,整體上與航飛數據對比的ME從0.93 mL/m3下降到0.51 mL/m3,改善幅度為45.16%,RMSE從2.50 mL/m3下降到2.34 mL/m3,改善幅度為6.40%。從圖7可以看出,與模擬實驗相比,同化實驗的垂直廓線更接近航飛觀測的數據,模擬實驗結果、同化實驗結果和CT2019B在2 km以下的垂直廓線與航飛觀測的偏差,明顯比2 km以上的偏差要大。

3.2.3與CT2019B的對比

當用CT2019B與模擬、同化實驗結果進行對比時,先將CT2019B的廓線濃度轉換成柱濃度,然后插值到GEOS-Chem設置的格點上,插值結果見圖8(a),具有明顯的季節變化。2018年4個季節和全年平均的模擬、同化實驗結果的XCO2與CT2019B XCO2的差異見圖8(b)和(c)。模擬、同化實驗結果的XCO2濃度顯著高于CT2019B的XCO2濃度,經過同化后,同化實驗結果的XCO2在北半球的中高緯度地區與CT2019B的XCO2分布是更接近的。模擬、同化實驗結果與CT2019B的統計指標見表6,其中與模擬實驗相比,同化的結果與CT2019B之間的ME從0.77 mL/m3下降到0.68 mL/m3,改善幅度為11.69%,而RMSE從0.93 mL/m3下降到0.75 mL/m3, 改善幅度為19.35%。

表6 模擬、同化實驗與CT2019B之間的誤差分析結果

圖8 2018年4個季節和年平均的CT2019B XCO2 以及模擬結果、同化結果和CT2019B XCO2差異

3.2.4實驗小結

與模擬實驗結果相比,同化了OCO-2衛星XCO2數據的同化實驗結果與3種觀測數據之間的偏差,都有明顯的改善,其中與TCCON、地面站和航飛觀測之間的平均誤差ME,分別改善了40.32%、42.25%和45.15%。將其與同化實驗結果相比,CT2019B與3種觀測數據之間的偏差更小,這是因為CT2019B本身就同化了Obspack里的數據,而同化實驗與這些觀測數據相互獨立,同化結果中不包含這些觀測數據的信息。

4 結 論

本文基于GEOS-Chem V12.6.3和OCO-2衛星XCO2數據,采用四維變分的方法構建了全球大氣CO2濃度同化系統。并利用有限差分法和獨立的TCCON、地面、航飛觀測對伴隨模塊和濃度同化系統進行驗證,結果見下。

1) 對于觀測算子、積云對流、行星邊界層和平流4個伴隨模塊,有限差分法計算出來的梯度與伴隨法計算出來的梯度一致,證明了伴隨模塊構造的正確性。

2) 在ME指標上,同化實驗結果與TCCON、地面和航飛觀測數據之間的ME分別為0.37 mL/m3、0.41 mL/m3和0.51 mL/m3,與對比的模擬實驗結果相比,分別改善了40.32%、42.25%和45.15%,在RMSE指標上,相比與模擬實驗結果分別改善了11.97%、2.40%和6.40%,表明同化OCO-2衛星數據能明顯提高對大氣CO2濃度估計的準確性。