[ Original Paper ]

Journal of Korean Society for Atmospheric Environment - Vol. 33, No. 6, pp.583-592

ISSN: 1598-7132 (Print) 2383-5346 (Online)

Print publication date 30 Dec 2017

Received 25 Oct 2017 Revised 10 Nov 2017 Accepted 10 Nov 2017

DOI: https://doi.org/10.5572/KOSAE.2017.33.6.583

도시 열환경 시뮬레이션을 위한 라그랑지안 열원 역추적 기법의 연구

김석철^* ; 이주성 ; 윤정임 ; 강종화 ; 김완수

㈜볼트시뮬레이션

Study On Lagrangian Heat Source Tracking Method for Urban Thermal Environment Simulations

Seogcheol Kim^* ; Joosung Lee ; Jeongim Yun ; Jonghwa Kang ; Wansoo Kim

Boolt Simulation, Inc

Correspondence to: * Tel : +82-(0)2-3477-1963, E-mail : sckim@boolt.co.kr

Abstract

A method is proposed for locating the heat sources from temperature observations, and its applicability is investigated for urban thermal environment simulations. A Lagrangian particle dispersion model, which is originally built for simulating the pollutants spread in the air, is exploited to identify the heat sources by transporting the Lagrangian heat particles backwards in time. The urban wind fields are estimated using a diagnostic meteorological model incorporating the morphological model for the urban canopy. The proposed method is tested for the horizontally homogeneous urban boundary layer problems. The effects of the turbulence levels and the computational time on the simulation are investigated.

Keywords:

Heat source tracking method, Urban temperature simulation, Morphological model, Lagrangian particle method

1. 서 론

전 세계적인 기후변화와 함께 도시열섬과 같은 도시개발에 따른 도시기후 변화가 도시민의 쾌적한 삶을 위협하는 요소로 부각되고 있는 가운데 최근 들어서는 폭염의 빈도와 강도, 지속시간이 증가하면서 폭염취약계층의 인명피해 위험성을 더욱 가중시키고 있다. 도시 열환경 개선을 위해서는 도시 열환경에 영향을 주는 미지의 요인들을 파악해 요인별 특성에 맞는 대책을 수립·시행해야 한다. 즉, 다양한 열원이 도시 열환경에 미치는 영향을 분석하고 예측하는 기술뿐만 아니라 열환경에 영향을 주는 미지의 요인들을 파악해 요인별 특성에 맞는 대책을 수립·시행해야 한다.

지금까지 수행된 도시 열환경 관련 국내·외 연구사례들을 살펴보면 열배출원 분포를 먼저 산정한 후 배출원으로부터 열이 대기 중으로 확산되어 가는 현상을 파악하는 데 초점을 두고 진행되어 왔다. 실례로 외국의 주요 도시들은 AWS 관측자료를 이용한 정확성이 향상된 도심 온도분포 시뮬레이션 모델들을 개발해 사용하고 있는데 대만 북부지역과 일본 내 주요 도시의 열환경 연구에서는 중규모 기상현상을 모의하기 위해 미국에서 개발된 WRF 모델이 이용되었으며 (Lin et al., 2016; Kusaka, 2012), 중유럽지역의 기후변화 연구에서는 상대적으로 협소한 영역에 대해서 정교한 분석이 가능한 CFD 모델인 ENVI-met이 사용되었다 (Huttner et al., 2009). 국내에서는 국립기상연구소가 베를린 공대와 공동으로 개발한 CAS (Climate Analysis Seoul) 모델과 한국건설기술연구원이 ㈜볼트시뮬레이션과 공동개발한 ACM (Air-City Model)이 있다.

지금까지 살펴 본 모델들 외에도 대부분의 연구들에서 열환경 변화를 분석하고 예측하는 정추적 모델들이 사용되고 있는데 이러한 모델들은 주어진 열원조건을 토대로 온도분포를 예측할 뿐 열원정보를 파악하는 데에는 한계가 있다. 따라서 도시 열환경 개선 연구뿐만 아니라 관련정책을 수립 및 시행하기 위해서는 측정된 열환경의 열원정보를 파악할 수 있는 역추적 기술이 필요한 실정이다.

도시 열환경 개선에 있어 열원의 정확한 파악은 한정된 재원으로 맞춤형 대책을 시행함으로써 비용효과적인 열환경 개선을 달성할 수 있도록 한다. 예컨대 특정지역의 고온현상이 건물 면으로부터의 발열에 의한 영향이 큰 것이 분명하다면, 옥상녹화와 쿨루프 (cool roof)와 같은 지붕면 발열에 대한 저감대책에 집중하는 것이 효과적일 것이다. 아울러, 열환경 정보로부터 열원을 추적해 나가는 역추적 기술은 하나의 대안으로서 현재 발생된 현상에 기초해 객관적이고도 과학적인 토대위에서 발생원인을 파악할 수 있을 뿐만 아니라 검증과정을 통해 정확도를 향상시킬 수 있다는 장점이 있다.

