Advanced Search
Article Contents

High-Order Statistics of Temperature Fluctuations in an Unstable Atmospheric Surface Layer over Grassland

doi: 10.1007/s00376-018-7248-x

  • Skewness (S) and kurtosis (K) of temperature in the surface layer over a grassland are investigated under unstable thermal stratifications. We find that both skewness and kurtosis generally obey Monin-Obukhov similarity theory and tend to be constant values (1.5 and 5.3, respectively) when the stability parameter z/L-2. Quantitative formulas of the similarity functions are proposed. The temperature probability density function (PDF) is close to Gaussian in near neutral stratification and non-Gaussian in unstable stratification. The influence of coherent motions on the PDF behavior is analyzed using the quadrant analysis technique. It shows that PDF behaviors are controlled by ejections and sweeps. The results also indicate that the PDF type of the ejections always follows a Gaussian distribution, while the PDF of the sweeps changes with stability.
    摘要: 本文研究了草原下垫面不稳定层结下温度脉动的高阶统计量特征。我们发现,温度脉动的偏度(S)和峰度(K)服从Monin-Obukhov相似理论,并给出了相似函数的定量表达式。在稳定度参数z/L-2时,S和K分别趋近于定值1.5和5.3。温度脉动的概率密度函数(PDF)在近中性层结下,接近正态分布;在不稳定层结下,为非正态分布。利用象限分析方法研究了相干运动对PDF的影响。结果表明,PDF主要受上冲和下扫运动控制,在上冲运动中,PDF始终服从正态分布,而在下扫运动中PDF随稳定度的变化而变化。
  • 加载中
  • Afanas'ev, A. L., V. A. Banakh, A. P. Rostov, 1999: Third- and fourth-order statistical moments of the turbulent fluctuations of wind velocity and temperature.Proc. SPIE Volume 3983, Sixth International Symposium on Atmospheric and Ocean Optics Tomsk, Russian Federation, SPIE, 377-383,
    Andreas E. L.,R. J. Hill, J. R. Gosz, D. I. Moore, W. D. Otto, and A. D. Sarma, 1998: Statistics of surface-layer turbulence over terrain with metre-scale heterogeneity.Bound.-Layer Meteor., 86 379-408,
    Antonia R. A.,A. J. Chambers, and E. F. Bradley, 1981: Temperature structure in the atmospheric surface layer.II. The budget of mean cube fluctuations. Bound.-Layer Meteor., 20 293-307,
    Babić, N., \vZ. Ve\vcenaj, S. F. J. De Wekker, 2015: Flux-variance similarity in complex terrain and its sensitivity to different methods of treating non-stationarity.Bound.-Layer Meteor., 159 123-145,
    Bisignano A.,L. Mortarini, and E. Ferrero, 2017: Evaluation of high-order concentration statistics in a dispersing plume.Physica A: Statistical Mechanics and Its Applications, 474, 115-126,
    Businger J. A.,1973. Turbulent transfer in the atmospheric surface layer.Workshop on Micrometeorology, D. A. Haugen, Ed., American Meteorological Society, 67- 100.
    De Bruin, H. A. R., B. J. J. M. Van Den Hurk, L. J. M. Kroon, 1999: On the temperature-humidity correlation and similarity.Bound.-Layer Meteor., 93 453-468,
    De Franceschi, M., D. Zardi, M. Tagliazucca, F. Tampieri, 2009: Analysis of second-order moments in surface layer turbulence in an Alpine valley.Quarte. J. Roy. Meteor. Soc., 135 1750-1765,
    Desjardins R. L.,J. I. MaCpherson, P. H. Schuepp, and F. Karanja, 1989: An evaluation of aircraft flux measurements of CO2, water vapor and sensible heat. Bound.-Layer Meteor., 47, 55-69,
    Effelsberg E.,N. Peters, 1983: A composite model for the conserved scalar PDF.Combustion and Flame, 50, 351-360,
    Emran M. S.,J. Schumacher, 2008: Fine-scale statistics of temperature and its derivatives in convective turbulence.J. Fluid Mech., 611 13-34,
    Foken T.,B. Wichura, 1996: Tools for quality assessment of surface-based flux measurements.Agricultural and Forest Meteorology, 78, 83-105,
    Foken T.,F. Wimmer, M. Mauder, C. Thomas, and C. Liebethal, 2006: Some aspects of the energy balance closure problem.Atmospheric Chemistry and Physics, 6 4395-4402,
    Garai A.,J. Kleissl, 2013: Interaction between coherent structures and surface temperature and its effect on ground heat flux in an unstably stratified boundary layer.Journal of Turbulence, 14, 1-23,
    Garai A.,E. Pardyjak, G.-J. Steeneveld, and J. Kleissl, 2013: Surface temperature and surface-layer turbulence in a convective boundary layer.Bound.-Layer Meteor., 148 51-72,
    Graf, A., Coauthors, 2010: Boundedness of turbulent temperature probability distributions,and their relation to the vertical profile in the convective boundary layer. Bound.-Layer Meteor., 134, 459-486,,
    Gryanik V. M.,J. Hartmann, 2002: A turbulence closure for the convective boundary layer based on a two-scale mass-flux approach. J. Atmos. Sci., 59, 2729-2744,<2729:ATCFTC>2.0.CO;2
    Jacobs A. F. G.,B. H. J. Van De Wiel, and A. A. M. Holtslag, 2001: Daily course of skewness and kurtosis within and above a crop canopy.Agricultural and Forest Meteorology, 110, 71-84,
    Kader B. A.,A. M. Yaglom, 1990: Mean fields and fluctuation moments in unstably stratified turbulent boundary layers.J. Fluid Mech., 212 637-662,
    Kaimal J. C.,J. J. Finnigan, 1994: Atmospheric Boundary Layer Flows: Their Structure and Measurement.Oxford University Press, 289 pp.
    Katsouvas G. D.,C. G. Helmis, and Q. Wang, 2007: Quadrant analysis of the scalar and momentum fluxes in the stable marine atmospheric surface layer.Bound.-Layer Meteor., 124 335-360,
    Katul G.,D. Poggi, D. Cava, and J. Finnigan, 2006: The relative importance of ejections and sweeps to momentum transfer in the atmospheric boundary layer.Bound.-Layer Meteor., 120 367-375,
    Katul G.,C. I. Hsieh, G. Kuhn, D. Ellsworth, and D. L. Nie, 1997: Turbulent eddy motion at the forest-atmosphere interface.J. Geophys. Res., 102 13 409-13 421,
    Klipp C. L.,L. Mahrt, 2004: Flux-gradient relationship,self-correlation and intermittency in the stable boundary layer. Quart. J. Roy. Meteor. Soc., 130, 2087-2103,
    Lenschow D. H.,J. Mann, and L. Kristensen, 1994: How long is long enough when measuring fluxes and other turbulence statistics? J.Atmos. Oceanic Technol., 11 661-673,<0661:HLILEW>2.0.CO;2
    Li D.,E. Bou-Zeid, 2011: Coherent structures and the dissimilarity of turbulent transport of momentum and scalars in the unstable atmospheric surface layer.Bound.-Layer Meteor., 140 243-262,
    Liang J. N.,L. Zhang, Y. Wang, X. J. Cao, Q. Zhang, H. B. Wang, and B. D. Zhang, 2014: Turbulence regimes and the validity of similarity theory in the stable boundary layer over complex terrain of the Loess Plateau,China. J. Geophys. Res., 119, 6009-6021,
    Liu L.,F. Hu, and X.-L. Cheng, 2011: Probability density functions of turbulent velocity and temperature fluctuations in the unstable atmospheric surface layer.J. Geophys. Res., 116 D12117,
    Mahrt L.,1991: Boundary-layer moisture regimes.Quart. J. Roy. Meteor. Soc., 117 151-176,
    Mahrt L.,2010: Computing turbulent fluxes near the surface: Needed improvement.Agricultural and Forest Meteorology, 150, 501-509,
    Moene A. F.,D. Schüttemeyer, and O. K. Hartogensis, 2006: Scalar similarity functions: The influence of surface heterogeneity and entrainment.Proc. 17th Symp. on Boundary Layers and Turbulence, San Diego, American Meteorological Society, 5. 1.
    Moeng C.-H.,P. P. Sullivan, 1994: A comparison of shear- and buoyancy-driven planetary boundary layer flows. J. Atmos. Sci., 51, 999-1022,<0999:ACOSAB>2.0.CO;2
    Monin A. S.,A. M. Yaglom, 1971: Statistical Fluid Mechanics: Mechanics of Turbulence, Vol. 1. Cambridge, Mass, 769 pp.,
    Nakagawa H.,I. Nezu, 1977: Prediction of the contributions to the reynolds stress from bursting events in open-channel flows.J. Fluid Mech., 80 99-128,
    Nironi C.,P. Salizzoni, M. Marro, P. Mejean, N. Grosjean, and L. Soulhac, 2015: Dispersion of a passive scalar fluctuating plume in a turbulent boundary layer.Part I: Velocity and concentration measurements. Bound.-Layer Meteor., 156 415-446,
    Price J. D.,2001: A study of probability distributions of boundary-layer humidity and associated errors in parametrized cloud-fraction.Quart. J. Roy. Meteor. Soc., 127 739-758,
    Quan L.,E. Ferrero, and F. Hu, 2012: Relating statistical moments and entropy in the stable boundary layer.Physica A: Statistical Mechanics and Its Applications, 391, 231-247,
    Quan L. H.,F. Hu, 2009: Relationship between turbulent flux and variance in the urban canopy.Meteor. Atmos. Phys., 104 29-36,
    Schopflocher T. P.,P. J. Sullivan, 2002: A mixture model for the PDF of a diffusing scalar in a turbulent flow.Atmos. Environ., 36 4405-4417,
    Schopflocher T. P.,P. J. Sullivan, 2005: The relationship between skewness and kurtosis of a diffusing scalar.Bound.-Layer Meteor., 115 341-358,
    Schotanus P.,F. T. M. Nieuwstadt, and H. A. R. De Bruin, 1983: Temperature measurement with a sonic anemometer and its application to heat and moisture fluxes.Bound.-Layer Meteor., 26 81-93,
    Stiperski I.,M. W. Rotach, 2016: On the measurement of turbulence over complex mountainous terrain.Bound.-Layer Meteor., 159 97-121,
    Stull R. B.,1988: An Introduction to Boundary Layer Meteorology.Kluwer Academic Publishers, 666 pp.
    Tillman J. E.,1972: The indirect determination of stability, heat and momentum fluxes in the atmospheric boundary layer from simple scalar variables during dry unstable conditions. J. Appl. Metero., 11, 783-792,<0783:TIDOSH>2.0.CO;2
    Träeumner, K., T. Damian, C. Stawiarski, A. Wieser, 2015: Turbulent structures and coherence in the atmospheric surface layer.Bound.-Layer Meteor., 154 1-25,
    Van De Boer, A., A. F. Moene, A. Graf, D. Schüettemeyer, C. Simmer, 2014: Detection of entrainment influences on surface-layer measurements and extension of Monin-Obukhov similarity theory.Bound.-Layer Meteor., 152 19-44,
    Varshney K.,K. Poddar, 2011: Experiments on integral length scale control in atmospheric boundary layer wind tunnel.Theor. Appl. Climatol., 106 127-137,
    Ve\vcenaj, \vZ., S. F. J. De Wekker, 2015: Determination of non-stationarity in the surface layer during the T-rex experiment.Quart. J. Roy. Meteor. Soc., 141 1560-1571,
    Vickers D.,L. Mahrt, 1997: Quality control and flux sampling problems for tower and aircraft data. Journal of Atmospheric and Oceanic Technology, 14, 512-526,<0512:QCAFSP>2.0.CO;2
    Wang L. L.,D. Li, Z. Q. Gao, T. Sun, X. F. Guo, and E. Bou-Zeid, 2014: Turbulent transport of momentum and scalars above an urban canopy.Bound.-Layer Meteor., 150 485-511,
    Wang S. Y.,Y. Zhang, S. H. Lü, H. P. Liu, and L. Y. Shang, 2013: Estimation of turbulent fluxes using the flux-variance method over an alpine meadow surface in the Eastern Tibetan Plateau.Adv. Atmos. Sci., 30 411-424,
    Wilczak J. M.,S. P. Oncley, and S. A. Stage, 2001: Sonic anemometer tilt correction algorithms.Bound.-Layer Meteor., 99 127-150,
    Wyngaard J. C.,2010: Turbulence in the Atmosphere.Cambridge University Press, 393 pp.
    Wyngaard J. C.,A. Sundararajan, 1979: The temperature skewness budget in the lower atmosphere and its implications for turbulence modeling,Turbulent Shear Flows I, F. Durst, B. E. Launder, F. W. Schmidt, and J. H. Whitelaw, Eds., Springer-Verlag, 319-326.
    Wyngaard J. C.,O. R. Cotè, and Y. Izumi, 1971: Local free convection, similarity, and the budgets of shear stress and heat flux. J. Atmos. Sci., 28, 1171-1182,<1171:LFCSAT>2.0.CO;2
    Young G. S.,1988: Convection in the atmospheric boundary layer.Earth-Science Reviews, 25, 179-198, 0012-8252(88)90020-7
    Young G. S.,D. A. R. Kristovich, M. R. Hjelmfelt, and R. C. Foster, 2002: Rolls, streets, waves, and more——A review of quasi-two-dimensional structures in the atmospheric boundary layer. Bull. Amer. Meteor. Soc., 83, 997-1001,<0997:RSWAMA>2.3.CO;2
    Yu B.,A. G. Chowdhury, and F. J. Masters, 2008: Hurricane wind power spectra, cospectra, and integral length scales. Bound.-Layer Meteor., 129, 411-430,
    Zuo J. Q.,J. P. Huang, J. M. Wang, W. Zhang, J. R. Bi, G. Y. Wang, W. J. Li, and P. J. Fu, 2009: Surface turbulent flux measurements over the Loess Plateau for a semi-arid climate change study.Adv. Atmos. Sci., 26 679-691,
  • [1] NIU Shengjie, ZHAO Lijuan, LU Chunsong, YANG Jun, WANG Jing, WANG Weiwei, 2012: Observational Evidence for the Monin-Obukhov Similarity under All Stability Conditions, ADVANCES IN ATMOSPHERIC SCIENCES, 29, 285-294.  doi: 10.1007/s00376-011-1112-6
    [2] QIAN Cheng, YAN Zhongwei, Zhaohua WU, FU Congbin, TU Kai, 2011: Trends in Temperature Extremes in Association with Weather-Intraseasonal Fluctuations in Eastern China, ADVANCES IN ATMOSPHERIC SCIENCES, 28, 297-309.  doi: 10.1007/s00376-010-9242-9
    [3] LIU Yonghe, FENG Jinming, MA Zhuguo, 2014: An Analysis of Historical and Future Temperature Fluctuations over China Based on CMIP5 Simulations, ADVANCES IN ATMOSPHERIC SCIENCES, 31, 457-467.  doi: 10.1007/s00376-013-3093-0
    [4] Wu Renguang, Chen Lieting, 1995: Interannual Fluctuations of Surface Air Temperature over North America and Its Relationship to the North Pacific SST Anomaly, ADVANCES IN ATMOSPHERIC SCIENCES, 12, 20-28.  doi: 10.1007/BF02661284
    [5] Yanting LIU, Yang ZHANG, Sen GU, Xiu-Qun YANG, Lujun ZHANG, 2023: A Cross-Seasonal Linkage between Arctic Sea Ice and Eurasian Summertime Temperature Fluctuations, ADVANCES IN ATMOSPHERIC SCIENCES, 40, 2195-2210.  doi: 10.1007/s00376-023-2313-5
    [6] Al-Jiboori M. H., Xu Yumao, Qian Yongfu, 2000: Local Similarity Relationships of Non-Dimensional Wind and Temperature Gradient in the Tower-Layer Atmosphere over Beijing City, ADVANCES IN ATMOSPHERIC SCIENCES, 17, 636-648.  doi: 10.1007/s00376-000-0025-6
    [7] S. S. Dugam, S. B. Kakade, 1995: Short-term Climatic Fluctuations in North Atlantic Oscillation and Frequency of Cyclonic Disturbances over North Indian Ocean and Northwest Pacific, ADVANCES IN ATMOSPHERIC SCIENCES, 12, 371-376.  doi: 10.1007/BF02656986
    [8] Liu Yangang, 1995: On the Generalized Theory of Atmospheric Particle Systems, ADVANCES IN ATMOSPHERIC SCIENCES, 12, 419-438.  doi: 10.1007/BF02657003
    [10] Zhe-Min TAN, Lili LEI, Yuqing WANG, Yinglong XU, Yi ZHANG, 2022: Typhoon Track, Intensity, and Structure: From Theory to Prediction, ADVANCES IN ATMOSPHERIC SCIENCES, 39, 1789-1799.  doi: 10.1007/s00376-022-2212-1
    [11] Gu Wei, Wu Rongsheng, 1992: The Theory Study of the Influence of the Topography on the Cold Frontal Motion, ADVANCES IN ATMOSPHERIC SCIENCES, 9, 167-172.  doi: 10.1007/BF02657507
    [13] GUO Xiaofeng, ZHANG Hongsheng, CAI Xuhui, KANG Ling, LI Wanbiao, DU Jinlin, 2007: Discrepancy and Applicability of Various Similarity Functions in Flux Calculations Under Stable Conditions, ADVANCES IN ATMOSPHERIC SCIENCES, 24, 644-654.  doi: 10.1007/s00376-007-0644-2
    [14] HU Dingzhu, TIAN Wenshou, XIE Fei, SHU Jianchuan, and Sandip DHOMSE, , 2014: Effects of Meridional Sea Surface Temperature Changes on Stratospheric Temperature and Circulation, ADVANCES IN ATMOSPHERIC SCIENCES, 31, 888-900.  doi: 10.1007/s00376-013-3152-6
    [15] HUANG Danqing, QIAN Yongfu, ZHU Jian, 2010: Trends of Temperature Extremes in China and its Relationship with Global temperature anomalies Relationship with Global Temperature Anomalies, ADVANCES IN ATMOSPHERIC SCIENCES, 27, 937-946.  doi: 10.1007/s00376-009-9085-4
    [17] Jiang Hao, Wang Keli, 2001: Analysis of the Surface Temperature on the Tibetan Plateau from Satellite, ADVANCES IN ATMOSPHERIC SCIENCES, 18, 1215-1223.  doi: 10.1007/s00376-001-0035-z
    [18] FAN Lijun, Deliang CHEN, FU Congbin, YAN Zhongwei, 2013: Statistical downscaling of summer temperature extremes in northern China, ADVANCES IN ATMOSPHERIC SCIENCES, 30, 1085-1095.  doi: 10.1007/s00376-012-2057-0
    [19] Philip JONES, 2016: The Reliability of Global and Hemispheric Surface Temperature Records, ADVANCES IN ATMOSPHERIC SCIENCES, 33, 269-282.  doi: 10.1007/s00376-015-5194-4
    [20] YE Xingnan, CHEN Tianyi, HU Dawei, YANG Xin, CHEN Jianmin, ZHANG Renyi, Alexei F. KHAKUZIV, WANG Lin, 2009: A Multifunctional HTDMA System with a Robust Temperature Control, ADVANCES IN ATMOSPHERIC SCIENCES, 26, 1235-1240.  doi: 10.1007/s00376-009-8134-3

