한반도 참나무 꽃가루 확산예측모델 개발

2.1 꽃가루 배출량

일중 꽃가루 배출량은 기상학적 요소인 바람과 상대습도 그리고 발원지 요소인 참나무 면적 등으로 구성할 수 있다(Zhang et al., 2014). 먼저 첫 번째 기상학적 요소에 의한 배출량은 로버스트 다중회귀식(robust multiple regression)을 이용한 방법이 적용되었다(NIMR, 2014). 다중회귀식에 적용된 요소별(적산온도, 풍속, 줄리안 데이, 강수량) 일최대 꽃가루 농도는 조건에 따라 최대농도가 일정한 형태를 보인다는 가정을 바탕으로 안정적인 꽃가루 관측자료가 존재하는 2007년 부터 2008년의 발생 빈도가 높은 3~6월의 전국 7개 지점(서울, 구리, 강릉, 대구, 광주, 부산, 제주) 자료를 이용하여 이동된 와이블(Weibull) 확률밀도함수 식(1)로 추정되었다.

여기서 이동된 와이블 분포의 모수는 척도(scale) 모수인 σ, 모양(shape) 모수인 c, 위치(location) 모수인 θ로 결정되며 각 요소별 결정 모수값을 Table 1에 나타내었다.

그 다음 단계인 일별 꽃가루 농도 추정을 위한 다중 회귀 모델은 종속변수로 요소별 최대 농도 추정에 사용된 로그 꽃가루 농도 중 1등급(NIMR, 2012)을 제외한 값을, 독립변수로 요소별 최대 농도 추정 값과 개별기상요소 값을 사용하여 개발하였다. 꽃가루 농도 1등급을 제외한 이유는 전체 관측 자료 중 1등급에 해당하는 비율이 80% 이상을 차지하고 있어 꽃가루 농도 범위가 너무 넓어지기 때문이다. 변수증감법으로 선택된 변수 중 다중공선성이 위배되지 않고 요소별 최대 농도 추정 값의 회귀계수가 양인 변수가 선택되었다. 선택된 독립변수로 최소절단제곱(LTS:least trimmed squares) 방법을 이용 로버스트 다중회귀식이 산정되었다(식 2).

여기서 Ta는 일평균 기온(^oC), hm는 일평균 상대습도(%), acc는 사례일까지 총 적산온도(^oC), ws는 일평균 풍속(m s^-1), rain은 일평균 강수량(mm), date는 줄리안 데이를 말하며 f()는 와이블 확률밀도함수이다.

두 번째로 꽃가루 발원지 요소는 산림이나 초목의 형태로 분포하며 이를 정확하게 산정하는 것은 불가능에 가깝다. 다양한 선행연구에서는 이들 꽃가루 발원지에 해당하는 정보를 해당 식생의 면적부분과 지표이용 자료 등을 이용하여 산정하고 있다. 본 연구에서는 우리나라 산림청 공간정보시스템에서 제공하는 임상도 자료(1:5000)를 GIS 프로그램을 이용하여 각 수종에 따른 면적을 수치모델 격자 크기(수평해상도 1.5 × 1.5 km)에 맞도록 면적비율로 계산하고, 이를 해당격자의 꽃가루 배출원으로 적용하였다. Figure 1은 각 격자별 우리나라 참나무과 수종 면적 비율을 나타낸 것으로 주로 산림이 우거진 강원산간지역과 지리산 일대에서 높은 값을 나타내고 있음을 알 수 있다. 이를 통해 꽃가루 배출량의 공간적 분포 모습을 알 수 있다. 따라서 본 연구에서 적용한 일 총 꽃가루 배출량(Pa, grains m^-3day^-1)은 식(3)과 같이 기상요소에 의한 잠재배출량에 발원지에 해당하는 참나무 수종 면적비율을 곱하여 정의하였다.

여기서 Oak는 로버스트 다중회귀식에 의한 꽃가루 배출량, γ는 참나무 꽃가루 수종 면적 비율을 나타낸다.

Table 1.

The parameters of the translated Weibull probability density functions for daily oak pollen.

Fig. 1.

Distribution of areal fraction of oak trees based on the FGIS (Forest Geographic Information System) data from the Korea Forest Service. Red circles indicate the pollen observation network.

