Determining Adequate Averaging Periods and Reference Coordinates for Eddy Covariance Measurements of Surface Heat and Water Vapor Fluxes over Mountainous Terrain

Two coordinate rotation approaches (double and planar-fit rotations) and no rotation, in association with averaging periods of 15 480 min, were applied to compute surface heat and water vapor fluxes using the eddy covariance approach. Measurements were conducted in an experimental watershed, the Lien-Hua-Chih (LHC) watershed, located in central Taiwan. For no rotation and double rotation approaches, an adequate averaging period of 15 or 30 min was suggested for better energy closure and small variations on energy closure fractions. For the planar-fit rotation approach, an adequate averaging period of 60 or 120 min was recommended, and a typical averaging period of 30 min is not superior to that of 60 or 120 min in terms of better energy closure and small variations on energy closure fractions. The Ogive function analysis revealed that the energy closure was improved with the increase of averaging time by capturing sensible heat fluxes at low-frequency ranges during certain midday hours at LHC site. Seasonal variations of daily energy closure fractions, high in dry season and low in wet season, were found to be associated with the surface dryness and strength of turbulent development. The mismatching of flux footprint areas among flux sensors was suggested as the cause of larger CF variations during the dry seasons as that indicated by the footprint analysis showing scattered source areas. During the wet season, the underestimation of turbulent fluxes by EC observations at the LHC site was attributed to weak turbulence developments as the source area identified by the footprint analysis was closer to the flux tower than those scattered in dry season.


