Estimation of Turbulent Sensible Heat and Momentum Fluxes over a Heterogeneous Urban Area Using a Large Aperture Scintillometer

Turbulent transfer of heat, moisture, and momentum continuously occurs at the interface between the atmosphere and underlying Earth surface. The atmospheric boundary layer (ABL) is formed as a consequence of the interactions on time scales of a few hours or less (Garratt, 1992; Foken, 2008). Therefore, accurate determination of the turbulent fluxes in the surface layer is of primary importance to understanding the spatial and temporal variability of atmospheric flow and thermal structure through the ABL, which sequentially influences the formation of the characteristic local climate and meteorology, as well as air quality.

Generally, eddy-covariance (EC) systems and scintillometers have been widely used to determine surface-layer turbulent fluxes. The EC system measures surface-layer turbulent fluxes on the order of hundreds of meters or less, which has been successfully used to provide the local-scale turbulent fluxes for various natural and urban surfaces (e.g., Offerle et al., 2006). The EC system directly measures atmospheric turbulent fluctuations, and thus the EC technique has great merit in producing accurate turbulent fluxes for relatively homogeneous flat surfaces. However, it is difficult to determine the sensor's location for measuring turbulent fluxes over highly heterogeneous surfaces, due to the uncertainty of surface source areas of turbulent fluxes observed at a particular location. Unlike the EC technique, the scintillometry technique uses Monin-Obukhov similarity theory to estimate area-averaged surface-layer turbulent fluxes. Thus, the turbulent fluxes can be determined at various spatial scales ranging from 100 m to a few kilometers, depending on the scintillometer type and optical path distance. For example, a small aperture scintillometer (SAS) uses a short optical path length (a few hundred meters), while a large aperture scintillometer (LAS) has a longer path length of up to 5 km.

Many field experiments based on the scintillometry technique have been conducted to determine surface-layer turbulent fluxes over natural surfaces (e.g., De Bruin et al., 1995, 2002; Chehbouni et al., 2000; Beyrich et al., 2002; Zeweldi et al., 2010; Evans et al., 2012) and urban areas (e.g., Kanda et al., 2002; Lagouarde et al., 2006; Pauscher, 2010; Ward et al., 2014), and showed the usefulness of LAS-derived turbulent fluxes through direct comparison with EC measurements (e.g., McAneney et al., 1995; Hartogensis et al., 2003; Lagouarde et al., 2006; Zeweldi et al., 2010; Liu et al., 2011, 2013; Xu et al., 2013; Ward et al., 2014). For example, (Lagouarde et al., 2006) showed the LAS-derived turbulent fluxes with a path length of 1.8 km compared well to the EC fluxes over the city of Marseille, France. (Ward et al., 2014) compared two LAS systems with different path lengths of 2.8 and 5.5 km and an EC instrument for surface-layer turbulent flux estimation over a residential area in Swindon, UK. These studies demonstrated the applicability of the LAS system over homogeneous residential areas in determining surface-layer turbulent fluxes in the urban environment.

Most urban areas at spatial scales of a few km or less are, morphologically, highly heterogeneous due to the diversity of buildings, natural/artificial land-use mix, and topography. Even though previous studies have successfully determined the surface-layer turbulent fluxes over homogeneous urban areas using the EC system and/or LAS system (e.g., Lagouarde et al., 2006; Ward et al., 2014), there are few studies concerning turbulent flux estimation over heterogeneous urban surfaces. As discussed, the EC technique is difficult to apply for highly heterogeneous surfaces, whereas the scintillation method using the LAS system can provide surface-layer turbulent fluxes integrated over the LAS-path distance, suffering less due to the surface heterogeneity than the EC system. In the present study, the aim was to investigate the ability of the LAS technique to estimate surface-layer turbulent fluxes over an urban area with a significant mix of different land-use types. First, surface-layer turbulent sensible heat and momentum fluxes in a highly heterogeneous urban area were estimated using the LAS system. Then, for verification of the LAS-derived turbulent fluxes, the results were compared to those of a land-surface process-resolving numerical model and radiative flux measurements. Three-dimensional LAS footprint modeling was also performed to determine the source areas of the estimated turbulent fluxes over the heterogeneous urban area.

The remainder of the paper is structured as follows: Section 2 briefly describes the scintillation methodology for estimating surface-layer turbulent fluxes and the 3D LAS footprint modeling. The experimental site, LAS instrumentation, and land-surface process-resolving modeling are described in section 3. Section 4 presents the LAS-derived turbulent fluxes and footprints over the heterogeneous urban area, and discusses the uncertainty of the LAS-derived turbulent fluxes through a series of sensitivity tests. A summary and conclusions follow in section 5.

