HTML
--> --> --> -->2.1. Experiments and instruments
The data are derived from the grass-covered boundary layer experiment carried out in the northeast of Xilin Hot (44.124°N, 116.29°E), Inner Mongolia,China.A 100-m tower stands over sufficiently flat and uniform grassland, at an altitude of 1159 m (Fig. 1a). The entire experiment began in May 2009 and ended in April 2010. Both slow- and fast-response instruments are installed on the tower (Fig. 1c). The instruments and their placement heights are shown in Table 1.Figure1. Basic information about the experiment: (a) topographic map (using ASTER Global Digital Elevation Model data) surrounding the observation site (indicated by the plus symbol); (b) photograph near the observation site during July 2009; (c) sketch of the 100-m tower and instruments.
Figure2. Wind roses at three levels for the analyzed period using vanes data: (a) 10 m; (b) 30 m; (c) 50 m. U represents the speed observed by vanes.
Three months (June to August 2009) of data were analyzed. During this period, large areas of ground were covered by dozens of centimeters of grass (Fig. 1b). The data used for turbulence calculation are derived from 10, 30 and 50 m. Although the height of 50 m is in the upper atmospheric surface layer, sometimes exceeding the surface layer, we still use data at this height for obtaining more near-ground information (Garai and Kleissl, 2013). Weak (<3 m s-1) and moderate (3-12 m s-1) intensity winds exist in all directions, with only small amounts of strong winds (>12 m s-1) (Fig. 2).
2
2.2. Data processing
We begin by defining some symbols used in this study. u, v and w are the longitudinal, lateral and vertical velocity components, while T is the air temperature. The turbulent fluctuations are computed following \(X'=X-\bar{X}\), where X is any variable, the overbar denotes the time averaging operator and the primed variable denotes the turbulent part.Before analysis, sonic temperature correction is performed to exclude the impact of water vapor using the relationship T s=T(1+0.51q), where T s is the sonic temperature and q is the specific humidity (Schotanus et al., 1983).
2.2.1. Averaging time
The choice of averaging time is a significant consideration when calculating the turbulence statistics. When a short averaging time is chosen, the statistics may exclude the contributions by turbulent motions in some scales. When an averaging time that is too long is chosen, the statistics may include the contributions by non-turbulent motions, such as mesoscale motions (Mahrt, 2010).
There are several methods that are often used in the selection of the averaging time. (Lenschow et al., 1994) discussed the selection of averaging time from the perspective of error analysis. They pointed out that the statistical error of estimating the moment decreases as the averaging time increases for the stationary stochastic process. When the averaging time is much larger than the integral time scale, the error is mainly determined by the random error, which represents the random scatter of sampling. Specially, for a Gaussian process, once the averaging time ta is chosen, the random error of the time mean of the nth-order central moment can be expressed by an(τ/ta)1/2, where τ is the integral time scale and an is a coefficient that increases with the increasing of order n. It indicates that a longer averaging time is needed to reduce the random error. However, the calculation of the integral time scale τ still depends on the averaging time (Yu et al., 2008; Varshney and Poddar, 2011), which makes this rule not so effective to determine the averaging time.
Another method of determining the averaging time, often used in micrometeorology, is the ogive test (Desjardins et al., 1989). This method is based on analysis of the second-order mixed moment of vertical velocity and temperature. It calculates the cumulative integral of the cross-spectral turbulent flux, expressed as \(\int_{\infty}^{f_0}Co_w'T'(f)df\), where Co is the cross spectrum and f0 is the cut frequency. The averaging time is acceptable when the integration is close to a constant at low frequencies. An interval of 30 minutes is usually determined using this method (Foken et al., 2006). This interval has also been used in some studies on high-order moments of turbulent fluctuations (Jacobs et al., 2001). In this study, we employ an averaging time of 30 minutes. Averaging times of 15 minutes and 1 hour were also used for comparison and no significant differences existed, as shown in section 3.1.
2.2.2. Quality control and coordinate rotation
Data quality control is performed as follows:
(1) Periods of relative humidity greater than 70% are removed to avoid the impact of potential precipitation.
(2) Spikes, amplitude resolution problems, dropouts and discontinuities are detected following (Vickers and Mahrt, 1997). A threshold of 5 times the standard deviation is chosen when detecting spikes in this study, instead of 3.5 times the standard deviation in the original literature to preserve more information in PDF tails for non-Gaussian distributions.
(3) Thresholds are given as (Vickers and Mahrt, 1997) and (Klipp and Mahrt, 2004) recommended. The absolute limit for horizontal wind is 30 m s-1 and for temperature it is -30°C to 50°C. The critical value for vertical gradients of mean velocity is 0.001 s-1.
Detrending and coordinate rotation are carried out after the above steps. We simply use the linear detrending approach (Kaimal and Finnigan, 1994). Coordinate rotation has no impact on calculation of high-order moments of temperature. But, it will influence calculation of the sensible heat flux and stability parameter. The most commonly used methods are the double rotation and planar fit methods (Wilczak et al., 2001). Before analysis, both methods were tested to calculate fluxes and no significant differences were found. The correlation coefficient of heat fluxes using these two methods exceeds 0.99 at each height. We only show the results with double rotation in the following part of the paper. Since linear detrending and double rotation are linear transformation methods, exchanging their operation orders does not affect the results.
2.2.3. Stationarity
Stationarity is a fundamental assumption in the statistical computation of turbulence and similarity theory analysis (Wyngaard, 2010). When the turbulence is stationary, the time-average statistics will be equal to the ensemble-average statistics, which is also called ergodicity. Strictly speaking, turbulence in the atmospheric boundary layer has always been in a non-stationary state, which makes it hard to select an appropriate criterion of stationarity. (Ve\vcenaj and De Wekker, 2015) compared several approaches of detecting non-stationarity and concluded that no single approach is most suitable for detecting non-stationarity of first- and second-order moments of atmospheric time series. They did, however, point out that the criterion proposed by (Foken and Wichura, 1996) performed better than other criteria in their investigation. This criterion detects the non-stationarity of second-order moments (\(\overline{u'u'} \), \(\overline{v'v'}\), \(\overline{w'w'}\), \(\overline{T'T'}\), \(\overline{u'w'}\), \(\overline{v'w'}\), \(\overline{w'T'}\)). Together with detrending, it can screen out the stationary time series effectively. Most intervals (more than 80% in our research) passing the test are still stationary when applying this approach to skewness and kurtosis.
After the above treatments, about 15% to 25% of the periods remain at each level (Table 2). We did not carry out an uncertainty assessment because it would have reduced the intervals that could be analyzed (Stiperski and Rotach, 2016). All the following analysis is based on intervals under unstable conditions.
-->
3.1. Similarity relationship of third- and fourth-order moments
As skewness and kurtosis are the normalized third- and fourth-order moments, we first investigate the behaviors of third- and fourth-order moments under unstable conditions.According to MOST, the scaled mean variable at the height of z in the surface layer is the only function of the dimensionless stability parameter ζ: \begin{equation} \label{eq1} \zeta=z/L , \ \ (1)\end{equation}
where L is the Obukhov length, defined by \begin{equation} \label{eq2} L=-\frac{\overline{\theta_{\rm v}}u_\ast^3}{\kappa g(\overline{w'\theta'_{\rm v}})_s} .\ \ (2) \end{equation} In Eq. (3), θ v is virtual potential temperature; \(u_\ast=(\overline{u'w'}_s^2+\overline{v'w'}_s^2)^{1/4}\) is the friction velocity; \(\overline{w'\theta'_v}\) is the kinematic heat flux; \(\kappa\) is the von Karman constant (0.4 in this article); and g is the gravitational acceleration (9.8 m s-2 in this article). The subscript s represents the ground value. In this article, we regard the value at 10 m as the near ground value. Virtual potential temperature used in this work can be computed as θ v≈ T v+gz/Cp=T s+gz/Cp, where T v is virtual temperature and close to T s (Stull, 1988) and Cp is the specific heat for air at constant pressure (1005 J kg-1 K-1 in this paper).
According to (Kader and Yaglom, 1990), high-order moments of temperature fluctuations in an unstably stratified boundary layer, where thermal production should be taken into account, follow \begin{equation} \label{eq3} \varphi_n=\frac{\overline{T'^n}}{T_\ast^n}=C_n(-\zeta)^{-n/3} .\ \ (3) \end{equation} In Eq. (3), φn represents a set of universal similarity functions. \(T_\ast=-\frac{(\overline{w'\theta'_v})_s{u_\ast}\) is the temperature scale in surface layer. The value of n should be greater than or equal to 2. Cn is the coefficient determined by experiment.
Another form often used can be written as (De Franceschi et al., 2009) \begin{equation} \label{eq4} \varphi_n=\frac{\overline{T'^n}}{T_\ast^n}=\alpha_n(1-\beta_n\zeta)^{-n/3} ,\ \ (4) \end{equation} where αn and βn are the corresponding coefficients.
When \(\zeta\to-\infty\), Eq. (4) is consistent with Eq. (3) and we have Cn=αnβn-n/3. When \(\zeta\to 0\), which means the thermal stratification is near neutral in condition, Eq. (4) tends to be a constant with a value of αn, while Eq. (3) tends toward infinity. We use the form of Eq. (4) to fit the similarity function of temperature fluctuations.
Figure 3 shows non-dimensional high-order moments versus the Monin-Obukhov stability parameter ζ. The fitting results using Eq. (4) are also shown in the figure. Under unstable conditions, high-order (at least fourth-order) moments of temperature obey MOST. Non-dimensional similarity functions of second-, third- and fourth-order moments are given as: \beginsubequations \begin{align} \label{eq5} \varphi_2&=15.5(1-58.8\zeta)^{-2/3} ;\ \ (5a)\\ \label{eq6} \varphi_3&=-10.5(1-8.2\zeta)^{-1} ;\ \ (5b)\\ \label{eq7} \varphi_4&=300.6(1-26.3\zeta)^{-4/3} .\ \ (5c) \end{align} \endsubequations Most existing studies focus only on the relationship between the second-order moment (often in the form of standard deviation) of the temperature and the stability parameter. When \(\zeta\to-\infty \), we have C2=α2β2-2/3=1.02, which is very close to other earlier studies (Wyngaard et al., 1971). Very few studies about third- and forth-order moments can be found within the framework of similarity theory. When \(\zeta\to-\infty\), we have C3=α3β3-1=1.3, which is smaller than the value given by Kader and Yaglom (C3=2 in their work). So far, according to our research, there is no literature on the accurate relationship between fourth-order moments and the stability parameter. The relationship between the fourth-order moments and the stability parameter shows C4=α4β4-4/3=3.8.
Figure3. Non-dimensional high-order moments of temperature as functions of the Monin-Obukhov stability parameter for 10 m (blue circles), 30 m (green circles) and 50 m (red circles): (a) second-order moment; (b) third-order moment; (c) fourth-order moment. Solid lines correspond to Eq. (5).
Figure4. Skewness (S) and kurtosis (K) as functions of the Monin-Obukhov stability parameter ζ for 10 m (blue circles), 30 m (green circles) and 50 m (red circles), using different averaging times: (a, b) 30 min; (c, d) 10 min; (e, f) 1 h. Solid lines are deduced from Eq. (9) and Eq. (10).
2
3.2. Similarity relationship of skewness and kurtosis
Under unstable conditions, we can express skewness (S) and kurtosis (K) as \begin{eqnarray} \label{eq8} S&=&\frac{\langle{T}'^3\rangle}{\langle{T}'^2\rangle^{3/2}}=-\frac{\langle{T}'^3\rangle/T_\ast^3}{\langle{T}'^2\rangle^{3/2}/(T_\ast^2)^{3/2}}= -\frac{\varphi_3}{\varphi_2^{3/2}} ;\ \ (6)\\ \label{eq9} K&=&\frac{\langle{T}'^4\rangle}{\langle{T}'^2\rangle^2}=\frac{\langle{T}'^4\rangle/T_\ast^4}{\langle{T}'^2\rangle^2/(T_\ast^2)^2}= \frac{\varphi_4}{\varphi_2^2} . \ \ (7)\end{eqnarray} Mathematical symbol \(\langle\ \rangle\) means ensemble average. There is a negative sign in Eq. (6) because T* is negative in unstable conditions. Equations (6) and (7) indicate that the behaviors of S and K are controlled by non-dimensional second-, third- and forth-order moments. If second-, third- and fourth-order moments obey MOST, it can be expected that S and K obey MOST.Figure 4 shows S and K as functions of the stability parameter. Though data points are scattered compared to those in Fig. 3, S and K generally obey MOST. When \(\zeta \to 0\), S approaches zero and K approaches three. This implies that the PDF of temperature tends to be Gaussian when thermal stratification is close to neutral. As ζ decreases, S and K increase. Both S and K seem to get near-constant values in strongly unstable condition (ζ<-2). That S approaches a constant in very unstable conditions has also been found by other studies. However, different constants are given by different researchers: 0.82 for (Andreas et al., 1998); 1 for (Wyngaard and Sundararajan, 1979); and 1.3 for Li and Bou-Zeid (2011, their Fig. 11). All these values seem to be smaller compared with our result, which approximately equals to 1.5. The constant for K in very unstable conditions (ζ<-2) approximately equals 5.3 in our investigation. Results using an averaging time of 10 minutes and 1 hour are also shown in Figs. 4c-f. No significant difference is found when using different averaging times.
One notable phenomenon is that lines deduced from Eq. (5) fail to predict an accurate relationship between S (or K) and ζ. The main reason is that the form of Eq. (4) is only an approximation. The coefficient and the form of Eq. (4) may change in different thermal stratification when buoyancy forces and dynamic forces change (Kader and Yaglom, 1990). The error is further amplified when Eq. (4) is used to derive skewness and kurtosis. In order to give a more precise function, we directly fit the S and K. Notice that the solid line in Fig. 4a describes the trend of the scattered points well. Thus, we slightly modify the form derived from Eqs. (5a) and (5b) and give the following form: \begin{equation} \label{eq10} S=\frac{a\zeta}{1-b\zeta} , \ \ (8)\end{equation} where a and b are coefficients determined by observation.
Figure5. Fitted results of (a) skewness and (b) kurtosis. Solid lines correspond to Eq. (13). The dotted line, from (Andreas et al., 1998), follows \(S=0.255\ln(-\zeta)+1.044\).
The solid line in Fig. 4b deviates from observed K considerably and a new form is needed to fit the data. A square relationship has always been found between S and K in scalars (Graf et al., 2010; Quan et al., 2012). A similar behavior is also found in our investigation. This inspired us to suggest that the following function may be a good approximation when describing the relationship between kurtosis and the stability parameter: \begin{equation} \label{eq11} K=\frac{c+d\zeta^2}{1+e\zeta^2} ,\ \ (9) \end{equation} where c, d and e are coefficients. We set c=3, as we find the kurtosis approaches 3 in neutral stability conditions.
Fitting Eqs. (11) and (12) to data, we have \begin{equation} \label{eq12} \left\{ \begin{array}{l} S=\dfrac{-7.1\zeta}{1-4.8\zeta}\\[3mm] K=\dfrac{3+9.5\zeta^2}{1+1.8\zeta^2} \end{array} \right..\ \ (10) \end{equation} The fitted curves are shown in Fig. 5, together with the fitted result used by (Andreas et al., 1998). The two forms of curves agree with each other on the trend, as shown in Fig. 5a. However, the linear form used by (Andreas et al., 1998) failed to predict the constants of skewness in near neutral and very unstable conditions.
Since skewness and kurtosis can be predicted by similarity theory, we speculate that the PDF of temperature also follows MOST. Figure 6 validates this idea. The PDF of normalized temperature fluctuations changes systematically with the stability parameter. In unstable conditions, the temperature PDF shows a non-Gaussian property. As the instability increases, the left tail of the temperature PDF becomes shorter and the right tail becomes longer, which has also been found by (Liu et al., 2011). The skewness (0.04, 0.32, 0.86, 1.22, 1.43 for each PDF) and kurtosis (3.01, 3.04, 3.45, 4.72, 5.11 for each PDF) increase with the increase in instability. Moreover, the peak value of the PDF increases as the instability increases. The shape of the PDF is consistent with experimental results carried out near a wall in the laboratory (Emran and Schumacher, 2008). We discuss the fine structures of these non-Gaussian PDFs in the next part.
Figure6. Normalized PDFs of temperature under different unstable conditions. The solid line is the Gaussian distribution curve. Data from 10 m, 30 m and 50 m are all used here.
2
3.3. Relationship between coherent motions and PDFs of turbulent temperature fluctuations
In this section we study the relationship between coherent motions and the fine structure of the temperature PDF. Quadrant analysis, which has been widely used in studying coherent motions (Katul et al., 2006; Katsouvas et al., 2007; Wang et al., 2013), is performed in this section. Under unstable thermal stratification, according to the sign of the fluctuations of temperature and vertical velocity, the movements can be divided into four quadrants:Quadrant 1: w'>0, T'>0, ejection
Quadrant 2: w'<0, T'>0, inward interaction
Quadrant 3: w'<0, T'<0, sweep
Quadrant 4: w'>0, T'<0, outward interaction
Based on this definition, the PDF of normalized temperature fluctuations can be decomposed into four parts: \begin{eqnarray} \label{eq13} f(T')&=&\int_1f(w',T')dw'+\int_2f(w',T')dw'\nonumber\\ &&+\int_3f(w',T')dw'+\int_4f(w',T')dw' , \ \ (11)\end{eqnarray} where, f(T') is the PDF of normalized temperature fluctuations and f(w',T') is the joint PDF of normalized vertical velocity and temperature; the subscript represents the respective quadrant.
Figure7. Quadrant compositions of PDFs under different unstable conditions: (a) -0.02<ζ <0; (b) -0.1<ζ <-0.02; (c) -0.5<ζ <-0.1; (d) -2<ζ <-0.5; (e) ζ <-2. Black dots represent entire PDFs. Dots with different colors represent contributors of each quadrant to PDFs.
Equation (11) decomposes the temperature PDF into contributions of ejections, inward interactions, sweeps and outward interactions. Figure 7 shows the result of PDF decomposition. For all conditions with different stabilities, the left tails of entire temperature PDFs coincide with those of sweeps' PDFs, while the right tails of entire temperature PDFs coincide with those of the ejections' PDFs. This indicates the left tail of the temperature PDF is mainly controlled by the sweep motions and the right tail of the temperature PDF is mainly controlled by ejection motions. Meanwhile, the peaks of sweeps' PDFs correspond to the peaks of entire temperature PDFs for ζ<-0.5. To illustrate how the temperature PDF changes systematically with stability, we analyze the influence of the stability parameter on the behaviors of ejections and sweeps. The contribution rate of each quadrant to temperature variance can be defined as: \begin{equation} \label{eq14} V_i =\iint_iT'^2f(w',T')dw'dT' , \ \ (12)\end{equation} where Vi denotes the contribution rate of the ith quadrant to variance (i=1,2,3,4). Since the variance of normalized temperature fluctuations is 1, the sum of Vi equals 1. Figure 8 shows the contribution rate of each quadrant. Under unstable thermal stratification, ejections and sweeps contribute most of the temperature variance (more than 70%). When the thermal stratification is near neutral (-0.02<ζ<0), the contributions of ejections and sweeps are almost equal. With the increase in instability, the contribution of ejections increases while the contribution of sweeps decreases. Since large variance relates to large deviation (or amplitude), the amplitude of the ejections increases with the increase in instability, and the amplitude of sweeps decreases with the increase in instability. This explains why the left tail of the PDF becomes shorter and the right tail becomes longer when the stratification becomes more unstable (as shown in Fig. 6). The PDFs of ejections and sweeps become more asymmetrical with the increase in instability. As a result, the skewness becomes larger.
Figure8. Contribution rate of each quadrant to temperature variance Vi as a function of the stability parameter, where the contribution rate Vi is calculated by summing up the squares of the normalized T' in the respective quadrants.
Figure9. PDF components of ejection and sweep motions under different unstable conditions. For clarity, dots are shifted up: 1 unit for -0.1<ζ <-0.02; 2 units for -0.5<ζ <-0.1; 3 units for -2<ζ <-0.5; 4 units for ζ <-2. The vertical dashed line represents T'=0. Solid lines are Gaussian and exponential curves (their functions are presented in Table 3).
We further study the time fraction of each quadrant, which is given by \begin{equation} \label{eq15} D_i=\frac{1}{t_a}\int_0^{t_a}{I_i(t)dt} , \ \ (13)\end{equation} where Di is the time fraction of the ith quadrant and Ii is an indicator function, expressed by \begin{equation} \label{eq16} I_i=\left\{ \begin{array}{l@{\quad}l} 1 & {\rm if} w',T' \hbox{is in quadrant} i\\ 0 & {\rm otherwise}. \end{array} \right.. \ \ (14)\end{equation}
We focus on D1 and D3, which are measures of time fractions of ejection and sweep events. Under the unstable condition, the time fraction of ejections almost doesn't change with the stability parameter and remains constant at 0.29 (Fig. 9). This value is consistent with other studies (Katul et al., 1997). The time fraction of sweeps is larger than that of ejections and increases as the instability increases. The larger the time fraction, the greater the probability of sweep events. This indicates the probability of sweep events is larger than that of ejection events. This explains why the peak of the PDF is controlled by sweeps and why the peak value becomes larger as the instability increases (for ζ<-0.5).
Although the contribution of ejections to variance becomes larger when the stratification becomes more unstable, the PDF of ejections remains Gaussian (Fig. 10 and Table 3).
Figure10. Time fraction of each quadrant Di as a function of the stability parameter. The dotted line represents the mean value (0.29) of the time fraction of the ejections for ζ <-0.02. The value of Di is calculated by dividing the number of T' in the respective quadrants by the total sample number.
In the unstable condition, rising warm air carries heat generated by solar heating from the surface during the ejection event. This warm air forms plumes or streaks, which are the main coherent structures in the surface layer (Moeng and Sullivan, 1994; Young et al., 2002; Garai and Kleissl, 2013; Tr?eumner et al., 2015). Warm air in plumes and streaks has similar statistical properties, and the probability distribution of temperature obeys a Gaussian distribution.
The PDF of sweeps changes with instability, being Gaussian in near neutral conditions and close to exponential in very unstable conditions (ζ<-2). The change of PDF type indicates that the property of the air mass may change with stability (Effelsberg and Peters, 1983). In near-neutral thermal stratification, shear forces play the main role and the shear instabilities occurs locally (Moeng and Sullivan, 1994). In this situation, the property of the air mass in sweep motions is similar to that in ejection motions and follows a Gaussian distribution. Whereas, in strongly unstable conditions, which often occur in a convective boundary layer, downdraughts with cold air from greater height may invade the surface layer (Tillman, 1972; Young, 1988; Graf et al., 2010). Cold air from high altitude differs from that in the surface layer and may cause the non-Gaussian behavior of sweeps. Some researchers have reported that this "high altitude" may reach the top of the boundary layer. For example, (Graf et al., 2010) pointed out that the lowest temperature in the entire boundary layer strongly affects the turbulent PDF, and (Mahrt, 1991) connected skewness with entrainment processes. These processes may also be the reasons that the skewness and kurtosis of temperature show scatter under the MOST framework, as shown in Fig. 5. Therefore, to analyze how such cold air influences the left tail of the PDF, more information on temperature in the mixed layer is needed.