InTRoDuCTIon
The eddy-covariance (EC) approach has been widely applied for studying exchanges of mass, energy, and momentum between terrestrial ecosystem and the atmosphere.This approach has been recognized as the most reliable technique for directly measuring surface fluxes over time scales ranging from hours to years.Variations in latent heat flux, sensible heat flux, or carbon flux over various land cover types, including cropland (Anderson et al. 1984), grassland (Verma et al. 1989), and forest (Verma et al. 1986;Valentini et al. 1991;Hollinger et al. 1994) have been studied over short periods due to the limitations of sensors in the early 1990s.Despite improvements to EC instruments (Moncrieff et al. 1997), the EC method still has significant difficulties in preserving energy closure.Energy closure uncertainties ranging from 10 -30% of available surface energy reported in numerous studies (Greco and Baldocchi 1996;Aubinet et al. 2000;Baldocchi et al. 2001;Wilson et al. 2002) and can seriously degrade the accuracy of long-term observations.Many factors in the field measurements, such as instrument setup, surface heterogeneity, flow stationarity, and the lack of a persistent footprint area, can affect EC energy closure for surface energy budget (Twine et al. 2000).Typical research issues of an averaging period (Finnigan et al. 2003), reference coordinates (Wilczak et al. 2001), data/instrument quality controls (Foken and Wichura 1996;Mauder et al. 2007b), and horizontal/vertical advections (Lee 1998;Aubinet et al. 2000;Baldocchi et al. 2000;Oncley et al. 2007) were investigated to improve the energy closure of EC measurements.Recent studies (Göckede et al. 2004(Göckede et al. , 2008;;Rebmann et al. 2005) suggested that flux footprint analysis should be included as one of the standard EC data processing to improve the accuracy of eddy covariance measurements.Regarding the energy closure issue, several studies (Mauder et al. 2007a;Foken 2008;Foken et al. 2010) indicated that the possible cause of underestimation of EC turbulent fluxes may due to the generation of secondary circulations at the site having the larger length scale of surface heterogeneity than the length scale of flux footprint, which are inherent not to be correctly measured by a signal tower.
To improve the accuracy of surface flux observations using the EC technique, one must investigate the uncertainties of different averaging periods and various reference coordinates at an individual flux observation site.Typically, using 30 min as the average interval when computing mean turbulence statistics is recommended for flat terrain (Stull 1988); however, a tall canopy or hilly terrain enhanced micro-scale atmospheric motions near the surface can contribute low frequency signal to turbulence and alter fluxes (Mahrt 1998).Therefore the 30 min averaging period may not be perfectly suitable to capture total flux for any site using the EC method.In recent years, some studies (Sakai et al. 2001;Finnigan et al. 2003;von Randow et al. 2004;Foken et al. 2006;Cava et al. 2008) suggested that an averaging period of 120 to 240 min is preferred for improving the energy closure of the surface fluxes over a complex terrain or rough surface.In this study, six averaging intervals, 15, 30, 60, 120, 240, and 480 min, were selected as the candidates of adequate averaging periods for surface fluxes calculations.
At a fixed point, only a local wind vector can be defined using the tower measurement; thus streamline-wise coordinates are preferred (Finnigan et al. 2003) because multiple sensors can rarely be aligned relative to each other with sufficient accuracy (Lee et al. 2004).The coordinate rotation has been often applied to fulfill the assumptions of 1D homogenous flow and to correct sensor tilt error due to instrument setup with procedures of aligning instrument coordinate to the mean wind direction and forcing the mean vertical velocity to zero in each averaging period (Tanner and Thurtell 1969;Wesely 1970;Kaimal and Finnigan 1994;Lee et al. 2004).However, a perfect site with a horizontal homogeneous condition is rare.In practice, most research sites are located in mountainous areas with mixed vegetation types.These characteristics are challenging for eddy flux measurements.Moreover, rotating the reference coordinate system to the streamwise over a complex terrain may generate errors due to sensor tilt, electronic offset, flow distortion, and adjective air flow, thereby increasing the uncertainties of turbulence flux measurements.Two coordinate rotation approaches, double rotation (DR) and planar-fit rotation (PFR), and no rotation (NR) with the designed averaging periods were investigated for determining adequate averaging periods and reference coordinates.The PFR was applied to consider wind variation in sloping terrain by calculating mean streamline on a plane from an ensemble observation period over 1 week or longer (Paw U et al. 2000;Wilczak et al. 2001;Chen and Li 2007;Yuan et al. 2007Yuan et al. , 2011;;Ono et al. 2008;Zuo et al. 2009).A nonzero mean vertical velocity may exist for individual averaging period, thus one can remove the systemic bias in turbulence flux calculations caused by tilt errors and instrument offset.
In this study, eddy covariance data, including wind velocity, water vapor, and temperature, collected in 2008 and 2009 were analyzed for improving energy closure fractions with adequate averaging periods and rotating coordinates.The effects of applying different averaging periods and reference coordinates on surface flux estimations over a mountainous terrain are investigated and identified by the Ogive function (Desjardins et al. 1989).Flux footprint analysis (Hsieh et al. 2000) was employed to provide logical explanations for the seasonal variations of energy closure fractions at the LHC site.

Site Description
This study was conducted in an experimental watershed, the Lien-Hua-Chih (LHC) watershed, located in central Taiwan (23°55'52''N, 120°53'39''E).Surface elevations vary between 700 and 800 m above sea level.A 22-m-high observation tower was built inside the No. 5 catchment.This study area is characterized by gentle rolling terrain with various tree species and a humid subtropical climate.The topographic contours and the tower position were shown in Fig. 1.Small streams dissect valleys and upland areas drain from the north to south.
In Taiwan, a typical dry season with less rainfall is from October to April and only the northeastern coastal areas of the island might have considerable precipitation (about 40% of annual rainfall) due to northeasterly winter monsoons.During the wet season (May -September), precipitation is mainly brought by the typhoons originating in the Western Pacific Ocean and South China Sea and partially from the southwesterly summer monsoon (Chen and Chen 2003).Average annual rainfall at the LHC study site is 2316.5 mm yr -1 of which 72% occurs during the wet season (Hwong 2010).The diurnal pattern of wind direction persistently changed from northeast (NE) counterclockwise to southwest (SW) in daytime and varied from southwest (SW) counterclockwise to eastward during the night.
The vegetation at the study site comprised warm-totemperate mountain rainforest of mixed evergreens and hardwoods, including Cryptocarya Chinensis, Engelhardtia Roxburghiana, Tutcheria Shinkoensis, and Helicia Formosana.Averaged canopy height was about 17 m.Some pteri-dophytes are commonly observed in lower canopy.The leaf area index measured by the LAI-2000 ranged at 2.5 -4.5 around the tower during different growing periods.The soil type is loam with an average bulk density of 1.29 g cm -3 and porosity of 0.5 for top soil within 1.0 m depth.

Instruments Setup and Data Collection
This study used an integrated hydrometeorology observation system consisting of a set of high-frequency instruments and several slow-response sensors mounted along the tower or buried in the ground as listed in Table 1.Air temperature and humidity profiles were measured at 5 m in-tervals from ground surface to the top of the tower at 20 m, where net radiation, wind speed/directions, rainfall and air pressure were measured.One soil heat flux plate was placed at one point -5 cm near the tower.Three Sentek capacitance probes for measuring soil moisture (θ s ) profile (at -10, -30, -50, -70, and -90 cm) were placed at three different locations in the catchment (see Fig. 1), and thermo couples for measuring soil temperature profile (at -10, -30, -50, -70, and -90 cm) was placed at one point near the flux tower.
Long-term meteorological data measured by slow-response instruments were recorded at 30 min intervals averaging from 1 min samplings.The EC system comprises a 3-D sonic anemometer (R. M. Young 81000) and a krypton  hygrometer (KH20) mounted with 30 cm separation on a mast at a height of 25 m, and a data logger (CR1000) with a 2G flash card for storing raw EC data.The output signals from the sonic anemometer were recorded in a digital format connected to the serial port of the CR1000 logger to avoid signal drift errors.The output signals from the KH20 hygrometer were recorded in an analog format.All highfrequency signals were recorded at a rate of 10 Hz for subsequent off-line flux calculations.A total of two years (2008 -2009) eddy data were analyzed in this study.

EC Data Processing
A computer program was written in FORTRAN for analyzing the EC data offline.This program comprised several options for calculating EC fluxes, including three coordinate rotation options, spiked data removal based on the criteria of 5 standard deviations in a fixed 5-min window, block average calculation, WPL corrections (Webb et al. 1980), integral turbulence characteristics test (ITC test, Foken and Wichura 1996), stationary test, Ogive function test, and footprint analysis (Hsieh et al. 2000).The ITC test is based on turbulent flux similarity assuming the nature turbulent fluxes was maintained the same statistic characters with friction velocity at different stability conditions.The coefficients of the ITC testing function used for momentum, sensible heat and latent heat under different stability were the same as those suggested in Foken et al. (2004).Detailed descriptions of coordinate rotations are given in the next section.Data collected on days with rainfalls were discarded.Only data of sunny days were remained for EC flux calculations for ensuring better data quality.A total of 71 days in the year 2008 and 70 days in the year 2009 are available for EC flux calculations.For year 2008, 24 and 47 days in wet and dry seasons, respectively, are available.For year 2009, 18 and 52 days in wet and dry seasons, respectively, are available.The reason of fewer days available in wet seasons is due to frequent summer convective rainfall in central Taiwan.Different block average intervals were conducted at 15, 30, 60, 120, 240, and 480 min for comparisons.

Coordinate Rotations
Two coordinate rotation approaches, DR and PFR, were applied for transformation of velocities measured by the EC system in this study.On the other hand, the third approach, NR, will not perform any coordinate rotation to observed wind velocities.The DR approach is the most commonly used rotation scheme (Tanner and Thurtell 1969;Wesely 1970;Kaimal and Finnigan 1994), which first sets the mean lateral wind speed to zero ( v 0 = ) by rotating the x-y plane about the z-axis to identify the first rotation angle -^h] and then sets mean vertical wind speed to zero ( w 0 = ) by rotating the new x-z plane about new y-axis to determine the second rotation angle [ tan u w -^h]; u, v, and w are longitudinal, lateral, and vertical wind velocities, respectively; subscripts m and 1 denote wind velocities originally measured and after the first rotation, respectively.The over bar indicates average value in a given block average interval.As a result, the DR approach aligns the new xaxis with the mean streamline and gives zero mean vertical velocity in each averaging period.
Assuming the mean flow field and sensor tilt are not changed over a prolonged experimental period, the PFR approach determines a regressed plane to correct the tilt angle and instrument offsets (Paw U et al. 2000;Wilczak et al. 2001).The regression plane is determined by the following equation .
wm is the measured mean vertical wind speed for each averaging period; b 0 , b 1 , and b 2 are the coefficients of the plane function obtained by multiple linear regression from long term mean wind filed observations.Here, we derived these parameters from 2008 dataset (71 days) and the value of b 0 , is -0.05 m s -1 ; b 1 is 0.0364; b 2 is -0.0469.This equation indicates that the contribution of mean horizontal wind velocity on mean vertical wind velocity is linear.Since the mean vertical velocity is zero over the entire observations, the subtraction of the residual mean vertical velocity (b 0 ) is required for each averaging period to ensure not to contribute to the Reynolds stress (Wilczak et al. 2001).
Based on the orthogonal properties of the u, v, and w directions, pitch angle tan b -^h can be derived.Wind velocities over this regressed plane (u p , v p , w p ) can be calculated by the transform matrix with the pitch and roll angles.Then, the third rotation angle tan v u p p 1 c = -^h, where u p and v p are mean longitudinal and lateral velocities over the regressed plane in each averaging period of the mean flow direction can be determined by setting the mean lateral wind velocity to zero.Unlike the DR approach directly determining the rotation angle (θ), the PFR approach removes the sensor tilt and instrument offset in advance then finds the rotation angle (γ) on a regression plane.