When electromagnetic radiation propagates through the atmosphere, its intensity is affected by the atmospheric turbulence, leading to fluctuations in the refractive index of the air (e.g., Frehlich and Ochs, 1990). The path-integrated structure parameter of the refractive index \(\langle C_n^2\rangle\) can be given (Wang et al., 1978) by \begin{equation} \langle C_n^2\rangle=\int_0^1C_n^2(r)W(r)dr , \end{equation} where r(=x/L _LAS) is the normalized distance along the LAS path (L _LAS) between transmitter and receiver, and C_n²(r) is the structure parameter at distance r. Here, x and L _LAS indicate distance from the transmitter and optical path length, respectively. W(r) is the LAS weighting function, which has a mathematical bell-shape curve with a maximum value in the middle of the optical path (e.g., Chehbouni et al., 2000; Meijninger and De Bruin, 2000; Timmermans et al., 2009; Geli et al., 2012): \begin{eqnarray} W(r)&=&16\pi K^2L_{\rm LAS}\int_0^\infty k{\O}_n(k)\times\nonumber\\ &&\sin^2\left[\dfrac{k^2L_{\rm LAS}r(1-r)}{2K}\right]\left[\dfrac{4J_1(x_1)J_1(x_2)}{x_1x_2}\right]^2dk ,\quad \end{eqnarray} where K(=2π/Λ) is the optical wavenumber, k is the turbulent spatial wavenumber, Ø_n(k)(=0.033k^-11/3) is the 3D Kolmogorov spectrum of the refractive index, J₁ is a Bessel function of the first kind for order one, with x₁=KDr/2 and x₂=[KD(1-r)]/2, where D is the aperture diameter of the LAS transmitter/receiver.

Scintillometers actually measure the intensity fluctuation (I) of the refractive index of the air ("scintillation") in terms of the structure parameter of the refractive index \(\langle C_n^2\rangle\), which can be written following (Tatarskii, 1961) and (Wang et al., 1978) as \begin{equation} \langle C_n^2\rangle=1.12\sigma_{lnI}^2D^{7/3}L_{\rm LAS}^{-3}, \end{equation} where σ_lnI² is the statistical variance of the logarithm of intensity fluctuation obtained from measured radiation intensity. Hereafter, \(\langle C_n^2\rangle\) is denoted as C_n² for simplicity.

Because the intensity fluctuation changes with air temperature and humidity, C_n² can be expressed as a combination of the structure parameters of temperature (C_T²), humidity (C_Q²), and their covariance fluctuations (C_TQ²). The influence by atmospheric pressure fluctuation can be negligible (Hill et al., 1980). The scintillation is predominantly influenced by temperature fluctuation for visible and infrared radiation (Wesely, 1976; Moene, 2003), leading to a simplified relation between C_n² (m^-2/3) and C_T² (K² m^-2/3) as follows: \begin{equation} C_T^2=C_n^2\left(\dfrac{T^2}{-0.78\times 10^{-6}P}\right)^2\left(1+\dfrac{0.031}{\beta}\right)^{-2} , \end{equation} where T is the air temperature (K), P is the atmospheric pressure (Pa), and β is the Bowen ratio, which represents a humidity correction in scintillation.