역추적 방법으로는 기상자료를 토대로 오염물질을 여러 입자로 나누고 각 덩어리의 궤적을 추출해 영향지역에 대한 오염원 정보를 파악하는 역추적 모델 (Backward Trajectory Model)이 흔히 활용된다. 이러한 역추적모델에 활용되는 대표적인 확산모델은 라그랑지안 확산모델 (Lagrangian particle dispersion model)로 HYSPLIT 4 (Draxler and Hess, 1998)와 FLEXTRA (Stohl, 1999) 등의 모델이 사용되고 있다.

역추적 기술은 대기오염물질 혹은 악취물질의 배출원을 추적하거나 혹은 화제발생 시 초기 시작점 등을 찾을 때 사용된다. 본 연구에서는 도시환경에서 열원에 대한 역추적 시뮬레이션 방법을 정립하고 이를 활용하여 기후변화 대응 및 도시 열환경 개선에 적용가능한 과학적 자료제공 방안을 모색하며 시사점을 도출하고자 한다. 이를 위해 본 연구에서는 연산속도가 빠른 모델시스템을 구성하고 이를 토대로 열원정보를 파악하기 위한 역추적 시뮬레이션 실험을 수행하였다.

2. 기상장 모델

유속분포와 온도를 계산하기 위해서는 Navier-Stokes 방정식 혹은 Euler 방정식을 이용한 수치해석모델들이 일반적으로 사용된다. 그러나 이러한 모델들을 넓은 도시규모에 적용할 경우 계산용량이 컴퓨터의 처리한도를 초과하거나 과도한 계산시간이 소요되는 한계가 있다. 도시규모 기상장을 계산하는 실용적인 방법으로 건물배치 형상통계에 기반한 morphological 모델이 제시된 바 있다 (Sutherland et al., 2004).

Morphological 모델은 CFD 기반 모델들과 비교해 정확도가 낮은 반면 수퍼컴이 아닌 PC에서도 운용이 가능하고 계산속도가 매우 빠르다는 장점이 있다. 본 논문에서는 morphological model을 토대로 도시공간에 대한 기상장 모델을 구성하였다. 여기서 도시공간은 건물 등의 도시구조물 효과가 바람장에 나타나는 공간범위로서, 시가화된 지표면에서 지상으로부터 수 백 m 고도 이하의 공간에 해당한다. 도시공간을 포함하는 넓은 주변 공간, 곧 지면으로부터 먼 상층과 도시주변 전원지역에 대해서는 기상진단모델을 사용하여 기상장을 계산하였다.

2. 1 풍속분포계산

경계층 내의 연직 풍속분포 함수로 널리 사용되는 다음 식을 사용하여 기상진단모델을 구성하였다.

u z = u * κ ln z L - ψ m z L + ψ m z 0 L,

(1)

여기서, κ는 von Kármán 상수, u^*는 마찰속도, z₀는 지표면 거칠기, L은 Monin-Obukhov 길이이며 z는 지면에서부터 측정한 연직고도이다. ψ_m은 대기안정도 효과를 반영하는 함수로 대기안정도에 따라, 곧 L의 부호에 따라 함수식 형태가 달라진다. 본 연구에서는 Panofsky and Dutton (1984), Ulden and Holtslag (1985)가 제시한 함수식을 적용하였다.

마찰속도 (u^*)와 Monin-Obukhov 길이 (L)는 기상변수이기는 하나 기상관측소에서 측정되는 물리량이 아니므로 기상 관측 값으로부터 계산되어야 하는데 그 과정은 다음과 같다. Monin-Obukhov 길이는 다음과 같이 정의된다.

L = - ρ c p T r e f u * 3 κ g H,

(2)

여기서, ρ는 공기밀도, c_p는 정압비열, g는 중력가속도, T_ref는 관측온도이며 H는 현열유속 (sensible heat flux)이다. 본 연구에서는 현열유속의 계산을 위하여 De Bruin and Holtslag (1982)의 경험식을 적용하였다. 식 (2)로부터 L을 계산하기 위해서는 u^*가 필요한데, u^* 또한 미지수이므로 다음과 같이 반복계산을 수행하였다. 먼저 식 (1)에서 ψ_m=0이라 설정하고 관측고도 (z_ref)와 관측풍속 (u_ref)을 적용해 u^*를 계산한다. 이후 식 (2)을 이용해 L값을 추정하고 이를 토대로 ψ_m의 형태를 결정한다. ψ_m이 포함된 식 (1)을 이용하여 u^*값을 다시 계산한다. 다시 계산된 u^*값을 식 (2)에 적용하여 L값을 다시 계산한다. 최초 계산된 u^*값과 초기의 L값을 다시 계산된 u^* 및 L과 비교하여 오차범위 이내로 일치하면 정답이므로 연산을 종료한다. 동일하지 않을 경우 계산된 u^*를 적용하여 L값을 다시 계산하는 것으로 전술한 과정을 반복한다. 이 알고리즘은 신속하게 정답으로 수렴되어 반복계산이 10회 이내에 종료되었다.

식 (1)의 풍속분포는 수목층 내부에서는 불합리하므로 적절한 조정이 필요하다. 본 연구에서는 Cionco (1965)에 의해 제시된 캐노피 풍속분포 모형을 적용하였다. Cionco 풍속분포는 수목층 상단에서 식 (1)에 의한 풍속분포로 연속적으로 이어지게 하였다. Cionco (1965)의 캐노피 모형은 도시건물군에 대한 것도 포함되어 있으나 본 연구에서는 수목 층과 같은 전원지역 캐노피에 대해서만 적용하였고 도시건물군에 대해서는 morphological 모델을 적용하였다.

