[ Article ]

Atmosphere - Vol. 26, No. 3, pp.473-482

ISSN: 1598-3560 (Print) 2288-3266 (Online)

Print publication date Sep 2016

Received 18 Jul 2016 Revised 06 Sep 2016 Accepted 13 Sep 2016

DOI: https://doi.org/10.14191/Atmos.2016.26.3.473

비균질 도시 지표에서 측정된 에디 공분산 난류 플럭스의 불확실성 분석: 좌표계 편향 영향

이두일 ; 이재형 ; 이상현^*

공주대학교 자연과학대학 대기과학과

Uncertainty Analysis of the Eddy-Covariance Turbulent Fluxes Measured over a Heterogeneous Urban Area: A Coordinate Tilt Impact

Doo-Il Lee ; Jae-Hyeong Lee ; Sang-Hyun Lee^*

Department of Atmospheric Science, Kongju National University, Gongju, Korea

Correspondence to: ^* Sang-Hyun Lee, Department of Atmospheric Science, Kongju National University, 56 Gongjudaehak-ro, Gongju-si, Chungcheongnam-do 32588, Korea Phone: +82-41-850-8526, Fax: +82-41-856-8527 E-mail: sanghyun@kongju.ac.kr

Abstract

An accurate determination of turbulent fluxes over an urban area is a challenging task due to its morphological diversity and associated flow complexity. In this study, an eddy covariance (EC) method is applied over a highly heterogeneous urban area in a small city (Gongju), South Korea to investigate the quantitative influence of ‘coordinate tilt’ in determining the turbulent fluxes of sensible heat, latent heat, momentum, and carbon dioxide mass. Two widely-used coordinate transform methods are adopted and applied to eight directional sections centered on the site to analyze a 1-year period EC measurement obtained from the urban site: double rotation (DR) and planar fit (PF) transform. The results show that mean streamline planes determined by the PF method are distinguished from the sections, representing morphological heterogeneity of the site. The sectional pitch angles determined by the DR method also compare well with those in the PF method. Both the PF and DR methods show large variabilities in the determined streamline planes at each directional section, implying that flow patterns may form in a complicate way due to the surface heterogeneity. Resulting relative differences of the turbulent fluxes, defined by (F_DR-F_PF)/F_DR, are found on average +13% in sensible heat flux, +21% in latent heat flux, +37% in momentum flux, and +26% in carbon dioxide mass flux, which are larger values than those reported previously for fairly homogeneous natural sites. The fractional differences depend significantly on wind direction, showing larger differences in northerly winds at the measurement site. It is also found that the relative fractional differences are negatively correlated with the mean wind speed at both stable/unstable atmospheric conditions. These results imply that EC turbulent fluxes determined over heterogeneous urban areas should be carefully interpreted with considering the uncertainty due to ‘coordinate tilt’ effect in their applications.

Keywords:

Coordinate transformation, double rotation, eddy covariance, heterogeneous urban surface, planar fit

1. 서 론

지표면과 대기 사이의 난류에 의한 열, 수분, 그리고 운동량 교환 과정의 이해는 대기 경계층의 발달과 국지적 미기후 특성을 규명하는 데 중요한 역할을 한다. 에디 공분산(eddy covariance) 방법은 대기 지표층(atmospheric surface layer)에서의 이들 난류 교환을 직접 측정하는 방법으로, 다양한 지표면에 적용되어 미기상학적 특성을 분석하는 데 널리 이용되고 있다(e.g., Baldocchi et al., 2001; Grimmond et al., 2004; Mestayer et al., 2005; Park et al., 2013). 이 방법은 측정된 고주파 바람 성분과 스칼라 성분(기온, 수증기, 이산화탄소)의 난류 공분산으로 대기 지표층 난류 플럭스를 산정한다(Baldocchi et al., 1988; Aubinet et al., 2012). 지형이 평편하고 지면 피복이 균질한 지역에서는 난류의 시간적인 정상성(stationarity)과 공간적인 균질성(homogeneity)을 쉽게 확보할 수 있고 연직 난류 교환이 수평 난류 교환에 비해 우세하여 에디 공분산 방법을 적용하기에 이상적인 조건이 된다(e.g., Zeweldi et al., 2010). 도시 지역의 경우에는 형태적 균질성이 확보되는 지역을 대상으로 에디 공분산 방법을 적용하여 대기-지표 난류 교환 특성 연구를 수행하고 있다(e.g., Grimmond et al., 2004; Velasco et al., 2009; Kotthaus and Grimmond, 2014).

에디 공분산 측정 자료를 이용하여 정확한 난류 플럭스를 산정하기 위해서는 엄격한 품질 관리 과정이 필요하다(Foken and Wichura, 1996; Vickers and Mahrt, 1997; Ruppert et al., 2006). 특히, 경사 지형이나 비균질 지면피복 지역에서 측정된 에디 공분산 자료는 난류 플럭스 산정 과정에서 평균 유선면을 이용한 좌표계 편향(coordinate tilt) 효과를 적절히 보정하는 것이 중요하다(e.g., Sun, 2007; Yuan et al., 2007, 2011; Lee et al., 2013; Li et al., 2013; Kim et al., 2014; Lee, 2015b; Ross and Grant, 2015). 유선 좌표계(natural wind coordinate)는 지형이나 지면 피복의 비균질성에 기인하는 국지 순환의 영향을 제거하기에 유용한 방법으로 제안되었다(Tanner and Thurtell, 1969). McMillen (1988)은 측기 좌표계를 기준으로 측정된 바람의 수평 및 연직 평균 성분 ( $v ̅$ 와 $w ̅$ )이 0이 되도록 하는 이중 회전(double rotation: DR) 좌표변환을 제시하였고, 비균질 자연 지표(숲 지역)에서 이 방법의 적용이 레이놀즈 응력, 난류 현열/잠열 플럭스의 산정 정확도 향상을 가져올 수 있음을 보였다. Wilczak et al. (2001)는 장기간 측정 자료로부터 평균 유선면을 찾는 좌표변환 방법 (planar fit: PF)을 제안하였으며, DR 방법의 과회전(over-rotation) 오류를 최소화하여 난류 플럭스를 현실적으로 산정하고자 하였다(Lee et al., 2004).