2.2 꽃가루 배출량 플럭스 모수화

꽃가루 배출식 식(3)을 통하여 하루 중 총 꽃가루 배출량 Pa를 얻을 수 있다. 하지만, 실제 대기중으로 배출되는 실시간 꽃가루 배출량은 그 때의 기상조건에 따라 달라질 수 있으며 하루 중 배출량 변화와도 관계가 있다. 이와 관련하여 많은 선행연구에서 Helbig et al. (2004)이 제안한 마찰속도(u*)에 근거한 모수화식을 많이 사용하고 있다(Sofiv et al., 2006; Vogel et al., 2008; Efstathiou et al., 2011; Zink et al., 2012; Zhang et al., 2014). 본 연구에서는 Zhang et al. (2014)이 제안한 마찰속도에 따른 모수화식을 도입하여 배출량 플럭스 모수화를 설계하였다.

매 시간별 각 격자에서 계산되는 배출량, E_p (grains m^-3 s^-1)는 일 총 배출량 Pa (grains m^-3)와 마찰속도u* (m s^-1), 풍속에 따른 무차원 변수 K_e로 구성된다.

여기서 H_s는 참나무 수종의 평균 캐노피 높이(9.6 m, Schuler and schlünzen, 2006), C는 일단위에서 초단위로 환산되는 단위환산 요소(86400)이다.

여기서 u*_t는 먼지의 풍동실험 자료에 의해서 구해진 회귀식 결과이며(Greeley and Iversen, 1985) ρ는 공기밀도(1.161 kg m^-3), ρ_p는 참나무 꽃가루 밀도(1058 kg m^-3), g는 중력가속도(9.8 m s^-2), d_p는 참나무 꽃가루 직경(31 μm)으로 구성되어 있다.

α는 저항계수로서 U₁₀ (10 m 고도 풍속, m s^-1)과 U_10e (경험상수, 2.9 m s^-1)의 비로써 나타낸다. 즉, 꽃가루 배출량 비산과 관련된 환경적 요소로 사용할 수 있다.

꽃가루 확산 과정에서 건성 침적과정은 중력 침강속도와 공기역학적 저항, 준 층류 저항값으로 계산되며(Wesely et al., 1977) 참나무 꽃가루 직경에 영향을 받는 중력 침강 속도를 식(8)에 나타내었다.