또한, 대기혼합고보다 높은 고도에서 식 (1)을 적용하는 것은 불합리하므로 대기혼합고보다 높은 고도의 풍속은 고도와 무관하게 일정한 것으로 가정하였다. 이러한 가정은 실제로 관측되는 현상과는 상당히 다를 수 있으나 본 연구의 관심현상이 지표면에 가까운 대기경계층 내부과정이므로 연구에 있어 큰 문제는 아니라고 판단된다.

난류량은 연직 온위분포로부터 추정할 수 있다. 본 연구에서는 Stull (1983)와 Cimorelli et al. (2004)에 의해 제안된 방법을 적용하여 난류량 연직분포를 계산하였다.

이상의 절차를 통하여 기상관측 지점, 예를 들면 AWS의 지면으로부터 상층에 이르는 모든 고도에서의 평균풍속과 난류량을 계산할 수 있다. 이러한 AWS 연직분포를 보간하여 시뮬레이션 대상공간의 모든 3차원 격자점에 대한 바람장 정보를 구성할 수 있다. AWS 연직분포를 3차원 격자점으로 보간하는 방식은 다양한데, 기존 연구 (Scire et al., 2000)에서는 AWS 풍속관측 고도의 수평 격자평면에서 보간하여 기상변수를 추정한 후에 식 (1)과 부합되도록 평균풍속과 난류량 연직분포를 계산하여 연직 격자점 데이터를 생성하였다. 이러한 방식은 상층풍속의 수평분포가 위치에 따라 지나치게 민감하게 변동되는 문제점이 있다고 생각된다. 계산결과를 살펴보면 수평거리가 불과 1~2 km 이격된 두 지점의 상층풍속이 수 내지 수십 배 차이가 발생하는 경우가 빈번하였다. 일반적으로 도시열섬이나 도시계획의 열환경 영향을 분석하기 위한 시뮬레이션의 수평범위는 최대 10~20 km 이내로 본 연구에서는 이 정도의 공간 범위 내에서 최상층 풍속은 동일한 것으로 가정하였다. 국내 최대도시인 서울특별시의 경우 동서방향으로 20 km, 남북방향으로 10 km 정도의 규모이다. 이 가정에 부합되도록 3차원 격자점으로 보간된 바람장 정보가 자동 조정되도록 설정하였다. 그 결과 본 연구에서 구성한 기상진단모델에서는 최종적으로 모든 격자점의 평균풍속과 난류량의 연직 분포가 식 (1)을 만족하면서 연직 최상층 격자점 풍속은 모두 동일하게 나타났다.

2. 2 도시건물효과

건물 등의 인공구조물이 밀집된 도시공간에서 기류는 구조물의 영향을 받게 된다. 이러한 도시건물 효과를 고려하기 위해 본 연구에서는 morphological 모델 (Kim et al., 2011)을 적용하였다. Morphological 모델에서는 건물형상을 개별적으로 고려하지 않고 일정 거리 이내의 모든 인접한 건물형상에 대한 통계량, 곧 morphological 인자를 이용하여 기류분포를 추정한다 (Bentham and Britter, 2003; Hanna and Britter, 2002; Macdonald, 2000; MacDonald et al., 1988). 본 연구에서는 도시공간을 100 m×100 m 격자 단위로 구분한 후에 각 격자별로 건물형상 통계량을 산출하였다. 건물형상 통계를 위해 설정한 격자크기는 도시 열환경 시뮬레이션에서 온도계산을 위해 설정한 격자 (10 m×10 m)의 100배 정도에 해당한다.

Morphological 모델에 의해 산출된 바람장은 건물효과가 나타나는 공간범위 이내, 곧 도시 캐노피 층 (urban canopy layer)에서만 유효하다. 본 연구에서는 도시지역의 경우 지상에서부터 고도 100 m에 이르는 공간을 morphological 모델 영역으로 설정하였다. 지상고도 100 m부터 200 m까지는 고도가 높아질수록 morphological 모델에 의한 바람장이 기상진단모델의 바람장으로 선형적으로 수렴하도록 설정하였다.

3. 라그랑지안 확산모델

열의 난류확산은 지표에서 가열된 공기덩어리들이 주변 공간으로 이동하는 현상이다. 라그랑지안 모델링은 공기덩어리 각각의 이동과정을 추적할 수 있는 모델링 기법이다. 기존 연구에서 열확산에 라그랑지안 모델링 기법을 적용해 온위가 보존된다는 가정으로부터 기온분포를 계산한 바 있다 (Kim and Yun, 2017). 본 연구에서도 동일한 기법을 적용하였다.

열확산에 라그랑지안 모델링 기법을 적용할 경우 열원별로 구분하여 확산경로를 추적할 수 있다. 이 점에 착안하여 본 연구에서는 열원으로부터 온도분포를 계산하거나 반대로 온도정보를 토대로 열원분포를 계산할 수 있는 정추적/역추적 라그랑지안 확산모델 알고리즘을 구현하였다.

3. 1 정추적 연산