Get Citation+


Share Article

Manuscript History

Manuscript received: 13 October 2017
Manuscript revised: 08 March 2018
Manuscript accepted: 21 March 2018
通讯作者: 陈斌,
  • 1. 

    沈阳化工大学材料科学与工程学院 沈阳 110142

  1. 本站搜索
  2. 百度学术搜索
  3. 万方数据库搜索
  4. CNKI搜索

High-Order Statistics of Temperature Fluctuations in an Unstable Atmospheric Surface Layer over Grassland

  • 1. State Key Laboratory of Atmospheric Boundary Layer Physics and Atmospheric Chemistry, Institute of Atmospheric Physics, Chinese Academy of Sciences, Beijing 100029, China
  • 2. University of Chinese Academy of Sciences, Beijing 100049, China

Abstract: Skewness (S) and kurtosis (K) of temperature in the surface layer over a grassland are investigated under unstable thermal stratifications. We find that both skewness and kurtosis generally obey Monin-Obukhov similarity theory and tend to be constant values (1.5 and 5.3, respectively) when the stability parameter z/L-2. Quantitative formulas of the similarity functions are proposed. The temperature probability density function (PDF) is close to Gaussian in near neutral stratification and non-Gaussian in unstable stratification. The influence of coherent motions on the PDF behavior is analyzed using the quadrant analysis technique. It shows that PDF behaviors are controlled by ejections and sweeps. The results also indicate that the PDF type of the ejections always follows a Gaussian distribution, while the PDF of the sweeps changes with stability.