SuRFACE EnERgy buDgET
The surface energy balance of forest canopy can be given as where R n is net radiation flux, LE is latent heat flux, H is sensible heat flux, ΔQ/Δt is the changing rate of heat storage in the surface layer of forest, G is ground heat flux, and each term in Eq. ( 1) has the unit of W m -2 .To account for the buoyancy effects on turbulent flux estimations, the WPL algorithm (Webb et al. 1980) was applied to correct the la-tent heat flux (LE).The changes of air temperatures, water vapor densities, and soil temperatures above the ground heat flux plate were considered for calculating the heat storage term as where a t is the density of air (kg m -3 ); c p is specific heat capacity of air (J kg -1 K -1 ); v t is the density of water vapor (kg m -3 ) λ v is the latent heat of vaporization (J kg -1 ); s t is the density of soil (kg m -3 ); c s is the heat capacity of soil (J kg -1 K -1 ) calculated using an empirical function [c s = θ s + 0.46 (1φ -X 0 ) + 0.6 X 0 , de Vries 1963; θ s is volumetric water content, φ is porosity, and X 0 is the percentage of organic matter]; Δz is the difference of measurement height (m); ΔT a /Δt, ΔT s /Δt, and Δ v t /Δt are change rates of air temperature (K s -1 ), soil temperature (K s -1 ), and water vapor density (kg m -3 s -1 ), respectively, and n is the number of canopy layers, dz is the thickness of soil layer above the heat flux plate (dz = 5 cm).Since air temperatures and relative humilities were measured at 0, 5, 10, 15, and 20 m, the value of n is 4 and Δz = 5 m.