라그랑지안 확산모델은 공간으로 퍼져나가는 열을 계산입자의 집합체로 표현하고 유속정보를 적용하여 각 계산입자의 움직임을 추적한다. 이때 유속정보, 곧 바람장은 기류를 통계적으로 표현하는 평균풍속 $u i ¯$ 와 난류성분 u_i′으로 나뉜다. 이 값들은 전술한 기상모델에서 제공된다.

계산입자의 이동은 평균풍속과 난류성분에 의한 합으로 계산된다. 평균풍속에 의한 이동은 입자에 가해지는 평균풍속을 적분하여 계산한다. 난류성분에 의한 이동과정은 Langevin 방정식에 의해 계산될 수 있는데 시간 dt 동안의 난류성분에 의한 계산입자의 속도변화는 다음 표현이 유용하다 (Kim, 2003).

d u i' = - u i' t L u i d t + C 0 ε d t ω i,

(3)

여기서, t_Lui, ε, ω_i는 각각 라그랑지안 시간척도 (Lagrangian time scale), 난류의 소산율 (dissipation rate), 난수이다. 난수의 확률밀도함수는 평균이 0이고 분산이 1인 정규분포이다. C₀는 무차원 상수인데 본 연구에서는 3을 적용하였다. 입자의 위치변화는 식 (2)를 적분하여 dt 시간 간격으로 계산된다. 이 모든 연산은 열배출과 함께 시작하여 물리적인 시간흐름과 동일한 순서대로 수행된다. 계산이 완료된 입자분포를 기반으로 온도분포를 ‘정추적’ 할 수 있다.

3. 2 역추적 연산

역추적 라그랑지안 입자확산모델은 계산입자의 과거 이동궤적을 추적한다. 이동궤적을 시간을 거슬러 끝까지 추적하면 배출원에 도달하게 되므로 계산입자의 현재위치를 토대로 배출원을 계산하게 된다. 계산입자의 현재 위치는, 곧 현재의 온도분포에 상응하므로 온도를 토대로 배출원을 추적하는 과정이다.

역추적 연산은 라그랑지안 입자모델에서 시간의 부호를 반대로 설정하면 되므로 매우 간단하게 구현할 수 있다. 본 논문에서는 바람벡터의 방향을 반대로 하여 시간의 부호를 변경하는 효과를 구현하였다.

d u i' = + u i' t L u i d t + C 0 ε d t ω i,

(4)

정추적 연산방정식 (3)과 역추적 방정식 (4)를 비교하면 수치연산의 관점에서 동일하다. 난류속도 성분 u_i가 양과 음의 값을 모두 가질 수 있는 물리량이기 때문이다. 따라서 역추적 모델링은 풍향이 정반대인 경우의 정추적 모델링이다. 그러므로 오일러리안 (Eulerian) 모델에서 시간부호를 변경하였을 때 발생할 수 있는 수치해석적 불안정성 문제는 역추적 라그랑지안 입자모델에서는 전혀 발생하지 않는다. 역추적을 진행하는 경우에도 라그랑지안 입자모델은 정추적과 마찬가지로 안정된 연산이 수행된다. 또한, 정추적 모델로써 역추적 모델링이 가능하기 때문에 역추적을 위한 별도의 코드를 만들 필요는 없다.

4. 시뮬레이션 실험 및 결과 고찰

전술한 morphological 기상진단모델과 라그랑지안 대기확산모델로부터 도시열환경에 대한 정추적/역추적 통합 시뮬레이션 시스템을 구성하였다. 시스템의 역추적 모델링 기능을 통해서 기온분포로부터 열원을 계산할 수 있다. 열원정보는 그 자체만으로 도시계획이나 도시열환경 관리 등의 목적으로 사용 가능하나 정추적 연산의 입력자료로 활용될 경우 더 정확한 기온모델링이 가능하다.

역추적을 통해 계산된 열원정보가 기존방식에 의한 것, 곧 열원모델을 통해 계산된 것보다 더 정확하다면 이를 입력자료로 적용하여 수행한 정추적 온도모델링의 정확도 또한 향상된다. 기존 열원모델을 통해 추정되는 열원정보, 특히 그 중에서도 열 배출량은 기온모델링 입력정보 가운데에서 불확실성이 높다. 운량과 같은 기존 열원모델의 주요 입력자료는 상시기상 관측과정에서 정확한 값이 확보되기가 어렵고, 기존 열원모델을 구성하는 이론 또한 현실에 대한 거친 근사식에 의존하는 경우가 많기 때문이다. 반면, 기온은 상시 관측되므로 이에 대한 역추적을 통하여 더 정확한 열원정보를 파악할 수 있다면 이를 토대로 정추적 연산을 수행하여 더욱 정확한 온도분포를 계산할 수 있다.

수평방향으로 균질한 대기경계층 조건에서 역추적 시뮬레이션 실험을 수행하였다. 대기경계층이 수평방향으로 균질하기 위해서는 지형이 평탄해야 하고 지표면 형상과 재질이 동일해야 한다. 이러한 실험 조건은 현실을 이상화한 상황이나 역추적 실험의 특징을 보다 쉽게 이해할 수 있기 때문에 연구조건으로 선정하였다.