摘要: 本文研究了草原下垫面不稳定层结下温度脉动的高阶统计量特征。我们发现,温度脉动的偏度(S)和峰度(K)服从Monin-Obukhov相似理论,并给出了相似函数的定量表达式。在稳定度参数z/L-2时,S和K分别趋近于定值1.5和5.3。温度脉动的概率密度函数(PDF)在近中性层结下,接近正态分布;在不稳定层结下,为非正态分布。利用象限分析方法研究了相干运动对PDF的影响。结果表明,PDF主要受上冲和下扫运动控制,在上冲运动中,PDF始终服从正态分布,而在下扫运动中PDF随稳定度的变化而变化。

1. Introduction
  • Research on probability density functions (PDFs) of scalars in the atmospheric boundary layer is mainly driven by boundary layer modeling and turbulence diffusion simulation. Theoretically, any order closure of the Reynolds equation can be performed with specific PDFs. Several forms of scalar PDF have been proposed, including beta distribution (Graf et al., 2010), pareto distribution (Schopflocher and Sullivan, 2002) and so on (Price, 2001). So far, no consensus has been reached on the scalar PDF due to its non-Gaussian property and complexity. Therefore, high-order moments of scalars attract more attention instead, because they can reflect the characteristics of probability density to a certain extent. Skewness and kurtosis, which are the third and fourth standardized moments, are most commonly used both in theoretical analysis (Schopflocher and Sullivan, 2005; Nironi et al., 2015) and modeling (Gryanik and Hartmann, 2002; Bisignano et al., 2017). Skewness is considered to have relevance to non-local transport, while kurtosis reveals microstructures and processes in turbulence (Quan et al., 2012).

    The scalar temperature has been widely investigated because of its importance in the surface layer. Research on temperature fluctuations have mainly focused on second-order moments, such as sensible heat flux (Quan and Hu, 2009; Zuo et al., 2009; Wang et al., 2013). So far, few works have reported on high-order moments of temperature fluctuations (Tillman, 1972; Antonia et al., 1981; Afanas'ev et al., 1999). Arguments on whether the behavior of temperature skewness can be predicted by Monin-Obukhov similarity theory (MOST) over horizontal homogeneity existed in the early days (Monin and Yaglom, 1971; Businger, 1973). Later, observation confirmed that MOST is applicable to temperature skewness (Tillman, 1972; Antonia et al., 1981; Andreas et al., 1998). Unfortunately, the relationship between skewness and stability obtained from these observations was based on a small amount of data, indicating the applicability of the above theory to be limited. Moreover, to the best of our knowledge, a quantitative relationship between kurtosis and stability has not yet been established, and to establish such a relationship is one of the main purposes of this paper. On the other hand, it has been widely recognized that the joint probability distribution of vertical velocity and temperature is closely related to coherent motion (Katul et al., 1997). Ejection and sweep motions are considered to be the two major types of coherent eddies. A third-order cumulative expansion of the joint probability distribution, which contains the information of skewness, is often used to predict the difference between flux contributions due to ejection and sweep motions (Nakagawa and Nezu, 1977; Katul et al., 1997). Moreover, how atmospheric stability affects coherent motion still attracts much attention (Wang et al., 2013), and an in-depth understanding of coherent motion will help to establish a refined PDF. All these issues are discussed.

    In this work, we investigate temperature fluctuations in the surface layer over grassland under the framework of similarity theory. There are two main purposes to the paper. One is to set up a quantitative relationship between the stability parameter and high-order statistics of temperature. The other is to investigate the influence of coherent motions on temperature's PDF. The paper is organized as follows: In section 2, the experimental design and data processing are described. In section 3.1, we discuss the applicability of MOST to third- and fourth- moments of temperature. Then, we derive the similarity functions of skewness and kurtosis in section 3.2. In section 3.3, the behavior of the fine structure of the PDF and its relationship with coherent motions is presented using the quadrant analysis technique. Finally, a conclusion is provided in section 4.

2. Data and processing
  • 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.

    Figure 1.  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.

    Figure 2.  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).

  • 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. Methods and results
  • 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 Cnnβ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 C22β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 C33β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 C44β4-4/3=3.8.

    Figure 3.  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).

    Figure 4.  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).

  • 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.

    Figure 5.  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.

    Figure 6.  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.

  • 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.

    Figure 7.  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.

    Figure 8.  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.

    Figure 9.  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).

    Figure 10.  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.