여기서, μ는 공기점성(1.78 × 10^-4 g cm^-1s^-1, C_c는 슬립 상관계수(100)이다.

습성침적 과정은 CMAQ의 RADM (Chang et al., 1987)을 채택하여 사용하였다. 시간별 꽃가루 농도는 대기 중에 배출된 꽃가루의 양과 수평 및 연직확산정도, 그리고 대기중 침적과정에 의해서 결정된다.

2.3 UM-CMAQ-Pollen 모델링 체계

본 연구에서 사용된 꽃가루 확산예측모델의 플랫폼은 CMAQ version 4.3에서 개발된 것이다. CMAQ 모델은 미국 Environmental Protection Agency (EPA)에서 개발된 3차원 광화학 수송 모델로 오존, 미세먼지, 가스상 · 입자상 대기오염물질, 산성강하, 시정 등을 모의할 수 있다.

CMAQ 모델은 주요 전처리 프로그램과 화학 · 수송모델(CMAQ Chemical Transport Model, 이하 CCTM)로 이루어져 있고, 전처리 프로그램은 기상입력자료 전처리를 위한 Meteorology-chemistry interface processor (MCIP), 광해리율 산출 프로그램인 Photolysis rate processor (JPROC), 초기 조건을 생성하는 Initial conditions processor (ICON)과 경계조건을 생성하는 Boundary conditions processor (BCON)으로 구성되어 있다(Fig. 2). UM-CMAQ-Pollen 모델은 기상 입력자료로 현업 국지모델(LDAPS, Local Data Assimilation and Prediction System)을 MCIP에 적용하고 이를 통해 일 총배출량 산정에 필요한 일 평균 풍속, 기온, 상대습도 변수를 계산한다. 그리고 기존 CCTM 프로세서에 화학반응 관련 프로세서는 제외하고 꽃가루 배출량 입력자료 처리 및 산정 모듈이 추가되어 시간별 꽃가루 농도가 계산되도록 수정되었다(Fig. 2).

Table 2에 UM-CMAQ-Pollen 모델의 구성 및 참나무 꽃가루의 물리적 특성을 나타내었다. 꽃가루의 침적과정에서 입자의 광학직경은 매우 중요한 변수로 일반적으로 참나무 꽃가루의 경우 20~40 μm 수준을 가진다. 본 연구에서는 Zhang et al. (2014)의 연구에서 활용된 31 μm 크기를 적용하였고 밀도는 1,058 kg m^-3 (Schuler and Schlünzen, 2006)가 고려되었다. 참나무 꽃가루는 비반응성으로 고려되었으며 에어로졸 중 거대입자(coarse particle)로 처리되었다. CMAQ 모델의 격자체계는 수평 358 × 439 (1.5 km 격자간격) 연직 24층으로 구성되었으며 모델의 최하층은 0.9987 σ 레벨로 약 39 m이다.

Fig. 2.

Flowchart of the regional UM-CMAQ-Pollen modeling framework with modified pollen emission and deposition processes shaded.

Table 2.

Description of the oak pollen model.

Fig. 3.

Hourly observed oak pollen counts and accumulated percentile during April-May 2014 for Guri, KMA, and Pocheon.

2.4 사례기간 선정

우리나라 꽃가루 관측은 1999년부터 버커드 샘플러(Burkard 7-day recording volumetric spore sampler) (Burkard, 2001)를 이용하여 일별 꽃가루 농도(grains m^-3) 관측이 시작되었다. 2014년에는 포천과 기상청, 구리 세 지점(Fig. 1 참조)에서 시간별 꽃가루 농도를 관측하였으며, 본 연구에서는 이들 시간별 자료가 존재하는 2014년 참나무 꽃가루 시즌을 대상으로 꽃가루 확산예측모델 결과를 관측자료와 비교, 분석 하였다.

Figure 3은 2014년 봄철 구리, 기상청, 포천 꽃가루 관측지점에서 관측된 참나무 꽃가루농도의 시간적 변화를 나타낸 것으로 4월 초부터 농도가 급격히 증가하기 시작하여 약 한 달간 높은 농도가 빈번히 나타나고 이후 점차 감소하여 5월 하순까지 약한 농도가 지속됨을 확인할 수 있다. 서울 기상청 지점은 4월 14일~20일 동안 꽃가루 농도가 높게 나타났고 수도권 관측지점 중 가장 짧은 기간에 고농도가 집중되었다. 이 기간 중 최대농도는 4월 15일(0600 LST) 기록된 4,600 grains m^-3이다. 구리지점은 4월 14일부터 5월 1일(18일 동안)까지 누적농도 10~90%가 나타났고 4월 16일(320 grains m^-3)과 18일(320 grains m^-3) 및 30일(390 grains m^-3)에 고농도가 기록되었다.

한편 포천지점은 약 2주 동안(4월 17일~4월 29일)에 꽃가루 비산이 집중되었고 4월 19일(0500 LST)에 2,547 grains m^-3의 최고농도가 기록되었다. 본 연구의 비산(계절)은 전체 누적 참나무 꽃가루 농도의 1~99%가 나타나는 시기로 정의하였다(Emberlin et al., 1994). 이 기준에 따른 2014년 참나무 꽃가루 비산계절은 4월 11일~5월 10일이며, 수치모델의 Spin up time 48시간을 고려하여 4월 13일~5월 10일 동안 UM-CMAQPollen 꽃가루 수치모의 결과와 관측값의 비교를 통해 꽃가루 확산예측모델을 평가하였다.

3.1 꽃가루 계절 농도

꽃가루 확산예측모델의 정확도에 있어서 배출량 산정부분은 매우 중요한 요소이다. Figure 4는 사례기간(4월 13일~5월 10일) 동안 꽃가루 관측지점의 평균 참나무 꽃가루 농도와 관측소 위치에 해당하는 꽃가루 확산예측모델 결과를 나타낸 것이다. 우리나라 꽃가루 관측소 지점 중 강릉은 관측자료의 오류로 분석대상에서 제외하였다.

2014년 참나무 꽃가루 계절 동안 관측된 대전의 평균 꽃가루 농도는 149 grains m^-3로 가장 높은 농도를 나타내었고 제주가 13 grains m^-3으로 가장 낮은 농도수준을 보였다. 꽃가루 확산예측모델에 의한 결과로 대전이 가장 높은 228 grains m^-3을 나타냈고 제주가 가장 낮은 14 grains m^-3 수준을 보였다. 꽃가루 관측농도가 비교적 높은 지역인 포천, 기상청, 전주, 대전 지점의 수치모의 결과가 관측값에 비해 약 1.7배 과대모의 하는 경향이 나타났으며 구리지점을 제외한 저농도 지역인 제주, 광주, 대전, 서울, 부산, 포항 지점들의 평균농도와 수치모델에 의한 계산값은 수치모델이 약 1.1배 과대모의 하는 모습을 보였다.

이와 같이 수치모델이 관측값보다 다소 과대모의하는 결과는 로버스트 다중회귀식에 의한 일 총배출량 산정방식에서 각 요소산출 모수화에 있어 최대농도 추정 방식을 적용하였기 때문이다. Zhang et al. (2014)의 연구에서는 참나무 꽃가루 수치모의 결과와 관측결과가 5배 이상 차이를 보였는데, 배출원에 해당하는 수종면적 정보의 불확실성, 꽃가루 침적속도의 과대 추정에 따른 장거리 수송 효과의 약화 등을 원인으로 제시하였고, 전반적인 배출량에 대한 많은 불확실성을 제기하였다. 본 연구결과에서 나타난 참나무 꽃가루 계절 동안 관측된 평균농도와 수치모의 결과의 정확도가 선행연구결과보다 다소 개선된 모습을 보이고 있어 수치모델에 적용된 참나무 꽃가루 배출량 산정방식이 유용한 것으로 판단된다.

Fig. 4.

Average concentrations of oak pollen by observation and simulation from 13 April to 10 May 2014.

3.2 꽃가루 일별 농도

앞서 참나무 꽃가루 계절 평균 농도 비교를 통해 모델에 적용된 배출량 산정식의 대략적 평가를 수행하였다. 보다 자세한 비교를 위해서 우리나라 꽃가루 관측망 자료 중 시간별 자료가 존재하는 구리, 포천, 기상청 세 지점과 가장 높은 평균농도를 보인 대전과 대구, 부산지역을 대상으로 사례기간 동안 관측값과 수치모델결과의 일별 평균농도를 비교하였다(Fig. 5).

최대 평균농도가 나타난 대전지점의 일별 꽃가루 농도는 4월 18일에 최대값이 관측되었으나, 이 기간 수치모델은 아주 낮은 농도를 모사하고 있었고, 비교적 낮은 농도를 유지했던 4월 28일과 29일에는 매우 높은 농도값을 모의하였다. 4월 28일과 29일은 약한 강수가 발생한 기간으로 수치모델의 고농도 경향은 대전뿐만 아니라 구리, 기상청, 포천지역에서도 동일하게 나타났다. 이러한 결과는 꽃가루 배출량 산정식(식 2)에서 강수 현상시 꽃가루 발생을 충분히 억제할 수 있도록 설계되지 않았기 때문이다. 즉, 강수시에도 풍속과 적산온도 변수 등에 의한 꽃가루 배출효과가 여전히 적용되고 있다. 물론 꽃가루 확산예측모델 내에서 고려되고 있는 강수율에 따른 습성침적(RADM) 과정에서 강수시 비산된 꽃가루를 제거하도록 설계되어 있지만, 강수율에 따른 함수로 실제 강수현상을 기상 수치모델(UM LDAPS)이 격자단위로 정확하게 재현하지 못하는 한계로 인해 오차가 크게 발생할 수 있는 여지가 있다. 따라서 향후 강수시에 꽃가루 배출을 억제할 수 있도록 산정식의 보완이 이루어져야 하며, 기상 수치모델의 예측 기상요소를 꽃가루 확산예측모델에 적용할 때 오차를 보완할 수 있는 방안을 마련해야 할 것으로 판단된다.

구리 지역의 경우 대체적으로 관측값에 비해 모델 값이 두드러지게 과대모의 하는 것으로 나타났는데, 이는 모델값의 과대모의라기보다 구리지점 관측소의 설치 위치에 의한 것으로 판단된다. 즉, 구리지점 관측값이 수도권 다른 지점(서울, 기상청, 포천)에 비해 낮은 값을 나타내는 것은 측정소 위치가 건물 옥상 조경지대에 위치하고 주변 건물에 의한 차폐의 영향이 크게 작용하기 때문이다.

대구와 부산의 경우는 비록 저농도 지역이지만, 농도값의 상승과 하락패턴에서 모델값이 관측값과 비교적 잘 일치하는 것을 알 수 있었다. 이것은 이 지역이 꽃가루 측정소 주변에 꽃가루 배출원에 해당하는 참나무 면적률이 비교적 작고 배출원에서 다소 떨어진 곳에 있기 때문인 것으로 판단된다. 즉, 오차가 크게 작용하는 측정소 주변 국지적 배출효과보다는 꽃가루 배출량 산정에 비교적 거시적으로 적용되는 바람에 의한 꽃가루 수송의 영향으로 나타난다고 볼 수 있다.

기상청과 포천지역 일별 꽃가루 농도 변화 패턴을 살펴보면 관측값에 비해 수치모의된 꽃가루 농도의 최대치가 나타나는 시간이 지연되는 현상을 보였다. 이와 관련하여 보다 자세히 분석하기 위해 기상청, 포천, 구리 세 지점에 대해서 참나무 꽃가루 농도에 대한 시간별 농도 변화를 살펴보았다.

3.3 꽃가루 시간별 농도

먼저 사례기간의 기상현상을 파악하기 위해 Fig. 6에 포천지역 시간별 꽃가루 농도와 바람정보를 나타내었다. 2014년 봄철 포천에서 측정된 참나무 꽃가루 농도는 크게 두 가지 패턴으로 구분될 수 있다. 첫 번째 패턴은 4월 27일 이전까지 발생한 참나무 꽃가루 농도 패턴으로 북서풍계열의 바람이 최대값을 가지는 1200 LST에 높은 농도가 나타났다. 두 번째는 4월 27일에서 30일 사이에 나타난 패턴으로 기압골의 영향으로 약한 강수와 더불어 동풍계열의 바람이 다소 강하게 나타남에 따라 강원도지역에서 배출된 꽃가루가 수송된 경우이다.

각 패턴사례일인 2014년 4월 21일과 29일에 대하여 시간에 따른 공간적 꽃가루 농도 분포를 Fig. 7에 나타냈다. 먼저 21일 사례를 살펴보면 경기도 일대 오전 0010 LST에는 약한 풍계로 꽃가루 농도가 거의 나타나지 않고 있다가 북서풍계 바람이 점차 발달하는 1300 LST에 들어서 참나무 배출원 지역을 중심으로 농도수준이 올라감을 알 수 있었다. 북서풍계열의 바람이 점차 강해지면서 1600 LST에서는 배출원 지역에서 풍하측으로 다소 확산되는 모습을 보였고, 1900 LST에 들어서 꽃가루 농도가 최대로 증가하는 모습을 확인할 수 있었다. 이는 수도권을 중심으로 꽃가루 농도 예측시, 북서풍계열의 바람이 불 때는 북한지역의 참나무 배출원에 대한 정보가 중요하다는 것을 시사한다. 29일 사례는 저기압이 우리나라 남부지역을 지나감에 따라 북동풍계열의 바람이 주요하게 작용한 것으로 꽃가루 농도 역시 바람을 따라 수송되어지는 모습이 뚜렷이 나타났다. 특히, 수도권 지역 꽃가루 농도 상승의 주요 원인은 가평일대의 참나무수림의 배출원에서 발생한 꽃가루가 바람을 타고 수송되기 때문이다.

보다 자세한 분석을 위하여 기상청, 포천, 구리 세지점에 대해서 참나무 꽃가루 농도에 대한 시간별 농도 변화를 살펴보았다(Fig. 8). 기상청 지점의 고농도 사례일인 4월 15일의 경우 최대농도가 0600 LST에 4,600 grains m^-3로 나타났지만 수치모델은 고농도현상을 모의하지 못하고 있었다. 두 번째 고농도를 보인 4월 19일의 경우, 관측값은 0600 LST에 1,760 grains m^-3으로 나타났지만 수치모델의 고농도는 1800 LST에 1,056 grains m^-3를 보여 고농도가 나타나는 시간이 12시간 지연된 모습으로 나타났다.

이러한 수치모델에 의한 고농도 꽃가루 발생의 시간지연 현상은 포천과 구리 지역에서도 약 5시간 지연된 모습으로 나타났는데 특히, 4월 21일에서 26일 사이에 뚜렷하게 나타났다. 이 기간의 기상상황은 주간에 약한 북서풍계열이 나타났고 야간에는 거의 무풍의 형태로 나타났으며 일 최고 22 ^oC를 보이는 맑은 날이었다(Fig. 6). 참나무에서 생산된 꽃가루의 양은 비산하기 전에 가장 많고 점차 그 양이 줄어들게 된다. 즉, 수종에서 발생하는 꽃가루 배출량은 비산시간이 늘어남에 따라 점차 줄어드는 효과가 있지만, 본 연구에서는 사용하고 있는 시간별 꽃가루 배출량 모수화식에서는 이를 반영하고 있지 못한 실정이다. 따라서 수치모델에 의한 꽃가루 농도의 최대치가 실제 관측값에 비해 지연되는 모습은 이러한 영향이 반영된 것으로 판단된다. 또한 관측 최대값이 나타난 시점에서 모델값이 저평가 되는 원인으로 북서풍계열의 바람은 경기 일대에서는 풍하측에 해당하며 북한지역의 참나무 배출원에서 발생한 꽃가루가 영향을 미칠 수 있다는 점이다. 본 연구에서는 남한지역 임상도 자료만을 사용하여 꽃가루 배출원으로 적용하였기 때문에 이 부분에서 북한지역 꽃가루 배출량 수송효과가 수치모델에 반영되지 못한 오류로 판단된다. 향후 북한지역 꽃가루 배출원 부분이 개선되어야 할 여지를 뚜렷이 보여주는 사례라 할 수 있다.

구리지역 4월 28일에서 29일의 경우 강수현상이 나타난 날로 실제 꽃가루 관측 농도값은 0에 가깝게 나타났지만 수치모델은 29일 최대농도 1,380 grains m^-3의 수준을 보여 많은 오차를 나타내고 있었다. 하지만 포천지역은 같은 사례일임에도 불구하고 비교적 관측값과 모델값이 유사하게 나타났다. 이러한 결과는 포천지역의 경우 구리 지점에 비해 강수현상이 28일 오전에 국한되어 실제 28일 오후부터 29일 동안에는 꽃가루 발생이 가능했었기 때문인 것으로 판단된다.

또한 수치모델결과에서는 기상청 지점의 강수에도 불구하고 높은 꽃가루농도가 모의되었는데, 이것은 격자점에 사용된 기상 입력 자료인 UM 기상장의 강수량 예보 정보가 실제 관측과 정확히 일치하지 않아 습성침적에 의한 제거효과가 제대로 반영되지 못한 것과 꽃가루 배출량 산정식(식 2)에서 강수시에 꽃가루 배출을 완전히 억제하지 않는 것에 기인한다. 즉, 실제 기상청 지점에서는 강수현상이 나타났지만, 기상모델이 예측한 기상청 지점에서는 강수현상이 없고 강한 풍속이 나타나 꽃가루 배출량 산정식에 의해 배출량이 높은 값으로 나타나기 때문이다. 따라서 앞서 일별 농도 분석결과에서 나타난 바와 같이 꽃가루 확산예측모델의 정확도 개선에 있어서 기상수치모델이 가질 수 있는 한계점을 이해하고 이를 보완하여 적용할 수 있는 방안이 추가되어야한다.

Fig. 5.

Observed and simulated daily concentration of oak pollen and precipitation during the flowering season in 2014 at Guri, Daejeon, KMA, Pocheon, Daegu, and Busan observation sites.

Fig. 6.

Hourly variation of oak pollen concentration and meteorological values at the Pocheon observational site for 21~30 April 2014.

Fig. 7a.

Simulated oak pollen concentration with wind vector overlays on 21 April 2014.

Fig. 7b.

Same as Fig. 7a except 29 April 2014.

Fig. 8.

Observed and simulated hourly oak pollen concentration during the flowering season in 2014 at the KMA, Pocheon, and Guri observational sites.