Provided that C_T² is estimated from scintillometer measurement, turbulent sensible heat and momentum fluxes can be derived based on surface-layer similarity relations (Wyngaard et al., 1971). First, the similarity relation that links C_T² to temperature scale (T_*) can be given as \begin{equation} \dfrac{C_T^2(z_{\rm eff}-d)^{2/3}}{T_*^2}=f_T\left(\dfrac{z_{\rm eff}-d}{L}\right)\quad (L<0) , \end{equation} where z _eff is the LAS effective height (m), d is the zero-plane displacement height (m), and f_T is a dimensionless universal function for the temperature structure function parameter. The L is the Obukhov length, defined by \begin{equation} L=\dfrac{u_*^2T}{kgT_*} , \end{equation} where k (=0.4) is the von Karman constant and g (=9.8 m s^-2) is the gravitational acceleration. The non-dimensional universal function for the temperature structure function parameter can be formulated for unstable atmospheric conditions following (Andreas, 1988): \begin{equation} f_T=4.9\left(1-6.1\dfrac{z_{\rm eff}}{L}\right)^{-2/3}. \end{equation} The friction velocity (u_*) is also estimated by the following similarity relation (e.g., Panofsky and Dutton, 1984; Foken, 2008): \begin{equation} \dfrac{u}{u_*}=\dfrac{1}{k}\left[\ln\left(\dfrac{z_{\rm eff}-d}{z_0}\right)-\psi_{\rm m}\left(\dfrac{z_{\rm eff}-d}{L}\right)+\psi_{\rm m}\left(\dfrac{z_0}{L}\right)\right], \end{equation} where U is the mean wind speed at z _eff, and ψ _m is the integrated universal function for momentum flux (e.g., Paulson, 1970; Businger et al., 1971). The integrated universal function used in this study for the unstable atmospheric surface layer is (H\"ogstr\"om, 1988) \begin{equation} \psi_{\rm m}=2\ln\left[\dfrac{1+x}{2}\right]+\ln\left[\dfrac{1+x^2}{2}\right]-2{\rm arctan}(x)+\dfrac{\pi}{2}\quad (L<0) , \end{equation} where \(x=\left[1-19.3\dfracz_\rm effL\right]^1/4\).

Given the similarity relations, the temperature scale (T_*), friction velocity (u_*), and the Obukhov length (L) can be solved iteratively. The turbulent sensible heat flux (Q _H; W m-2) and momentum flux (τ; kg m^-1 s^-2) can be calculated as follows: \beginsubequations \begin{align} Q_{\rm H}&=-\rho c_pu_*T_* ,\\ \tau&=\rho u_*^2 , \end{align} \endsubequations where ρ is the air density (kg m^-3) and c_p is the specific heat of the air at a constant pressure (J kg^-1 K^-1).

Additional parameters required for estimation of the turbulent fluxes are the LAS effective height (z _eff), surface aerodynamic roughness length (z₀), and zero-plane displacement height (d). First, the LAS effective height can be estimated iteratively from the following relation (Hartogensis et al., 2003): \begin{equation} z_{\rm eff}^{-2/3}f_T\left(\dfrac{z_{\rm eff}}{L}\right)=\int_0^1z(r)^{-2/3}f_T\left(\dfrac{z(r)}{L}\right)W(r)dr , \end{equation} where z(r) is the optical beam height at distance r. This relation indicates that z _eff is dependent on atmospheric stability, and thus it was estimated at every time increment (30 min) of the available data in this study.

Many flux footprint models have been suggested based on analytical and numerical approaches and applied in the interpretation of EC flux measurements (e.g., Schuepp et al., 1990; Horst and Weil, 1992, 1994; Schmid, 1994; Horst, 1999; Hsieh et al., 2000; Kormann and Meixner, 2001; Kljun et al., 2004). The 3D LAS footprint can be calculated as an extension of the flux footprint model; that is, a combination of the flux footprint calculated at each segment location along the optical path and the LAS spatial weighting function (e.g., Hoedjes et al., 2002; Meijninger et al., 2002; Hartogensis et al., 2003). The turbulent flux footprint relation used for the LAS footprint modeling in this study is briefly introduced below.

Figure 1.  (a) Aerial map of the experimental site (source: Google Maps). The LAS (large aperture scintillometer) beam path (yellow line) and the location of the meteorological site (red star) are shown. The LAS optical path length between the receiver (A) and transmitter (B) is 2.1 km. (b) Topography along the LAS beam path (shaded), the LAS beam height (dashed line), and the LAS weighting function (solid line with dots).

The turbulent flux F(x,y,z_m) measured at height z_m can be related to the footprint function f(x,y,z_m) and the spatial distribution of the surface fluxes F₀(x',y',z'=0) as follows (Horst and Weil, 1992): \begin{equation} F(x,y,z_m)\!=\!\int_{-\infty}^\infty\int_{-\infty}^x F_0(x',y',z'\!=\!0)f(x-x',y-y',z_m)dx'dy'\,. \end{equation} Here, the flux footprint function (units: m^-2) is determined by the product of the crosswind distribution function D_y(x,y) (units: m^-1) and the crosswind-integrated footprint function f_y(x,z_m) (units: m^-1) (van Ulden, 1978; Horst and Weil, 1992) as \begin{equation} f(x,y,z_m)=D_y(x,y)f_y(x,z_m) , \end{equation} where x is the downwind distance from the measurement location, and y is the crosswind distance from the centerline. The dispersion to the crosswind direction is assumed to be symmetric using a Gaussian distribution function (Pasquill, 1974): \begin{equation} D_y(x,y)=\dfrac{1}{\sqrt{2\pi}\sigma_y}\exp\left(-\dfrac{y^2}{2\sigma_y^2}\right) , \end{equation} where σ_y is the standard deviation of the plume in the y dimension, which relies on atmospheric turbulence intensity, upwind distance (x), and mean wind speed [u(x)] as σ_vx/u(x) for a short-range dispersion (Schmid, 1994). Here, the standard deviation of lateral wind fluctuations, σ_v(=1.9u_*), was calculated following (Panofsky and Dutton, 1984). On the other hand, the footprint model by (Hsieh et al., 2000) was used to calculate the crosswind-integrated footprint: \begin{equation} f_y(x,z_m)=\dfrac{1}{k^2x^2}Dz_u^P|L|^{1-P}\exp\left(-\dfrac{1}{k^2x}Dz_u^P|L|^{1-P}\right), \end{equation} where D and P are specific constants that depend on atmospheric stability, and z_u is a length scale defined as \begin{equation} z_u=z_m\left[\ln\left(\dfrac{z_m}{z_0}\right)-1+\dfrac{z_0}{z_m}\right] . \end{equation} The constant parameters are D=0.28 and P=0.59 for unstable conditions (Hsieh et al., 2000).

The field experiment was conducted in a small urban area in Gongju, South Korea, where there is a population of 117 000 people. The city is located in a basin surrounded by extensive agricultural and hilly regions. It is also located inland, about 50 km west of the Yellow Sea. Figure 1a shows the experimental site within the city. The urbanized area is largely composed of single (1-3 story) residential houses, and a river (named "Geum River") crosses over the city. Relatively new buildings are located alongside the northern stretches of the river, showing the directions of the city's sprawl. The vegetation fraction within the built-up region ranges from 5% to 10%, and hilly mountains surround the city in the south and east directions. The measurement site was highly heterogeneous in terms of its land cover and topography, which is a typical morphological configuration that can be found in many cities of South Korea.

The LAS instruments were deployed at two buildings of the Kongju National University, located within the city, to investigate the ability in estimating turbulent sensible heat and momentum fluxes over the heterogeneous urban area (Fig. 1a). The receiver (marked A in Fig. 1a) and transmitter (marked B in Fig. 1a) were installed on the rooftops of buildings that are 49 m and 48 m above mean sea level, and the optical path extends 2.1 km from southwest to northeast, crossing the Geum River (A-B in Fig. 1a). This configuration of the LAS instruments measured area-averaged turbulent sensible heat and momentum fluxes that were representative of the heterogeneous urban area (residential area, bare ground, vegetated area, water body, etc.).

The LAS system used was the LAS MkII manufactured by Kipp & Zonen, with an aperture diameter of D=0.149 m. The light source of the LAS MkII transmitter operates at a near-infrared wavelength of 850 nm, at which the scintillations measured by C_n² are primarily influenced by turbulent temperature fluctuations. Based on (Kleissl et al., 2010), the threshold C_n² value for signal saturation was estimated as 4.06× 10^-14 m^-2/3 for our LAS experimental configuration. The LAS sampled C_n² at a 1 Hz frequency, which were averaged to 15 s intervals for turbulent flux estimation. The air temperature, wind speed, and atmospheric pressure were measured at the same temporal frequency from separate meteorological sensors, but which were connected to the LAS system at the same height of the receiver. In addition, a 7 m meteorological tower on the rooftop of an eight-story building was also operated at the location of 320 m north of the receiver along the optical path. The meteorological data measured from the meteorological tower were used for the calculation of the LAS-derived turbulent fluxes.

The field measurements were conducted for 29 days during 5-19 November and 5-18 December 2013, in which adverse meteorological conditions (e.g., fog) affected the measurements for many days. In this study, two sequential clear days on 11 and 12 November 2013 were selected to investigate the capability of the LAS in estimating the turbulent fluxes over the heterogeneous surfaces. Furthermore, since the LAS cannot discriminate atmospheric stability from the measured C_n² values, the daytime period (0800-1700 LST) data were analyzed for estimation of the LAS-derived turbulent fluxes in unstable atmospheric conditions.

Figure 1b shows the topographical variation along the LAS beam path and the optical beam height relative to the ground. From Eq. (11), the LAS effective heights (z _eff) were estimated as 25.1-27.2 m during the daytime periods of 11 and 12 November, which satisfied a "saturation condition" required for reasonable estimation of turbulent fluxes using the LAS system (Kipp & Zonen, 2012). The Bowen ratio (β) was initially assigned as 1. The effective aerodynamic roughness length and displacement height were assigned as z₀=0.3 m and d=2.1 m, respectively. The aerodynamic parameters were determined following the general rule of z₀≈ 0.1h and d≈0.7h (Garratt, 1992; Grimmond and Oke, 1999), in which h(∼ 3 m) indicates the mean height of obstacles along the beam path. Generally, the determination of surface aerodynamic parameters can be obtained by applying a morphological method (e.g., Raupach, 1994; Macdonald et al., 1998) and/or surface-layer similarity relation using EC measurements taken in the atmospheric surface layer for homogeneous surfaces (e.g., Grimmond and Oke, 1999; Roth, 2000). However, accurate estimation of the surface aerodynamic parameters is a very difficult task for heterogeneous surfaces due to intrinsic limitations of the methodologies. No standard method that determines the surface aerodynamic parameters for heterogeneous surfaces exists, which inevitably leads to uncertainty in the LAS-derived surface-layer turbulent fluxes. In this study, a series of sensitivity tests were performed to quantify the influence of errors in the surface input parameters on the LAS-derived turbulent fluxes.

The ability of the LAS equipment used in this study in estimating surface-layer turbulent fluxes has been validated against the EC system for homogeneous land surfaces (e.g., Wilson et al., 2013). Due to the difference of the footprints of the EC and LAS-derived turbulent fluxes, direct comparison of the two turbulent fluxes is difficult for highly heterogeneous surfaces. In order to verify the LAS-derived turbulent fluxes over the heterogeneous field experimental area, land-surface process-resolving modeling was used in the present study. Specifically, we conducted very high resolution urban numerical modeling using the Weather Research and Forecasting (WRF) model (Skamarock et al., 2008). Given a proper experimental configuration, the WRF model has the capability to explicitly represent significant surface land-use heterogeneity over the field experiment region and to resolve associated land-surface physical (radiative/dynamic/hydrological) processes.

The WRF model was configured with five nested domains, covering East Asia with a 16.2 km grid spacing in domain 1 and nesting down to the Gongju region with a 200 m grid spacing in domain 5 (Fig. 2). To realistically represent the significant heterogeneity of land-use and topography over the experimental region, very high resolution static datasets of the Moderate Resolution Imaging Spectroradiometer (MODIS)-based land-use (Kang et al., 2010; spatial resolution of 7.5") and Shuttle Radar Topography Mission (SRTM) topography (http://www2.jpl.nasa.gov/srtm/index.html; spatial resolution of 3") were implemented as input data of the WRF model. Physical processes applied in the simulation included: Goddard shortwave radiation (Chou and Suarez, 1994), Rapid Radiative Transfer Model (RRTM) longwave radiation (Mlawer et al., 1997), Lin cloud microphysics (Lin et al., 1983), and a turbulent kinetic energy (TKE)-based turbulence scheme (Janji\'c, 2002). Land-surface processes were parameterized following the Noah land surface model (Chen and Dudhia, 2001) with a single-layer urban canopy model (Kusaka et al., 2001) for urban patches.

The WRF simulation was conducted starting from 0000 UTC 9 November 2013 to 0000 UTC 13 November 2013. In order to realistically reproduce surface-layer meteorology and relevant turbulent fluxes, observed wind fields from the meteorological site were nudged during the simulation period. The WRF model calculates surface-layer turbulent sensible heat and momentum fluxes at every grid mesh. The WRF-simulated turbulent fluxes with a 30 min time interval were extracted along the LAS beam path, which were multiplied by the LAS weighting function for comparison with the LAS-derived fluxes. In this study, a WRF model sensitivity simulation was also conducted, by replacing the Geum River with a natural surface type ("shrubland") to further interpret the LAS-derived turbulent fluxes (discussed in section 4.2).

Figure 2.  (a) WRF (Weather Research and Forecasting) model domains for land-surface process-resolving simulation. Five domains were configured with a grid nesting ratio of 3. (b) Spatial distribution of land-use types in the innermost domain 5. Urban patches were classified to a low density residential area in the urban canopy parameterization used in this study.

Figure 3.  Temporal variation of (a) daytime wind speed and direction and (b) air temperature measured at the meteorological tower on 11 and 12 November 2013.

The weather conditions on 11 and 12 November 2013 were fairly good, with little influence by fog and clouds. As shown in Fig. 3, the prevailing winds were northeasterly or northwesterly during the daytime on both days, except for the morning hours on 12 November. The wind speed changed within the range of 0.5-4 m s^-1 during the daytime, being relatively low in the morning hours (0.5-2 m s^-1) and relatively high in the afternoon (2-4 m s^-1), while the air temperature showed a typical temporal variation for both the days, ranging from 0^°C to 9^°C.

Based on the meteorological conditions (Table 1), 3D LAS footprint modeling was conducted to interpret the LAS-derived turbulent fluxes. In this study, the flux footprints were calculated every 10 m in discrete locations along the LAS optical path, and then the estimates were multiplied by the LAS weighting function to model the 3D LAS footprints. As shown in Fig. 4, the calculated LAS footprints clearly revealed the source regions of the LAS-derived turbulent fluxes over the heterogeneous urban area. The source regions of the LAS-derived fluxes depended highly on the angle between the wind direction and the LAS beam path, as well as the meteorological conditions (e.g., atmospheric stability, wind speed, friction velocity). When the wind direction was almost parallel to the LAS optical path (Figs. 4a and c), the calculated footprint indicated that 90% of the daytime turbulent fluxes were originated from areas within a width of 100-300 m along the optical path (0.2-0.6 km²). It is also clear that the width of influence became narrower as the atmospheric stability decreased (more unstable). On the other hand, the daytime turbulent fluxes of 90% were attributed to larger areas, expanding up to about 1 km (∼2 km²) upwind from the LAS optical path, when the wind direction was almost perpendicular to the optical path (Figs. 4b and d). The results also show that the LAS-derived turbulent fluxes were the least influenced by sources near the LAS transmitter/receiver, as discussed by (Chehbouni et al., 2000).

Figure 4.  Representative 3D LAS (large aperture scintillometer) footprints for different time periods on 11 and 12 November 2013. Contour lines (20%, 40%, 60%, 80%, and 90%) represent the cumulative source attribution of the LAS-derived turbulent fluxes (%). The thick arrows indicate mean wind direction over the time period.

By virtue of the 3D footprint analysis, it was found that the LAS system, operating with an optical path length of 2.1 km and effective height of about 26 m, can represent the area-averaged turbulent fluxes over this heterogeneous urban area on a micro-α scale, i.e., 0.2-2 km. The results imply that the LAS's spatial-integration capability will be very useful to evaluate modeled surface-layer turbulent fluxes over significantly heterogeneous surfaces on the spatial scale.

Figure 5.  Temporal variation of LAS (large aperture scintillometer)-derived and WRF (Weather Research and Forecasting) model-simulated momentum fluxes over the heterogeneous urban area on 11 and 12 November 2013. "WRF" and "WRF (no river)" indicate the base simulation result and sensitivity simulation result (in which the Geum River was replaced with "shrubland"), respectively. The vertical bars indicate the temporal standard deviation values of the LAS-derived momentum fluxes at 30 min time intervals.

Figure 6.  Temporal variation of LAS-derived and WRF (Weather Research and Forecasting) model-simulated sensible heat fluxes over the heterogeneous urban area on 11 and 12 November 2013. The vertical bars indicate the temporal standard deviation values of the LAS-derived sensible heat fluxes at 30 min time intervals.

The measured C_n² ranged from 0.1 to 3.5× 10^-14 m^-2/3 and 0.1 to 2.8× 10^-14 m^-2/3 during the daytime periods (0800-1700 LST) on 11 and 12 November, respectively, thus satisfying the saturation-free condition described in section 3.2. The LAS-derived turbulent sensible heat and momentum fluxes were calculated every 15 s from the C_n² values following the procedure described in section 2, and then the fluxes were averaged over 30 min intervals. Figure 5 shows that the LAS-derived momentum fluxes (τ) on 11 and 12 November 2013, which were well correlated with the mean wind speed profiles shown in Fig. 3, ranged from 0.05 to 0.35 N m^-2 during the daytime. The land-surface process-resolving WRF simulation also compared well to the estimation, indicating the LAS estimation is able to represent the spatial-integrated momentum fluxes over heterogeneous surfaces. The ratio of friction velocity to mean wind speed (u_*/U) varied from 0.14 in the early morning to 0.11 at later times on 11 November, and from 0.22 during morning to 0.10 in the afternoon on 12 November. The u_*/U difference between the morning and the afternoon is attributed to the temporal variation of atmospheric stability, following Eq. (8).

Figure 6 shows that the LAS-derived daytime sensible heat flux (Q _H) varied between 25 and 160 W m^-2 on 11 November, and between 20 and 120 W m^-2 on 12 November. The LAS-derived Q _H quantities were different in magnitude and temporal variation between the two days, with a dramatic reduction at around 1400 LST 12 November. Temporal fluctuations in the LAS-derived turbulent fluxes reflected the temporal variability of the measured scintillation intensity for the daytime periods. The calculated temporal variability (1σ) in Q _H was estimated to be as much as about 21% relative to the 30 min averaged values on 11 November, and 22% on 12 November. The temporal variability in τ and its footprint dependency were not as large as in Q _H (Fig. 5) because the momentum fluxes were not directly associated with the scintillation intensity.

The LAS-derived heat fluxes were in good agreement with the land-surface process-resolving simulation results, especially during morning hours (0800-1200 LST). Further comparison to the net radiative fluxes measured at the meteorological site shows that the estimated daily mean Q _H (107 and 72 W m^-2) corresponded to about 49% and 37% of the net radiation of 219 and 194 W m^-2 on 11 and 12 November 2013, respectively. Despite the intrinsic difficulty in the quantitative verification of the LAS-derived turbulent fluxes over heterogeneous surfaces, the comparison results of the LAS-derived Q _H to the WRF simulation and the flux ratio analysis show that the scintillation method can reasonably estimate area-averaged surface-layer turbulent fluxes.

The distinctive features of the LAS-derived Q _H temporal variation are that the peak Q _H (especially on 11 November) occurred earlier in the day and that morning Q _H was enhanced compared to that in the afternoon (asymmetric distribution). Generally, a peak Q _H value in the surface-layer energy balance is found at 1300-1500 LST on clear days (e.g., Lee and Park, 2008; Lee, 2011; Loridan et al., 2013). In addition, the maximum values tend to lag in time over urban areas compared to natural surfaces (Lee and Baik, 2011). To further investigate the features in the LAS-derived Q _H, another WRF simulation was conducted by replacing the Geum River with a natural surface type ("shrubland"). The WRF simulation revealed that the peak time lag and the morning enhancement in the LAS-derived Q _H can be partly attributed to the existence of the river, where significantly larger Q _H values (∼200 W m^-2) were found in morning hours (0800-1000 LST) due to the large temperature difference between the water body (∼14.0^°C) and the overlying atmosphere in the cold season.

The LAS system produced surface-layer turbulent fluxes over this heterogeneous urban area within a reasonable uncertainty range (1σ). In this section, the results from several sensitivity tests are discussed to further quantify the influence of input parameters (surface aerodynamic parameters and Bowen ratio) and surface-layer similarity function for the temperature structure function on the LAS-derived turbulent fluxes.

The surface aerodynamic parameters (z₀ and d) could be estimated by an average of the surface parameters of each homogeneous patch consisting of heterogeneous source areas, as described in section 2.1. However, in practice it is difficult to determine the surface parameters for highly heterogeneous surfaces because of the complexity of the morphology of buildings/obstacles and the topography. The uncertainty in the surface aerodynamic parameters also influences the determination of the LAS effective height (z _eff). As a realistic ballpark guess of the surface parameters for the study area (e.g., Garratt, 1992; Foken, 2008), an aerodynamic roughness length of 0.1 m and 1.0 m and an LAS effective height of about 20% (21 and 31 m) from the base estimation were applied in the sensitivity tests. Besides the surface aerodynamic parameters, the influence of the Bowen ratio was also considered in the sensitivity tests, with values of β=0.5 and β=2.0.

Figure 7.  Sensitivity of LAS-derived momentum fluxes to the Bowen ratio (0.5 and 2.0), aerodynamic roughness length (0.1 m and 1.0 m), and LAS effective height (21 m and 31 m). Base estimation was made using a Bowen ratio of 1, aerodynamic roughness length of 0.3 m, and LAS effective height of 26.1 m.

Figure 8.  Sensitivity of LAS-derived sensible heat fluxes to the Bowen ratio (0.5 and 2.0), aerodynamic roughness length (0.1 m and 1.0 m), and LAS effective height (21 m and 31 m).

Figure 7 shows that the estimated τ was very sensitive to z₀, leading to a reduction of about -34% for z₀=0.3 m but an increase of about +80% for z₀=1.0 m, compared to the base estimation on both days. The influence of the Bowen ratio on the LAS-derived τ was negligible, while the change in z _eff led to momentum-flux difference within 5%. On the other hand, as shown in Fig. 8, the estimated Q _H was more variable to the change in z _eff, ranging between 17% on a daytime average. The change in z₀ led to the flux difference of between -6% and 14% on 11 November, which was lower than that of between -4% and 9% on 12 November. The difference in β led to the smallest (<4%) changes in Q _H. Overall, the sensitivity tests indicate that the uncertainty in the surface parameters may affect the LAS-derived Q _H within 20%.

To quantify the uncertainty by the choice of the temperature structure function parameter, the LAS turbulent fluxes were estimated using four different non-dimensional similarity functions for the temperature structure function parameter. When calculating turbulent sensible heat flux using the scintillation method, previous studies have used various non-dimensional similarity functions (e.g., Meijninger and De Bruin, 2000; Lagouarde et al., 2006; Evans et al., 2012; Geli et al., 2012). Along with the surface-layer structure function suggested by (Andreas, 1988), the function by (De Bruin et al., 1993) has been popularly used, which is given for unstable conditions (z/L<0) as follows: \begin{equation} f_T=4.9\left(1-9.0\dfrac{z}{L}\right)^{-2/3} . \end{equation} (Foken and Kretschmer, 1990) and (Thiermann and Grassl, 1992) proposed different functional forms as follows: \begin{equation} f_T=\left(\dfrac{0.95}{k}\right)^2\left(1-11.6\dfrac{z}{L}\right)^{-1/2} , \end{equation} where k=0.4, and \begin{equation} f_T=6.34\left[1-7\dfrac{z}{L}+75\left(\dfrac{z}{L}^2\right)\right]^{-1/3} , \end{equation} respectively. By comparing the non-dimensional similarity functions for unstable atmospheric conditions, it was found that the function by (Andreas, 1988) had higher f_T than that by (De Bruin et al., 1993), except for near neutral conditions (-z/L<0.01), while it had relatively lower f_T than the functions proposed by (Foken and Kretschmer, 1990) and (Thiermann and Grassl, 1992) for a wide range of unstable atmospheric conditions (Fig. 9).

As shown in Fig. 10, the Q _H from the Andreas function was lower by ∼14% (range: 8%-17%) in daytime average than that by the De Bruin function, but larger by about 9% and 10% than those from the functions of (Foken and Kretschmer, 1990) and (Thiermann and Grassl, 1992), respectively. The estimated Q _H quantities were well correlated in time with each other, with a correlation coefficient of 0.99.

Figure 9.  Comparison of non-dimensional similarity functions for the temperature structure function parameter for unstable conditions as a function of atmospheric stability.

Figure 10.  Comparison of LAS (large aperture scintillometer)-derived sensible heat fluxes using different non-dimensional similarity functions for the temperature structure function parameter. Plotted are the 30 min averaged daytime sensible heat fluxes for the two clear days of 11 and 12 November 2013. Q _{H_AN} on the x-axis indicates the sensible heat fluxes derived using the function by Andreas (1988).

The influence of the temperature similarity functions on τ were less than 2% in daytime mean flux (not shown). It is worth noting that (Lagouarde et al., 2006) reported the EC-derived Q _H roughly ranged between the LAS-derived Q _H using the De Bruin and Andreas functions for a reasonably homogeneous urban area. Overall, the results indicate that the uncertainty caused by the non-dimensional similarity function for the temperature structure function parameter is within 15% in the LAS-derived Q _H.

This study tested the ability of the LAS technique in determining surface-layer turbulent sensible heat and momentum fluxes over a heterogeneous urban area, in comparison with land-surface process-resolving numerical modeling. The LAS system was deployed and operated over a highly heterogeneous urban area (Gongju, South Korea) during the cold season with an optical path length of 2.1 km. Three-dimensional LAS footprint modeling was also performed to interpret the LAS-derived turbulent fluxes over the heterogeneous surfaces. The WRF model, which was configured to explicitly represent the land-use heterogeneity over the experimental site and associated physical processes, was introduced for verification of the LAS-derived turbulent fluxes.

The LAS-derived turbulent fluxes (Q _H and τ) showed reasonable temporal variation during the daytime on the two clear days of 11 and 12 November 2013 when compared to the results of the land-surface process-resolving modeling and radiative flux measurements. A series of sensitivity tests showed that the overall uncertainty in the daytime Q_H estimate had 20%-30% as an upper limit: within 20% in surface aerodynamic parameters and 15% in the non-dimensional similarity function for the temperature structure function parameter. In addition, the estimated error in τ was very sensitive to the determination of the aerodynamic roughness length. It is noted that only two days of data were used in the analysis, which may also have contributed uncertainties to some extent. On the other hand, the footprint modeling analysis showed that the LAS-derived surface-layer turbulent fluxes estimated with the 2.1 km optical path distance were contributed by source areas within 0.2-2 km² ("micro-α scale") depending on local meteorological conditions. Unlike the flux footprint modeling of the EC system, the LAS footprint (or source area) was significantly influenced by an incident angle of wind to the LAS beam path, as well as by atmospheric conditions such as atmospheric stability and wind speed.

Recently, mesoscale meteorological and environmental models have undergone rapid development by the inclusion of urban canopy models (e.g., Masson, 2000; Kusaka et al., 2001; Martilli et al., 2002; Lee and Park, 2008; Oleson et al., 2008; Porson et al., 2010; Lee, 2011; Ryu et al., 2011) for better representation of urban land-surface processes (e.g., Lee et al., 2011). Along with numerical model development, quantitative validation of surface-layer turbulent fluxes modeled by 3D mesoscale meteorological models is very important to improving forecast skill, as well as understanding boundary layer structures. Even though previous studies have successfully determined the surface-layer turbulent fluxes over various homogeneous surfaces using the EC system and/or LAS system (e.g., Lagouarde et al., 2006; Zeweldi et al., 2010; Evans et al., 2012; Ward et al., 2014), few studies have provided estimations over heterogeneous urban surfaces. The present study demonstrates that the LAS technique can provide realistic surface-layer turbulent fluxes over highly heterogeneous urban areas within a reasonable uncertainty range. In addition, it shows that LAS footprint modeling can provide valuable insight for model-measurement comparisons by the interpretation of sources' contributions to the LAS estimates. Furthermore, this study indicates that LAS-derived turbulent fluxes can be useful for the verification of mesoscale meteorological and environmental modeling conducted over heterogeneous surfaces.