The Energy Closure Fraction
The energy closure fraction (CF) is commonly used to evaluate the quality of EC data by examining surface energy conservation (Aubinet et al. 2000;Foken et al. 2004;Barr et al. 2006).The CF value of surface energy budget is given as In this study, the CF values were used to help determine the adequate coordinate rotation approaches and averaging periods required for calculating latent heat and sensible heat fluxes over mountainous terrain.

Averaging Period Evaluation
Desjardins et al. (1989) introduced a cumulative spectrum density function to investigate turbulent fluxes characteristics for airborne eddy covariance platform.This function, named Ogive function, was proposed as a test to determine that whether all low frequency parts of turbulent eddies are included to calculate EC fluxes (Foken et al. 2004(Foken et al. , 2006)).The Ogive function is calculated as the cumulative integral of the co-spectrum starting with the highest frequencies and is defined as where Co w, x is the co-spectrum of a turbulent flux, w is vertical wind velocity, x is the horizontal wind component or scalar considered, f 1 is the Nyquist frequency, and f 0 = (1/t) is the lowest resolvable frequency.Here, t is the duration of the time series (s).The Ogive function of all frequencies equals the kinematic flux of the corresponding time series to be examined.The Ogive functions also show an asymptotic shape from highest toward the lowest frequencies, suggesting that all flux-carrying scales are constrained within the sampling duration.This study calculated each Ogive curve by a 2-hour fixed duration within a diurnal course from 00:00 -24:00 local time with the representative sunny days collected during the study period (i.e., 71 days in 2008 and 70 days in 2009).

