DBPIA-NURIMEDIA

Transcription

1 Atmosphere. Korean Meteorological Society Vol. 22, No. 4 (2012) pp 연구논문 앙상블칼만필터기반탄소추적시스템의아시아지역탄소순환진단에의적용 김진웅 1) 김현미 1), * 조천호 2) 1) 연세대학교대기과학과, 대기예측성및자료동화연구실 2) 국립기상연구소기후연구과 ( 접수 : 2012 년 10 월 12 일, 수정 : 2012 년 11 월 15 일, 게재확정일 : 2012 년 11 월 20 일 ) Application of Carbon Tracking System based on Ensemble Kalman Filter on the Diagnosis of Carbon Cycle in Asia JinWoong Kim 1), Hyun Mee Kim 1), *, and Chun-Ho Cho 2) 1) Atmospheric Predictability and Data Assimilation Laboratory, Department of Atmospheric Sciences, Yonsei University, Seoul , Korea 2) Climate Research Division, National Institute of Meteorological Research, Korea Meteorological Administration, Seoul , Korea (Received: 12 October 2012, Revised: 15 November 2012, Accepted: 20 November 2012) Abstract CO 2 is the most important trace gas related to climate change. Therefore, understanding surface carbon sources and sinks is important when seeking to estimate the impact of CO 2 on the environment and climate. CarbonTracker, developed by NOAA, is an inverse modeling system that estimates surface carbon fluxes using an ensemble Kalman filter with atmospheric CO 2 measurements as a constraint. In this study, to investigate the capability of CarbonTracker as an analysis tool for estimating surface carbon fluxes in Asia, an experiment with a nesting domain centered in Asia is performed. In general, the results show that setting a nesting domain centered in Asia region enables detailed estimations of surface carbon fluxes in Asia. From a rank histogram, the prior ensemble spread verified at observational sites located in Asia is well represented with a relatively flat rank histogram. The posterior flux in the Eurasian Boreal and Eurasian Temperate regions is well analyzed with proper seasonal cycles and amplitudes. On the other hand, in tropical regions of Asia, the posterior flux does not differ greatly from the prior flux due to fewer CO 2 observations. The root mean square error of the model CO 2 calculated by the posterior flux is less than the model CO 2 calculated by the prior flux, implying that CarbonTracker based on the ensemble Kalman filter works appropriately for the Asia region. Key words: Ensemble Kalman filter, CarbonTracker, data assimilation, carbon cycle, inverse modeling 1. 서론 *Corresponding Author: Hyun Mee Kim, Department of Atmospheric Sciences, Yonsei University, 50 Yonsei-ro, Seodaemun-gu, Seoul , Republic of Korea. Phone : , Fax : khm@yonsei.ac.kr 이산화탄소는지구온난화에가장큰영향을주는기체로알려져있으며, 산업혁명이후인류의활동에의해농도가꾸준히증가하고있는추세이다. 지구온난화에미치는이산화탄소의영향을정확히진 415

2 416 앙상블칼만필터기반탄소추적시스템의아시아지역탄소순환진단에의적용 단하기위해서는이산화탄소의흡수및배출량에대한정확한정보가필요하나, 현재이산화탄소의흡수및배출에대한정량적인정보는부족한상황이다. 이산화탄소는대기중에서안정한기체이기때문에지표면에서의탄소순환에의해대기중의양이결정되고, 대기중의이산화탄소농도분포를알기위해선지표면에서의탄소플럭스흡수및방출량을정확하게알고있어야한다. 지표면의탄소플럭스는대기중의이산화탄소관측자료를구속요건으로수송모델의역모델링을통해추정할수있다 (Enting, 2002). 기상수치모델에서는초기조건의불확실성을줄여좀더정확한초기조건을얻기위해서배경값과관측값을적절히조합하는자료동화과정을거친다 (Park and Kim, 2010). 이러한자료동화방법들이화학추적물질의분포와지표면탄소플럭스를추정하기위해서도이용되어왔다. Enting (2002) 은추정하고자하는지표면탄소플럭스와관측자료사이의관계에대한선형방정식을한번에계산하는 Bayesian synthesis 방법을이용하였으나, 이방식은추정하고자하는지표면탄소플럭스의시공간적해상도와사용가능한관측자료가증가함에따라역행렬을계산하는데필요한계산비용이급격히증가하는문제가발생한다. 이러한문제를해결하기위해 4 차원변분자료동화방법 (four dimensional variational data assimilation method, 4D-Var; Baker et al., 2006; Chevallier et al., 2009a, 2009b; Engelen et al., 2009) 를이용하거나앙상블칼만필터 (Ensemble Kalman Filter, EnKF; Peters et al., 2005, 2007, 2010; Feng et al., 2009; Kang et al., 2011) 를이용하여지표면탄소플럭스를추정할수있다. 이런방법들을이용하면, 지표면탄소플럭스를추정할때, 좀더다양한관측자료 (e.g., 위성에서관측된복사값 ) 를사용하거나시공간적해상도가높은관측자료를이용하는것이가능해진다. Evensen (1994) 이 EnKF 의이론적배경을제시한이래로 EnKF 는다양한지구물리학분야에서성공적으로사용되고있다 (Houtekamer and Mitchell, 2001; Anderson, 2001; Bishop et al., 2001; Whitaker and Hamill, 2002; Ott et al., 2004, Zupanski, 2005; Hunt et al., 2007). 미국해양대기청 (National Oceanic and Atmospheric Administration: NOAA) 산하지구시스템연구소 (Earth System Research Laboratory: ESRL) 에서는앙상블제곱근칼만필터 (Ensemble Square Root Kalman Filter, EnSRF; Whitaker and Hamill, 2002) 를기반으로한 SEAT-A (System for Ensemble Assimilation of Tracers in the Atmosphere) 를개발하였으며, 모의된이산화탄소관측자료를이용해고해상도의지표면탄소플럭스를적은계산비용으로추정할수있음을보인바있는데 (Peters et al., 2005), 이는앙상블칼 만필터를지표면탄소플럭스추정에이용한최초의연구이다. NOAA 에서는앞서개발된 SEAT-A 를개선시켜탄소추적시스템 (CarbonTracker) 을개발하였으며, 이시스템과실제지표면관측자료를이용해 2000 년부터 2005 년까지의북미대륙의지표면탄소플럭스에대한연구를수행하였다 (Peters et al., 2007). Peters et al. (2010) 은유럽지역에특별히적용된탄소추적시스템의유럽버전을개발하였으며, 이를이용해 2000 년부터 2007 년까지유럽지역을중심으로한지표면탄소플럭스를분석하였다. 하지만현재까지앙상블칼만필터방법이나 4 차원변분자료동화방법을적용하여아시아지역중심의지표면탄소플럭스를고해상도로추정한연구는보고된바없다. 따라서, 이연구에서는미국해양대기청에서개발한탄소추적시스템 2008 년버전 (CT2008) 을이용하여, 아시아지역을중심으로탄소추적시스템을시험운영하고, 2000 년부터 2006 년까지의결과자료분석을통해탄소추적시스템이아시아지역에서잘동작하는지살펴보려한다. 2 장에서는탄소추적시스템의이론적배경및사례실험을설명하고, 3 장에서는실험운영된이산화탄소농도의앙상블스프레드, 최적화된탄소플럭스, 그리고이에기반한탄소몰농도결과를검증한다. 마지막으로결론및요약을 4 장에제시하였다. 2. 방법 2.1. 앙상블칼만필터탄소추적시스템에서사용하는자료동화방법은 Whitaker and Hamill (2002) 의앙상블제곱근칼만필터 (EnSRF) 이다. 이방법은관측자료를연속적으로동화하는방식으로분석값을계산할때한단계에서하나의관측자료만을이용하며, 여러관측자료를계산하기위해서는관측개수만큼반복계산이필요하다. 이장에서는탄소추적시스템에서사용되는 EnSRF 의알고리즘에대해간략하게설명하며, EnSRF 에대한자세한수식은 Park and Kim (2010) 을참고하기바란다. EnSRF 에서는분석값을구할때평균값과앙상블섭동을따로업데이트하며구하는식은 (2.1), (2.2) 와같이표현된다 : x t a x i 'a = = x t b K y o b + ( H( x t )) x i 'b 'b k H( x i ) (2.1) (2.2) 여기서 x 는상태벡터로서자료동화과정에서계산하고자하는값으로, 탄소추적시스템에서는지표면플 한국기상학회대기제 22 권 4 호 (2012)

