Spatiotemporal characteristics of atmospheric turbulence over China estimated using operational high-resolution soundings

Large-scale in situ observations are sorely lacking, leading to poor understanding of nationwide atmospheric turbulence over China. Nevertheless, high-resolution soundings have become available starting in 2011, providing a unique opportunity to investigate turbulence across China. Here, we calculated the mean turbulence dissipation rate (ϵ) from radiosonde measurements across China for the period 2011–2018 using Thorpe analysis. The atmospheric layers that had stronger turbulence indicated by larger ϵ generally came with larger Thorpe length but with smaller Brunt–Väisälä frequency. Overall, the clear-air ϵ in the free atmosphere exhibited large spatial variability with a ‘south-high north-low’ pattern. Large clear-air ϵ values were observed in both the lower stratosphere (LS) and upper troposphere (UT), especially over the Tibetan Plateau (TP) and its neighboring regions with complex terrain likely due to large-amplitude mountain waves. Particularly, less frequent but more intense clear-air turbulence was observed in both lower troposphere (LT) and UT over the TP, while more frequent, less intense clear-air turbulence was found in northern China. The all-sky turbulence considering the moist-saturation effects was much stronger in the troposphere, notably in southern China where convective clouds and precipitation oftentimes dominated. In the vertical direction, the altitude of peak clear-air ϵ in the troposphere was found to decrease poleward, broadly consistent with the meridional gradient of tropopause height in the Northern Hemisphere. A double-peak mode stood out for the profiles of clear-air ϵ at midlatitudes to the north of 30° N in winter: one peak was at altitudes of 15–18 km, and another at altitudes of 5–8 km. The strong shear instabilities around the westerly jet stream could account for the vertical bimodal structures. The seasonality of ϵ was also pronounced, reaching maxima in summer and minima in winter. Our results may help understand and avoid clear-air turbulence, as related to aviation safety among other issues.


