基于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浓度估计的准确性。