단일 열원에서 배출된 열을 일정시간 정추적한 후 다시 같은 시간 동안 역추적하여 열원을 재구성하여 보았다. 2차원 실험공간을 x축과 z축으로 균질한 격자로 나누었다. 이때 x는 풍하방향이며 z는 연직방향이다. 좌표 원점 (0, 0) 위치에 발열체인 12 m (4층) 높이의 건물을 위치시켰다. 태양복사열에 의해 가열된 건물벽면과 건물 내 에너지 사용으로 인해 건물에서 발생한 현열은 주변 공기로 유입된다. 지상고도 10 m (z = 10)에서 풍속은 2 m/s이다. Morphological 기상진단모델에 의해 계산된 바람장의 연직분포는 그림 1과 같다. 여기서, u는 수평방향 풍속이고, σ_u와 σ_ω는 각각 수평 및 연직방향 난류속도 성분크기이며, t_Lu와 t_Lω는 각각 수평 및 연직방향 라그랑지안 시간척도이다.

Fig. 1.

The vertical profiles of meteorological variables by the morphological diagnostic flow model: u is the horizontal wind speed, σu and σω are the horizontal and vertical (z) turbulence, ε represents the dissipation rate, tLu and tLω represent the horizontal and vertical Lagrangian time scales.

Morphological 파라미터는 수평면에서 균질한데 평균 건물높이는 12 m, 배제두께 (displacement thickness)는 9 m, 지면거칠기 (z₀)는 1.5 m, 풍차단비율은 0.4, 평균건폐율은 0.3이다. 여기서 풍차단비율은 건물에 의해 차단되는 기류단면적을 morphological grid (본 연구에서 100 m×100 m) 면적으로 나눈 값이다. 지면의 평균적 알베도와 보웬비는 각각 0.16, 2.0이며, 태양고도각은 6월 10일 14시 서울에 대한 값을 적용하였는데 약 66°이다. 이러한 조건하에서 열배출모델을 통해 계산된 지표면 현열유속 (sensible heat flux)은 353 W/m²이다. 기상입력조건으로부터 산정된 마찰속도 (u^*)는 0.5 m/s이고, Monin-Obukhov 길이 (L)는 -40 m로 대기는 불안정한 상태이다. Morphological 기상진단모델에서 u^*와 L을 구하는 반복계산은 7번 만에 상대오차 0.1% 이내로 수렴하였다. 대기혼합고는 1,800 m이다.

건물 배출열로 인한 초기 온도분포는 건물위치를 중심으로 그림 2와 같이 반지름 12 m의 반원모양으로 단순화하여 모델링하였다. 초기 온도장은 라그랑지안 확산모델링에서 계산입자의 초기조건을 설정하기 위한 것으로 열원형상과 열방출조건 등에 의해 결정된다. 시뮬레이션 코드에서 초기 온도장으로부터 열원정보를 구성할 수 있으므로 초기 온도장은, 곧 열원분포에 상응한다. 라그랑지안 확산모델링에서 확산프로세스에 대한 계산은 초기장을 형성하는 입자위치에서 시작되므로 역추적을 통해 초기 온도장을 찾는 것은, 곧 열원정보를 찾는 것에 해당한다.

Fig. 2.

Initial temperature profile, an idealization of an 12 m height building structure as heat source located in the origin: air temperature within the semicircle is uniformly 1°C higher than the surrounding values.

그림 2는 초기온도 분포이다. 여기서 표현된 온도는 실제 온도와 배경온도의 차이이다. 반원 내부에서 온도는 외부보다 1℃ 더 높다. 라그랑지안 입자모델을 적용하기 위해 반원의 초기 온도장 내부는 계산입자로 균질하게 채웠다. 각 입자들은 먼저 30초 동안 정추적 연산에 의해 확산되도록 하였고 그리고 다시 최초시간까지 역추적 연산에 의해 돌아오도록 하였다. 역추적 결과가 초기온도 분포상태에 얼마나 근접한 지 살펴보았다. 계산입자는 3,653개를 사용하였다.

그림 3은 30초 정추적에 이어 30초 역추적한 결과를 제시한 것으로 모델링 시간간격 (dt)은 0.3초로 설정하였다. 확산되었던 입자들이 다시 중심점 (0, 0)으로 돌아와 온도분포가 그림 1의 초기상태와 비슷한 모양으로 형성된 것을 관찰할 수 있다.

Fig. 3.

Results of the backward dispersion modeling started from the state after the particles having been spreaded for 30 s: (a) particle positions and (b) the temperature profile. The solid line represents the boundary of the initial high temperature region.

역추적된 온도분포 형태는 시간간격 (dt)에 영향을 받지 않는다. 시간간격을 0.03, 0.003, 0.0003초로 변경하면서 실험을 반복하여도 결과는 그림 3과 거의 동일하였다. 이는 dt=0.3초가 현재의 시뮬레이션 실험에 적절한 값이기 때문이다. 시간간격 dt는 물리적인 의미를 지니지 않는 계산변수이므로 모델링 결과에 영향을 주지 않는다. 반면, 정추적 확산시간은 역추적 결과에 큰 영향을 준다.

