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

지금까지 수행된 도시 열환경 관련 국내·외 연구사례들을 살펴보면 열배출원 분포를 먼저 산정한 후 배출원으로부터 열이 대기 중으로 확산되어 가는 현상을 파악하는 데 초점을 두고 진행되어 왔다. 실례로 외국의 주요 도시들은 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) 등의 모델이 사용되고 있다.

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

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

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

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

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

전술한 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, t_Lu and t_Lω 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 m²/s³ and (b) for ε=10^-6 m²/s³.

이러한 현상은 입자 움직임 계산에서 확률적 인자의 영향정도와 관련되어 있다. 역추적 시간이 길어지거나 바람장의 소산율이 높을수록 입자이동 계산에서 무작위 확률의 기여도가 더 커진다. 식 (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)의 사례에 대해서도 명확한 해답을 구하는 것이 가능해진다.

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

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

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

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

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

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

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