도시 지표에서 에디 공분산 방법을 이용하여 대기 지표층 난류 플럭스를 산출하는 경우 좌표변환에 의한 영향을 정량적으로 파악하는 것은 대상 지표면에서의 난류 교환의 정량적 이해뿐 만 아니라 독립된 플랫폼(e.g., 수치 모형, 위성, 섬광계)에서 얻어진 동일 물리량에 대한 정량적 상호 검증을 위해서도 매우 중요하다. 본 연구에서는 비균질 도시 지표로 대표되는 지역에서 측정된 에디 공분산 자료를 이용하여 회전 좌표변환이 난류 현열/잠열/운동량/이산화탄소 플럭스 산정에 미치는 영향을 정량적으로 평가하고자 한다. 이를 위해 에디 공분산 방법에서 널리 활용되고 있는 DR 방법과 PF 좌표변환 방법을 적용하여 난류 플럭스를 산정하고, 이를 통해 비균질 도시 지표에서 좌표계 편향 효과가 난류 플럭스 산정에 미치는 불확실성을 파악할 것이다.

2. 회전 좌표변환 방법

2.1 DR 좌표변환(McMillen, 1988)

측기 좌표계는 동-서(x), 남-북(y), 그리고 수평면에 대한 연직(z) 방향으로 구성되는 직교 좌표계(orthogonal coordinate)를 따른다. DR 회전 좌표변환은 난류 플럭스 분석 시간 구간(매 30분)에서 평균 수평 유선면 방향과 측기 좌표계의 x축이 이루는 각(yaw angle: γ)과 평균 유선면과 측기 좌표의 z축이 이루는 각(pitch angle: α)을 이용하여 측정된 바람 성분을 순서대로 회전 좌표변환을 수행한다(Tanner and Thurtell, 1969; McMillen, 1988). 1차 회전 좌표변환은 γ를 이용하여 x축이 평균 유선면에 나란하도록 조정한다( $v ̅$ =0):

u 1 v 1 w 1 = cos γ sin γ 0 - sin γ cos γ 0 001 u 0 v 0 w 0,

(1)

여기서 u₀, v₀, w₀는 측정된 바람 성분이며, u₁, v₁, w₁는 회전 좌표변환 후 바람 성분이다. γ는 매 난류 플럭스 분석 시간 구간에서 평균된 수평 바람 성분으로부터 다음의 식을 이용하여 계산할 수 있다:

γ = ⁡ tan - 1 v 0 ̅ u 0 ̅ .

(2)

2차 회전 좌표변환은 1차 회전 좌표변환을 수행한 바람 성분을 α를 이용하여 측기 좌표의 z축을 평균 유선면에 수직인 방향으로 조정한다( $w ̅$ =0):

u 2 v 2 w 2 = cos α 0 sin α 010 - sin α 0 cos α u 1 v 1 w 1,

(3)

여기서 u₂, v₂, w₂는 2차 회전 후 바람 성분을 나타내며, α는 다음과 같이 계산한다:

α = ⁡ tan - 1 w 1 ̅ u 1 ̅ .

(4)

2.2 PF 좌표변환(Wilczak et al., 2001)

PF 좌표변환은 측정된 3차원 바람 성분으로부터 평균 유선면을 결정한 후, 연직축을 평균 유선면에 수직이 되도록 회전 좌표변환을 수행하고 수평면은 평균 유선면과 평행하게 다시 수평 회전 좌표변환을 수행한다(Wilczak et al., 2001). 평균 유선면은 지점의 흐름 특성을 잘 반영할 수 있도록 수 주 이상의 장기간 자료를 이용하여 다음의 다중 선형 회귀식을 적용하여 결정한다:

w 0 ̅ = b 0 + b 1 u 0 ̅ + b 2 v 0 ̅,

(5)

여기서 $u 0 ̅, v 0 ̅, w 0 ̅$ 는 측정된 평균 바람 성분이며, b₀, b₁, 는 다중 선형 회귀식의 계수로 연직 바람 성분의 측기 변위(instrumental offset)와 연직 바람 성분에 대한 수평 바람 성분들의 기여를 의미한다. b₁과 b₂로부터 평균 유선면의 기울기 α (pitch angle)와 β (roll angle)는 다음과 같이 계산된다:

α = ⁡ sin - 1 - b 1 b 12 + b 22 + 1,

(6)

β = ⁡ tan - 1 b 2 .

(7)

1차 회전 좌표변환은 계산된 α와 β를 이용하여 다음의 행렬식을 따라 수행한다:

u 1 v 1 w 1 = cos α sinαsinβ - sin αcosβ 0 cosβ sinβ sin α - cosαsinβ cos αcosβ u 0 v 0 w 0 - b 0,

(8)

여기서 연직 바람 성분에 b₀를 적용하여 측기 변위를 보정한다.

2차 회전 좌표변환은 1차 회전된 바람성분의 수평 성분이 유선 좌표계(natural wind coordinate)를 따르도록 조정한다( $v ̅ 2$ =0):

u 2 v 2 w 2 = cos γ sinγ 0 - sinγ cosγ 0 001 u 1 v 1 w 1,

(9)

γ = ⁡ tan - 1 v 1 ̅ u 1 ̅ .

(10)

여기서 u₂, v₂, w₂는 PF 회전 좌표변환이 수행된 바람 성분을 의미하며, 2차 회전 좌표변환은 매 분석 구간에서 반복 수행한다. 마지막으로 플럭스 계산 전에 분석 구간별 연직 바람 성분( $w ̅ 2$ )이 0이 되도록 조정하여 난류 플럭스 계산에 잔류 평균 연직 성분이 영향을 미치지 않도록 한다. 본 연구에서 잔류 $w ̅$ (평균 0.001 m s^-1, 최대 ± 0.3m s^-1)의 제거는 난류 플럭스 산정에 거의 영향을 미치지 않았다(현열: ~0.25W m^-2 이내).

3. 난류 관측 및 자료 처리

3.1 측정 지역과 관측 기기

비균질 도시 지역에서의 난류 측정을 위해 국립공주대학교(공주시) 내 건물 옥상(36.4717^oN, 127.1426^oE)에 높이 7 m의 기상 타워(~33 m a.g.l.)를 설치하여, 2013년 10월부터 연속적으로 난류 측정을 수행하고 있다(Fig. 1a). 기상 타워에서는 3차원 초음파 풍속계(CSAT3, Campbell Scientific Inc., USA)와 적외선 기체분석기(LI-7500A, LI-COR Inc., USA)를 이용하여 3차원 바람 성분(u, v, w), 수증기(ρ_H2O)와 이산화탄소(ρ_CO₂)의 질량 농도를 20 Hz로 측정한다. 이와 함께 온도/습도계(HMP155, Väisäla, Finland), 복사계(CNR4, Campbell Scientific Inc., USA), 강수량계(WDR-205, Wedean, South Korea)를 이용하여 주요 기상 요소들을 연속적으로 측정하고, 자료 저장 장치(CR3000, Campbell Scientific Inc., USA)를 통해 저장한다. 측정지점은 주변에 아파트 지역, 소규모 상업 지역, 도로, 자연 식생 지역이 다양하게 분포하고 있고, 방위에 따라 지면피복 구성이 뚜렷이 구별되는 비균질 도시 지역에 해당한다(Fig. 1a). 북쪽 지역은 건물과 도로와 같은 인공 지표 지역이 높은 비율을 차지하고 있고, 남쪽 지역은 자연 지표의 비율이 상대적으로 높게 구성되어 있다(Table 1).

Fig. 1.

(a) Spatial photographic view of a micrometeorological tower located in the Kongju National University (Source: Google Earth). Topography is contoured with an interval of 10 m using the digital elevation data (version 2.0) from National Geographic Information Institute (NGII), Ministry of Land, Infrastructure and Transport (MLIT) of Korea. Eight directional sections (black lines) are divided by an interval of 45o (S1: 0~45o, S2: 45~90o, S3: 90~135o, S4: 135~180o, S5: 180~225o, S6: 225~270o, S7: 270~315o, S8: 315~360o). (b) Sectional distributions of topography of four paths crossing the center of measurement site (dashed lines in Fig. 1a). The measurement location is denoted by an open circle.

Table 1.

Land-cover fractions for eight directional sections within a radius of 500 m.

3.2 품질 관리와 자료 분석

기상 타워에서 측정된 에디 공분산 자료를 통한 난류 플럭스 분석을 위해 일반적으로 널리 적용되는 품질 관리 기준을 따랐다(Vickers and Mahrt, 1997; Hong and Kim, 2002; Foken et al., 2004). 기상 요소에 따라 물리적 허용 한계값(|U| ≤ 30 m s^-1, |w| ≤ 10 m s^-1, |T| ≤ 40^oC, ρ_H2O ≤ 50 g m^-3, ρ_CO₂ ≤ 1000 mg m^-3)을 설정하여 범위를 벗어나는 값은 제거하였으며, 스파이크(spike) 제거는 |x_i - $x i ̅$ | > D × δ_i를 따라 수행하였다(Schmid et al., 2000). 여기서 $x i ̅$ , δ_i, D는 각각 시계열의 정해진 분석 구간의 평균값, 표준 편차, 그리고 필터 계수를 나타내며, 필터 계수는 D = 3.5 + 0.3k로 반복 회수(k)에 따라 증가시키면서 3회 반복 수행하였다. 분석 구간은 5분(Vickers and Mahrt, 1997)과 15분(Schmid et al., 2000)을 각각 적용한 결과, 측정 지점의 경우 분석 구간 5분을 적용하는 경우에 보다 효율적인 스파이크 제거가 가능하였으며, 이를 품질 검사 과정에 적용하였다.