그림 4는 정추적 확산시간을 30초에서 60초로 늘인 경우의 역추적 결과이다. 나머지 조건은 그림 3과 같다. 그림 3과 비교해보면 되돌아온 입자들이 원점을 중심으로 더 넓게 퍼져 있다. 역추적 시간이 길어질수록 역추적된 온도분포는 초기와 비교하여 더 많이 달라진다. 더 긴 시간을 거슬러 역추적을 할수록 역추적된 형태로부터 직관적으로 초기조건, 곧 열원형태를 파악하는 것이 더 어려워진다.

Fig. 4.

Results of the backward dispersion modeling started from the state after the particles having been spreaded for 60 s: (a) particle positions and (b) the temperature profile. The solid line represents the boundary of the initial high temperature region.

역추적 실험에서 또 다른 중요한 시뮬레이션 환경 변수는 소산율 (ε)이었다. 그림 3과 동일한 조건을 유지한 채 ε만 변경한 시뮬레이션의 결과를 그림 5에 제시하였다. 그림 5(a)는 고도와 상관없이 ε값이 10^-2 m²/s³ 으로 더 큰 경우이고, 그림 5(b)는 ε값이 10^-6 m²/s³으로 더 작은 경우이다. ε값이 작아질수록 역추적된 입자들이 더 집중된 분포를 보였는데 ε값이 커지면 역추적된 입자위치가 넓게 산개하여 초기열원의 형태가 잘 드러나지 않는다는 문제점이 발견되었다. 실제 대기에서 발생하는 소산율은 10^-2 m²/s³ 수준이 많기 때문에 본 연구의 역추적 기법이 현실적으로 활용될 수 있기 위해서 이는 반드시 해결되어야 할 문제이다.

Fig. 5.

Results of the backward dispersion modeling started from the state after the particles having been spreaded for 30 s: (a) temperature profiles (a) for ε=10-2 m2/s3 and (b) for ε=10-6 m2/s3.

이러한 현상은 입자 움직임 계산에서 확률적 인자의 영향정도와 관련되어 있다. 역추적 시간이 길어지거나 바람장의 소산율이 높을수록 입자이동 계산에서 무작위 확률의 기여도가 더 커진다. 식 (3)과 식 (4)로 각각 표현된 정추적 및 역추적 연산에서 우측 두 번째 항은 무작위로 입자의 위치를 변화시키는데 이 항은 정추적에서 역추적으로 전환되어도 감소되지 않기 때문에 누적효과는 총 연산시간 (정추적시간과 역추적 시간을 합친 값)에 정비례하여 증가한다. 반면, 우측 첫 번째 항은 물리적 시간의 흐름에 따라 정확히 반대의 부호로 증가 혹은 감소하기 때문에 역추적에 따른 누적효과가 영으로 수렴한다.

그림 6은 소산율이 10^-8으로 매우 낮게 설정된 실험환경에서 정추적 20분에 이어 역추적 20분을 수행한 결과이다. 그림 3과 비교하여 40배의 긴 시간에 대하여 역추적을 수행하였음에도 불구하고 초기온도 분포와 유사하게 복귀된다. 그러나 이 경우에도 역추적 시간이 1시간을 넘으면 확률적 효과가 지배적이 되므로 역추적 결과로부터 초기온도 분포가 잘 구별되지 않는 것으로 나타났다.

Fig. 6.

Temperature field of the backward dispersion modeling started from the state after the particles having been spreaded for 1,200 s. The solid line represents the boundary of the initial high temperature region.

계산입자가 무작위로 움직인다면 현재 위치에서 그 다음 위치를 정확히 예측할 수 없고 마찬가지로 그 이전의 계산위치도 알 수가 없다. 라그랑지안 입자모델에서 계산입자의 이동은 결정론적인 부분과 무작위적인 부분의 합이다. 입자의 움직임에서 무작위적인 움직임이 차지하는 비중이 클수록 입자의 다음시간 위치를 예측하는 것과 이전 시간의 위치로 복귀하는 것이 점점 더 불확실해진다. 시뮬레이션 실험의 역추적 결과에는 이러한 특성이 나타난다. 결론적으로 역추적을 통해서 정확하게 복귀할 수 있는 입자의 위치는 무작위적인 움직임을 제외한 결정론적인 움직임에 의한 것만 가능하다. 소산율이 클수록 혹은 역추적시간이 길수록 복귀된 온도분포가 초기온도 분포와 비교하여 더 크게 달라지는 까닭이 여기에 있다.

이러한 한계는 속도장이 모델링 시에 시간평균과 분산, 곧 난류량으로 축약되면서 많은 정보가 유실되기 때문에 발생하는 현상이다. 속도장은 모든 시간과 위치에서 값이 정의되는 연속함수이나 모델링에는 평균풍속과 난류량만 활용된다. 속도장 정보가 유실되기 때문에 속도장에 따라 이동하는 계산입자의 과거 위치와 미래 위치를 정확하게 결정할 수 없다.