Footprint Analysis
A hybrid analytical footprint model (Hsieh et al. 2000) was used for footprint estimation.The footprint function (f t ) is defined as where D and P are similarity constants, k = 0.4 is the von Karman constant, x is the horizontal coordinate, L is the Obukhov length, and z u is a length scale defined as a function of the measurement height (z m ) and roughness height (z 0 ), z u = z m [ln (z m /z 0 ) -1 + (z 0 /z m )].The peak footprint contribution location (x m ) can be derived as which permits explicit estimation of the peak footprint contribution location as a function of atmospheric stability and z u .Hourly flux was directly allocated into the x m and the "footprint density" was further evaluated as the cumulative flux divided to the total measured flux at 5 m × 5 m spacing over the domain of interest.

Micrometeorology
The distribution of the wind speed and direction in the day and at night are presented in Fig. 2. The wind direction was mainly northwest and northeast during the day and mainly south and northeast at night.The daytime mean horizontal wind speed was in an order of 1.0 -2.0 m s -1 .However, most of the nighttime wind had smaller magnitude, in an order of 0.5 m s -1 , than that in the daytime.Figure 3 depicts the variation of daily averaged surface soil moisture (from 0 to 20 cm) from three selected locations and indicates notable seasonal variations at the LHC site during the study period (71 days in 2008 and 70 days in 2009).During dry seasons, the value of soil moisture can be lower than 0.3 as those of the 1 st -90 th Julian days and the 300 th Julian day to next year.The wet season soil moisture can maintain values above 0.3.diminished at 150 W m -2 due to a limitation of available soil moisture.

Variations of Daily Energy Closure Fractions by Different Averaging Periods and Coordinate Rotations
Figure 5 depicts the variations of daily energy closure fractions by using different coordinate rotations and averag-ing periods for dry season, wet season, and all seasons.For all three different coordinate rotation approaches, the CF values gradually increased when applying longer averaging periods from 15 -60 min.This agreed with the findings in recent studies (e.g., Sakai et al. 2001;Malhi et al. 2002;Finnigan et al. 2003;von Randow et al. 2004;Barr et al. 2006;Cava et al. 2008).No significant enhancement of CF values can be obtained by further increasing the averaging  The CF values in the dry season were higher than those in the wet season, which contradicts recent studies (Turnipseed et al. 2002;Barr et al. 2006;Tanaka et al. 2008) showing the average CF may reach the unity or higher in wet season.The reason for this seasonal trend will be detailed in section 3.5.
Among three coordinate rotation approaches investigated at this study, only the PFR approach can maintain relatively stable CF values when the averaging periods exceeding 60 min while others show significant decrease of CF values as the averaging period becomes greater than 60 min for both the dry and wet seasons.The calculated pitch and roll angles both are less than 3° (α = -2.08°;β = -2.69°)which is very small compared with the topographic variation.However, the PFR approach did correct the tilt offset because the instrument setup had a positive wind speed offset.Table 2 shows the computed daily averaged H and LE fluxes with various averaging periods and coordinate rotation approaches for all seasons.
For longer averaging periods of 120 -480 min, all three approaches introduced large variability in turbulent flux estimations (i.e., higher standard deviations of CF in Table 2) due to the stationary assumption cannot be preserved as the increase of averaging periods.Among three coordinate rotation approaches investigated at the LHC site, the PFR approach has relatively small variation compared with the other two approaches as shown in Table 2.The PFR approach obtains the CF results with stable variation for averaging periods from 15 to 60 min, while other approaches enlarge variation on CF as increasing averaging period.In addition, the PFR approach obtains nearly a constant CF value for large averaging periods, which is very different from other approaches (NR and DR).
The adequate combination of averaging periods and coordinate rotations for EC flux measurements should assure surface energy balance and capture dominant scales of turbulent fluxes.By increasing averaging periods, the CF can be improved as the low-frequency component of turbulent fluxes can be included for flux calculations (detailed in the next section).As numbers and trends of CF values given in Table 2 and Fig. 5, respectively, only the PFR approach presents a consistent CF trend toward the unity with relatively smaller variations, while both NR and DR approaches show downward trends with larger variations (see Table 2 and Fig. 5).Besides, CF values of both NR and DR in dry season (Fig. 5a) are much higher than those of PFR in dry season causing a better compensation of overall CF values (all seasons) toward the unity as shown in Fig. 5c, while three approaches have similar lower CF values in wet season.
Therefore, it is better to recommend several possible combinations for each rotation method at the LHC site.As shown in Table 2, 15 or 30 min is an adequate averaging period for the NR and DR methods with a CF value close to the unity and small standard deviations.For the PFR method, 60 or 120 min is recommended and a typical 30 min may not be good choice due to low CF values.For all three rotation approaches an averaging period of 480 min is not Table 2. Average daily integrated latent heat (LE) and sensible heat (H) fluxes with different averaging periods and coordinate rotations for all seasons, the (std) stands for the standard deviation of daily energy closure fraction (CF), and B is the Bowen ratio for whole available data from 2008 -2009 (71 + 70 days).
Note: (1) Available Energy: R n -G -ΔQ/Δt = 9.822, units in (MJ m -2 ); ( 2) CF values marked with star were the recommended averaging periods for each rotation.recommended due to bad energy closure and higher variations of CF values.The cause of the underestimation using typical averaging period with PFR approach will be further examined by the Ogive function on kinematic sensible heat and latent heat fluxes in next section.