본 연구에는 기상 타워에서 측정된 1년간(2014년 1월 1일~12월 31일)의 에디 공분산 자료를 분석 대상으로 하였다. 대상 자료 중 강수일 자료(89일, 24%)와 품질 관리 과정을 통해 불량 자료로 판별된 자료(7%)를 분석에서 제외하였으며, 비현실적인 회전 좌표변환 효과를 제거하기 위해 풍속이 0.5 m s^-1 이하인 자료(16%)도 분석에서 제외하였다(Lee et al., 2004; Lee, 2015b). 측정 지역의 비균질 지표 특성을 구분하기 위하여 기상 타워를 중심으로 풍향에 따라 8 방위로 나누어(Fig. 1) 분석을 수행하였으며, 두 회전 좌표변환 방법을 적용하여 각 방위별 난류 현열/잠열/운동량/이산화탄소 플럭스를 매 30분 간격으로 선형 경향 제거(linear detrending) 방법을 적용하여 계산하였다. 분석 기간에 대해 블록 평균 제거(mean removal) 방법은 선형 경향 제거 방법에 비해 평균 ~5% 높은 현열 플럭스 값을 나타내었다. 난류 현열/잠열/이산화탄소 플럭스는 WPL (Webb-Pearman-Leuning) 밀도 섭동(density fluctuation) 효과를 보정하여 산정하였다(Webb et al., 1980; Tsukamoto and Kondo, 2008). 운동량 플럭스(τ)는 초음파 풍속계에서 측정된 3차원 바람 성분(u, v, w)을 이용하여 x와 y 방향 공분산( $u'w' ̅$ 와 $v'w' ̅$ )을 계산한 후 이들 성분의 벡터 합으로 산정하였다( $τ = ρ u'w' ̅ 2 + v'w' ̅ 2$ ). 대기 안정도는 회전 변환 후 산정한 난류 현열속과 운동량 속을 이용하여 계산한 오브코브 길이(Obukhov length: L)를 적용하여 구분하였다.

4. 결 과

4.1 측정 지점의 평균 유선면 특성

분석 기간 동안 측정 지역의 풍향별 발생 빈도는 남풍 계열(S4, S5, S6)에서는 16~20%의 높은 빈도를, 북풍 계열(S1, S2, S8)에서는 10% 이내의 낮은 빈도를 나타내었다(Fig. 2). 최저 빈도는 S2 영역으로 분석 기간 중 약 2.3%를 차지하였으며, 최대 빈도는 S4 영역으로 19.6%를 나타내었다. 경사진 지형이나 국지 순환이 발달하는 지역의 경우 특정 풍향의 바람이 주/야간에 지배적으로 나타나는 특징을 보이는데(e.g., Yuan et al., 2007), 측정 지점에서는 방위에 따른 바람의 발생 빈도가 특정 방향에서 두드러지게 나타나지는 않았다(Fig. 2). 측정 지점의 풍향별 발생 빈도는 주/야간 시간대에 유사하게 나타나고 있어, 남동쪽(S3과 S4)에 위치한 낮은 산(120 m)이나 남서쪽(S5와 S6)에 위치한 금강에 의한 국지 순환 발달 가능성은 크지 않음을 알 수 있다.

Fig. 2.

Occurrence frequency of eight directional sections during a year period of 2014. Daytime and nighttime are defined as 2200-0600 LST and 0900-1700 LST, respectively. The total number of each section is 664 (S1), 217 (S2), 1198 (S3), 1819 (S4), 1498 (S5), 1774 (S6), 1244 (S7), 864 (S8).

Figure 3은 DR과 PF 좌표변환 방법에 적용된 평균 유선면의 기울기 각도(tilt angle)를 나타내며, 각 방위별 측정 지점의 지면 고도 경사각을 함께 제시하였다. 지면 고도 경사각(α_topo = tan^-1(Δz/Δx))은 측정 지점을 중심으로 50 m 간격으로 계산한 후, 이를 풍상쪽 300 m 구간의 평균으로 계산하였다. Δx와 Δz는 각각 두 지점 사이의 수평 거리와 지면 고도의 차이를 의미한다. Table 2는 PF 방법에서 각 방위별로 계산된 다중 선형 회귀식의 계수를 보여준다. 평균 α 값(DR과 PF 방법)은 방위별로 -3~14^o 범위에서 뚜렷한 차이를 보인다. 또한 β 값(PF 방법)은 측정 지점에서 -2~11^o 범위를 나타낸다. DR 방법과 PF 방법에서 결정된 평균 α 값은 S1 구간(PF: -3.0^o, DR: 2.5^o)에서 상대적으로 큰 차이를 보이나 대체로 전 방위에서 잘 일치한다. 하지만, DR 방법에서 각 방위별 평균 α 값에 비해 표준 편차가 크게 나타나며, PF 방법에서는 다중 선형 회귀식의 결정 계수(R²)가 대체로 낮은 값을 나타낸다. 이러한 경향은 남풍 계열에 비해 북풍 계열에서 상대적으로 뚜렷하게 나타난다(Table 2). 평균 α_topo는 전 방위에서 -0.4~7^o의 변화를 보이며, 측정 지점을 기준으로 동쪽과 남쪽 방위에서 상대적으로 큰 지면 고도 경사각을 가진다. S6와 S7 영역에서는 다른 방위와 달리 평균 유선면의 기울기와 지면 고도 경사각의 크기가 다소 큰 차이(~10^o)를 보이고 있는데, 이는 이들 풍향에서 흐름은 측정 타워 주변의 국지 지형에 의한 영향이 큰 것으로 판단된다. 흥미롭게도 이들 영역은 PF 방법에서 평균 유선면 회귀식의 결정 계수가 높은 영역과 대체로 잘 일치한다. 복잡한 지면 피복과 건물들로 구성된 도시 지역에서 바람 패턴은 수평적으로 균질한 지역과 달리 복잡한 3차원 패턴을 보이며, 플럭스 발자국(flux footprint)은 풍향에 매우 민감한 변화를 보인다(Järvi et al., 2009). 비균질 도시 지역의 특성을 가지는 본 연구의 측정 지점에서 동일한 방위에서도 유선면의 기울기 각도의 변동이 크게 나타나는 특징은 바람 패턴이 3차원적이며 매우 복잡한 패턴이 나타남을 시사한다.