그런데 역추적 연산을 통해서 얻고자 하는 정보는 개별 계산입자에 대한 초기위치가 아니다. 다수 입자를 동시에 역추적하면 통계량이 산출된다. 다수입자의 역추적을 통해 초기분포에 대한 확률적 기대치가 계산될 수 있다. 역추적 시간에 따라 위치추정 불확실성이 얼마나 증가되는 지는 충분한 수의 계산입자를 적용한 역추적 결과로부터 파악될 수 있다. 다수입자를 사용한 정추적의 경우에는 최종분포에 대한 확률적 기대치가 계산되는데 입자가 충분히 많을 경우 이 기대치는 온도분포로 수렴되는 것과 같다.

그림 4 또는 그림 5(a)의 사례는 전술한 방식의 역추적만으로는 명확한 해답을 얻을 수 없는 문제에 해당한다. 즉, 의미있는 정확도의 역추적이 가능한 시간이 매우 짧다. 역추적을 통하여 추정된 분포가 초기분포에 비해서 매우 넓게 펼쳐져 있기 때문이다. 지금까지 살펴본 바와 같이 이는 문제의 본질로서 역추적 연산이 완벽하게 정확히 수행된다고 할지라도 해결될 수가 없다. 이상의 시뮬레이션 실험에서는 한 가지 중요한 정보가 활용되지 않았는데, 곧 역추적 연산을 적용한 계산입자가 그림 1의 초기 온도장을 형성했던 모든 계산입자라는 사실을 활용하지 않았다.

정추적과 역추적 연산에 사용한 계산입자를 각각 정추적 입자, 역추적 입자로 구분하면 보다 명확한 논의가 가능하다. 이상의 시뮬레이션 실험에서는 정추적 입자를 역추적 입자로 적용하였다. 그 결과 정추적 계산입자가 있는 위치에 대해서만 역추적 연산을 수행한 반면 정추적 입자가 없는 공간에 대해서는 역추적 연산을 수행하지 않았다. 정추적 입자가 없는 공간에 대한 역추적은 정추적 입자가 없는 초기 공간분포를 확률적으로 제시하게 된다. 따라서 정추적 입자가 없는 공간에 대한 역추적을 추가로 활용할 경우 그림 4와 그림 5(a)의 사례에 대해서도 명확한 해답을 구하는 것이 가능해진다.

정추적 입자가 없는 공간으로부터의 역추적은 역추적 입자가 지닌 열 에너지량을 충분히 큰 음수로 설정하는 등의 방식을 통해 현재 알고리즘에 대한 수정 없이도 가능하다. 충분히 큰 음의 에너지가 할당된 계산입자가 양의 에너지를 지닌 계산입자를 만나면 모델 알고리즘에 포함된 에너지 연산과정에서 양의 에너지 계산입자는 소멸되므로 역추적 시간이 진행되면서 계산입자의 수가 감소함에 따라 연산부하가 점차 줄어들게 된다. 그러나 역추적 초기조건을 구성하는 입자배열 형태에 따라 매우 많은 수의 음의 에너지를 지닌 역추적 입자가 소요될 수 있으므로 연산시간이 크게 증가될 것으로 예상된다. 정추적 입자가 도달하지 않는 공간을 대상으로 보다 효율적인 역추적 알고리즘에 대한 연구가 필요하다.

5. 결 론

도시열환경 모델링 시스템을 구성하고 시뮬레이션 실험을 통하여 열원 역추적 문제를 분석하였다. 도시열환경 모델링 시스템은 기상진단모델과 라그랑지안 입자모델을 결합하여 구성하였다. 기상진단모델은 기상관측자료를 토대로 도시공간의 바람장을 구성하기 위한 것으로 도시건물 캐노피를 고려하기 위해 morphological 모델을 적용하였다. 라그랑지안 입자모델은 정추적과 역추적 연산을 단일 코드 내에서 통합적으로 구현하여 열원으로부터 온도를 계산하거나 반대로 온도로부터 열원을 계산할 수 있도록 하였다.

수평적으로 균질한 도시 대기경계층에 대한 역추적 실험을 수행하였다. 초기 온도장으로부터 일정시간 동안 정추적을 시행하고, 다시 같은 시간을 거슬러 역추적을 시행한 후 결과를 초기온도장과 비교하였다. 결과를 요약하면 다음과 같다.

첫째, 역추적 시간이 길어질수록 또는 바람장의 소산율이 클수록 역추적을 통해 추정할 수 있는 초기 온도장의 불확실성은 증가한다.

둘째, 위의 결과는 역추적 문제에 내재한 특성으로 역추적 알고리즘의 정확도와는 무관하다. 평균풍속과 난류량으로 바람장이 축약됨에 따라 정보가 소실되어 나타난 결과로서 단일입자의 역추적 과정에서 이러 한 불확실성을 제거할 수 있는 역추적 알고리즘은 없다.

셋째, 역추적에 충분히 많은 수의 계산입자를 동시에 적용하고 아울러 정추적 입자가 도달하지 못하는 공간에 대해서도 역추적을 시행한다면 역추적을 통해 초기 온도장을 정교하게 추정하는 것이 가능하다.

본 논문에서 다룬 열원 역추적 모델은 도시 열환경 개선을 위한 다양한 활동에 과학적이고 객관적인 열원정보를 제공하는 도구가 될 수 있다. 그러나 시작단계에 불과하므로 열원 역추적 연산알고리즘을 효율적으로 최적화하여야 하고 실제 사례에 대해 검증하는 추가적인 연구가 필요하다.