Effects of Increasing Averaging Period Examined by the ogive Function
In order to investigate the enhancement of CF values with an increased averaging period on EC flux calculations, the Ogive function was applied to analyze the turbulent fluxes captured at the local time scale.Figures 6 and 7 show the Ogive curves depicted by every 2-hour average at the local time for kinematic latent heat and kinematic sensible heat fluxes, respectively, with the PFR approach.With an averaging time of 30 min, a significant increase of kinematic fluxes (i.e., rising of the Ogive curves for midday hours) can be observed in both figures, which means most turbulent eddies can be sufficiently resolved.However, further rising of the Ogive curves for the kinematic sensible heat fluxes in midday hours, especially from 10:00 to 14:00 local time, can be detected after the 30 min integration time as shown in Fig. 7. Unlike the response of sensible heat flux to larger averaging time, the further increase of the Ogive function for the latent heat flux mainly observed at early morning hours (i.e., 08:00 to 10:00 local time) as shown in Fig. 6.As analyzed herein, the enhancement of CF values with 60 min averaging time was mainly contributed by the improvement of sensible heat flux captured at midday hours.The current finding was in good agreement with the previous study of Sun et al. (2006) indicting that a half-hour averaging period might underestimate low-frequency transport.Moreover, we found that the improvement of sensible heat flux captured during midday hours was the major contribution of improvement in energy closure with the extended averaging time of 60 min.
In Twine et al. (2000), several methods were proposed for dealing with the residual energy to balance the surface energy budget equation.For example, the residual energy can be allotted to sensible heat flux and latent heat flux by preserving a fixed Bowen ratio, estimated based on original EC fluxes, to achieve the surface energy closure.However, the fixed Bowen ratio approach would not be appropriate because of the different contribution of low frequency range to flux between sensible and latent heat fluxes, as different increasing trends at low frequency range shown in Figs. 6  and 7.As a result, by maintaining a fixed Bowen ratio to obtain surface energy closure may overestimate the latent heat flux and underestimate the sensible heat flux.To sufficiently capture characteristics of kinematic sensible and latent heat fluxes, it is suggested that an averaging time of 60 min should be adequate for calculating turbulence fluxes at the LHC site.