4. Summary and conclusions
  • The high-order statistics of temperature fluctuations are of great significance for the parameterization schemes in modeling and for basic theory of atmospheric turbulence. Observational data over flat and uniform grassland are used to study the skewness and kurtosis of temperature under the framework of MOST. The temperature skewness and kurtosis generally follow the similarity theory with slight scatter. When the thermal stratification is close to neutral, the temperature follows a near Gaussian distribution. As the instability becomes larger, both the skewness and kurtosis increase and approach constant values (1.5 for skewness and 5.3 for kurtosis) in very unstable conditions (ζ<-2). Similarity functions of skewness and kurtosis derived directly from higher-order moment similarity functions show considerable deviations compared with observation. Therefore, we present new forms of similarity functions of skewness and kurtosis. In the unstable condition, the relationship between skewness and the stability parameter follows S=-7.1ζ/(1-4.8ζ), and that between kurtosis and the stability parameter follows K=(3+9.5ζ)/(1+1.8ζ2).

    The fine structure of the temperature PDF under different stratifications is also analyzed, using the quadrant analysis technique. Results indicate a close relationship between the PDF and ejections and sweeps. In unstable conditions, sweeps control the left tail, while ejections control the right tail, of the PDF. The difference in statistical characteristics between ejections and sweeps results in the asymmetry of the PDF. Increasing the duration of sweeps causes the peak of the PDF to become larger with instability. Refined analysis of the temperature PDF reveals that the temperature fluctuations in the ejection movements obey a Gaussian distribution and that the PDF of temperature in the sweep movements varies with the stability. It should be noted that both "statistical factors", such as stationarity (Liang et al., 2014; Babić et al., 2015), and "process factors", such as advection and entrainment (Mahrt, 1991; De Bruin et al., 1999; Moene et al., 2006; Van De Boer et al., 2014), will cause departure from MOST. How to distinguish and quantify these factors is an interesting topic and needs more combined observations of the surface layer and mixed layer.

    This research was supported by the National Key Research and Development Program of China (2017YFC0209605) and the National Natural Science Foundation of China (Grant Nos. 11472272, 41605010 and 41675012).




    DownLoad:  Full-Size Img  PowerPoint