3 김진웅 김현미 조천호 417 럭스를직접추정하는대신 scaling factor를사용해간접적으로추정하는데, 이에대해서는다음장에서자세히설명한다. a는분석값, b는초기추정값, 그리고 t는해당시간을나타낸다. 는앙상블평균을, ' 은앙상블섭동을나타내며, y o 는관측값벡터이다. H는관측연산자 (observation operator) 로서상태변수를모델공간에서관측공간으로변환하는역할을한다. 탄소추적시스템에서는 TM5 수송모델이이에해당하며, 이모델에대한자세한설명은 Peters et al. (2004) 와 Krol et al. (2005) 를참고하기바란다. 그리고 K와 k 는각각칼만게인행렬 (Kalman gain matrix) 및감소된칼만게인 (reduced Kalman gain) 이며이를계산하는식은 (2.3), (2.4) 와같다. b T b T 1 K = ( P t H ) ( HP t H + R) k = K α P t b (2.3) (2.4) 여기서는배경오차공분산, R은관측오차공분산이이며, 관측지점별로관측오차공분산값이정해져있다. 그리고 α는칼만게인행렬에곱해지는스칼라값으로하나의관측자료가분석단계에서사용될때마다변하는값으로아래식으로계산된다 : R 1 α = b T HP t H + R (2.5) 이렇게평균값과앙상블섭동을독립적으로계산 a 함으로써분석오차공분산 P t 의과소평가를방지할 b T 수있다 (Whitaker and Hamill, 2002). 이때 P t H 와 b T HP t H 를계산하기위해서는아래의식을이용한다 : PH T ( x' N 1 1, x' 2,, x' N ) ( H( x' 1 ), H( x' 2 ),, H( x' N )) T (2.6) HPH T ( H( x' N 1 1 ), H( x' 2 ),, H( x' N )) ( H( x' 1 ), H( x' 2 ),, H( x' N )) T (2.7) 여기서 N은앙상블의개수이다. 연속적인 EnSRF 자료동화에서하나의관측값을사용해분석값을계산한후, 다음단계의관측값을사용할때이전단계에서업데이트된정보를반영하기위해계산된 CO 2 농도값인 H( x' b i ) 를수송모델을이용해다시계산해야한다. 하나의관측자료를사용할때마다수송모델을이용해이값을다시계산하는것은매우비효율적이기때문에, 다음식 (2.8), (2.9) 를이용해, 각관측자료의위치및시간에대응 b 하는, 모델에의해계산된이산화탄소값인 H( x t ) 의평균값과앙상블섭동값을업데이트한다 : b b H( x t ) m H( x t ) m H m K y o b = + ( H( x t )) H( x' b t ) m = H( x' b t ) m H m k ( y o H( x' b t )) (2.8) (2.9) 여기서 m은자료동화에사용된각관측자료에해당하는값이며, H m 은업데이트를하기위해도입된행렬이다. H m K을계산하기위해서는앞서칼만게인 b T 행렬을계산하기위한식 (2.3) 의 P t H 항을식 (2.7) 을이용해계산하면간단하게구할수있으며, H m k 또한같은방법으로계산할수있다. 이때 H m 을따로구현할필요는없으며, 식 (2.6) 과 (2.7) 을이용해각관측자료를이용한분석시필요한업데이트된이산화탄소농도값을복잡한과정없이계산할수있다. 자료동화과정중한분석시점에서, 식 (2.1) 부터 (2.7) 까지는관측자료의개수만큼계산이필요하고, 식 (2.8) 과 (2.9) 은관측자료개수의제곱만큼의계산이반복적으로이루어진다. 기존의수치예보모델에적용된자료동화방법과는달리탄소추적시스템에서는상태벡터가현재분석하고자하는시점의자료뿐만아니라이전시점에해당하는정보들도포함하고있다. 일반적으로 Bayesian synthesis를이용한방법에서는분석하고자하는상태벡터가모든시점의정보를포함하고있는데, 이는이산화탄소가대기중으로확산되는데시간이걸려관측값과상태벡터사이의관계가오래지속되기때문이다 (Peters et al., 2005). 따라서이러한 지연 을자료동화단계에서고려할필요가있으며, 이것은자료동화과정에서분석되는상태벡터에현재단계의정보뿐만아니라이전단계의정보를포함하는방식으로구현할수있다. 본연구에서사용한지연시간은 Peters et al. (2007) 과 Peters et al. (2010) 에서사용된것과같은 5주이며, 이것은현재관측된이산화탄소농도값이 4주전의지표면플럭스추정값에영향을준다는의미이다. EnSRF에서는한정된앙상블개수를사용할때발생하는표본추출오차 (sampling error) 를줄이고배경오차공분산의과소추정으로인한필터발산 (filter divergence) 을방지하기위해공분산국지화 (covariance localization) 방법이사용된다 (Houtekamer and Mitchell, 2001). 하지만탄소추적시스템에서최적화되는 scaling factor 사이에는물리적인거리를정의할수없기때문에, scaling factor의앙상블섭동으로부터각 scaling factor 사이의상관관계를계산한후, 이값들에대한유의성검증을하여정해진 t-value를넘은 scaling factor에대한칼만게인의성분을 0으로하는방식으 Atmosphere, Vol. 22, No. 4. (2012)

4 418 앙상블칼만필터기반탄소추적시스템의아시아지역탄소순환진단에의적용 로국지화를적용한다. 그리고국지화를적용할때는 Marine Boundary Layer (MBL) 가아닌지역의관측자료에대해서만적용하는데, 이것은 MBL 지역에서관측된자료는먼지역의플럭스정보를포함하고있다고간주하기때문이다 (Peters et al., 2007) 탄소추적시스템 탄소추적시스템운영체계탄소추적시스템은전구에서관측된이산화탄소관측량과상응하는지표면탄소플럭스를계산하기위한시스템으로서크게 4 개의부분으로이루어져있는데, 각각플럭스모듈, 관측자료, 이산화탄소수송계산을위한 TM5, 앙상블칼만필터기반의자료동화시스템이다 (Fig. 1). 4 개의모듈을포함한탄소추적 시스템의알고리즘이 Fig. 2 에제시되어있다. 플럭스모듈은총 4 개의지표플럭스로구성되어있는데, 각각 (1) 생권플럭스 ; Carnegie Ames Stanford Approach (CASA) 모델로계산한월간격플럭스를기반으로하는 3 시간간격의지표순생태계교환량 (net ecosystem exchange, NEE) (van der Werf et al., 2006), (2) Jacobson et al. (2007) 로부터얻어진대기 - 해양플럭스교환자료로구성된 3 시간간격의해양플럭스, (3) Carbon Dioxide Information and Analysis Center (CDIAC) 자료와 EDGAR 자료를이용해계산된, 인류활동으로인해발생한화석연료플럭스, 그리고 (4) CASA 모델을기반으로계산된 Global Fire Emissions Database version 2 (GFEDv2) 의화재플럭스 (van der Werf et al., 2006) 이다. 이러한플럭스자료들은 1 o 1 o 해상도로구성되어있으며, Peters et al. (2007) 에서처럼다음과같은식을통해각모델격자에서의시간단위별플럭스값을계산하는데이용된다 : Fxyt (,, ) = λ r F bio ( xyt,, ) + λ r F ocn ( xyt,, ) + F ff ( xyt,, ) + F fire ( xyt,, ) (2.10) Fig. 1. Schematic of components of the CarbonTracker, which are composed of flux module, observation database, TM5 model, and EnKF data assimilation system. 여기서 F bio, F ocn, F ff, 그리고 F fire 는각각초기추정된생권, 해양, 화석연료및화재플럭스이다. λ r 는생태지역 (ecoregion) 에해당하는 1주일단위의 scaling factor로서육지생태계와해양플럭스에만적용되는변수이다. 탄소추적시스템에서는전구를 240 개의생태지역으로구분하며, 자료동화과정에서는지표플럭스가아닌 scaling factor가대신최적화된다. 이 scaling factor의국지화는위의 2.1절에서논의된 Fig. 2. Schematic of the CarbonTracker algorithm. 한국기상학회대기제 22 권 4 호 (2012)

5 김진웅 김현미 조천호 419 바와같이한다. 탄소추적시스템에서사용되는 scaling factor λ 는특정생태지역에대한 1주일단위의진폭값이다. 이때생태지역의구분은 Olson et al. (1985) 의 29개범주를 19개로축소한것과 TransCom 연구 (Gurney et al., 2002) 에서사용한지역범주에서대륙에해당하는 11개의곱, 그리고 30개의해양범주의합으로이루어져있다. 이때각육지 Transcom 지역범주마다 19개의모든 Olson 지역범주가존재하는것은아니기때문에실제로사용되는 scaling factor는 156개이다. 이러한추가적인구속요건을사용함으로써밀도가낮은관측망에서오는정보를모델영역으로효과적으로전달할수있다 (Peters et al., 2010). 식 (2.10) 에서계산된플럭스자료를이용해, TM5 수송모델 (Krol et al., 2005) 이양방형둥지격자 (two way nested grid) 를사용하여대기중이산화탄소의시공간적분포를계산한다. TM5 모델은 off-line 수송모델로서 ECMWF의현업예보모델의기상장을이용해이산화탄소의수송을계산한다. 또한이과정에서자료동화시계산된이산화탄소의농도를관측자료와비교하기위해관측자료의위치와시간에해당하는모델값이별도로계산된다 탄소추적시스템예보단계앞서 2.1 장에서소개한앙상블제곱근칼만필터는자료동화의분석단계에서사용되는알고리즘이다. 다음분석단계에서사용될배경값을계산하기위해서는수치예보에서처럼예보단계가필요한데, 탄소추적시스템에서는다음과같은식으로다음분석단계에서사용할 scaling factor 의배경값을계산한다 : a a b ( λ λ t 2 + λ t 1 +λ p ) t = (2.11) 여기서는현재단계에서의배경값이며, λ t 2 와 a λ t 1 는 2주전과 1주전에해당하는분석값이다. λ p 는미리정해진고정값으로서현재는 1로정해져있다. 이런가정을사용하는이유는, 분석하고자하는지표면탄소플럭스의수보다관측자료의수가적어서관측자료수와의균형을맞춰주기위한추가적인구속요건이필요하기때문이다. 또한를 1로정 λ t b λ p a 해줌으로써새로운관측정보가없을때에는미리정해진초기추정값으로 scaling factor 가회귀할수있도록한다 (Peters et al., 2007) 탄소추적시스템의관측자료탄소추적시스템에서사용되는관측자료를측정하는지점은전구에흩어져있으며, 북미와유럽지역에상대적으로많은관측지점이존재한다. 본연구에서사용되는관측자료들은 NOAA, Commonwealth Scientific and Industrial Research Organization (CSIRO), Environment Canada (EC), National Center for Atmospheric Research (NCAR), 그리고 Lawrence Berkeley National Laboratory (LBNL) 에서관측한자료들이다 (Masarie et al., 2011). 이산화탄소관측자료들은연속적으로관측된자료와비연속적으로관측된자료들이있는데, 탄소추적시스템에서자료동화에이용될때는 LST 사이관측값의평균을사용한다. 이는현재수송모델의능력으로는정확한이산화탄소의농도를계산할수없으며, 특히밤시간대지표근처의강한이산화탄소농도의연직경도를모의할수없기때문이다 (Masarie et al., 2011) 사례및실험설계 Table 1 에실험의주요설정들이제시되어있다. 시험운영에사용된탄소추적시스템은 CT2008 버전이며, 격자를아시아지역에설정하였다 (Fig. 3). 실험에서급격한해상도변화의영향을줄이기위해해상도가 3 o 2 o 인 domain 2 가 domain 3 밖에설정되었다. 실행기간은 2000 년 1 월 1 일부터 2006 년 12 월 31 일까지이며, 설정된지연값은 5 주이다. Peters et al. (2007) 에따르면북미지역에 1 o 1 o 해상도의둥지격자를설정하였을때플럭스패턴이잘구분되었고, 지연값을 5 주이상으로해서실험을수행하면북미지역의플럭스신호가너무많이퍼져버리고상대적으로작은신호들이수송모델의연직교환 ( 확산및대류 ) 방정식에민감해져서큰수송오차가생길수있다고한다. 아시아지역에서적당한지연값에대한연구는아직까지이루어지지않았으므로, 이연구에서는북미지역연구와마찬가지로 5 주의지연값을사용 Table 1. Experimental framework of running the CarbonTracker. (domain 1) 6 o 4 o Globe Model domain (domain 2) 3 o 2 o 2 o S~70 o N, 30 o ~168 o E (domain 3) 1 o 1 o 14 o ~54 o N, 66 o ~153 o E Period 2000/01/01 ~ 2006/12/31 Weeks of lag 5 The number of ensemble member 150 Atmosphere, Vol. 22, No. 4. (2012)

6 420 앙상블칼만필터기반탄소추적시스템의아시아지역탄소순환진단에의적용 해아직적어, 일반적으로 1 주일에 10 개이하의관측자료가있다. 3. 결과 Fig. 3. Domain 2 and 3 of the experimental domain. 하였다. 앙상블멤버수는 150 개이다. 실험에사용된관측값의개수가 Fig. 4 에나타나있다. 전체적인관측값의개수는시간에따라증가하지만, 플럭스최적화에사용된관측값의개수는 1 주일에약 개정도이다. 최초사이클에해당하는 2000 년 1 월 1 일부터 2000 년 1 월 7 일사이에는특별히관측값이 364 개가사용되었는데, 이는최초의사이클에서최적화에사용된관측값은현재사이클부터 4 주후에해당하는관측값까지를모두포함했기때문이다. 관측자료중일부는최적화에사용되지않으며 ( 항공관측자료 ), 최적화에사용되지않는자료가관측된지역에대해서도 TM5 모델은실험된값을제공하므로그런지역에서모델값을관측값과비교해탄소추적시스템의성능을검증할수있다. 아시아지역에서관측된자료의개수는전체관측개수에비 3.1. Rank histogram Rank histogram (Hamill, 2001) 은앙상블스프레드가적절하게표현되었는지살펴볼수있는척도로사용된다. Rank histogram 에서 adjusted missing rate 을 (Hou et al., 2001; Meng and Zhang, 2008) 살펴봄으로써앙상블스프레드가적절하게편평한지살펴볼수있다. 이값은 rank histogram 이이상적으로편평할때의양극값의합 (e.g., 150 개의앙상블에대해선 1.32%) 과실제로나타난양극값의차이로계산할수있으며, 10% 이하의값을보일때적절하다고간주한다. Fig. 5 는전구지표관측지점에대해서전체분석기간및여름과겨울에대해계산된 rank histogram 이다. 관측오차를고려하지않았을때는양극값이큰 U 모양을보인다 (Figs. 5(a), (c), (e)). 이것은앙상블이작게모의되어관측값이앙상블의표현범위밖에존재한다는것을의미한다. 반면에 Fig. 4. Number of observations used at each cycle in the experiment. Solid line indicates the number of observations used in assimilation and dashed line indicates the total number of observations including aircraft measurements, and thick solid line indicates the number of observations measured in Asia. Observations not used in assimilation (e.g., aircraft measurements) are able to be used to evaluate the performance of the CarbonTracker. Fig. 5. Rank histogram of prior ensembles accumulated for: (a, b) all seasons, (c, d) summer, and (e, f) winter, verified with observations over the globe. The observational error is incorporated into the rank histograms in b, d, and f. The dashed lines denote the ideal value for the rank histograms (e.g., 100/151 = 0.66%). The numbers in the top right of each panel represent the adjusted missing rate. 한국기상학회대기제 22 권 4 호 (2012)

7 김진웅 김현미 조천호 421 Hamill (2001) 이제시한방법에따라각앙상블값에관측오차와동일한잡음을섞어주면양극값이줄어드는모습을보인다 (Figs. 5(b), (d), (f)). 이때관측오차는각관측지점별로미리지정된값을사용하는데, 관측지점의위치에따라 7 개의범주로구분된값이지정되어있다. 이값들은수송모델의각관측지점에대한이산화탄소모의능력의기대치에기반하여주관적으로정한것이다 (Peters et al., 2005). 관측오차를고려함으로써앙상블스프레드가모두좀더평평한값을보이고, adjusted missing rate 도 10% 이하로적은값을보인다. Rank histogram 에서관측오차를고려하는이유는, 앙상블예측과관측이비슷한지를나타내는또다른척도인 innovation covariance 가관측오차가예보오차와독립적이라는가정하에서는배경오차공분산과관측오차공분산의합과유사하기때문이다 (Dee, 1995; Houtekamer et al., 2005; Whitaker et al., 2008). 즉앙상블예측이예보오차를포함하고있는것처럼관측도참값이아니며관측오차를포함하고있기때문에, 그두오차를포함하여앙상블예측과관측을비교하는것이앙상블스프레드가적정한지판단하는데필요하다. Fig. 6. Rank histogram of prior ensembles accumulated for: (a, b) all seasons, (c, d) summer, and (e, f) winter, verified with observations over Asia. The observational error is incorporated into the rank histograms in b, d, and f. The dashed lines denote the ideal value for the rank histograms (e.g., 100/151 = 0.66%). The numbers in the top right of each panel represent the adjusted missing rate. 계절별로앙상블스프레드의분포경향을살펴보면여름철에는앙상블스프레드가양의바이어스를보인다. 이는여름철에이산화탄소농도가관측값보다크게모의되는경향이있기때문이며, 초기추정플럭스가북반구여름철지표면탄소플럭스의흡수량을과소평가해서수송모델에서대기중이산화탄소농도를실제보다크게모의했기때문으로보인다. 반면에겨울철에는돔형태의스프레드가나타났는데, 이것은관측보다앙상블스프레드가크게표현됐다는것을나타낸다. 실험의목적이아시아지역중심의탄소추적시스템시험운영이기때문에, 아시아지역의관측지점에대해서 rank histogram 을살펴보았다 (Fig. 6). Fig. 5 와마찬가지로관측오차를고려하지않았을때에는여전히큰 U 모양의스프레드가나타난다 (Figs. 6(a), (c), (e)). 그렇지만관측오차를고려했을경우전구에대해서계산하였을때보다 adjust missing rate 이더크게줄어 10% 이하의값을나타낸다 (Figs. 6(b), (d), (f)). 따라서이지역에대한앙상블스프레드는전구보다더평평하게모의되었다고할수있다. 또한전구관측시점에대해서계산했을때보다양극값을제외하고도더평평한스프레드를보이고있다. 따라서, 아시아지역의관측지점에대해서는탄소추적시스템이더적절한앙상블스프레드를모의했다고할수있다 최적화된지면탄소플럭스 Table 2 는최적화된지표면탄소플럭스를 TransCom 지역구분중, Eurasian Boreal (EB), Eurasian Temperate (ET), 그리고 Tropical Asia (TA) 지역별로합산하여연도별로표현한평균값과표준편차값이다. 이값은최적화된식생플럭스뿐만아니라화석연료및화재플럭스값도포함된것이다. 다른지역에비해 ET 지역은시간에따라플럭스배출량이증가하는추세를보이는데, 이는이지역에많은인구가살고있고주요공업국이위치하여화석연료의사용량이계속증가하기때문이다. 반면에 EB 지역은화석연료에의한탄소배출량이적고침엽수림의식생작용에의해지표면에서탄소의흡수가이루어져평균지표면플럭스값이음의값을나타낸다. TA 지역은전반적으로방출이이루어지는것으로나타났다. 연도별로살펴보면 2003 년에 EB 지역의흡수량이제일적었으며 2005 년에 ET 지역에서방출량이제일컸다 년에는유럽지역에큰가뭄이발생하여여름철식생에의한탄소의흡수량이크게줄었는데 (Cias et al., 2005; Peters et al., 2010), 이러한가뭄의영향이 EB 지역의식생플럭스에도영향을미쳤기때문인것으로보인다. Atmosphere, Vol. 22, No. 4. (2012)

8 422 앙상블칼만필터기반탄소추적시스템의아시아지역탄소순환진단에의적용 Table 2. Annual flux results [PgC/yr] on TransCom regions with variance. Region Eurasian Boreal -0.97± ± ± ± ± ± ±0.10 Eurasian Temperate 1.67± ± ± ± ± ± ±0.067 Tropical Asia 0.36± ± ± ± ± ± ±0.027 Fig. 7 은 Table 2 에서나타낸지역과동일한지역에대한최적화된플럭스의시계열자료이다. EB 와 ET 지역은중위도지역에위치한지역이기때문에초기추정플럭스에서계절패턴이모사되어있으며, 최적화이후에도이패턴이잘유지되면서흡수및방출량의크기변화만있었다 (Figs. 7(a), (b)) 년부터 2005 년까지겨울철을제외하고최적화된플럭스값은초기추정플럭스값보다낮은값을보이고있으며, 이것은겨울철에는초기추정보다플럭스방출 Fig. 7. Time series of flux [PgC/yr] on TransCom region of (a) Eurasian Boreal (EB), (b) Eurasian Temperate (ET), and (c) Tropical Asia (TA). Black line represents the prior flux, red line represents the posterior flux, blue shaded represents the uncertainty of flux. 이더이루어지고여름철에는초기추정보다플럭스흡수가더이루어진다는것을의미한다. EB 지역에서는이러한패턴이잘나타나있으며 (Fig. 7(a)), ET 지역도인위적인배출량의증가추세를제외한다면이러한패턴이잘나타난다 (Fig. 7(b)). 따라서여름철에는초기추정플럭스가이들지역의식생에의한플럭스흡수량을과소추정하였음을알수있다. EB 지역에서는봄철에초기추정플럭스의피크값이보이는데, 최적화된후에는이피크값이제거되는것을확인할수있으며, 여름철에는특정시점에서초기추정플럭스의두배가넘는최적화된플럭스양이계산되기도한다 (Fig. 7(a)). 따라서 EB 지역에서관측자료의자료동화과정을통해최적화된플럭스는 bottom up 방식으로생성된초기추정플럭스와크게다를수있음을알수있다. ET 지역은 Table 2 에서나타난것처럼시간에따라플럭스값이증가하는추세가잘나타나있다 (Fig. 7(b)). 반면에 TA 지역에서는중위도지역과는달리뚜렷한계절패턴이나타나지않고있으며, 초기추정플럭스와최적화된플럭스의차이가크게나타나지않았다 (Fig. 7(c)). 이는관측자료가부족하여적도지방의정확한초기추정플럭스의산정이아직까지이루어지지않기때문이다. 그리고이지역의지표면탄소플럭스는중위도지역의지표면탄소플럭스방출량과의연관성도낮아, 중위도지역관측자료를통한지표면탄소플럭스추정도어려운상태이다 (Gurney et al., 2002). 이로인해 Peters et al. (2007) 에따르면적도지방의지표면탄소플럭스는인버젼연구별로계산된플럭스값의차이가크게나타난다. 따라서 TA 지역의정확한지표면탄소플럭스계산을위해선추가적인관측과연구가필요할것으로사료된다. Fig. 8 은 2000 년의초기추정및최적화된플럭스평균값의공간적인분포를나타낸다. 그림에표시된지역은수송모델의둥지격자지역에해당한다. 그림에서살펴볼수있듯이 EB 지역에서는최적화된플럭스의흡수량이초기추정플럭스보다크게나타나고있다. 반면에 ET 지역중한국과중국및일본지역에서는큰양의플럭스값을보이고있는데, 이는이지역들에서화석연료에의한탄소방출량이많기때문이다. 식 (2.10) 에서따라, 식생및해양플럭스만최적화되기때문에화석연료및화재플럭스의값은 한국기상학회대기제 22 권 4 호 (2012)

9 김진웅 김현미 조천호 423 Fig. 8. (a) Prior and (b) posterior flux (gcm 2yr 1) with 1o 1o horizontal resolution in These fluxes include biosphere, ocean, fossil fuel, and fires fluxes. Fig. 9. Correlation coefficients of (a) prior fluxes and (b) posterior fluxes, between Eurasian Boreal and other TransCom region. 최적화에 의한 변화가 없다. 최적화된 이후 EB 지역 에서 흡수가 많이 일어나지만 동북아 인구 밀집 지역 에서 화석 연료에 의한 탄소 배출이 많아 아시아 지 역의 총 탄소 수지는 방출 값이 된다 (Table 2). TA 지역에서는 Fig. 7(c)에서 초기 추정 및 최적화된 플 럭스의 시간적 차이가 크지 않았던 것처럼 공간적 차 이도 크게 나타나지 않았다. Fig. 9는 최적화되기 전 후의 EB 지역의 플럭스 시 계열과 다른 Transcom 지역 구분의 플럭스 시계열과의 상관관계를 2000년부터 2006년까지 나타낸 것이다. 전 체적으로 EB와 ET 지역의 플럭스 패턴이 서로 유사 했지만, 다른 지역과의 상관관계에는 차이가 발생했 다 (Figs. 9(a), (b)). 초기 추정 플럭스 값은 EB 지역 과 특정 중위도 지역 (North American Boreal, North American Temperate, Europe) 간의 상관관계가 높았 으며, 그 외 다른 지역과의 상관관계는 낮거나 반대 로 나타났다 (Fig. 9(a)). 특히 열대 지방에서 상관관 계가 낮게 나타나는데, 이것은 Fig. 7(c)에 보인 바와 같이 열대 지역의 플럭스가 중위도 지역과는 달리 계 절별 패턴이 존재하지 않기 때문이다. 최적화된 플럭 스로 계산한 결과를 살펴보면, 여전히 북반구 중위도 지역과의 상관 관계가 높게 나타나고 있으며, 그 외 나머지 지역들과는 상관관계가 더 낮아졌음을 확인할 수 있다 (Fig. 9(b)). Fig. 10은 ET 지역의 플럭스와 다른 지역의 플럭스 간의 상관관계를 나타낸 것이다. 전반적으로 살펴볼 때 EB 지역과 마찬가지로 ET 지역도 다른 북반구 중 위도 지역과 상관 관계가 높게 나타났다. 그러나 EB 지역과 달리 ET 지역의 플럭스는 남반구 Temperate 지 역들 (South Pacific Temperate, South Atlantic Temperate, South Indian Temperate)의 플럭스 및 Northern Ocean지역의 플럭스와의 상관관계가 좀 더 높게 나 타났으며, South American Temperate와 Southern Africa 지역의 플럭스와는 음의 상관관계가 크게 나타 났다. ET 지역의 초기 추정 플럭스에서는 North Atlantic Temperate, Atlantic Tropics, Southern Ocean, 그리고 Indian Tropical 지역의 플럭스와 큰 음의 상관 관계가 나타났으나 (Fig. 10(a)), 최적화된 플럭스와의 상관관계는 크게 줄어든 모습이 나타났다 (Fig. 10(b)) 모의된 이산화탄소 농도 이 장에서는 실험에서 계산된 관측 지점별 이산화 Atmosphere, Vol. 22, No. 4. (2012)

10 424 앙상블 칼만 필터 기반 탄소추적시스템의 아시아 지역 탄소 순환 진단에의 적용 Fig. 10. Correlation coefficients of (a) prior fluxes and (b) posterior fluxes, between Eurasian Temperate and other TransCom region. 탄소 몰 농도를 관측 값과 비교한 결과를 살펴본다. Fig. 11은 초기 추정 플럭스로 계산된 이산화탄소의 RMSE와 분석된 플럭스로 계산된 이산화탄소의 RMSE 값을 나타낸다. 이때 분석된 플럭스로 계산된 이산화 탄소 자료의 RMSE는 자료동화에서 사용된 관측자료 와 자료동화에서 사용되지 않은 항공 관측자료에 대 해 따로 계산되었다. 전반적으로 자료동화과정을 통해 최적화된 플럭스 로 계산된 이산화탄소의 RMSE 값이 초기 추정 플럭 스로 계산된 이산화탄소의 RMSE 값보다 더 작았다. 이로서 자료동화과정을 통해 좀 더 정확한 지표면 탄 소 플럭스 추정이 이루어졌다고 할 수 있다. 계절별 로 결과를 살펴보면 여름철이 겨울철보다 RMSE 값 이 크다. 아시아 지역에 대한 관측자료만을 살펴 보 면 이 지역에서는 전구에 대해 계산했을 때 보다 더 큰 RMSE 값이 나타났으며, 이는 초기 추정 플럭스 및 최적화된 플럭스로 계산된 이산화탄소 값 모두에 해당한다 (Fig. 11(b)). 이 결과는 이 지역의 초기 추 정 플럭스 인벤토리의 정확성을 향상시킬 필요가 있 다는 것을 나타내며, 이것은 정확한 식생 플럭스 산 정 및 현실적인 화석 연료 사용량 계산을 통해 이루 한국기상학회 대기 제22권 4호 (2012) Fig. 11. Times series of RMSE calculated over (a) global observational site and (b) Asia observational site. Black line represents observation minus model CO2 by the prior scaling factor, red line represents surface observation minus model CO2 transported by TM5 after applying the posterior scaling factor, and green line represents aircraft observation minus model CO2 transported by TM5 after applying the posterior scaling factor. 어질 수 있다. 또한 아시아 지역의 관측자료에 대해 서는 RMSE의 계절별 패턴이 전구에 대해 계산했을 때보다 약하게 나타났는데, 이는 부족한 샘플 수의 영 향일 수 있어 이 지역의 관측자료가 좀 더 필요하다 는 것을 의미한다. Peters et al. (2007) 이 지적한 바와 같이 수송 모델 이 가지고 있는 연직 수송 능력의 한계로 인해 대류 권 자유 대기에서 모의된 이산화 탄소의 농도에 오차 및 바이어스가 발생할 수 있다. Fig. 11(a)에서 분석된 플럭스를 이용해 TM5 모델로 계산된 이산화탄소와 항공 관측자료간의 RMSE는 계산된 이산화탄소와 지 표 관측자료간의 RMSE 보다 대부분 더 적고, 그 두 RMSE 간의 계절별 패턴도 일치하는 모습을 보여주 고 있기 때문에, 본 실험에서 이산화탄소의 연직 모 의는 적절하게 된 것으로 보여진다. 4. 결론 및 요약 이 연구에서는 탄소의 흡수 및 배출량 산정을 위해 개발된 앙상블 칼만 필터 자료동화 방법을 이용한 탄 소추적시스템의 자료동화 알고리즘, 운영 체계, 및 예 보 단계에 대해 소개하고, 아시아 지역을 중심으로 탄

11 김진웅 김현미 조천호 425 소추적시스템을시험운영하였으며, 실험결과를분석하였다. 미국 NOAA ESRL 에서개발한탄소추적시스템은 4 개의플럭스모듈, 이산화탄소몰비율관측값, TM5 모델, 앙상블칼만필터자료동화시스템으로구성되어있다. 각플럭스모듈에서미리정해진초기추정값에 scaling factor 를곱하여이산화탄소방출및흡수량을정하고이값을 TM5 모델을이용해수송시켜전구에대한이산화탄소몰비율을얻는다. 이때의모델값은관측이이루어진장소와시간에대응된값이며, 이값과관측값을이용해앙상블제곱근칼만필터방법으로 scaling factor 를최적화한다. 이때공분산지역화를통해분석하려는 scaling factor 와상관관계가낮은관측값은 scaling factor 의최적화에사용하지않는다. 최적화된규모인자를초기추정플럭스에곱하면분석된플럭스가되고, 이값을이용해 TM5 모델을실행하면분석된이산화탄소모델값을얻을수있다. 시험운영을위해 2000 년부터 2006 년까지아시아지역을중심으로하는둥지격자를설정하여탄소추적시스템을이용한실험을수행하였고, 초기추정앙상블의 rank histogram, 지역별플럭스의연평균값과시계열, EB 및 ET 지역과다른 TransCom 지역간의상관관계, 모의된이산화탄소농도의 RMSE 분석을하였다. 그결과 rank histogram 의 adjusted missing rate 이전구보다는아시아지역만고려하였을때더적어, 아시아지역에대해앙상블스프레드가더적절하게표현됐다는것을알수있었다. 계절별로살펴보았을때는여름철에양의 bias 가나타났는데, 이는초기추정플럭스가식생에의한흡수량을과소평가했기때문이다. 반면에겨울철에는관측값의분포에비해앙상블스프레드가크게모의된것으로보인다. 최적화된플럭스결과를살펴보면 EB 지역에서여름철에이산화탄소의흡수량이크게나타났고이로인해연평균플럭스값은음의값을나타나낸다는것을알수있었으며, ET 지역에서도여름철에식생에의한이산화탄소의흡수가컸지만화석연료방출량에의해그효과가상쇄되었고, 연평균으로따지면전체적으로방출이더크다는것을알수있었다. 반면에 TA 지역은계절별변화가뚜렷하지않았으며중위도지역에위치한 EB 및 ET 지역과의초기추정및최적화된지표면플럭스패턴과의상관관계도낮게나타나, 이지역의정확한플럭스추정을위해서는추가적인관측및연구가필요하다는것을알수있었다 년에대해초기추정및최적화된플럭스평균값의공간적인분포를살펴보면, 최적화된이후 EB 지역에서흡수가더많아짐을알수있었다. 실험된이산화탄소농도의 RMSE 를살펴보았을때 탄소추적시스템의자료동화과정을통해더적은오차를갖는이산화탄소몰농도를계산할수있었고, 이것은탄소추적시스템이관측자료와좀더일치하는지표면탄소플럭스를추정할수있다는것을의미한다. 계절별로살펴보았을때겨울철보다여름철에 RMSE 값이더컸으며, 시간에따라 RMSE 가증가하는추세를보였다. 아시아지역만살펴보면이지역에서도자료동화과정을통해좀더정확한플럭스값을추정할수있었으나앞으로추가적인관측자료가더필요할것으로판단된다. 탄소추적시스템분석과아시아지역을중심으로하는시험운영을통해앞으로한반도를포함한아시아지역뿐만아니라전구에대한이산화탄소배출및흡수량을연구할수있는기반이마련되었다. 앞으로탄소추적시스템을아시아지역에적합하게개선하기위한추가적인연구가수행될예정이다. 감사의글 이연구는기상청기후변화감시 예측및국가정책지원강화사업 (CATER ) 의지원으로수행되었습니다. 참고문헌 Anderson, J. L., 2001: An ensemble adjustment Kalman filter for data assimilation. Mon. Wea. Rev., 129, Baker, D. F., S. C. Doney, and D. S. Schimel, 2006: Variational data assimilation for atmospheric CO 2. Tellus, Ser. B, 58, Bishop, C. H., B. J. Etherton, and S. J. Majumdar, 2001: Adaptive sampling with the ensemble transform Kalman filter. Part I: Theoretical aspects. Mon. Wea. Rev., 129, Chevallier, F., R. J. Engelen, C. Carouge, T. J. Conway, P. Peylin, C. Pickett-Heaps, M. Ramonet, P. J. Rayner, and I. Xueref-Remy, 2009a: AIRS-based versus flaskbased estimation of carbon surface fluxes. J. Geophys. Res., 114, D20303, doi: /2009jd , F., S. Maksyutov, P. Bousquet, F.-M. Bréon, R. Saito, Y. Yoshida, and T. Yokota, 2009b: On the accuracy of the CO 2 surface fluxes to be estimated from the GOSAT observations. Geophys. Res. Lett., 36, L19807, doi: /2009gl Ciais, P. and Coauthors, 2005: Europe-wide reduction in primary productivity caused by the heat and drought in Atmosphere, Vol. 22, No. 4. (2012)

12 426 앙상블칼만필터기반탄소추적시스템의아시아지역탄소순환진단에의적용 Nature, 437, Dee, D. P., 1995: On-line estimation of error covariance parameters for atmospheric data assimilation. Mon. Wea. Rev., 123, Engelen, R. J., S. Serrar, and F. Chevallier, 2009: Fourdimensional data assimilation of atmospheric CO 2 using AIRS observations. J. Geophys. Res., 114, D03303, doi: /2008jd Enting, I., 2002: Inverse Problems in Atmospheric Constituent Transport. Cambridege Univ. Press, New York, doi: /cbo Evensen, G., 1994: Sequential data assimilation with a nonlinear quasi-geostrophic model using Monte Carlo methods to forecast error statistics. J. Geophys. Res., 99, Feng, L., P. I. Palmer, H. Bosch, and S. Dance, 2009: Estimating surface CO 2 fluxes from space-born CO 2 dry air mole fraction observations using an ensemble Kalman Filter. Atmos. Chem. Phys., 9, Gurney, K. R. and Coauthors, 2002: Towards robust regional estimates of CO 2 sources and sinks using atmospheric transport models. Nature, 415, Hamill, T. M. 2001: Interpretation of rank histograms for verifying ensemble forecasts. Mon. Wea. Rev., 129, Hou, D., E. Kalnay, and K. K. Droegemeier, 2001: Objective verification of the SAMEX'98 ensemble forecasts. Mon. Wea. Rev., 128, Houtekamer, P. L. and H. L. Mitchell, 2001: A sequential ensemble Kalman filter for atmospheric data assimilation. Mon. Wea. Rev., 129, , H. L. Mitchell, G. Pellerin, M. Buehner, M. Charron, L. Speak, and B. Hansen, 2005: Atmospheric data assimilation with an ensemble Kalman filter: Results with real observations. Mon. Wea. Rev., 133, Hunt, B. R., E. Kostelich, and I. Szunyogh, 2007: Efficient data assimilation for spatiotemporal chaos: A local ensemble transform Kalman filter. Physica D, 230, Kang, J.-S., E. Kalnay, J. Liu, I. Fung, T. Miyoshi, and K. Ide, 2011: Variable localization in an ensemble Kalman filter: Application to the carbon cycle data assimilation. J. Geophys. Res., 116, D09110, doi: /2010jd Krol, M., S. Houweling, B. Bregman, M. Van Den Broek, A. Segers, P. Van Velthoven, W. Peters, F. Dentener, and P. Bergamaschi, 2005: The two-way nested global chemistry-transport zoom model TM5: Algorithm and applications. Atmos. Chem. Phys., 5, Masarie, K. A., G. Pétron, A. Andrews, L. Bruhwiler, T. J. Conway, A. R. Jacobson, J. B. Miller, P. P. Tans, D. E. Worthy, and W. Peters, 2011: Impact of CO 2 measurement bias on CarbonTracker surface flux estimates. J. Geophys. Res., 116, D17305, doi: / 2011JD Meng, Z. and F. Zhang, 2008: Test of an ensemble Kalman filter for mesoscale and regional-scale data assimilation. Part IV: performance over a warm-season month of June 2003, Mon. Wea. Rev., 136, Olson, S. C., J. A. Watts, and L. J. Allison: Major World Ecosystem Complexes Ranked by Carbon in Live Vegetation: A Database. NDP017, Carbon Dioxide Info. Analy. Cent., Oak Ridge Nat. Lab., Oak Rideger, Tenn., Ott, E., B. R. Hunt, I. Szunyogh, A. V. Zimin, E. J. Kostelich, M. Corazza, E. Kalnay, D. J. Patil, and J. A. Yorke, 2004: Estimating the state of large spatiotemporally chaotic systems. Phys. Lett. A, 330, Park, J. I. and H. M. Kim, 2010: Typhoon Wukong, 2006: Prediction Based on The Ensemble Kalman Filter and Ensemble Sensitivity Analysis (in Korean with English abstract). Atmosphere, 20, Peters, W. and Coauthors, 2004: Toward regional-scale modeling using the two-way nested global model TM5: Characterizatino of transport using SF 6. J. Geophys. Res., 109, D19314, doi: / 2004JD and Coauthors, 2010: Seven years of recent European net terrestrial carbon dioxide exchange constrained by atmospheric observations. Glob. Change Bio., 16, , doi: /j x and Coauthors, 2007: An atmospheric perspective on North American carbon dioxide exchange: CarbonTracker. Proc. Nat. Acad. Sci. U.S.A., 104, , J. B. Miller, J. S. Whitaker, A. S. Denning, A. Hirsch, M. C. Krol, D. Zupanski, L. Bruhwiler, and P. P. Tans, 2005: An ensemble data assimilation system to estimate CO 2 surface fluxes from atmospheric trace 한국기상학회대기제 22 권 4 호 (2012)

13 김진웅 김현미 조천호 427 gas observations. J. Geophys. Res., 110, D24304, doi: /2005jd van der Werf, G. R., J. T. Randerson, L. Giglio, G. J. Collatz, P. S. Kasibhatla, and A. F. Arellano Jr. 2006: Interannual variability of global biomass burning emissions from 1997 to Atmos. Chem. Phys., 6, Whitaker, J. S. and T.M. Hamill, 2002: Ensemble Data Assimilation without Perturbed Observations. Mon. Wea. Rev., 130, ,, X. Wei, Y. Song, and Z. Toth, 2008: Ensemble data assimilation with the NCEP global forecast system. Mon. Wea. Rev., 136, Zupanski, M., 2005: Maximum likelihood ensemble filter: Theoretical aspects. Mon. Wea. Rev., 133, Atmosphere, Vol. 22, No. 4. (2012)