Table 2.

Regression coefficients in Eq. (5) and determination coefficients (R2) for eight directional sections. N indicates the number of data at each section. α and β denote the pitch and roll angles in Eqs. (6)-(7), respectively.

Fig. 3.

Coordinate tilt angles (pitch and roll angles) in planar fit (PF) method and double rotation (DR) method for eight directional sections. Vertical bars denote standard deviation values of pitch angle found in DR method.

4.2 산정 난류 플럭스 비교

Figure 4는 두 회전 좌표변환을 적용하여 계산된 난류 현열/잠열/운동량/이산화탄소 플럭스를 보여준다. 여기서는 방위별 평균 유선면의 특성을 고려하여(Fig. 2) 북풍 계열(S1, S2, S7, S8), 남동풍 계열(S3, S4), 남서풍 계열(S5, S6)의 세 그룹으로 구분하여 비교하였다. 북풍 계열은 남풍 계열에 비해 상대적으로 큰 β 값을 나타내며, 남서풍 계열은 남동풍 계열에 비해 큰 α 값을 나타내는 특징을 가진다. 전체적으로 DR 방법을 적용하여 계산된 난류 플럭스가 PF 방법을 적용했을 때에 비해 평균적으로 높은 값을 나타내며(현열: 13%, 잠열: 21%, 운동량: 37%, CO₂: 26%), 스칼라(현열/잠열/CO₂) 플럭스에 비해 운동량 플럭스의 민감도가 높게 나타난다. 이러한 차이는 남풍 계열에 비해 북풍 계열에서 상대적으로 크게 나타난다(Table 3). 이는 북풍 계열에서 두 방법의 유선면에 따른 회전 각도의 차이가 크게 나타나는 것과 일치한다(Fig. 3). 두 방법간 평균 유선면(α와 β)이 유사한 남풍 계열의 경우에도 산정된 플럭스의 차이가 크게 나타나는 이유는 측정지점 주변 흐름의 복잡성에 기인하여 DR 방법에서 산정한 α의 표준편차 값이 크기 때문이다. 또한, 대상지역에서 두 회전 좌표변환에 의해 산정된 난류 플럭스의 차이는 상대적으로 균질한 자연 지표에서 보고되고 있는 결과에 비해 상대적으로 큰 값을 나타내는 특징이 있다(e.g., Yuan et al., 2007; Li et al., 2013; Shimizu, 2015).

Fig. 4.

Comparison of the PF method against the DR method in turbulent fluxes of (a) sensible heat flux (W m-2), (b) latent heat flux (W m-2), (c) momentum flux (kg m-1 s-2), and (d) CO2 flux (mg m-2 s-1). Colored symbols denote ‘Northerly’ (S1, S2, S7, S8) (black circle), ‘Southeasterly’ (S3, S4) (red triangle), and ‘Southwesterly’ (S5, S6) (blue square). The regression relations for the whole data are given with determination coefficients.

Table 3.

Regression coefficients of a linear regression function (y = ax + b) between turbulent fluxes by the PF and DR methods. R2 values are determination coefficients. Directional sections are classified into ‘Northerly’ (S1, S2, S7, S8), ‘Southeasterly’ (S3, S4), and ‘Southwesterly’ (S5, S6).

Figure 5는 두 회전 좌표변환 방법에 의해 산정된 난류 플럭스의 상대편차를 대기 안정도 조건에 따라 구분하여 비교하고 있다. 대기 안정도는 0.1 < z/L ≤ 10을 안정 조건으로, -10 ≤ z/L < -0.1 을 불안정 조건으로 구분하였다. DR 방법과 PF 방법을 적용한 후 계산된 L 값이 동일한 안정도 조건을 나타내는 자료만을 분석에 사용하였으며, 전체 분석 자료 중 불안정과 안정 조건의 발생 빈도는 각각 63%와 9%를 차지하였다. PF 방법에 의해 산정된 현열 플럭스는 DR 방법에 비해 불안정 조건에서 최소 3% (S6)부터 최대 19% (S7과 S8)까지 과소 산정하였다. 안정 조건에서의 두 방법간의 상대적 편차는 불안정 조건에 비해 상대적으로 크게(최대 44% (S7)) 나타나며, 남서풍 계열인 S6에서는 3%로 과대 산정을 보였다(Fig. 5a). 잠열 플럭스는 현열 플럭스에서와 유사하게 불안정과 안정 조건에서 각각 2~28%와 1~42%의 범위에서 과소 산정하였다(Fig. 5b). 운동량 플럭스는 불안정 조건과 안정 조건에서 각각 15~44%와 2~41%로 과소 산정하였다(Fig. 5c). CO₂ 플럭스는 불안정 조건에서 3~26%, 안정 조건에서 -3~41% 과대/과소 산정하였다(Fig. 5d). 스칼라(현열/잠열/CO₂) 플럭스는 불안정 조건에 비해 안정 조건에서 대체로 큰 상대편차를 나타내었으며, 운동량 플럭스의 상대편차는 대기 안정도에 상대적으로 덜 민감하였다. 측정지점에서 두 회전 좌표변환 방법을 적용하여 산정한 난류 플럭스는 안정과 불안정 조건 모두에서 북서풍 계열(S7과 S8)에서 상대편차가 크게 나타나고 남서풍 계열인 S6에서 가장 낮게 나타나, 풍향에 따른 상대편차의 의존성이 높게 나타나는 특징을 보였다. 두 방법을 통해 산정한 난류 플럭스의 상대편차는 안정과 불안정 조건 모두에서 풍속이 증가함에 따라 대체로 감소하는 경향을 나타내었다(Fig. 6).