Variation of Energy Closure Fraction
The issue of surface energy closure in EC measurements has been discussed in many studies (e.g., Turnipseed et al. 2002;Barr et al. 2006;Foken et al. 2006;Tanaka et al. 2008).Lacking sufficient turbulent mixing may contribute underestimation of surface fluxes accounting for 5 -10% of net available energy and causes the variation of CF values.Otherwise, as pointed out by Foken (2008), surface heterogeneity may cause large scale advective transport unable to be captured by the EC system to achieve surface energy closure.
Figure 8 shows the hourly CF values with the PFR approach for 15, 30, and 60 min averaging periods at different friction velocities {u u w v w As expected, the surface energy closure cannot be preserved at insufficient turbulent mixing (i.e., u * < 0.35 m s -1 ), no matter what averaging time used for calculating EC fluxes.However, results of using a 60 min averaging period persistently have higher CF values than those of using 15 or 30 min averaging periods when u * < 0.35 m s -1 .The enhancement of energy closure by increasing averaging time from 15 to 60 min can be attributed to the improvement of energy closure at low u * because CF values of the three averaging periods are similar when u * > 0.35 m s -1 .Seasonal variations of daily energy CF was depicted in Fig. 9, indicating low CF values in the wet season and high CF values in the dry season.As shown in Fig. 3, the soil moisture also has significant seasonal characteristics, high in the wet season and low in the dry season at the LHC site.In order to investigate the effects of surface heterogeneity on seasonal CF variation, the soil moisture was classified into four classes (Class 1: θ s < 0.175, Class 2: 0.175 E θ s < 0.275, Class 3: 0.275 E θ s < 0.375, Class 4: θ s > 0.375) in association with CF, wind speed, and net radiation variations (Fig. 10).As shown in Fig. 10a, a decreasing trend of CF values was found as the increase of soil moistures, and higher CF values with larger variations were obtained at low soil moistures (Classes 1 and 2) than those at high soil moistures (Classes 3 and 4).However, the increase of mean net radiation was consistent with the increase of soil moistures (Fig. 10c).The tendency correlation between mean wind speeds and soil moistures was less significant (Fig. 10b).Slightly higher variations of mean wind speed at low soil moistures than those at high soil moistures did attract our attention.A footprint analysis was applied to investigate the effect of wind speed variations on CF calculations.
Wind climatology and meteorological conditions are important to understand the representativeness of measured fluxes at the domain interest and influences of surrounding land use patterns (Göckede et al. 2004;Rebmann et al. 2005); thus flux footprint analysis is often recommended as a tool for flux quality check and quality assurance (Göckede et al. 2008).An explicit footprint estimation model (Hsieh et al. 2000) was applied to investigate the flux source area.Figure 11 depicts the spatial distributions of peak source locations of measured fluxes at 4 different soil moisture classes as aforementioned.The footprint density for each grid cell (5 m spacing) was calculated as the percentage of grid cumulative LE flux to the total LE flux in each soil moisture class.The spatial distributions of footprint density at dry surface state (Classes 1 and 2) were more scattered than those at wet surface state (Classes 3 and 4).Some peak source locations in a dry surface state are more than 50 m away from the flux tower (Figs.11a and b), which enhances the influence of surface heterogeneity on calculating EC fluxes.Unlike scattered source areas of EC fluxes (H and LE), net radiation was measured around the tower across the seasons.Scattered source areas at dry state contribute the cause of high CF variations at dry seasons.In addition, Fig. 10b also shows a lower wind speed compared with the dry state, which implies that the weak turbulent strength could be the cause of the underestimation of CF during wet seasons.
Significant decrease of the CF values at friction velocity of < 0.3 m s -1 indicates the underestimation of EC fluxes under weak turbulent condition.Footprint analysis shows the source area was close to the tower in wet season further support the logical explanation of lower CF values in wet season due to weaker turbulent condition as compared to that of dry seasons.Two main reasons for the seasonal CF variations can be summarized according to above discussions.One is the mismatch between available energy and turbulent fluxes and another is the strength of turbulence development during different seasons.However, more sophisticated footprint models should be applied to quantitatively support this conclusion.Other energy components, such as ground heat flux and heat storage, have less significant effects on CF variations due to their small magnitudes (Fig. 4).

ConCluSIonS
Among three coordinate rotation approaches investigated in this study, the PFR approach has persistent CF val-ues with the lowest variations for averaging periods from 15 -60 min.Based on the frequency analysis of different eddy sizes contributing to sensible heat fluxes and latent heat fluxes, a steep rising of the Ogive functions as increasing averaging windows can be observed mainly before the averaging period of 30 min.By further increasing the averaging period, the enhancement of capturing turbulent kinematic energy was attributed by the sensible heat fluxes during midday hours, while the latent heat flux was insignificant.It was suggested that using a fixed Bowen ratio to achieve energy closure or to fill flux gaps may introduce systematic bias on estimating surface fluxes.Several combinations for each rotation method were suggested for applying eddy covariance approach at this study site.For both NR and DR approaches, the adequate averaging period was 15 or 30 min.However, using an averaging period of 60 or 120 min was recommended for applying the PFR approach and the typical averaging period of 30 min may underestimate the surface turbulent fluxes.For all three rotation approaches, an averaging period of 480 min was not recommended.Another interesting observation was that the seasonal CF values were scattered in dry season and less scattered in wet season.Such seasonal CF variations were correlated with wind speed/direction changes and seasonal soil moisture states at the LHC site.The mismatching of flux footprint areas among flux sensors was suggested as the causes of larger CF variations during the dry seasons as that indicated by the footprint analysis.In addition, the possible cause of underestimating the turbulent fluxes by EC observations was attributed to weak turbulence developments during the wet seasons.However, more sophisticated footprint investigations or numerical eddy simulations are required to reveal the clear relationship between large scale transport and surface fluxes over complex terrain.Takagi, Prof. Shufen Sun, and an anonymous reviewer are greatly appreciated.

Fig. 1 .
Fig. 1.The topographic contour map of the LHC study site.The flux tower (red-circle) coordinates are 23°55'52''N, 120°53'39''E.The star symbols indicate locations of the three soil moisture probes; the arrows indicate the flow directions of nearby intermittent streams.

Figures
Figures 4a and b show the diurnal course of hourly average R n , G, ΔQ/Δt, H, and LE during the wet and dry

Fig. 4 .
Fig. 4. The diurnal cycle of hourly averaged R n , G, ΔQ/Δt, H, and LE at the LHC site.The bar indicates the standard deviation (NR with 60 min averaging period for H and LE fluxes): (a) dry season, (b) wet season.

Fig. 5 .
Fig. 5. Variations of daily energy closure fraction computed with different averaging periods and rotation approaches.The solid blue line with a circle is for the Planar-Fit rotation.The dashed and dotted red line with a square is for the double rotation.The dashed green line with a delta is without any rotation; the bar indicates the standard deviation of daily energy closure for each run: (a) dry season, (b) wet season, and (c) both dry and wet seasons.
-480 min for all rotation approaches in both dry and wet seasons.

Fig. 6 .
Fig.6.The Ogive curves of kinematic latent heat fluxes for every 2 hours local time (PFR with 2 hr averaging period).

Fig. 7 .
Fig. 7.The Ogive curves of kinematic sensible heat fluxes for every 2 hours local time (PFR with 2 hr averaging period).

Fig. 8 .
Fig. 8.The hourly energy closure fraction with PFR at different friction velocity classes.The dashed and dotted line with a circle is the 15 min averaging period.The dashed line with a square is the 30 min averaging period, and the solid line with a delta is the 60 min averaging period.Data are binned at a friction velocity interval of 0.05 m s -1 .

Fig. 10 .
Fig. 10.Relationship between volumetric soil moisture and daily (a) energy closure fraction, (b) wind speed, and (c) net radiation.The symbols represent the binned average for each soil moisture class.The bars show the standard deviation of variables at each soil moisture class.The soil moistures are classified into four classes: Class 1 is less than 0.175; Class 2 is between 0.175 and 0.275; Class 3 is between 0.275 and 0.375; and Class 4 is greater than 0.375 (PFR with 60 min averaging period).

Fig. 11 .
Fig. 11.Spatial distributions of maximum footprint locations for different soil moisture classes: (a) Class 1, (b) Class 2, (c) Class 3, and (d) Class 4. The flux tower is located at the origin; colors give the footprint density values in percentage (%), with the same color scales used for all graphs.The maximum footprint locations are parameterized with atmospheric stability and a length scale along the mean wind direction after Hsieh et al. (2000) (PFR with 60 min averaging period).

Table 1 .
Descriptions of models, sampling rates and measurement heights of instruments installed at the LHC site.