Journal of Korean Society for Atmospheric Environment CMAQ-pollen 모델을이용한참나무꽃가루확산고해상도수치모의및검증 31 Vol. 33, No. 1, February 2017, pp. 31-44 https://doi.org/10.5572/kosae.2017.33.1.031 p-issn 1598-7132, e-issn 2383-5346 CMAQ-pollen 모델을이용한참나무꽃가루확산고해상도수치모의및검증 A High-resolution Numerical Simulation and Evaluation of Oak Pollen Dispersion Using the CMAQ-pollen Model 오인보 김규랑 1) 방진희 임윤규 2) 조창범 1) 오재원 3) 김양호 4) 황미경 5), * 울산대학교의과대학환경보건센터, 1) 국립기상과학원응용기상연구과, 2) 국립기상과학원환경기상연구과, 3) 한양대학교구리병원소아청소년과, 4) 울산대학교의과대학직업환경의학교실, 5) 부산대학교대기환경과학과 (2016년 12월 6일접수, 2017년 1월 4일수정, 2017년 1월 25일채택 ) Inbo Oh, Kyu Rang Kim 1), Jin-Hee Bang, Yun-Kyu Lim 2), Changbum Cho 1), Jae-Won Oh 3), Yangho Kim 4) and Mi-Kyoung Hwang 5), * Environmental Health Center, University of Ulsan College of Medicine 1) Applied Meteorology Research Division, National Institute of Meteorological Research 2) Environmental Meteorology Research Division, National Institute of Meteorological Research 3) Department of Pediatrics, Hanyang University Guri Hospital, Hanyang University College of Medicine, Guri 4) Department of Occupational & Environmental Medicine, University of Ulsan College of Medicine 5) Department of Atmospheric Sciences, Pusan National University (Received 6 December 2016, revised 4 January 2017, accepted 25 January 2017) Abstract The aim of this study is to evaluate the accuracy and variability of the oak pollen concentrations over the Seoul metropolitan region (SMR) simulated by the Community Multiscale Air Quality (CMAQ)-based pollen dispersion model, which is the CMAQ-pollen model integrated with the improved oak pollen emission model (PEM-oak). The PEM-oak model developed is based on hourly emission flux parameterization that includes the effects of plantspecific release, meteorological adjustment, and diurnal variations of oak pollen concentrations. A 33 day-run for oak pollen simulation was conducted by the CMAQ-pollen model with a 3 km spatial resolution for the SMR during the 2014 spring pollen season. Modeled concentrations were evaluated against the hourly measurements at three Burkard sampling sites. Temporal variations of oak concentrations were largely well represented by the model, but the quantitative difference between simulations and measurements was found to be significant in some periods. The model results also showed that large variations in oak pollen concentrations existed in time and space and high concentrations in the SMR were closely associated with the regional transport under strong wind *Corresponding author. Tel : +82-(0)51-583-2652, E-mail : hmk1001@pusan.ac.kr J. Korean Soc. Atmos. Environ., Vol. 33, No. 1, 2017
32 오인보 김규랑 방진희 임윤규 조창범 오재원 김양호 황미경 condition. This study showed the effective application of the CMAQ-pollen modeling system to simulate oak pollen concentration in the SMR. Our results could be helpful in providing information on allergenic pollen exposure. Further efforts are needed to further understand the oak pollen release characteristics such as interannual variation of the oak pollen productivity and its spatio-temporal flowering timing. Key words : Oak pollen, Seoul metropolitan region, CMAQ-pollen model, Pollen emission model 1. 서론꽃가루는알레르기질환을일으키는중요한공기중항원 (aeroallergens) 으로잘알려져있다 (Sofiev et al., 2013). 최근에는기후변화와연계된꽃가루비산량의증가및알레르기유발성 (allergenictiy) 의변화 (Sofiev et al., 2013; Singer et al., 2005; Beggs, 2004) 로인해알레르기질환관리에중요한대기환경인자로부각되고있다. 아울러국내 외알레르기질환유병인구의증가경향과 (Healthcare Bigdata Hub, 2014; D Amato et al., 2007; Asher et al., 2006) 꽃가루항원에대한높은감작률은 (Oh et al., 2009a, 2009b; Kim et al., 1999) 시 공간적꽃가루분포의정확한이해및예측의필요성을시사한다. 알레르기꽃가루의시 공간적변동과예측에관한연구들은유럽을중심으로 2000년전후로시작되었다 (Helbig et al., 2004; Campbell et al, 1999; Corden and Millington, 1999; Cabezudo et al., 1997; Norris-Hill, 1995). 최근에는제한된꽃가루모니터링자료와통계모델의한계를보완하기위해꽃가루수치모델링방법이구체적으로제안되었고 (Efstathiou et al., 2011; D Amato et al., 2007; Sofiev et al., 2006), 이후배출량추정식의개선과해상도높은입력자료의적용을통해수치모델을이용한꽃가루농도예측의신뢰성이높아지고있다. Sofiev et al. (2006) 은 SILAM (System for Integrated modelling of Atmospheric composition) 모델을개발하여서유럽전역을대상으로꽃가루장거리수송을모의하였고이후 SILAM에꽃가루비산시기결정모델을추가하고자작나무꽃가루배출모수화식을개선하였다 (Sofiev et al., 2013). Efstathiou et al. (2011) 은수정된 Helbig et al. (2004) 모수화식을 CMAQ (Community Multi-scale Air Quality) 모델에적용하여 2002년미국뉴욕지역을대상으로꽃가루모델링을수행하였다. Zink et al. (2013) 은개화후꽃가루의대 기유입전단계인 pollen reservoir 단계 ( 꽃이맺혀있는단계 ) 를도입하여배출량산정모수화식을개선하고기존 COSMO-ART (COnsortium for Small-scale MOdelling - Aerosols and Reactive Trace gases) 모델의정확도를평가하였다. 국내에서는꽃가루농도예측과관련하여기상인자와의관련성에근거한통계모델개발연구가최근있었고 (Oh, 2009; Park et al., 2008), 이를기반으로기상청에서는참나무, 소나무, 잡초류를대상으로꽃가루농도위험지수예보를시행중에있다 (KMA, 2016). 하지만통계모델의지역적제한성과농도값에대한불확실성이크고해상도높은정량적예측에는한계가있다. 최근 Oh et al. (2012) 이울산지역을대상으로 WRF (Weather Research and Forecasting Model)-CMAQ을이용하여참나무꽃가루모델링을최초로수행하고꽃가루시 공간적변화를전반적으로재현하였다. 이를근간으로수도권지역을대상으로한참나무꽃가루수치모델링시스템을구축하고개선하는연구가수행되었다 (National Institute of Meteorological Sciences, 2015, 2014, 2013). 또한 Lim et al. (2015) 는로버스트다중회귀식을이용하여꽃가루배출량을추정하고, 전국기반 UM (Unified Model)-CMAQ 모델을이용한참나무꽃가루모델링을수행하여국내예보를위한기초연구가되었다. 이러한국내 외연구들에서최근시도되고있는꽃가루수치모델링은알레르기질환관리와예방에중요한잠재적노출정보를제공해줄수있다는측면에서의미가크다. 기후변화와연계해서도미래의꽃가루비산기간, 비산량변화를예측하고대응할수있는중요한부분이기도하다 (Duhl et al., 2013; Zhang et al., 2013). 하지만성공적인수치모델링에있어꽃가루배출량추정의불확실성, 물리적특성과침적과정의정확한고려, 식생입력자료의한계극복을위해서지속적인연구가필요할것이다. 한국대기환경학회지제 33 권제 1 호
CMAQ-pollen 모델을이용한참나무꽃가루확산고해상도수치모의및검증 33 이연구에서는우리나라의대표적인수종이자알레르기성 (allergenicity) 이높은참나무꽃가루를 (Oh, 2007; Lee et al., 2006) 모의할수있는수정된배출량모델과 CMAQ 모델을제시하고, 수도권지역봄철사례기간참나무꽃가루확산에대한해상도높은모델링을수행하였다. 아울러대상영역내모니터링된실측꽃가루농도와의비교를통해모델링결과를검증하고향후참나무꽃가루수치예보의가능성을평가하였다. 2. 연구방법 2. 1 참나무꽃가루수치모델링시스템꽃가루수치모델링은크게꽃가루배출량추정과대기로배출된꽃가루의확산이류과정을포함한다. 배출량추정과정은신뢰성있는꽃가루농도예측을위해매우중요한부분으로식생 ( 발생원 ) 의분포, 꽃가루생성량과생물학적계절변화, 주요기상인자들과의관계식등을통해추정가능하다. 특히기상인자들과의관계는꽃가루배출량의단기적시 공간적변화를계산하는데있어중요하다. Helbig et al. (2004) 은배출량추정에마찰속도를비롯한기온, 습도, 풍속을고려하였고, 이는이후여러꽃가루대기확산모델연구 (Efstathiou et al., 2011; Vogel et al., 2008; Sofiev et al., 2006) 에적용및응용되었다. 최근에는 Siljamo et al. (2013) 이자작나무꽃가루대상으로민감도모델링을통해풍향, 강수, 습도가배출량추정에중요한기상인자임을제시하였다. 그림 1은이연구에서제시하는참나무꽃가루수치모델링시스템의구성과수치모의흐름도이다. 시스템은크게참나무꽃가루배출량추정모델 (PEM-oak) 과기상입력자료를생성하는 WRF 모델, 그리고참나무꽃가루의대기중이류, 확산, 침적과정을시 공간적으로계산할수있도록수정된 CMAQ 모델인 CMAQpollen 모델로구성된다. 먼저 PEM-oak는 Helbig et al. (2004) 과 Efstathiou et al. (2011) 연구에서제시된모수화식을근간으로, 모니터링된참나무꽃가루의분포특성 ( 계절 / 일중변화 ) 과기상요소와의관계를개선하여고안된배출량모델이다. 영역별배출량산정은참나무분포와필요한기상요소를입력하여격자별로계산된다. 여기서참나무분포는수치임상도 (Korea Forest Service, FGIS) 를이용하여격자별참나무면적으로산출되며, 격자별기온, 습도, 풍속, 마찰속도, 엽면적지수 (LAI, Leaf Area Index) 는 WRF 모델링결과로부터입력된다. CMAQ-pollen 모델은대기중으로배출된참나무 Fig. 1. Schematic flowchart depicting the modified WRF-CMAQ modeling system linked with the PEM-oak model for oak pollen simulations. J. Korean Soc. Atmos. Environ., Vol. 33, No. 1, 2017
34 오인보 김규랑 방진희 임윤규 조창범 오재원 김양호 황미경 꽃가루의이류, 확산, 침적과정을모의할수있도록, 본연구에서기존 CMAQ 모델에관련모듈과 I/O 프로세스를추가 / 수정한모델이다. WRF 모델결과가입력되어 MCIP (Meteorology-Chemistry Processor) 전처리과정을통해산출된기상인자, PEM-oak 모델에서계산된참나무꽃가루배출량, 영역별꽃가루농도의초기 / 경계조건등이입력자료로요구된다. 참나무꽃가루의이류, 확산과정은수정된 CCTM module에서처리되어매시간별격자별참나무꽃가루농도가계산된다. 2. 2 대상영역, 측정자료, 모델링기간참나무꽃가루모델링대상영역은그림 2와같이구성된다. 한반도영역 (D1, 67 79, 9-km 격자간격 ) 의경우대기정책모델링지원시스템 (Clean Air Policy Modeling System: CAPMOS) 에서제시하는국가표준 도메인과동일하다. 수도권영역 (D2, 64 58, 3-km 격자간격 ) 은서울중심의수도권전체를포함하며주변산림지역의지역규모꽃가루수송영향을적절히고려하기위해강원도일부를비롯한인근산림지역을포함하여구성하였다. 수도권영역안에는기상청에서운영하는세개의꽃가루모니터링지점 ( 서울 : P1, 구리 : P2, 포천 : P3) 이위치하며배출량모델개발과 CMAQpollen 결과검증에이용되었다. 참나무꽃가루모니터링자료는버커드샘플러 (Burkard 7-day recording volumetric spore sampler) (Burkard, 2001) 를이용하여채집된시료를검경하여얻어진농도 ( 단위 : grains m -3 ) 자료이다. PEM-oak 모델개발과관련해서는, 6년간 (2009~2014년) P1과 P2 지점에서모니터링된일평균농도자료와 2014년세지점의 1시간평균농도자료를이용하였다. 40 N 30 N (a) 120 E 130 E (b) Oak (%) 75~ 75 70 65 60 55 50 45 40 35 30 25 20 15 10 5 0 (c) Height (m) 1000 900 800 700 600 500 400 300 200 100 0 Fig. 2. (a) Map showing 9-, 3-km grid nests (D1 and D2) used for CMAQ-pollen modeling. (b) Density (%) of oak trees from the FGIS digital forest type map. (c) Topographic feature of the Seoul and its surrounding area. The read dots indicate pollen monitoring sites (P1: Seoul, P2: Guri, P3: Pocheon) and the blue dots indicate meteorological monitoring sites (M1: KMA-ASOS, M2: Guri-AWS M3: Pocheon-AWS). 한국대기환경학회지제 33 권제 1 호
CMAQ-pollen 모델을이용한참나무꽃가루확산고해상도수치모의및검증 35 농도자료기반참나무꽃가루비산계절 (oak pollen season) 을결정하고, 이기간은 PEM-oak 모델비산계수및가중치추정과모델링기간설정에이용되었다. 본연구에서는 Emberlin et al. (1994) 이제시한방법으로 ( 전체누적꽃가루농도의 1% 이상시점부터 99% 시점까지 ) 비산계절을정의하였고이에따라 2014년비산계절기간은 P1 지점은 4월 12일 ~5월 10일, P2 지점은 4월 11일 ~5월 10일, P3 지점은 4월 15일 ~5월 9 일로각각결정되었다. 모델링기간은세지점의참나무꽃가루비산계절기간을포함하도록선정하여, WRF 모델링은 2014년 4월 7일 (0000 UTC) 부터 5월 11일 (0000 UTC) 로 35일간을대상으로수행되었고, CMAQ-pollen 모델링은하루의초기값적응기간을포함하여 4월 9일 (0000 UTC) 부터 5월 11일 (0000 UTC) 까지총 33일간이루어졌다. 2. 3 WRF 모델링 WRF 모델링을통해얻어지는기상인자는 PEM-oak 모델과 CMAQ-pollen 모델의입력자료가된다. 모델의초기및경계조건은기상청에서제공하는 6시간간격의 KLAPS (Korea Local Analysis and Prediction System) 5 km 분석장과 GFS (Global Forecast System) 1도토양수분자료를이용하였다. 연직층은 100 hpa 이하 44개 sigma layer로구성되었으며최하층고도는약 32 m이다. 모델링은단방향네스팅 (one-way nesting) 방법을통해수행되었다. 물리옵션은모든모델링영역에대해, Microphysics는 WSM6 scheme, Long wave radiation은 RRTM scheme, Shortwave radiation은 Dudhia scheme, PBL physics는 Yonsei University scheme, Land surface physics는 Noah Land-Surface Model을각각적용하였다. 대상지역의현실적인지형과지표면상태를 WRF 모델링에반영하기위하여 SRTM (Shuttle Radar Topography Mission) 3초 ( 약 90 m) 지형자료와국내환경부에서제작한환경공간정보서비스 EGIS (Environmental Geographic Information System, 2007) 의중분류 ( 해상도 5 m급, 22개항목 ) 토지피복자료를적용하였다. 2. 4 참나무꽃가루배출량모델입력자료그림 3은 PEM-oak 모델의모수화식 (3장에구체적으로제시 ) 에필요한입력자료, 격자별참나무화분배 출량추정과정을간략히보여준다. 모수화식의연중비산계수산정을위해 6년간 (2009~2014년) P1과 P2 에서모니터링된일평균꽃가루농도자료를분석하였다. 또한세지점모두에서 2014년모니터링된 1시간평균농도자료를이용하여기상조절항과가중치항 ( 이후설명 ) 의계수를만들었다. PEM-oak 모델의입력자료로격자별참나무배출강도의대리인자인참나무림의분포면적이요구된다. 고해상도 5차수치임상도 (Korea Forest Service, FGIS) 자료 ( 축적 1 : 5,000, 한반도영역 ) 에서 GIS 수종별속성정보를처리하여참나무분포를얻었다. 이자료를 MIMS (Multi-scale Integrated Modeling System) Spatial Allocator에입력하여각격자별참나무림면적을산출하였다 ( 그림 2(b) 참조 ). 추가적으로요구되는격자별기상인자 ( 지상 2 m 기온, 10 m 풍속, 2 m 비습 ) 와엽면적지수 (LAI), 마찰속도정보는 WRF 결과에서얻을수있다. 여기서엽면적지수는토지피복도에의존하며, EGIS 토지피복자료를통해현실적인값을산출할수있도록하였다. 2. 5 CMAQ-pollen 모델 CMAQ-pollen 모델은대기중으로배출된참나무꽃가루의수송 ( 이류 / 확산 ), 침적과정을모의할수있도록기존 CMAQ 모델 (ver. 5.0.1) 에서관련모듈과 I/O 프로세스를수정한모델이다. 입력자료로는 WRF 모델결과로부터 MCIP 전처리과정을통해산출된기상인자, PEM_oak 모델에서계산된참나무꽃가루배출량, 영역별꽃가루농도의초기 / 경계조건등이요구된다. 참나무꽃가루의수송과정은수정된 CCTM module ( 수평이류 : hyamo, 연직이류 : vwrf, 수평확산 : multiscale, 연직확산 : acm2, 에어로졸 : aero5, 침적 : m3dry) 에서처리되어시 공간적격자별참나무꽃가루농도가계산되어진다. 각시간대의꽃가루농도계산은수평이류 (x y방향 ), 연직이류, 수평확산 (x y방향 ), 연직확산, 에어로졸관련모듈이순차적으로처리되어이루어진다. 참나무꽃가루의물리적특성을고려하기위해수정된연직확산및에어로졸모듈에서밀도 (ρ) 를 1,058 kg m -3 (Zhang et al., 2013) 로입력하였다. 참나무꽃가루의광학직경 (D) 은 27.7 μm로 Coarse particle로처리되며 (Schueler and Schlunzen, 2006), 식 1은 CMAQ 내에어로졸건성침적속도를계산하는기본식으로 (Bink- J. Korean Soc. Atmos. Environ., Vol. 33, No. 1, 2017
36 오인보 김규랑 방진희 임윤규 조창범 오재원 김양호 황미경 Fig. 3. Flowchart for calculating the gridded oak pollen emissions using the PEM-oak model. owski and Shankar, 1995) 위의수치를적용하여참나 무꽃가루에대한건성침적속도 v d 가계산된다. g ρ p v d =[ D2] Cc (1) 18v ρ air ------ ( ------ ) 여기서, v=kinematic viscosity of air ρ p =density of particles ρ air =density of air D =particle diameter Cc=linearized slip correction 모델링수행에있어참나무꽃가루농도의초기장은모든영역에서 0 grains m -3 으로입력하여생성시켰다. 경계조건농도값은한반도영역에대해서는모든방향으로 0 grains m -3 으로처리하였고수도권영역과서울영역은각각상위도메인으로부터계산된격자농도값으로만들었다. 초기값과경계값이현실의수치와분명차이가있으나현재선행연구나참나무꽃가루의배경농도에대한정보가없어정확한값의고려가어렵다. 량이다르며기상조건또는식물학적리듬에따라달라질수있다. 또한수종에따라특징있는계절및일변화를가지며, 기온, 습도, 풍속등의기상조건에영향을크게받는다 (Oh et al., 2013). PEM-oak 모델의기본모수화방정식을식 2에나타내었다. 이는 Efstathiou et al. (2011), Helbig et al. (2004), Oh et al. (2012) 을바탕으로수정된모수화식이다. F p =c* c e K e u* w f(hr) (2) 여기서참나무꽃가루의시간별방출되는플럭스 (F p, grains m -2 s -1 ) 는수종의특성농도 (characteristic concentration, c*), 개화관련식물고유특성을반영하는인자 ( 일종의비산계수 ) (plant-specific factor, c e ), 기상조절인자 (meteorological adjustment factor, K e ), 마찰속도 (friction velocity, u*) 에의존한다. 추가적으로일중꽃가루배출및비산과관련된가중치 (w f(hr) ) 에영향을받게된다. 먼저수종의특성농도 c* 항은식 3과같이정의된다 (Efstathiou et al., 2011; Helbing et al., 2004). 3. 참나무꽃가루배출량추정꽃가루배출량은기상학적조건및이와관련된식물학적주기 (Mullins and Emberlin, 1997) 와밀접한관련을가지고있다. 수종의종류에따라꽃가루총생성 P c*=( ------------- q ) (3) LAI C h 여기서, P q 는연중단위면적당참나무꽃가루총생산량 (seasonal total oak pollen production) 으로, 본연구에서는연간참나무한그루당평균꽃가루개수 한국대기환경학회지제 33 권제 1 호
CMAQ-pollen 모델을이용한참나무꽃가루확산고해상도수치모의및검증 37 ce 0.10 0.08 0.06 0.04 0.02 0.00 Oak pollen Ce 0 5 10 15 20 25 30 35 Number of days from onset day (day) 250 200 150 100 50 Oak pollen concentration (grains m -3 ) Fig. 4. Variation of oak pollen concentrations from the onset day and its production rate estimated by the log-normal distribution during oak pollen season, 2009~2014. 2.6 10 10 grains tree -1 yr -1 (Jato et al., 2007), 참나무밀도 3,390 tree ha -1 (Saito et al., 2006) 을이용하여 P q 를 8.814 10 9 grains m -2 로입력하였다. 엽면적지수 LAI는기상모델 WRF에서계산된값을이용하였고 C h 는참나무의 canopy height로 5 m를입력하였다. 식 3을통해얻어진 c* 값과식 2의다른변수 ( 이후 3.1절부터설명 ) 들을고려하여매시간별대기로방출되는꽃가루플럭스를얻을수있다. 모수화식을통해격자별 / 매시간별참나무꽃가루배출량이계산되면, 이는 CMAQ-pollen 모델의연직최하층에입력된다. 3. 1 비산계수 (c e ) 산정비산계수 c e 는개화와관련하여식물고유의특징을반영하는무차원계수 (dimensionless plant-specific factor) 로서계절적변화에따른꽃가루배출량을결정하는중요한인자이다. 먼저 P1과 P2 지점에서모니터링된일평균농도자료를이용하여참나무꽃가루비산계절 (oak pollen season, 6년간평균 35일 ) 의비산시작일을결정하였다. 그림 4의막대그래프는 P1과 P2 두지점의비산시작일을기준으로이후의평균참나무꽃가루농도의일별변화를보여주는것으로, 비산시작일시작후급격히농도가증가하고약 4~5일경에최고값이나타난후점차감소하는패턴을가진다. 식 4는이러한참나무꽃가루농도변화를근거로최적의 fitting curve를 log-normal distribution을가정하여추정된모수화식이다. P(d)=143.2808 exp(-0.5(ln(d/4.7538)/0.8481) 2 ) (4) 여기서, d는비산시작일로부터의일수 (number of day from pollen onset day) 이며, P(d) 는추정된일별꽃가루농도 (daily potential pollen concentration during pollen season) 이다. 비산시작일이예측되면식 4에서일별 P(d) 값이계산된다. 식 4를통해추정된값과모니터링값과의상관도는 r =0.88로높게나타났다. 최종적으로일별 c e 는 oak pollen season 기간의누적농도 (P total ) 에대한일별 P(d) 의비로나타낼수있다 ( 식 5). P(d) N c e (d) = -------- with P total = P(d) (5) P total d=1 계산된 c e 값의변화는그림 4의실선으로나타나있다. 참나무꽃가루비산시기동안 c e 값의변화가실제꽃가루농도의변화경향을잘설명함을확인할수있다. 하지만꽃의개화로인한꽃가루의배출은대기중꽃가루농도와다른개념으로 c e 의정확한계산을위해서는순림지역에서의장기간모니터링정보가필요하다. 3. 2 기상조절항 (K e ) 산정참나무꽃가루배출량산정에있어기상학적조절항 (meteorological adjustment factor, K e ) 은기상모델로부터얻어지는기온 (T), 풍속 (ws), 상대습도 (RH) 가인자 J. Korean Soc. Atmos. Environ., Vol. 33, No. 1, 2017
38 오인보 김규랑 방진희 임윤규 조창범 오재원 김양호 황미경 로고려되어계산되는항으로, 꽃가루의비산에대한 저항을고려한식으로구성된다 (Oh et al., 2012; Helbig et al., 2004) ( 식 6). K e 항의 T te, ws te, RH te 는각각기온, 풍속, 상대습도에대한임계치이다. 이들값은수종에 따라달라지며꽃가루비산과의관련성을통해특정값이입력된다. 여기서기온및풍속이임계치이상으로증가할경우 K e 값은증가하는반면, 상대습도는값이커질수록저항으로작용하여 K e 값을감소시킨다. 이연구에서는 K e 모수화식을국립기상연구소 (2014) 의결과를바탕으로임계치와가중치 (c 1, c 2, c 3 ) 를설정을하였다 ( 식 6). 기온, 풍속, 습도모두꽃가루비산에대한저항을의미하는 α항의분모항으로고려하였다. 가중치의경우, 세개의기상인자와참나무꽃가루농도와의상대적인관련성을근거로 c 1 을 0.5, c 2 를 2, c 3 를 1로설정하였다. K e =1-α (α<1) (6) =0 (α 1) 3 여기서, α=----------------------------------- T ws RH te c 1 -----+c 2 ------+c 3 ------- T te ws te RH T te =temperature threshold values (8 ) ws te =wind speed threshold values (2.5 m s -1 ) RH te =relative humidity threshold values (90%) c 1, c 2, c 3 =plant-specific constants to weigh the influence of the corresponding meteorological resistance (c 1 =0.5, c 2 =2, c 3 =1) 3. 3 일중비산가중치 (w f(hr) ) 산정 PEM-oak 모델의 K e 항에서높은기온과낮은상대습도값이입력되면참나무꽃가루배출량이많아진다. 즉일반적으로기온상승과습도감소가나타나는오후의경우배출량이증가하게된다. 하지만참나무의꽃가루농도는오전 ~ 정오경에최고농도가나타나며, 경우에따라오후에 2차정점을보이기도한다 (Corden and Millington, 1999; Kapyla, 1984). 이러한경향을모델이적절히반영토록하기위해 1시간간격으로모니터링된대기중농도자료로부터가중치 (w f (hr) ) 를추정하여모델에고려하였다. 이를위해 2014년 P1, P2, P3 세지점에서측정된시간별참나무꽃가루농도와기상조건을함께분석하였 다. 대상일선정은지점별꽃가루집중비산기간 ( 참나무꽃가루누적농도의 10~90% 에해당되는기간 ) 중강수가없는날, 임의의특정풍계가 20시간이상지속적으로존재하는 6일로하였다. 이러한조건은꽃가루농도측정값이참나무산림지역에서얻어진것이아니므로수송효과와관련한바람변화에따른꽃가루농도변화를최대한배제하기위해서이다. 선정된 6일 (P1: 4 월 15, 19일 ; P2: 4월 19일 ; P3: 4월 19, 20, 29일 ) 의시간별평균농도를이용하여최적의 fitting curve를가우시안분포를가정하여식 7을추정하였다. w f(hr) (t)= 0.0763 exp(-((t-5.4775) 2 /(2 9.2140 2 ))) (7) 여기서, 시간 (t, hour (LST)) 별추정된값과모니터링값의일중비율의상관도는 r =0.87 (p<0.01) 이상으로꽃가루농도의일중변화를잘설명하는것으로평가되었다. 4. 참나무꽃가루고해상도수치모의결과와검증 4. 1 WRF 모델결과검증및참나무꽃가루배출분포 CMAQ-pllen 모델링수행에앞서 WRF 모델링결과를검증하고, 격자별참나무분포면적및기상모델결과를입력하여수행한 PEM-oak 모델결과를제시하였다. WRF 수행결과검증은모델링기간 (2014년 4월 9일 5월 11일 ) 동안계산된기온 (2 m) 과풍속 (10 m) 을수도권 ASOS 기상관측지점 ( 서울, 인천, 수원, 구리, 포천 ) 의값들과비교하였다. 표 1에는통계지표 MBE (mean bias error), MAGE (mean absolute gross error), RMSE (root mean squared error), IOA (index of agreement) 를이용한검증결과들을제시하였다. 전반적으로 WRF 모델이기온과풍속의시간별변화를잘재현하는모습이다. 기온은모델치가관측치와비교해다소과소평가 ( - 2.6 0.58 ) 되었으나정확도는 0.9 이상의높은값을나타내었다. 풍속은모델이약 1 m s -1 과대평가 ( 수원제외 ) 하였으나풍속 RMSE, IOA가모든지점에서신뢰구간에포함되었다. 그림 5는참나무꽃가루농도가높게나타난 2014년 한국대기환경학회지제 33 권제 1 호
CMAQ-pollen 모델을이용한참나무꽃가루확산고해상도수치모의및검증 39 Table 1. Evaluation of WRF modeling results: 2-m temperature and 10-m wind speed. Meteorological parameters Temperature (2 m agl.) Wind speed (10 m agl.) Statistical index (Confidence range 1) ) ASOS/AWS sites for evaluation Seoul (108) Inchoun (112) Suwon (119) Pocheon (504) Guri (569) MBE (±0.5 ) - 1.24-2.61 0.58-0.17-0.58 MAGE (<2 ) 1.68 3.05 1.45 1.64 1.42 IOA( 0.8) 0.93 0.77 0.94 0.95 0.95 MBE (±0.5 m/s) 0.29 1.05-0.22 1.15 1.23 RMSE (<2 m/s) 1.36 1.72 1.02 1.99 1.84 IOA( 0.6) 0.82 0.74 0.81 0.70 0.69 1) Confidence range is included in the EPA Draft Guidance on meteorological model evaluation (2009) 0300 LST 19 April 0600 LST 1200 LST 1800 LST (g/s) 10.00 8.00 7.00 6.00 5.00 4.00 3.00 2.00 1.00 0.50 0.00 Fig. 5. Horizontal distributions of oak pollen emission fluxes (Domain 2) calculated with the PEM-oak in April 19, 2014. 4 월 19 일수도권영역주요시간대의참나무꽃가루배 출량분포를보여준다. 그림 2 에서제시된참나무밀도 의공간적분포와비교해보면, 배출량의정도가참나 무분포에크게의존함을확인할수있다. 또한서울외 곽경기도지역에상대적으로높은배출량이계산되었고, 도시역을중심으로참나무꽃가루배출이거의없 J. Korean Soc. Atmos. Environ., Vol. 33, No. 1, 2017
40 오인보 김규랑 방진희 임윤규 조창범 오재원 김양호 황미경 Oak pollen concentration (grains m -3 ) 4500 3000 2000 1500 1000 500 0 2000 1500 1000 500 3000 0 1500 1000 500 P1 P2 P3 Measurement CMAQ-pollen 0 4/12 4/14 4/16 4/18 4/20 4/22 4/24 4/26 4/28 4/30 5/02 5/04 5/06 5/08 5/10 Date (LST) Fig. 6. Comparison between modeled and measured oak pollen concentrations for the 12 April-10 May, 2014. 는반면북한산과도봉산을중심으로높은배출량이존재함을알수있다. 참나무꽃가루배출량의시간적변화는기상요소, 일중가중치등이복합적으로고려되어다양한분포를확인할수있다. 0300 LST부터꽃가루배출량이점차증가하여오후 (1200~1500 LST) 에최대값을나타내고다시감소하는모습이다. 4. 2 CMAQ-pollen 모델링결과및검증 CMAQ-pollen 모델이재현한참나무꽃가루농도 (4 월 12일 0100 LST~5월 11일 0000 LST) 를모니터링값과의비교를통해시간적변화의경향과정량적인차이를평가하였다. 그림 6에는모델링영역내에위치한세지점 ( 그림 1의 P1~P3) 에서모니터링된참나무꽃가루의시간별농도변화와함께각모니터링지점들위치에해당하는 CMAQ-pollen 수치모의농도를제시하였다. 세지점에서모니터링된꽃가루농도는정량적인차이가존재하지만 4월초부터증가하기시작하여 4월중후반에최고치가나타나고이후약한달간증감을반복하면서점차감소하는경향을보인다. 특히서울에위치한 P1 지점은 4월 14일 ~20일기간동안급격한농도상승현상과고농도빈도가상대적으로뚜렷함을 볼수있다. 높은농도가집중되어나타난 4월 14일부터 20일까지늦은저녁과새벽에일최고농도가나타나는특징이있고최고값은 15일 0600 LST에 4,600 grains m -3 로기록되었다. 구리에위치한 P2 지점은 P1 지점농도와비교해농도수준이낮고큰증감의변화를볼수없었는데, Lim et al. (2015) 는이러한원인으로측정소주변환경영향 ( 건물차폐영향 ) 을제시된바있다. 고농도는 4월 16일과 18일에 328 grains m -3, 4월 30일에 390 grains m -3 로기록되었다. 포천의 P3 지점역시 P1, P2와는다소차이가있는변화를보였지만 4월 19일의고농도패턴은 P1 지점과유사하다. 최고값은 4월 19일 0500 LST에 2,547 grains m -3 로기록되었고고농도발생시점은 P1 지점과유사하게주로새벽시간에많았다. 모델링기간세측정지점에해당되는격자의시간별모델값 (n =2,088) 을대표적인통계지표를이용하여해당측정값과비교검증한결과 r =0.74, MB=44.9, NMB (%)=80.0%, NME (%)=130.8%, IOA=0.85로나타났다. 이러한통계검증수치들은모델의성공적인수행으로충분히설명하기어렵지만, 이연구에서제시한 CMAQ-pollen 모델이참나무꽃가루농도의전반적인변화를예측할수있음을보여준다. 그림 6에서 한국대기환경학회지제 33 권제 1 호
CMAQ-pollen 모델을 이용한 참나무 꽃가루 확산 고해상도 수치모의 및 검증 41 Reference vectors 0.5 0300 LST 0600 LST 1200 LST 1800 LST 7 (grain m-3) Fig. 7. Horizontal distributions of modeled oak pollen concentrations and wind vectors (Domain 2) in April 19, 2014. 알 수 있듯이 모델을 통해 재현된 참나무 꽃가루 농도 는 측정된 농도수준과 다소의 차이는 있지만 전반적인 그림 7은 사례일 주요 시간대 CMAQ-pollen과 WRF 모델 결과로 참나무 꽃가루 농도와 지상바람 (10 m)의 변화의 패턴은 유사하다. 또한 P1 지점의 경우 CMAQ- 수평 분포를 보여준다. 전반적으로 참나무 꽃가루의 배 잘 재현되었음이 확인된다. 한편 18일의 경우, P1과 P3 가 높게 나타나며 배출량과 유사하게 공간적 농도차가 지점에서 모델결과가 과대평가 되었는데, 이는 오전 큼을 알 수 있다. 특히 서울의 북쪽 인근지역과 수도권 중의 기온 과대모의 (P1, P3 지점 기준, 최대 + 3.9 ) 북동쪽 외곽 강원도 지역의 고농도 현상이 뚜렷하며 pollen 수치모의를 통해 15일, 19일의 고농도 현상이 로 꽃가루 배출량 증가가 주된 원인으로 분석된다. 출량이 많은 지역 (그림 5 참고)과 인근에 대기 중 농도 바람에 의한 수송현상을 확인할 수 있다. 19일 새벽과 수치모델의 적용은 해상도 높은 공간적 정보를 제공 이른 아침에 나타난 비교적 강한 풍속조건과 남동내지 하고 고농도 및 농도변화의 원인을 이해하는 데 중요한 동풍계열의 바람의 조건들은 산림지역 꽃가루 비산을 역할을 할 수 있다. 앞선 모델링 결과 검증에서 고농도 강화하고 (식 6 참고) 측정지점인 P1 (1,807 grains m-3, 일 측정값의 변화를 비교적 잘 재현한 4월 19일을 대상 0600 KST)과 P3 (2,547 grains m-3, 0500 KST)에서의 으로 주요 시간대별 (0300, 0600, 1200, 1800 LST) 모의 고농도 원인이 될 수 있다. 정확한 고농도 원인을 파악 된 참나무 꽃가루 농도의 공간적 분포를 파악하였다. 하기 위해서는 지상 고농도의 원인이 될 수 있는 대기 J. Korean Soc. Atmos. Environ., Vol. 33, No. 1, 2017
42 오인보 김규랑 방진희 임윤규 조창범 오재원 김양호 황미경 경계층에서의혼합과정 (Raynor et al., 1974) 을분석하고산림지역에서비산된참나무꽃가루의공간적수송현상을보다면밀히분석할필요가있다. 5. 결론이연구에서는알레르기질환과관련한대표적인바이오에어로졸인참나무꽃가루를수치모의할수있는배출량모델 (PEM-oak) 과수정된 CMAQ 모델 (CMAQpollen) 을제시하였다. 이모델들을이용하여봄철약한달간 (2014년 4월 9일 5월 11일 ) 수도권지역을대상으로고해상도수치모의를수행하고측정값과의검증을통해모델의적용가능성을평가하였다. PEM-oak 모델에는고해상도참나무림분포정보와함께참나무꽃가루생성의식물계절학적특성을반영된비산계수및기상의영향을적절히고려하였고, 측정자료를근거로일중비산의가중치를추정하여배출량의시간적변화에대한정보의신뢰성을높였다. 그리고 PEM-oak 모델을이용하여한반도-수도권영역의참나무꽃가루배출량을산정하였고참나무꽃가루의수송, 침적과정들을처리할수있게수정된 CMAQpollen 모델을이용하여모델링을수행하였다. 그결과모델이동일기간영역내세지점 ( 서울, 구리, 포천 ) 에서측정된농도의변화패턴을전반적으로잘재현함을보여주었고, 향후모델개선을통한참나무꽃가루수치모의의가능성이높게평가되었다. 하지만단기적고농도현상을모델이충분히재현하지못하는한계가있었고, 이는배출량산정의불확실성이주된원인으로평가되며모델입력자료의해상도와정확도를함께향상시키는연구가반드시필요함을보여주었다. 현재산림현장자료와선행연구들이부족해참나무꽃가루배출량추정의불확실성은크다. 대기오염물질과는달리꽃가루의대기중배출은기후조건과관련하여매년, 계절적특성이강하게반영된다. 즉정확한배출량추정을위해참나무순림에서의꽃가루및기상인자집중측정이지속적으로필요하며, 이는비산시작과연간배출량, 계절적변화등식물계절학적특성을모델이적절히반영할수있는중요한자료가될것이다. 이와함께지속적인모수화식개선을통해신뢰성있는꽃가루배출량산정이가능할것이다. 또한참 나무림의정확한분포정보생성과경계조건값의불확실성해소 ( 예로북한지역배출량고려 ) 역시모델링예측결과의정확도를높일수있을것이며, 참나무꽃가루관련물리적특성과관련한정보개선, 재비산과정고려등도수치모델개선연구에중요한향후연구과제이다. 이연구에서시도한참나무꽃가루수치모델링은봄철꽃가루알레르기질환예방과관리에있어중요한실외환경노출정보를제공해줄수있다는측면에서매우의미있는연구라사료된다. 아울러기후변화와알레르기질환자의증가추세와대응하여미래의꽃가루농도변화를예측할수있는기초연구로서활용가치가있다. 향후꽃가루모니터링지점의확대와모델개선연구를통해참나무꽃가루에대한정량적농도예측연구의고도화가필요할것이다. 아울러참나무외에국내넓게분포하는소나무나알레르기성이큰삼나무등과같은수목들의꽃가루에대한모델링역시, 상기언급된배출량추정의정확성문제가해결된다면수종분포에대한공간정보가활용가능함으로충분히수치모의가가능하리라사료된다. 또한이연구에서제시한배출량모델의추가적인개선과모델의검증연구가이루어진다면대기질예보와연계하여현업에활용할수있는꽃가루수치예보의도입역시가능할것이다. 감사의글본연구는국립기상과학원응용기상기술개발연구와환경부환경보건센터재원에의해이루어졌습니다. References Asher, M.I., S. Montefort, B. BjorLSTen, C.K. Lai, D.P. Strachan, S.K. Weiland, H. Williams, and the ISAAC Phase Three Study Group (2006) Worldwide time trends in the prevalence of symptoms of asthma, allergic rhinoconjunctivitis, and eczema in childhood: ISAAC phases one and three repeat multicountry cross-sectional surveys, Lancet, 368, 733-743. 한국대기환경학회지제 33 권제 1 호
CMAQ-pollen 모델을이용한참나무꽃가루확산고해상도수치모의및검증 43 Begges, P.J. (2004) Impacts of climate change on aeroallergens: past and future, Clinical & Experimental Allergy, 34, 1507-1513. Binkowski, F.S. and Shankar, U. (1995) The regional particulate matter model: 1. Model description and preliminary results. Journal of Geophysical Research: Atmospheres (1984-2012), 100(D12), 26191-26209. Burkard (2001) 7-day recording volumetric spore trap. (http:// www.burkard.co.uk/7dayst.htm accessed November 2016) Cabezudo, B., M. Recio, J.M. SanchezLaulhe, M.D. Trigo, F.J. Toro, and F. Polvorinos (1997) Atmospheric transportation of marihuana pollen from North Africa to the southwest of Europe. Atmospheric Environment, 31, 3323-3328. Campbell, I.D., K. McDonald, M.D. Flannigan, and J. Kringayark (1999) Long-distance transport of pollen into the Arctic, Nature, 399, 29-30. Corden, J. and W. Millington (1999) A study of Quercus pollen in the Derby area, UK, Aerobiologia, 15(1), 29-37. D amato, G., L. Cecchi, S. Bonini, C. Nunes, I. Annesi Maesano, H. Behrendt, G. Liccardi, T. Popov, and P. Van Cauwenberge (2007) Allergenic pollen and pollen allergy in Europe, Allergy, 62(9), 976-990. Duhl, T.R., R. Zhang, A. Guenther, S.H. Chung, M.T. Salam, J.M. House, and E. Salathe (2013) The Simulator of the Timing and Magnitude of Pollen Season (STaMPS) model: a pollen production model for regional emission and transport modeling, Geoscientific Model Development Discussions, 6(2), 2325-2368. Efstathiou, C., S. Isukapalli, and P. Georgopoulos (2011) A mechanistic modeling system for estimating largescale emissions and transport of pollen and coallergens, Atmospheric Environment, 45, 2260-2276. Emberlin, J., S. Jones, J. Bailey, E. Caulton, J. Corden, S. Dubbels, J. Evans, N. Mcdonagh, W. Millington, J. Mullin, R. Russel, and T. Spencer (1994) Variation in the start of the grass pollen season at selected sites in the United Kingdom 1987-1992, Grana, 33(2), 94-99. Environmental Geographic Information Service from Ministry of Environment (2007) https://egis.me.go.kr/ [accessed November 2016]. Healthcare Bigdata Hub (2014) http://opendata.hira.or.kr/ [accessed November 2016]. Helbig, N., B. Vogel, H. Vogel, and F. Fiedler (2004) Numerical modelling of pollen dispersion on the regional scale, Aerobiologia, 3, 3-19. Jato, V., F.J. Rodríguez-Rajo, and M.J. Aira (2007) Use of Quercus ilex subsp. ballota phenological and pollen-production data for interpreting Quercus pollen curves, Aerobiologia, 23(2), 91-105. Käpylä, M. (1984) Diurnal variation of tree pollen in the air in Finland, Grana, 23(3), 167-176. KMA (Korea Meteorological Adminstration) (2016) http:// www.kma.go.kr/ [accessed November. 2016]. Kim, M.K. and S.W. Oh (1999) Change of causative inhalant allergens in respiratory allergic patients in Chungbuk district, Allergy, Asthma & Immunology Research, 19(5), 696-702. (in Korean with English abstract) Lee, C.-S., W.-K. Lee, J.-H. Yoon, and C.-C. Song (2006) Distribution pattern of pinus densiflora and quercus spp. stand in Korea using spatial statistics and GIS, Journal of Korean Forest Society, 95(6), 663-671. (in Korean with English abstract) Lim, Y.-K., K.R. Kim, C. Cho, M. Kim, H. Choi, M.J. Han, I. Oh, and B. Kim (2015) Development of a oak pollen emission and transport modeling framework in South Korea, Atmosphere, 25(2), 221-233. (in Korean with English abstract) Mullins, J. and Emberlin, J. (1997) Sampling pollens, Journal of Aerosol Science, 28(3), 365-370. National Institute of Meteorological Research (2015) Development of applied biometeorology model (III): Meteorological characteristics analysis affected to high airborn oak pollen concentration and improvement of forecasting model. (in Korean with English abstract) National Institute of Meteorological Research (2014) Advanced Research on Applied Meteorology (II) Advanced Research on Bio-Meteorology: improvement for pollen dispersion modeling system using WRF- CMAQ. (in Korean with English abstract) National Institute of Meteorological Research (2013) Development of applied biometeorology model (I) Development of basic technology for pollen forecast system. (in Korean with English abstract) Norris-Hill, J. (1995) The modelling of daily Poaceae pollen concentrations, Grana, 34(3), 182-188. Oh, I., Y. Kim, K.-R. Choi, and J.H. Lee (2013) Relationship between Pollen Concentration and Meteorological, J. Korean Soc. Atmos. Environ., Vol. 33, No. 1, 2017
44 오인보 김규랑 방진희 임윤규 조창범 오재원 김양호 황미경 Journal of Korean Society for Atmospheric Environment, 29(6), 780-788. (in Korean with English abstract) Oh, I., Y. Kim, K.R. Choi, M. Suzuki, and J. Lee (2012) Pollen simulations in a coastal urban area of Ulsan, Korea: Preliminary results using WRF-CMAQ model, Proceeding of the 13th International Palynological Congress and 9th International Organization of Palaeobotany Conference, Tokyo, Japan, Paper No SS28-O05. 118. Oh, J.-W. (2009) Development of Pollen Concentration Prediction Models, Journal of the Korean Medical Association, 52(6), 579-591. (in Korean with English abstract) Oh, J.-W., I.J. Kang, S.W. Kim, M.H. Kook, B.S. Kim, J.T. Cheong, and H.B. Lee (2009a) The Association between the Concentration of Pollen and Outbreak of Pollinosis in Childhood, Allergy Asthma & Respiratory Disease, 19(1), 4-11. (in Korean with English abstract) Oh, J.-W. (2007) Characteristics of allergic pollen and the pollen amount was recently changed in Korea, Korean Korean Journal of Asthma, Allergy and Clinical Immunology, 27(1), 1-7. (in Korean with English abstract) Oh, Y.-C. H.-A Kim, I.-J. Kang, J.-T. Cheong, S.-W. Kim, M.-H. Kook, B.-S. Kim, H.-B. Lee, and J.-W. Oh (2009b) Evaluation of the relationship between pollen count and the outbreak of allergic diseases, Allergy Asthma & Respiratory Disease, 19(4) 354-364. (in Korean with English abstract) Park, K.-J., H.-A. Kim, K.R. Kim, J.-W. Oh, S.-Y. Lee, and Y.-J. Choi (2008) Characteristics of Regional Distribution of Pollen Concentration in Korean Peninsula, Korean Journal of Agricultural and Forest Meteorology, 10(4), 167-176. (in Korean with English abstract) Raynor, G.S., E.C. Ogden, and J.V. Hayes (1970) Dispersion and deposition of ragweed pollen from experimental sources, Journal of Applied Meteorology, 9(6), 885-895. Saito, H., M. Ohkubo, and J. Kunitomo (2006) Pollen production of a Quercus phillyraeoides stand in Shodoshima Island, Kagawa [in Japanese], Japanese Journal of Palynology, 52(1), 47-52. Schueler, S. and K.H. Schlunzen (2006) Modeling of oak pollen dispersal on the landscape level with a mesoscale atmospheric model, Environmental Modeling & Assessment, 11, 179-194. Siljamo, P., M. Sofiev, E. Filatova, Ł. Grewling, S. Jäger, E. Khoreva, and J. Kukkonen (2013) A numerical model of birch pollen emission and dispersion in the atmosphere. Model evaluation and sensitivity analysis, International Journal of Biometeorology, 57(1), 125-136. Singer, B.D., L.H. Ziska, D.A. Frenz, D.E. Gebhard, and J.G. Straka (2005) Research note: Increasing Amb a 1 content in common ragweed (Ambrosia artemisiifolia) pollen as a function of rising atmospheric CO 2 concentration, Functional Plant Biology, 32(7), 667-670. Sofiev, M., P. Siljamo, H. Ranta, and A. Rantio-Lehtimaki (2006) Towards numerical forecasting of longrange air transport of birch pollen: theoretical considerations and a feasibility study, International Journal of Biometeorology, 50, 392-402. Sofiev, M. and K-C. Bergmann (2013) Allergenic pollen: A review of the production, release, distribution and health impacts. Edited by L. Cecchi, 1 st Ed., Springer Netherlands. 1-7. Sofiev, M., P. Siljamo, H. Ranta, T. Linkosalo, S. Jaeger, A. Rasmussen, A. Rantio-Lehtimaki, E. Severova, and J. Kukkonen (2013) A numerical model of birch pollen emission and dispersion in the atmosphere. Description of the emission module, International Journal of Biometeorology, 57, 45-58. Vogel, H., A. Pauling, and B. Vogel (2008) Numerical simulation of birch pollen dispersion with an operational weather forecast system, International Journal of Biometeorology, 52, 805-814. Zhang, R., T. Duhl, M.T. Salam, J.M. House, R.C. Flagan, E.L. Avol, F.D. Gilliland, A. Guenther, S.H. Chung, B.K. Lamb, and T.M. VanReken (2013) Development of a regional-scale pollen emission and transport modeling framework for investigating the impact of climate change on allergic airway disease, Biogeosciences Discuss, 10, 3977-4023. Zink, K., A. Pauling, M.W. Rotach, H. Vogel, P. Kaufmann, and B. Clot (2013) EMPOL 1.0: a new parameterization of pollen emission in numerical weather prediction models, Geoscientific Model Development, 6(6), 1961-1975. 한국대기환경학회지제 33 권제 1 호