Fig. 5.

Fractional differences (= (FPF - FDR)/FDR) of the turbulent fluxes determined by the PF and DR methods for eight directional sections. The fractional difference values are first calculated at every 30-min intervals, and then the values are averaged for each section. Median flux values calculated from the DR method are presented for atmospheric unstable and stable conditions, respectively.

Fig. 6.

Fractional differences (= (FPF - FDR)/FDR) of the turbulent fluxes determined by the PF and DR methods as a function of wind speed class for (a) unstable and (b) stable condition.

5. 요약 및 결론

도시 지표와 대기 사이의 운동량, 에너지, 질량의 난류 교환 과정의 정량적 이해는 도시 미기후 분석뿐만 아니라 도시 대기 모델링의 예측 검증을 위해서도 필수적인 요소이다. 에디 공분산 방법은 대기 지표층 난류 교환을 정량적으로 측정하기 위해 널리 활용되어 오고 있는 방법이며, 측정 지점에서의 경사진 지형이나 국지 순환에 의한 오류를 보정하기 위해 품질관리 과정에서 대부분 회전 좌표변환 과정을 거친다. 여러 연구자들에 의해 좌표변환 방법들 간의 차이에 대한 연구가 있었으나(e.g., Wilczak et al., 2001; Sun, 2007; Yuan et al., 2007), 대부분 산림지역에 대해 주로 수행되었다. 본 연구에서는 비균질 도시 지표면에 설치된 기상 타워에서 측정된 2014년 1년간의 에디 공분산 자료를 이용하여 회전 좌표변환이 난류 현열/잠열/운동량/이산화탄소 플럭스 산정에 미치는 영향을 분석하였다. 이를 위해 널리 활용되는 DR 방법과 PF 회전 좌표변환 방법을 적용하였으며, 비균질 지표면의 특성에 따른 차이를 분석하기 위해 측정 지점을 중심으로 8 방위로 구분하여 분석을 수행하였다.

측정 지점에 대해 PF 방법에 의한 평균 유선면은 방위별로(α = -3.0~13.8^o와 β = -2.0~10.9^o)를 나타내었으며, DR 방법에 의한 평균 유선면의 기울기는 α = 2.2~14.5^o로 PF 방법과 전 방위에 대해 대체로 잘 일치하였다. PF 방법에 의해 산정된 난류 플럭스는 DR 방법에 의해 산정된 값에 비해 평균적으로 현열 플럭스는 13%, 잠열 플럭스는 21%, 운동량 플럭스는 37%, 그리고 CO₂ 플럭스는 26% 과소 산정되었다. 두 방법에 의해 산정된 난류 플럭스의 편차는 풍향에 따라 큰 차이를 보였으며, 평균 유선면의 차이가 상대적으로 큰 북풍 계열(S1, S2, S7, S8)에서 남풍 계열에 비해 높은 편차를 나타내었다. 평균 유선면의 차이가 크지 않은 풍향에 대해서도 산정 플럭스의 큰 편차가 나타나는데, 이는 비균질 도시 지표에서 유선면의 시/공간적 변동성이 큰 것에 기인하는 것으로 에디 공분산 난류 플럭스 산정 과정의 주요한 오차 요소로 평가할 수 있다. 두 방법간 산정 플럭스의 상대편차는 안정과 불안정 조건 모두에서 평균 풍속이 증가할수록 대체로 감소하는 경향을 보였다.

본 연구의 결과는 에디 공분산 방법을 이용하여 비균질 도시 지표에서의 난류 플럭스를 산정하는 경우 회전 좌표변환 과정의 산정 오차가 균질한 지표면에서 평가된 값들에 비해 현저하게 크게 나타날 수 있음을 보여준다. 이러한 EC 방법에 의한 난류 플럭스의 산정 불확실성은 에디 공분산 난류 플럭스를 이용하여 미기후 특성을 분석하거나 혹은 수치 모형이나 원격 측정 플랫폼 등 다양한 형태로 얻어진 지표면 난류 플럭스(e.g., Lee, 2015a; Lee et al., 2015)와 상호 비교할 경우 측정 오차를 정량적으로 포함하여 분석을 수행하여야 함을 시사한다.

Acknowledgments

논문의 개선을 위해 좋은 의견을 주신 두 분의 심사위원께 감사 드립니다. 이 연구는 기상청 기상산업지원 및 활용기술개발사업(과제번호: KMIPA2015-5043)의 지원으로 수행되었습니다.

References