Introduction
Turbulence and mixing in the free atmosphere are fundamental physical processes determining the cross-tropopause exchange of mass and energy (Jaeger and Sprenger 2007). It is known that turbulence is ubiquitous and is involved in a myriad of complicated physical processes, including atmospheric pollutant dispersion, cloud microphysics, entrainment, the jet stream, and synoptic weather (Bodenschatz et al 2010). For instance, the structure of turbulence in the LT has important effects on the diffusion of momentum and atmospheric pollutants near the ground (Fritts et al 2012). Furthermore, turbulence is found to dramatically influence the exchanges of mass and energy between the troposphere and stratosphere (Dutta et al 2009). This in turn affects passenger comfort and safety in aviation especially when aircraft vortices are dramatically modulated (Gerz et al 2005). Therefore, an improved understanding of the distribution of turbulence and its underlying formation mechanism is essential for better analyses of atmospheric pollution, forecasting weather, and climate predictions.
To date, the large-scale spatiotemporal characteristics of turbulence has yet to be fully understood, due to the lack of sufficient in situ detection instruments in the atmosphere. As such, there is much room for improvement in turbulence-related parameterizations in global weather and climate prediction models (Arnfield 2003). A wide spectrum of probes attached to aircraft (Shi et al 2015, Wen et al 2015, balloons (Dutta et al 2009, Zhang et al 2012, and ground-based radar (Luce et al 2002, Kantha and Hocking 2011, Fritts et al 2012 have been used to derive turbulence. The radar probing of atmospheric turbulence is fraught with uncertainties because of the absence of simultaneous in situ measurements (Whiteway et al 2003), but it has potential advantages over in situ measurements, especially in the stratosphere (Clayson and Kantha 2008). On a side note, both aircraft-and ground-based radar suffer from limited spatial coverage. By comparison, high-resolution soundings (HRS) from the worldwide radiosonde network provide a unique dataset for obtaining the large-scale distribution of turbulence, although the dataset contains measurement uncertainties in the stratosphere due to solar radiation and the sensors (Sun et al 2010).
Recently, a host of pioneering work (e.g. Luce et al 2002, Clayson andKantha 2008) applied Thorpe analysis (Thorpe 1977) to the estimation of free atmosphere turbulence using HRS, which showed reasonable and promising results, even though Thorpe analysis was previously used to estimate local turbulence in the deep ocean. However, false turbulent layers were oftentimes inferred due to instrument noise when using Thorpe analysis. To address this issue, statistical methods were developed to discriminate between background stability and instrument noise (Wilson et al 2010(Wilson et al , 2011, which was able to dramatically reduce or eliminate the aforementioned uncertainties (Schneider et al 2015, Sunilkumar et al 2015. The turbulence results from Thorpe analysis were generally consistent with those from other instruments or methods, including scanning radar (Kantha and Hocking 2011, Wilson et al 2014, Li et al 2016a, the Leibniz Institute of Turbulence Observations in the Stratosphere instrument (Schneider et al 2015), and numerical simulations (Balsley et al 2008, Fritts et al 2016. Consequently, climatological features of atmospheric turbulence have been recorded across the world from HRS measurements, including those over the Bay of Bengal and the Arabian Sea during the 2006 pre-monsoon season (Alappattu and Kunhikrishnan 2010), at Gadanki and Trivandrum in the Indian subcontinent (Sunilkumar et  nationwide turbulence has not been reported over China. Therefore, one of the objectives in this study is to investigate the turbulence in the free atmosphere over China using long-term HRS measurements. The rest of this paper is structured as follows: section 2 describes the data and method used; section 3 examines the spatiotemporal variability of turbulence in detail, and discusses the potential factors influencing turbulence followed by the key findings summarized in section 4.

High-resolution radiosonde measurements
The dataset used for estimating turbulence comprised HRS measurements from the China Radiosonde observational Network in China operationally run by the China Meteorological Administration (Guo et al 2016). Since the beginning of 2011, the network has been expanded to 120 operational radiosonde sites (figure S1 (available online at stacks.iop.org/ERL/16/ 054050/mmedia)), which provides profiles of temperature, relative humidity, pressure, and wind at 1 s intervals (approximately 5-8 m in the vertical) twice per day: 0800 and 2000 Beijing time (BJT = UTC (Universal Time Coordinated) + 8). The HRS data for the period from 1 January 2011 to 31 December 2018 were used and comprised 788 162 radiosonde profiles. Figure S2 shows the temporal variation of the number of radiosonde stations that were available for estimation of turbulence for each month of the study period. It can be seen that more than 115 stations have valid radiosonde measurements for all months, and all the 120 stations have observations for as long as 18 months. The maximum height that the majority of the balloons (97.5%) reach is higher than 16 km AGL all year round, with the peak height of higher than 18 km in summer (figure S3), which indicates sufficient sounding samples in the height ranges analyzed here. Note that the HRS data with a maximum altitude less than 10 km were discarded. To obtain coherent turbulence dissipation rate profiles, the original soundings were resampled because the vertical ascent velocities of balloons vary greatly. To this end, the original soundings were interpolated with an equally spaced resolution of 10 m using the cubic spline interpolation method.

Determination of turbulence in the free atmosphere using Thorpe analysis
Thorpe analysis has been widely used to calculate turbulence dissipation rate in the free atmosphere from HRS, which differs greatly by clear-sky and cloudy conditions. For the calculation of clear-air turbulence in the free atmosphere, the effect of water vapor could be negligible (Wilson et al 2013). In this case, the potential temperature θ (K) is a conservative variable during reversible adiabatic process, which is calculated from the following formula: where R is the gas constant for dry air, c p is the specific heat of dry air at constant pressure, T is the temperature, and P is the atmospheric pressure. Both T and P were directly measured by HRS. Then, a series of θ values were sorted to obtain a monotonically increasing profile, resulting in a profile of sorted θ (i.e. θ * ). For instance, suppose that the sample at height z n needs be moved to height z m in order to create a stable layer; the resulting displacement d = z m − z n is known as the Thorpe displacement (L D ), and its root-mean-square value within an inversion is called the Thorpe length (L T ). Typically speaking, inversion is considered to occur when the cumulative L D simultaneously satisfies the following formulae: where k is less than n (Wilson et al 2010). As illustrated by Dillon (1982), the fundamental hypothesis of the Thorpe method lies in the L T being linearly proportional to the Ozmidov scale (L O ), which is expressed as: where c is a constant, which has been widely discussed in previous studies (Dillon 1982, Wijesekera et al 1993, Wesson and Gregg 1994, Sunilkumar et al 2015, Zhang et al 2019b. Meanwhile, the Ozmidov scale can be expressed as: Here ε is the turbulence dissipation rate and N represents Brunt-Väisälä frequency. Combining equations (4) and (5), one can derive ε as: where C K equals c 2 (Caldwell 1983, Fer et al 2004, and N is Brunt-Väisälä frequency under unsaturated atmospheric conditions (Wilson et al 2013), which can be calculated using the following formula: where θ * represents sorted potential temperature, g represents the gravitational acceleration, and z represents the height. However, when the moist-saturation occurs, the potential temperature is not conserved due to the latent-heat release in clouds, so the moist-saturation significantly affects the magnitude of the Thorpe scale and the resultant turbulence dissipation rate. To consider the moist-saturation effects, Thorpe analysis should be conducted using re-constructed potential temperature, which is detailed in the supplementary materials.
Notably, the profiles of relative humidity (RH) measured by radiosondes are used to determine whether the atmosphere is saturated (i.e. the presence of clouds) according to the methods as described by Zhang et al (2010). For a given atmospheric profile, the RH values within the altitudes grid from ith to jth are equal to or greater than RHmin, and at the same time there exists a kth (i < k < j) altitude grid with RH(k) ⩾ RHmax (table 1). In this case, this profile is deemed as cloudy or moisture saturated profile. Consequently, we got the spatial distribution of the frequency of the atmospheric profiles containing clouds for the study period from 2011 to 2018 ( figure  S4). On average, 31.8% of the soundings occurred in the presence of clouds, which was roughly congruent with cloud fraction in China as observed from geostationary images (Chen et al 2018).
Specifically, the value of C K does contain some uncertainty, which has been widely discussed in ocean-related research (Dillon 1982, Wesson and Gregg 1994, Moum 1996. Also, recent numerical attempts indicated that C K varied with the temporal evolution of turbulence and tended to increase with time (Fritts et al 2016). Nevertheless, C K was found to exhibit a pronounced lognormal distribution (Wijesekera et al 1993, Mater et al 2015, Sunilkumar et al 2015. As such, the ensemble mean of C K is generally seen as a near-constant value, which makes it possible for the Thorpe analysis and thus for the estimate of turbulence from a statistical perspective. This has been further confirmed by a recent study based on Monte Carlo simulation (Zhang et al 2019b), which showed that C K remained near-constant and its optimal prediction was 0.54. Note that the following analysis is based on the assumption of C K = 0.54, unless noted otherwise.
On the basis of the above-mentioned procedures, figure S5 illustrated the profiles of potential temperature, Thorp length, Brunt-Väisälä frequency, and turbulence dissipation rate at Xisha (16.83 • N, 112.33 • E) at 00 UTC on 11 June 2017, which was under clear-sky conditions. Meanwhile, a case under all-sky conditions (considering moisture saturation effect) was shown in figure S6, which is based on the sounding at Beijing at 11 UTC on 15 June 2011. Interestingly, the enhanced turbulence was found in the cloud layers detected by the RH threshold methods as described above. According to Thorpe analysis, the values of L T can be grouped into three subsets: L T > 0, L T = 0 and L T = invalid. Among them, L T = 0 means there is no overturing, which also suggests that there exists no turbulence, and L T = invalid indicates noise-induced turbulence. Meanwhile, for any profiles of L T , turbulence is deemed to occur only when L T > 0 (Zhang et al 2019a), and the occurrence rate of turbulence (P T ) shown in figures S7 and 8 only represents the ratio of positive L T cases to total cases. Instrument noise cannot be ignored when using soundings, especially in weakly stratified layers. True inversions are distinguished from the false turbulence layers induced by instrument noise using the method detailed in the supplementary material. The measurement uncertainties associated with the soundings in China were previously found to increase with altitude, leading to large uncertainties in the calculation of turbulence in the upper stratosphere (Bian et al 2011). Convection-induced turbulence dominates in the planetary boundary layer (PBL), which is generally below 3 km above ground level (AGL, Guo et al 2019). Given the large biases in estimating ε in the PBL (Scotti 2015), our analysis of turbulence is limited to the free atmosphere ranging from 3 to 20 km AGL for all 120 radiosonde sites in China. The turbulence in the free atmosphere was previously found to exhibit significant vertical stratification and reach a peak near the tropopause (Alappattu and Kunhikrishnan 2010, Ko et al 2019, Zhang et al 2019b. To better describe the features in the vertical, we focused on three height intervals: the lower tropospher (LT), the Uper troposphere (UT), and the lower stratosphere (LS). LT refers to all the samples at the heights ranging from 3 km AGL to half the tropopause height, while UT refers to those samples from half the tropopause height to the tropopause itself. In light of the large variability of tropopause height with latitude and season (Chen et al 2019b), we calculated the tropopause height from HRS, which were then normalized relative to their climatological means. Note that the cold point tropopause height (CPT-H) was used here, which was defined as the altitude with the minimum temperature in the vertical temperature profile (Highwood and Hoskins 1998). Figure 1 shows the horizontal distribution of ε of clear-air turbulence in four seasons over the period 2011-2018. By and large, the mean turbulence in the UT was the strongest, followed by the LT, and the weakest occurred in the LS. In the LT, it was worthy noting that the Tibetan Plateau (TP), characterized by complex terrain, had the largest ε throughout the year, with peak values in spring March-April-May (MAM) and summer June-July-August (JJA). This is generally similar to the highamplitude turbulence observed over western parts of the United States where Rocky Mountains dominate (Ko et al 2019). Conversely, the weakest turbulence over northeast China occurred in winter December-January-February (DJF), and the relatively weaker turbulence over south China tended to occur in summer. In terms of the frequency of occurrence of clearair turbulence in the LT ( figure S7(a)), it was relatively high over northern China compared with southern China and the TP, irrespective of season. A closer look revealed that there was a low frequency of turbulence over the TP (figure S7(a)), which had a large ε, indicating that the turbulence intensity was strong therein.

Horizontal distribution characteristics of atmospheric turbulence
In the UT, the annual mean ε of clear-air turbulence exhibited a 'south-high north-low' pattern with the smallest meridional difference observed in summer. Nevertheless, the peak nationwide mean ε of 2.34 × 10 -4 m 2 s -3 occurred in summer, almost the same magnitude as in the United States (Ko et al 2019, Zhang et al 2019b. Strikingly, the strongest turbulence was observed over the TP in the UT. Overall, the spatial pattern of turbulence frequency ( figure S7(b)) matched the horizontal distribution of ε, except over the eastern TP, where it was significantly lower than over other areas in the UT. The frequency of turbulence also exhibited patterns of seasonal variability similar to that of the horizontal distribution of ε, which was large over northern China in summer. In the LS, the spatial distribution of ε appeared to be more heterogeneous than in the troposphere, albeit with much smaller magnitudes. And ε exhibited the same spatial 'south-high north-low' pattern as that observed in the UT.
Both the intensity (ε) and frequency of clear-air turbulence in the troposphere were much stronger than in the LS (figure 1). This statistical distribution of log 10 ε broadly resembled the HRS-derived distributions reported in previous studies that used the same Thorpe analysis method in the U.S. and Norway (Li et al 2016a, Ko et al 2019, Zhang et al 2019b with smaller ε values, which may be due to the fact that only the atmospheric profiles without clouds were considered in our results. Nevertheless, the magnitude of ε in our results was larger than that obtained from the instruments attached to commercial aircraft which is on the order of 10 -5 m 2 s -3 (Cho and Lindborg 2001, Frehlich and Sharman 2010, Sharman et al 2014, likely because the principles used to estimate turbulence intensity differ considerably. Importantly, the clear-air turbulence in the troposphere over the TP was less frequently detected (figure S7) but with larger magnitude (figures 1 and 2), which could be due to the dominant roles of the breaking of mountain waves therein. Moreover, the 'south-high north-low' spatial pattern of turbulence intensity agreed well with the strong meridional gradient observed for CPT-H (figure 2), corroborating the notion of the dominance of turbulence in the troposphere.
In terms of all-sky turbulence when considering the effect of moisture (figure 3), it was on average found much stronger than the clear-air turbulence in the troposphere (figure 1), irrespective of season. In the LT, much stronger turbulence was found in South China in which the presence of cloud dominated ( figure S4). In summer, the largest ε value was observed in both the LT and the UT, and interestingly, no significant difference could be observed in turbulence intensity between the South and North China. This indicates that moisture saturation effects could be associated with the enhanced turbulence intensity. Also noteworthy is the larger bias in humidity measurements in the presence of clouds by radiosonde (Bian et al 2011, Zhang et al 2018b, only the soundings without saturated layers are used here so that clear-air turbulence is analyzed in the following paragraphs, unless noted otherwise. Figure 4 illustrates the features of the percent distributions of log 10 ε, L T , and N from the eight years of HRS from 120 radiosonde sites sampled within the LT, UT, and LS. The plot of log 10 ε in the UP approximates a Gaussian distribution with a positive skew ( figure 4(a)). The mean values of log 10 ε were −4. 35, −4.12, and −4.82 in the LT, UT, and stratosphere, respectively. This indicates that turbulence intensity was high in the UT, which could also be seen in the results for L T calculations in figure 4(b), although L T was more broadly distributed in the troposphere. As shown in figure 4(c), much smaller Brunt-Väisälä frequency (N) was found in the UT and LT compared with in the LS, indicating lower stability in both the UT and the LT where larger L T was coincidently observed ( figure 4(b)). The statistical results for N were similar to those reported in previous studies (Bellenger et al 2017, Ko et al 2019. The positive skew of the distribution of log 10 ε at altitudes of 8-16 km is likely related to the twomode structure of N. Given that turbulence varies significantly with latitude, we also investigated the seasonal variation in turbulence at three latitude belts: 15 • -30 • N, 30 • -40 • N, and 40 • -50 • N. Figure 5 shows the vertical profiles of the seasonal mean log 10 ε averaged over 120 sites in China and fitted using a cubic spline approximation. Overall, ε exhibited significant seasonal variability, being strongest in summer, followed by spring and autumn, and weakest in winter.

Vertical distribution of atmospheric turbulence
The vertical profiles of ε were found to vary dramatically with latitude in most seasons (figure 5). The height with the strongest turbulence in the mid-latitudes varied seasonally, especially between 30 • N and 40 • N and the maximum turbulence was observed at an altitude of 10-13 km in summer and 6-8 km in other seasons ( figure 5(b)). Regarding the annual cycle of ε, at lower-and mid-latitudes (15 • -40 • N), the strongest ε tended to occur at a height of 12 km in summer, whereas the weakest ε occurred in winter. This seasonality agreed well with that in the United States (Ko et al 2019). The vertical distribution of P T is shown in figure S7. Below the tropopause, P T was relatively higher throughout the year, which was consistent with the stronger average turbulence intensity observed below the tropopause (figure 3). On average, the strongest turbulence intensity was observed in summer, while the weakest turbulence was observed in winter (figure 5(d)), highly consistent with the seasonality of P T (figure S8). As shown in figure S8, the peak P T in summer mainly occurred in the UT, in which synoptic-scale dynamic forcings such as the jet stream tended to frequently occur (Chen et al 2019a).

Potential factors influencing the distribution of atmospheric turbulence
As shown in figure S9, the seasonal propagation of the Western Pacific Subtropical High (WPSH) and East Asian monsoon over space was closely associated with the occurrence of strong convection in summer over eastern China, consistent with previous studies (Luo et al 2013, Li et al 2016b. Coincidently, strong turbulence was observed in summer over China, notably over southeastern China (figure 1). It was known that convectively induced turbulence mainly originated from unstable upper tropospheric thunderstorm outflow (Trier and Sharman 2009). In this case, strong vertical wind shear tended to occur in the outflow regions of Mesoscale Convective Systems frequently occurring over eastern China in summer (Chen et al 2019a), which in turn reduced the Richardson number until it reached a critical threshold value and turbulence was generated (Storer et al 2019). Simultaneously, strong static instability, which was related to the adiabatic cooling in the convective updrafts, facilitated the subsequent formation of turbulence. In this study, we found that the ε under all-sky conditions (i.e. the cases considering the moist-saturation effect) was significantly larger than that of clear-air turbulence, especially in the LT over southern China where the convection dominated, especially in summer. It is reasonable to argue that moist convection made a great contribution to the frequent and strong turbulence observed in southern China. Additionally, the propagating WPSH was closely linked to the summertime convection over southern China, which in turn resulted in the breaking of gravity waves, thereby providing favorable conditions for the occurrence of turbulence.
On the other hand, the high-magnitude turbulence observed in the troposphere over the TP and its neighboring regions (figures 1 and 2) could be likely related to orographically-induced convection (Liu et al 2009, Chen and Bordoni 2014, Guo et al 2014 and mountain-wave breaking (Nastrom and Fritts 1992, Ralph et al 1997, Worthington 1998, Wolff and Sharman 2008, Ko et al 2019, owing to the extremely complex terrain (Wu et al 2007). Additionally, the roles of orographic gravity wave were well recognized in generating turbulence, especially under the conditions of terrain-induced convective instability and unfavorable lower boundary structures (Storer et al 2019). Among others, wind shear was an important factor in the breaking of orographic gravity waves, given its primary contribution to shear instability (Fritts et al 1996, Wurtele et al 1996, Epifanio and Qian 2008. More interestingly, a double-peak mode stood out for the wintertime profiles of turbulence dissipation rates at midlatitudes (figures 5(b) and (c)): one peak was at altitudes of 15-18 km, and another at altitudes of 5-8 km. A less pronounced bimodal structure of turbulent dissipation rate could also be found in spring and autumn to the north of 30 • N. As shown in figure 6, the core of westerly jet stream was at an altitude of about 12 km, irrespective of seasons. Another noteworthy thing is that the zonal wind speed at the jet stream core in winter exceeded 60 m s −1 , much larger than that other seasons. Nevertheless, there existed extremely dense contour lines in winter (figure 6(d)), which meant a large vertical gradient of zonal wind speed occurring above and under the jet stream core. This indicated strong vertical wind shear occurring at least two places in the vertical. As such, these strong shear instabilities around the jet stream could well account for the vertical

Concluding remarks
We examined the spatiotemporal distribution of turbulence dissipation rate (ε) for all sky conditions over China for the 2011-2018 period based on HRS measurements using Thorpe analysis. To our knowledge, the HRS-derived estimates of ε was a first characterization of turbulence in the free atmosphere across China on a national scale, which exhibited large spatial, seasonal, and terrain dependence.
By and large, the ε of clear-air turbulence in the free atmosphere showed a pronounced 'south-high north-low' pattern, which was opposite to the spatial pattern of occurrence frequency of turbulence. In particular, relatively less frequent, but more intense turbulence is observed in both LT and UT over the TP, and more frequent, less intense turbulence in northern China, which provides a key reference for aviation safety on a large scale. Another hot spot of ε was found along the coast of South China where convective clouds and precipitation oftentimes dominated, especially in summer (Zhai et al 2005, Chen et al 2019a. Meanwhile, the turbulence that considered the moist-saturation effects (i.e. under all-sky conditions) was significantly stronger in the troposphere, especially in southern China where clouds appeared more frequently. The former largeamplitude ε was likely due to large-amplitude mountain wave breaking around the TP and the shear instability around the jet stream over northern China, whereas the latter was probably owing to the frequent convection. In terms of vertical variation, strong ε tended to occur in the UT, and much weaker values were found in the LT and the LS, irrespective of season and location. Overall, the vertical distribution of ε exhibited a double-peak pattern, notably at midlatitudes in winter in midlatitude region of China. The strong shear instabilities above and under the westerly jet stream seemed accountable for the vertical bimodal structures.
Moreover, turbulence exhibited a significant seasonality, with maximum intensity and frequency in summer in the UT, which may be related to the frequent convective cloud that is also a main source of the turbulence (Storer et al 2019). Because China is located in the East Asian monsoon region, convection in China has a prominent seasonal cycle with the prevailing monsoon and the seasonal activities of the WPSH. The effect of the convection on turbulence requires more in-depth research. The height of peak ε in the troposphere decreased poleward, consistent with the meridional gradient of CPT-H in the Northern Hemisphere.
As extensively reviewed by (Stohl et al 2003), the extratropical stratosphere-troposphere exchange (STE) generally resulted in laminae in ozone profiles (Dobson 1973) and filamentary structures of water vapor on satellite images (Wirth et al 1997). Turbulence has been well recognized to be one of the crucial factors in the processes of transport, both laterally and vertically, that models failed to reproduce (Cristofanelli et al 2003, Meloen et al 2003. Observations of the intensity and frequency of atmospheric turbulence are needed. Meanwhile, aviation turbulence is one of the drivers behind the high-impact weather events that contribute to the aviation accidents (Gultepe et al 2019), which accounts for over 70% of weather-related accidents at cruising levels of commercial jet airliners. Pilot report (PIREP) is commonly used to report turbulence in flight, but it is subjective and highly dependent on cruise height and flight route. Consequently, PIREP was not proposed to be used for research or forecasting applications (Schultz andPolitovich 1992, Kelsch andWharton 1996). Therefore, the free atmospheric turbulence obtained here provides a key reference for turbulencerelated researches over large scale, and more importantly helps gaining a panoramic picture of atmospheric turbulence in China. Overall, our findings have great implications for the STE of energy and mass and for aviation safety.

Data availability statement
All turbulence data that support the findings of this study are available from the corresponding author upon reasonable request (Dr J Guo, Email: jpguocams@gmail.com). The climatological mean data of turbulence dissipation rate in China are deposited at the 'Figshare' data repository at https:// doi.org/10.6084/m9.figshare.12564071.v1. Part of the radiosonde dataset that used to produce turbulence can also be made available by contacting the National Meteorological Information Center of China (http://data.cma.cn/en/?r=data/ index&cid=565b6087618affa1).
The data that support the findings of this study are openly available at the following URL/DOI: https:// doi.org/10.6084/m9.figshare.12564071.v1.