Acknowledgments

본 연구는 국토교통부/국토교통과학기술진흥원의 지원으로 수행되었습니다 (과제번호 17AUDP-B10240603).

References

Bentham, T., Britter, R., (2003), Spatially averaged flow within obstacle arrays, Atmospheric Environment, 37(15), p2037-2043. [https://doi.org/10.1016/s1352-2310(03)00123-7]
Cimorelli, A.J., Perry, S.G., Venkatram, A., Weil, J.C., Paine, R.J., Wilson, R.B., Brode, R.W., (2005), AERMOD: A dispersion model for industrial source applications. Part I: General model formulation and boundary layer characterization, Journal of Applied Meteorology, 44(5), p682-693. [https://doi.org/10.1175/jam2227.1]
De Bruin, H.A.R., Holtslag, A.A.M., (1982), A simple parameterization of the surface fluxes of sensible and latent heat during daytime compared with the Penman- Monteith concept, Journal of Applied Meteorology, 21(11), p1610-1621. [https://doi.org/10.1175/1520-0450(1982)021<1610:aspots>2.0.co;2]
Draxler, R.R., Hess, G.D., (1998), An overview of the HYSPLIT_4 modelling system for trajectories, Australian Meteorological Magazine, 47(4), p295-308.
Hanna, S.R., Britter, R.E., (2002), Wind flow and vapor cloud dispersion at industrial sites, American Institute of Chemical Engineers, New York, p47-117.
Huttner, S., Bruse, M., Dostal, P., Katzschner, A., (2009, June), Strategies for mitigating thermal heat stress in Central European cities: The project KLIMES, In The Seventh International Conference on Urban Climate (Vol. 29).
Kim, B., Lee, C., Joo, S., Ryu, K., Kim, S., You, D., Shim, W., (2011), Estimation of roughness parameters within sparse urban-like obstacle arrays, Boundary-Layer Meteorology, 139, p457-485. [https://doi.org/10.1007/s10546-011-9590-8]
Kim, S., Yun, J., (2017), Study on Urban temperature prediction method using Lagrangian particle dispersion model, Journal of Korean Society for Atmospheric Environment, 33(1), p45-53, (in Korean with English abstract). [https://doi.org/10.5572/kosae.2017.33.1.045]
Kim, S., (2003), Evaluation of one-particle stochastic Lagrangian models in horizontally homogeneous neutrallystratified atmospheric surface layer, Journal of Korean Society for Atmospheric Environment, 19(4), p397-414, (in Korean with English abstract).
Kusaka, H., Hara, M., Takane, Y., (2012), Urban climate projection by the WRF model at 3-km horizontal grid increment: dynamical downscaling and predicting heat stress in the 2070’s August for Tokyo, Osaka, and Nagoya metropolises, Journal of the Meteorological Society of Japan. Ser. II, 90, p47-63.
Lettau, H., (1969), Note on aerodynamic roughness-parameter estimation on the basis of roughness-element description, Journal of Applied Meteorology, 8(5), p828-832. [https://doi.org/10.1175/1520-0450(1969)008<0828:noarpe>2.0.co;2]
Lin, C.Y., Su, C.J., Kusaka, H., Akimoto, Y., Sheng, Y.F., Huang, Jr. C., Hsu, H.H., (2016), Impact of an improved WRF urban canopy model on diurnal air temperature simulation over northern Taiwan, Atmospheric Chemistry and Physics, 16(3), p1809-1822. [https://doi.org/10.5194/acp-16-1809-2016]
Macdonald, R.W., (2000), Modelling the mean velocity profile in the urban canopy layer, Boundary-Layer Meteorology, 97(1), p25-45. [https://doi.org/10.1023/a:1002785830512]
Panofsky, H.A., Dutton, J.A., (1984), Atmospheric turbulence: models and methods for engineering applications, Wiley.
Scire, J.S., Robe, F.R., Fernau, M.E., Yamartino, R.J., (2000), A user’s guide for the CALMET meteorological model (version5), p332.
Stohl, A., (1999), The FLEXTRA Trajectory Model Version 3.0, User’s Guide, p41.
Stull, R.B., (1983), A heat-flux history length scale for the nocturnal boundary layer, Tellus A: Dynamic Meteorology and Oceanography, 35(3), p219-230. [https://doi.org/10.3402/tellusa.v35i3.11435]
Sutherland, J., Peet, A.H., Soulsby, R.L., (2004), Evaluating the performance of morphological models, Coastal Engineering, 51(8), p917-939. [https://doi.org/10.1016/j.coastaleng.2004.07.015]
Van Ulden, A.P., Holtslag, A.A.M., (1985), Estimation of atmospheric boundary layer parameters for diffusion applications, Journal of Climate and Applied Meteorology, 24(11), p1196-1207. [https://doi.org/10.1175/1520-0450(1985)024<1196:eoablp>2.0.co;2]
Yi, C.Y., Eum, J.H., Choi, Y.J., Kim, K.R., Scherer, D., Fehrenbach, U., Kim, G.H., (2011), Development of Climate Analysis Seoul (CAS) maps based on landuse and meteorological model, Journal of the Korean Association of Geographic Information Studies, 14(1), p12-25, (in Korean with English abstract).