Aubinet, M., T. Vesala, and D. Papale Eds, (2012), Eddy covariance: a practical guide to measurement and data analysis, Springer, p438.
Baldocchi, D. D., B. B. Hincks, and T. P. Meyers, (1998), Measuring biosphere-atmosphere exchange of biologically related gasses with micrometeorological methods, Ecology, 69, p1331-1340. [https://doi.org/10.2307/1941631]
Baldocchi, D. D., and Coauthors , (2001), FLUXNET: a new tool to study the temporal and spatial variability of ecosystem- scale carbon dioxide, water vapor, and energy flux densities, Bull. Amer. Meteor. Soc, 82, p2415-2434. [https://doi.org/10.1175/1520-0477(2001)082<2415:fantts>2.3.co;2]
Foken, T., and B. Wichura, (1996), Tools for quality assessment of surface-based flux measurements, Agric. For. Meteor, 78, p83-105. [https://doi.org/10.1016/0168-1923(95)02248-1]
Foken, T., M. Göckede, M. Mauder, L. Mahrt, B. Amiro, and W. Munger, (2004), Post-field data quality control, Handbook of Micrometeorology, W. Lee, W. Massman, and B. E. Law Eds, Kluwer Academic Publishers, p181-208.
Grimmond, C. B. S., J. A. Salmond, T. R. Oke, B. Offerle, and A. Lemonsu, (2004), Flux and turbulence measurements at a densely built-up site in Marseille: Heat, mass (water and carbon dioxide), and momentum, J. Geophys. Res, 109, pD24. [https://doi.org/10.1029/2004jd004936]
Hong, J., and J. Kim, (2002), On processing raw data from micrometeorological field experiments, Korean J. Agri. Forest Meteor, 4, p119-126.
Järvi, L., Ü Rannik, I. Mammarella, A. Sogachev, P. P. Aalto, P. Keronen, E. Siivola, M. Kulmala, and T. Vesala, (2009), Annual particle flux observations over a heterogeneous urban area, Atmos. Chem. Phys, 9, p7847-7856. [https://doi.org/10.5194/acp-9-7847-2009]
Kim, S., Y.-H. Lee, K. R. Kim, and Y.-S. Park, (2014), Analysis of surface energy balance closure over heterogeneous surfaces, Asia-Pac. J. Atmos. Sci, 50, p553-565. [https://doi.org/10.1007/s13143-014-0045-2]
Kotthaus, S., and C. S. B. Grimmond, (2014), Energy exchange in a dense urban environment-Part I: Temporal variability of long-term observations in central London, Urban Climate, 10, p261-280. [https://doi.org/10.1016/j.uclim.2013.10.002]
Lee, S.-H., (2015a), LAS-derived determination of surface-layer sensible heat flux over a heterogeneous urban area, Atmosphere, 25, p193-203, (in Korean with English abstract). [https://doi.org/10.14191/Atmos.2015.25.2.193]
Lee, S.-H., J.-H. Lee, and B.-Y. Kim, (2015), Estimation of turbulent sensible heat and momentum fluxes over a heterogeneous urban area using a large aperture scintillometer, Adv. Atmos. Sci, 32, p1092-1105. [https://doi.org/10.1007/s00376-015-4236-2]
Lee, X., W. J. Massman, and B. Law Eds, (2004), Handbook of micrometeorology: a guide for surface flux measurement and analysis, Kluwer, p250.
Lee, Y.-H., (2015b), Characteristics of wind direction shear and momentum fluxes within roughness sublayer over sloping terrain, Atmosphere, 25, p591-600, (in Korean with English abstract). [https://doi.org/10.14191/Atmos.2015.25.4.591]
Lee, Y.-H., B. Lee, K. Kahng, S.-J. Kim, and S.-O. Hong, (2013), Quality control and characteristic of eddy covariance data in the region of Nakdong river, Atmosphere, 23, p307-320, (in Korean with English abstract). [https://doi.org/10.14191/Atmos.2013.23.3.307]
Li, M., W. Babel, K. Tanaka, and T. Foken, (2013), Note on the application of planar-fit rotation for non-omnidirectional sonic anemometers, Atmos. Meas. Tech, 6, p221-229. [https://doi.org/10.5194/amt-6-221-2013]
McMillen, R. T., (1988), An eddy correlation technique with extended applicability to non-simple terrain, Bound.-Layer Meteor, 43, p231-245. [https://doi.org/10.1007/BF00128405]
Mestayer, P., and Coauthors , (2005), The urban boundary layer field experiment over Marseille UBL/CLU-ESCOMPTE: experimental set-up and first results, Bound.-Layer Meteor, 114, p315-365. [https://doi.org/10.1007/s10546-004-9241-4]
Park, S.-J., T.-J. Choi, and S.-J. Kim, (2013), Heat flux variations over sea ice observed at the coastal area of the Sejong station, Antarctica, Asia-Pac. J. Atmos. Sci, 49, p443-450. [https://doi.org/10.1007/s13143-013-0040-z]
Ross, A. N., and E. R. Grant, (2015), A new continuous planar fit method for calculating fluxes in complex forested terrain, Atmos. Sci. Let, 16, p445-452. [https://doi.org/10.1002/asl.580]
Ruppert, J., M. Mauder, C. Thomas, and J. Lüers, (2006), Innovative gapfilling strategy for annual sums of CO₂ net ecosystem exchange, Agric. Forest Meteorol, 138, p5-18. [https://doi.org/10.1016/j.agrformet.2006.03.003]
Schmid, H. P., C. S. B. Grimmond, F. Cropley, B. Offerle, and H.-H. Su, (2000), Measurements of CO₂ and energy fluxes over a mixed hardwood forest in the mid-western United States, Agric. Forest Meteorol, 103, p357-374. [https://doi.org/10.1016/S0168-1923(00)00140-4]
Shimizu, T., (2015), Effect of coordinate rotation systems on calculated fluxes over a forest in complex terrain: a comprehensive comparison, Bound.-Layer Meteor, 156, p277-301. [https://doi.org/10.1007/s10546-015-0027-7]
Sun, J., (2007), Tilt corrections over complex terrain and their implication for CO₂ transport, Bound.-Layer Meteor, 124, p143-159. [https://doi.org/10.1007/s10546-007-9186-5]
Tanner, C. B., and G. W. Thurtell, (1969), Anemoclinometer measurements of Reynolds stress and heat transport in the atmospheric surface layer, Research and Development Technical Report ECOM 66-G22-F to the US Army Electronics Command, p200.
Tsukamoto, O., and F. Kondo, (2008), Experimental validation of the Webb correction for CO₂ flux with an open-path gas analyzer, 18th Symposium on Boundary Layers and Turbulence, Stockholm University, p5.
Velasco, E., and Coauthors , (2009), Eddy covariance flux measurements of pollutant gases in urban Mexico City, Atmos. Chem. Phys, 9, p7325-7342. [https://doi.org/10.5194/acp-9-7325-2009]
Vickers, D., and L. Mahrt, (1997), Quality control and flux sampling problems for tower and aircraft data, J. Atmos. Oceanic Technol, 14, p512-526. [https://doi.org/10.1175/1520-0426(1997)014<0512:QCAFSP>2.0.CO;2]
Webb, E. K., G. I. Pearman, and R. Leuning, (1980), Correction of flux measurement for density effects due to heat and water vapour transfer, Quart. J. Roy. Meteor. Soc, 106, p85-100. [https://doi.org/10.1002/qj.49710644707]
Wilczak, J. M., S. P. Oncley, and S. A. Sage, (2001), Sonic anemometer tilt correction algorithms, Bound.-Layer Meteor, 99, p127-150. [https://doi.org/10.1023/A:1018966204465]
Yuan, R., M. Kang, S. Park, J. Hong, D. Lee, and J. Kim, (2007), The effect of coordinate rotation on the eddy covariance flux estimation in a hilly KoFlux forest catchment, Korean J. Agri. Forest Meteor, 9, p100-108. [https://doi.org/10.5532/KJAFM.2007.9.2.100]
Yuan, R., M. Kang, S.-B. Park, J. Hong, D. Lee, and J. Kim, (2011), Expansion of the planar-fit method to estimate flux over complex terrain, Meteorol. Atmos. Phys, 110, p123-133. [https://doi.org/10.1007/s00703-010-0113-9]
Zeweldi, D. A., M. Gebremichael, J. Wang, T. Sammis, J. Kleissl, and D. Miller, (2010), Intercomparison of sensible heat flux from large aperture scintillometer and eddy covariance methods: field experiment over a heterogeneous semi-arid region, Bound.-Layer Meteor, 135, p151-159. [https://doi.org/10.1007/s10546-009-9460-9]

Section	Residential	Commercial/Industrial/Public	Transportation	Agricultural	Forest	Bare land
S1	9.8	11.0	41.2	22.6	14.8	0.6
S2	11.2	9.6	58.9	3.9	15.3	1.1
S3	1.5	10.3	25.0	5.2	57.2	0.8
S4	1.0	9.4	31.9	0.8	56.9	0.0
S5	0.0	23.6	33.5	0.1	35.5	7.3
S6	2.3	25.4	32.0	0.3	35.0	5.0
S7	4.8	21.2	34.1	21.0	18.8	0.1
S8	11.0	15.4	46.3	8.9	14.7	3.7

Section	N	b_o	b₁	b₂	R²	α	β
S1	664	0.07	–0.05	-0.06	0.45	-3.0	-3.5
S2	217	0.05	-0.04	–0.04	0.06	2.3	–2.0
S3	1198	0.05	-0.04	-0.00	0.33	2.0	-0.2
S4	1819	0.04	-0.04	-0.00	0.19	2.1	-0.2
S5	1498	0.04	-0.19	-0.01	0.89	10.7	-0.8
S6	1774	0.01	-0.25	-0.00	0.94	13.8	–0.3
S7	1244	0.12	-0.15	-0.19	0.54	8.5	10.9
S8	864	0.11	-0.06	-0.08	0.43	3.7	-4.8

	Sensible heat flux (W m^-2)			Latent heat flux (W m^-2)			Momentum flux (kg m^-1 s^-2)			CO₂ flux (mg m^-2 s^-1)
	a	b	R²	a	b	R²	a	b	R²	a	b	R²
Total	0.87	1.06	0.96	0.79	1.23	0.83	0.63	0.01	0.88	0.74	0.05	0.80
Northerly	0.76	2.76	0.95	0.68	1.26	0.81	0.57	0.01	0.82	0.63	0.05	0.76
Southeasterly	0.91	0.25	0.99	0.84	0.74	0.87	0.73	0.00	0.92	0.85	0.01	0.90
Southwesterly	0.91	0.82	0.98	0.83	2.16	0.83	0.63	0.01	0.89	0.79	0.08	0.81