Occurrence frequency of subcritical Richardson numbers assessed by global high-resolution radiosonde and ERA5 reanalysis

Kelvin–Helmholtz instability (KHI) is most likely to be the primary source for clear-air turbulence, which is of importance in pollution transfer and diffusion and aircraft safety. It is indicated by the critical value of the dimensionless Richardson (Ri) number, which is predicted to be $\mathrm{1}/\mathrm{4}$ from linear stability analysis. However, Ri is fairly sensitive to the vertical resolution of the dataset; a higher resolution systematically leads to a finer structure. The study aims to evaluate the performance of ERA5 reanalysis in determining the spatial–temporal variabilities in subcritical Ri by comparing it against a near-global high-resolution radiosonde dataset during the years 2017 to 2022 and further highlights the global climatology and dynamical environment of subcritical Ri. Overall, the occurrence frequency of $\mathit{\text{Ri}}<\mathrm{1}/\mathrm{4}$ is inevitably underestimated by the ERA5 reanalysis over all climate zones at all heights from the near-ground atmosphere up to 30 km, compared to radiosonde, due directly to the severe underestimation in wind shears. Otherwise, the occurrence frequency of Ri<1 in ERA5 is climatologically consistent with that from $\mathit{\text{Ri}}<\mathrm{1}/\mathrm{4}$ in radiosondes in the free troposphere, especially over the midlatitude and subtropics in the Northern Hemisphere and Southern Hemisphere. Therefore, we argue that the threshold value of Ri could be approximated as 1 rather than $\mathrm{1}/\mathrm{4}$ when using ERA5-based Ri as a proxy for KHI. The occurrence frequency of subcritical Ri revealed by both datasets exhibits significant seasonal cycles over all climate zones. In addition, it is positively correlated with the standard derivation of orography at low levels and is exceptionally strong over the Niño-3 region at heights of 6–13 km. Furthermore, a high occurrence of subcritical Ri would likely be accompanied by strong wind speeds and intensive orographic or non-orographic gravity waves.

1 Introduction

Kelvin–Helmholtz instability (KHI) is a common phenomenon in the atmospheric boundary layer and the free atmosphere (Muschinski and Wode, 1998), and its wavelengths and depths span a wide range of scales throughout the atmosphere, varying from a few meters or less to tens of kilometers (Fritts et al., 2011). It contributes to vertical mixing of heat, momentum, and constituents, and it acts to limit the maximum shears, just to name a few (Fritts et al., 2011). KHI and gravity wave (GW) breaking are the most recognized instabilities in stably stratified flows (Fritts and Rastogi, 1985). KHI arises preferentially from micro- and mesoscale wind shear intensification, with the maximal occurrence frequency near synoptic-scale upper-level frontal zones near jet streams, with mountain waves, and above the tops of severe thunderstorms (North et al., 2014). Large wind shears are commonly associated with regions where stability changes rapidly (e.g., near the top of the boundary layer) and with the large wind gradient in jet streams (Grasmick and Geerts, 2020). In a changing climate, wind shears in the North Atlantic upper-level jet stream could be increased (Lee et al., 2019), which may increase clear-air turbulence at cruise altitudes. In turn, KHI can reduce wind shears and alter tracer gradients where turbulence and mixing are most intense (Fritts et al., 2022).

KHI influences depend on the spatial scales at which they lead to turbulence (Fritts et al., 2022). Turbulence is by far the most common cause of serious damage to aircraft (Williams and Joshi, 2013). Convective instability, shear instability, KHI, and GW breaking are known to be the major sources for turbulence (Sharman et al., 2012; Ko et al., 2019; Ko and Chun, 2022; Lazarus et al., 2021). KHI requires a sufficiently large Reynolds number and a Richardson (Ri) number sufficiently below $\mathrm{1}/\mathrm{4}$ to enable KHI formation and subsequent secondary instability leading to turbulence (Fritts et al., 2022). Ri is not a good guide to instability character in general, and $\mathit{\text{Ri}}>\mathrm{1}/\mathrm{4}$ does not assure flow stability for superpositions of mean and GW motions. Despite these caveats, $\mathit{\text{Ri}}<\mathrm{1}/\mathrm{4}$ does provide a reasonable guide to expected local KHI structure in cases where clear Kelvin–Helmholtz billows arise, according to the simulation in the mesosphere and lower-thermosphere region (Fritts et al., 2014). Values of Ri close to zero favor strong instability, deep billows, and relatively intense turbulence, whereas values of Ri closer to $\mathrm{1}/\mathrm{4}$ favor weak instability and shallow billows (Fritts et al., 2011). The Richardson number criterion can be applied as a turbulence diagnostic in numerical model outputs (e.g., Sharman and Pearson, 2017), and it has been used as such in climatological studies on the occurrence of clear-air turbulence (Jaeger and Sprenger, 2007). Kunkel et al. (2019) include a brief discussion on the capability of European Centre for Medium-Range Weather Forecasts (ECMWF) models based on case studies to resolve subcritical Richardson numbers and argue that the threshold value of Ri (Ri_t) taken as 1 might be a good proxy for observed KHI. A very recent study by Lee et al. (2023) also sets Ri_t from 0–1 in their climatology on the upper-troposphere and lower-stratosphere turbulence diagnostics. Moreover, Zhang et al. (2022b) show that over half of turbulence exists below Ri<1 when the environment is beneficial for the development of turbulence.

Turbulent mixing is of crucial importance to mass, energy, momentum transfer, the dispersion of pollutants, and stratosphere–troposphere exchange. In numerical models, turbulent dissipation rate, turbulent diffusivity, and other parameters representing turbulent mixing efficiency are the most basic parameters which need to be accurately parameterized to evaluate the impact of turbulence effect on matter and energy distribution (Gavrilov et al., 2005). However, due to the intermittent nature of turbulence, it is generally not resolved in (global) numerical weather prediction models, even at nowadays common or state-of-the-art horizontal resolutions of the order of tens of kilometers (Sandu et al., 2019), and it presents a challenge for both observation and numerical modeling (Sharman et al., 2012; Homeyer et al., 2014; Plougonven and Zhang, 2014). For this reason, the indices of turbulence, such as large wind shears, small Ri, and the negative squared Brunt–Väisälä frequency, could comprise a great tool to characterize turbulence (Jaeger and Sprenger, 2007).

The Richardson number is estimated by the finite differences across thin layers and is quite sensitive to the vertical resolution of measurements (Haack et al., 2014). Thus, a proper estimation of Ri requires a high-resolution measurement of temperature and wind speed. The near-globally distributed radiosonde site offers a unique opportunity to investigate the climatology of subcritical Ri occurrence frequency. An overview of subcritical Ri occurrence by using near-global high-resolution (10 m) radiosonde data was presented in Zhang et al. (2022b), and a close association between subcritical Ri occurrence frequency and turbulence fraction has been found. However, the global climatology characteristic of subcritical Ri remains most unclear, especially over oceans where the radiosonde network has a poor coverage.

By comparison, ERA5 global reanalysis can provide a seamless coverage of temperature and wind. It is the latest generation of the ECMWF atmospheric reanalysis and is based on the state-of-the-art Integrated Forecasting System (IFS) Cy41r2 (Hersbach et al., 2020; Gu et al., 2023). Its predecessor, ERA-Interim, was found in particular to have wind shear a factor of 2–3 lower than the measurement based on high-resolution radiosondes (Houchi et al., 2010). Moreover, results show that whatever the location and the geophysical conditions considered, biases between ERA-Interim and balloon wind measurements increase as a function of altitude (Duruisseau et al., 2017). Recent studies have suggested that the structure of and variability in the trade winds in the lower troposphere are reasonably reproduced in the ERA5 reanalysis based on the EUREC4A field campaign (Savazzi et al., 2022). However, a similar comparison between ERA5 and high-resolution radiosonde data across a near-global area has largely been undetermined. The proper estimation of wind shear and Brunt–Väisälä frequency is essential for the determination of Ri.

Thus, our objectives are as follows: (1) evaluate the performance of ERA5 at different heights and climate zones in estimating wind shear and small Richardson number occurrence frequencies, in comparison with a large high-resolution radiosonde dataset spanning the years from 2017 to 2022. (2) Based on the validation and comparison results, we pose the question of how to use ERA5 for subcritical Ri estimation. (3) Evaluate the global climatology of subcritical Ri occurrence based on versatile measurements and model products. (4) Evaluate the dynamic environment (GWs and mean flow) of subcritical Ri. To this end, this analysis is organized as follows. Section 2 shows the data and methods used. Section 3 represents the climatological variation in subcritical Ri and its comparison with radiosonde. Section 4 ends with a summary.

2 Data and methods

2.1 High-resolution radiosonde dataset

As described in Guo et al. (2021) and Zhang et al. (2022b), a high-vertical-resolution radiosonde (HVRRS) dataset gained from several organizations was adopted, spanning January 2017 to October 2022 over a total of 5.8 years. The organizations include the China Meteorological Administration (CMA), the US National Oceanic and Atmospheric Administration (NOAA), the Global Climate Observing System (GCOS) Reference Upper-Air Network (GRUAN), the Centre for Environmental Data Analysis of the United Kingdom (CEDA), the University of Wyoming, the Deutscher Wetterdienst, and ECMWF. In total, around 0.95 million radiosonde profiles from 434 radiosonde stations released at regular synoptic times of 00:00 and 12:00 UTC were collected to determine the value of Ri. These profiles were sampled at 0.5 or 1 Hz, corresponding to a vertical resolution of approximately 10 or 5 m. Thus, all the profiles were evenly interpreted to 10 m resolution in the vertical by applying a cubic spline interpolation. In addition, the sounding with the burst height lower than 10 km above ground level (a.g.l.) was directly discarded for further study. Meteorological variables, including temperature and wind speed, were prepared for the Ri estimation.

One of the shortages of radiosonde measurements is its inadequate concentration over the polar and ocean regions (Xia et al., 2021). The geographical distribution of the total profile number of each radiosonde station is demonstrated in Fig. S1 in the Supplement. The radiosoundings released over Europe, the United States, and Australia have good geographical coverages and time durations. Over some islands in oceans (e.g., the Pacific Ocean), there are dozens of stations that can provide high-resolution measurements. In the polar regions, there are around 30 stations.

2.2 ERA5 reanalysis and the collocation procedure

ERA5 is the latest version of ECMWF meteorological reanalysis, benefiting from a decade of developments in model physics, core dynamics, and data assimilation (Hersbach et al., 2020). The wind and temperature fields are modeled by the ERA5 reanalysis on a spatial resolution of 0.25^∘ latitude–longitude and a temporal resolution of 1 h. The reanalysis has 137 model levels, giving a vertical resolution of approximately 300 m in the middle and upper troposphere. The vertical resolution of ERA5 is illustrated in Fig. S2 in the Supplement. Compared to ERA5, the HVRRS does not provide global seamless observations. Thus, the collocation procedure between reanalysis and the HVRRS goes as follows: (1) the matched grid of ERA5 reanalysis is the nearest neighbor of the radiosonde station; (2) the regular synoptic start time of the radiosonde and reanalysis needs to be kept exactly the same; (3) the model level of reanalysis that follows a hybrid sigma–pressure coordinate is converted into geopotential height to match the HVRRS.

In addition, the standard deviations of orography (SDORs) and the gravity wave dissipation due to the effects of stress associated with unresolved valleys, hills, and mountains in ERA5 reanalysis are extracted.

The relative error between HVRRS-based and ERA5-based quantities is estimated by the ratio of deviations between HVRRS- and ERA5-derived quantities to the HVRRS one.

2.3 The occurrence frequency of subcritical Ri and its uncertainty

Based on a linear theory, the threshold Ri (Ri_t) defines the boundary where the airflow changes from stability to turbulence, and it is usually suggested to be $\mathrm{1}/\mathrm{4}$ (Haack et al., 2014). Ri is formulated as

$$\begin{array}{}\text{(1)}& \mathit{\text{Ri}}={\stackrel{\mathrm{‾}}{N}}^{\mathrm{2}}/{\stackrel{\mathrm{‾}}{S}}^{\mathrm{2}},\end{array}$$

where N is the Brunt–Väisälä frequency ($\sqrt{\frac{g}{\mathit{\theta }}\frac{\mathrm{d}\mathit{\theta }}{\mathrm{d}Z}}$), S is the vertical shear of horizontal wind ($\sqrt{(\frac{\mathrm{d}U}{\mathrm{d}Z}{)}^{\mathrm{2}}+(\frac{\mathrm{d}V}{\mathrm{d}Z}{)}^{\mathrm{2}}}$), and the overbar denotes a moving average in a 200 m bin to eliminate the influence of measurement noises and small-scale fluctuations such as turbulence and small-scale waves. Therefore, the wind shear and Brunt–Väisälä frequency are computed at 10 m resolution, and then those estimates are averaged over 200 m (20 points) and squared. More exactly, the averaged parameter at altitude i can be represented as $\stackrel{\mathrm{‾}}{A}\left(i\right)=\frac{\mathrm{1}}{n}{\sum }_{j=i-\mathrm{10}}^{i+\mathrm{10}}A\left(j\right)$, where A denotes wind shear or Brunt–Väisälä frequency and n is the number of vertical bins. In addition, horizontal winds measured under radiosonde at the scale of a few tens of meters are affected by the chaotic movements of the gondola due to the pendulum and to the balloon's own movements (Ingleby et al., 2022). However, it is hard to quantify those movements in the present study.

The Richardson number calculated from Eq. (1) depends on the vertical resolution of the underlying data, as well as on the averaging interval. Ultimately, this influences the estimated occurrence frequency for subcritical Richardson numbers as a proxy for KHI. We resample the HVRRS data to 50 and 100 m and range the length scale of the overbar from 100 to 500 m to diagnose the uncertainties raised by the length scale of segments and the vertical resolution of the dataset. As indicated in Fig. 1, at the same length scale of the overbar, a sparser vertical grid inevitably leads to a lower occurrence frequency of subcritical Ri. For instance, at the length scale set to 100 m, the occurrence frequency of $\mathit{\text{Ri}}<\mathrm{1}/\mathrm{4}$ at 0–2 km above sea level (a.s.l.) decreases from 22 % when vertical resolution is equal to 10 m to 16 % for a vertical resolution of 50 m. Moreover, a longer length scale of the segment generally yields a smaller occurrence frequency. For example, when the vertical resolution of the radiosonde is equal to 10 m, the occurrence frequency at 10–15 km a.s.l. decreases from 9 % when the length scale of segment equals 100 m to 1 % when it equals 500 m. It is interesting to note that the occurrence frequency with a vertical resolution of 50 m and a segment interval of 100 m is a bit larger than that with a vertical resolution of 10 m and a segment of 200 m, possibly implying that a shorter segment interval could be expected for a sparser vertical resolution.

Figure 1 The averaged occurrence frequencies of $\mathit{\text{Ri}}<\mathrm{1}/\mathrm{4}$ at heights of 0–2 (a) and 10–15 km a.s.l. (b), with vertical resolutions ranging from 10 to 100 m and moving-point numbers increasing from 0 to 50. The error bars correspond to the standard deviation. The metrics are counted based on all radiosonde profiles during the years 2017–2022.

2.4 Gravity wave energy

The GW energy is extracted based on the broad spectral method, according to Wang and Geller (2003). In this method, zonal wind (u), meridional wind (v), and temperature (T) consist of background states (u₀, v₀, and T₀) that are determined by applying a second-order polynomial fit (Zhang et al., 2022a) and perturbations. Therefore, total perturbations are derived as

$$\begin{array}{}\text{(2)}& (u^{\prime },v^{\prime },T^{\prime })=(u,v,T)-(u_{\mathrm{0}},v_{\mathrm{0}},T_{\mathrm{0}}).\end{array}$$

The perturbations include measurement noises, Kelvin–Helmholtz waves, GWs, and planetary waves. Only the perturbations with vertical wavelengths of 0.3–6.9 km are considered GWs (Wang and Geller, 2003). The mean vertical wavelength of GWs is about 2 km (Wang et al., 2005), and therefore, the lowermost threshold of 0.3 km could have limited influence on the GW energy. However, the retrieval of the largest wavelength is not well determined, which is acknowledged as the radiosonde's “observational filter” (Alexander, 1998). By applying this band-pass filter, the average gravity wave kinetic energy per unit mass (energy density) and the average potential energy density can be expressed as

$$\begin{array}{}\text{(3)}& {\displaystyle }{E}_{\mathrm{k}}={\displaystyle \frac{\mathrm{1}}{\mathrm{2}}}\left[\stackrel{\mathrm{‾}}{u^{\prime \mathrm{2}}}+\stackrel{\mathrm{‾}}{v^{\prime \mathrm{2}}}\right],\text{(4)}& {\displaystyle }{E}_{\mathrm{p}}={\displaystyle \frac{\mathrm{1}}{\mathrm{2}}}{\displaystyle \frac{g^{\mathrm{2}}\stackrel{\mathrm{‾}}{{\widehat{T}}^{\prime \mathrm{2}}}}{N^{\mathrm{2}}}},\end{array}$$

where g is the gravitational constant; ${\widehat{T}}^{\prime }=T^{\prime }/\stackrel{\mathrm{‾}}{T}$ is the normalized temperature perturbation; and the overbar indicates an averaging over the tropospheric segment, which is chosen as 2–8.9 km for all regions expect the polar region and is selected as 2–7.4 km in the polar region (Wang and Geller, 2003). Eventually, the total GW energy E_t is the sum of E_k and E_p.

3 Results and discussions

3.1 Comparisons of wind shear between HVRRS and ERA5 reanalysis

The variations in vertical shear of horizontal wind speed and the squared Brunt–Väisälä frequency entirely determine the Ri magnitude. Figure 2 provides an overview of the spatial distribution of wind shear at heights of 0–2 and 10–15 km a.g.l. obtained from the HVRRS and ERA5 reanalysis. HVRRS-based wind shears are taken from Eq. (1), with a vertical resolution of 10 m. The ERA5-derived shear at heights of 0–2 km a.g.l. demonstrates a strong spatial variation, and it is clearly dependent on underlying terrains and latitudes (Fig. 2a). For example, large values can most likely be observed along the coastline, which could be attributed to the prevailing sea-breeze circulation. As compared to the HVRRS, these shears are slightly underestimated by 3.30 $\mathrm{m}{\mathrm{s}}^{-\mathrm{1}}{\mathrm{km}}^{-\mathrm{1}}$, based on all sounding measurements (Fig. 2b). Nevertheless, a close association between averaged ERA5-retrieved shears and HVRRS-determined shears can be noticed in terms of geospatial distribution, with a correlation coefficient of 0.48 (Fig. 2b).

Figure 2 The spatial distribution of mean wind shear in ERA5 reanalysis at heights of 0–2 (a) and 10–15 km a.g.l. (c). The overlaid colored circles represent the result in the HVRRS at the same height levels. Each data point represents a vertically averaged value of the wind shear at one radiosonde station during the whole study period. Density plots (b, d) show the correlation between wind shears in HVRRS and ERA5 reanalysis. The ERA5-derived wind shears are spatially and temporally collocated with those of the HVRRS. In addition, the red lines represent a least-squares linear regression, and the star superscripts indicate that values are statistically significant (p<0.05).

It is noteworthy that shears in the ERA5 reanalysis at heights of 10–15 km a.g.l. are substantially underestimated compared to the HVRRS, especially at middle latitudes, with a mean absolute error for all stations of about 8 $\mathrm{m}{\mathrm{s}}^{-\mathrm{1}}{\mathrm{km}}^{-\mathrm{1}}$ (Table 1). The underestimation could partly be due to the coarse vertical resolution (around 300 m) in the ERA5 reanalysis in this height interval. Nevertheless, the geographical distribution of the ERA5 shear still exhibits a significant positive correlation with the HVRRS shear, with a correlation coefficient of 0.44 (Fig. 2d).

Table 1 Comparisons of mean wind shears between HVRRS and ERA5 reanalysis at heights of 0–2 (a), 10–15 (b), and 20–25 km a.g.l. (c).

Figure 3 Monthly mean wind shears during the years 2017–2022 in HVRRS (a–g) and ERA5 reanalysis (h–n) in different climate zones. The ERA5-derived wind shears are spatially and temporally collocated with those of the HVRRS. NH: Northern Hemisphere; SH: Southern Hemisphere.

Following Houchi et al. (2010), the monthly averaged shears over seven typical climate zones are separately investigated (Fig. 3) and are defined as follows: polar (70–90^∘ N or S), midlatitudes (40–70^∘ N or S), subtropics (20–40^∘ N or S), and tropics (20^∘ S–20^∘ N). Over the polar regions in the Northern Hemisphere and Southern Hemisphere, HVRRS-based shears are exceptionally strong in the lower stratosphere compared to those in the troposphere (Fig. 3a and g), which could be attributed to the stratospheric polar jet. However, the similar altitude variation can hardly be found in ERA5-based shears, which are dramatically underestimated by around 12 $\mathrm{m}{\mathrm{s}}^{-\mathrm{1}}{\mathrm{km}}^{-\mathrm{1}}$ in the lower stratosphere (Fig. 3h and n, also seen in Table 1). The results in midlatitudes reach a similar conclusion (Fig. 3b, f, i, and m). In the subtropical region, HVRRS-based shears are consistently strong at heights of 16–21 km a.g.l., just above the subtropical jet stream (Fig. 3c and e). However, in the ERA5 reanalysis, the region with consistently strong shears can be noticed at approximately 16 km a.g.l. (Fig. 3j and l), which is about 3 km lower than that in the HVRRS. One possible reason might be that the model fails to resolve the further increasing shear in the lower stratosphere. In the tropics, the signature of quasi-biennial oscillation (QBO) can be identified in the lower stratosphere (Fig. 3d and k).

The comparison between HVRRS-based and ERA5-based shears in three typical regimes is tabulated in Table 1. These metrics highlight that ERA5-based shears are underestimated by approximately 3.92, 7.65, and 11.99 $\mathrm{m}{\mathrm{s}}^{-\mathrm{1}}{\mathrm{km}}^{-\mathrm{1}}$ at heights of 0–2, 10–15, and 20–25 km a.g.l., respectively, which are roughly consistent with Houchi et al. (2010).

By comparison, the ERA5-acquired N² averaged over four height intervals (e.g., 0–5, 5–10, 10–15, 15–20 km a.g.l.) is reliably estimated, with a relative error of around 11 %, as illustrated in Fig. S3 in the Supplement. This finding indicates that the ERA5 reanalysis can properly present the static stability of the background atmosphere, but it is not properly coincident with radiosonde in terms of the small-scale variability in dynamical structures. Due to a lack of global measurement of the fine structure of the upper-air wind, however, the accuracy of ERA5-resolved shears is hard to globally validate.

3.2 Occurrence frequency of $\mathit{\text{Ri}}<\mathrm{1}/\mathrm{4}$ in HVRRS and ERA5 reanalysis

As a prominent example, the monthly occurrence frequency of $\mathit{\text{Ri}}<\mathrm{1}/\mathrm{4}$ over the Corpus Christi station (27.77^∘ N, 97.5^∘ W) from January 2017 to October 2022 is illustrated in Fig. 4. As a result, the monthly occurrence rate of $\mathit{\text{Ri}}<\mathrm{1}/\mathrm{4}$ in the low troposphere determined from the HVRRS is lower than the ERA5-based one, with mean values of around 10.6 % and 16.9 %, respectively. In the lowermost 2 km, the vertical resolution of ERA5 reanalysis is less than 200 m, and it is less than the moving segment interval in Eq. (1). The high occurrence frequency in the low troposphere could be likely related to the negative or small N². Especially during the daytime, the planetary boundary layer (PBL) is well mixed due to strong turbulence induced by uprising thermals (Song et al., 2018). In addition, an obvious seasonal cycle of occurrence frequencies is revealed by the HVRRS in the middle and upper troposphere and has a maximum in the winter (December–January–February) and spring (March–April–May) seasons, which is consistent with the finding in Zhang et al. (2019). In the vicinity of jet streams, the occurrence frequency of $\mathit{\text{Ri}}<\mathrm{1}/\mathrm{4}$ is generally enhanced by large wind shears. However, the ERA5 reanalysis does not provide such a seasonal cycle pattern, and the occurrence frequency of $\mathit{\text{Ri}}<\mathrm{1}/\mathrm{4}$ is significantly underestimated by around 8 % (Fig. 4b), which could be attributed to the underestimation in wind shears. In the lower stratosphere, both the HVRRS and the ERA5 reanalysis provide a low estimation of occurrence frequencies, with a value of around 1 %.

Figure 4 The monthly occurrence frequency of $\mathit{\text{Ri}}<\mathrm{1}/\mathrm{4}$ at the Corpus Christi station (27.77^∘ N, 97.5^∘ W) in HVRRS (a) and ERA5 reanalysis (b). Note that the contour curves in (a) and (b) concern the mean horizontal wind speed and that the ERA5-derived quantities are spatially and temporally collocated with those of the HVRRS.

Furthermore, on a large geographical scale the occurrence frequency of $\mathit{\text{Ri}}<\mathrm{1}/\mathrm{4}$ retrieved by ERA5 reanalysis is remarkably underestimated in the free atmosphere, as compared to the HVRRS. The annual variation in the occurrence frequency of $\mathit{\text{Ri}}<\mathrm{1}/\mathrm{4}$ over seven climate zones at 10–15 km a.g.l. indicated by HVRRS and ERA5 reanalysis is further demonstrated in Fig. 5. It is clearly seen that the occurrence frequency of $\mathit{\text{Ri}}<\mathrm{1}/\mathrm{4}$ provided by ERA5 reanalysis is underestimated in all months, over all climate zones, possibly implying that, in the free atmosphere, the threshold value of $\mathrm{1}/\mathrm{4}$ in Eq. (1) is too small for the ERA5 reanalysis to capture the occurrence of KHI.

Figure 5 The annual cycles of the occurrence frequency of $\mathit{\text{Ri}}<\mathrm{1}/\mathrm{4}$ in different climate zones at 10–15 km a.g.l. The red and blue boxes represent the frequencies in HVRRS and ERA5 reanalysis, respectively. The ERA5-derived Ri is spatially and temporally collocated with that of the HVRRS. NH, Northern Hemisphere; SH, Southern Hemisphere.

However, the ERA5 reanalysis data are non-uniformly sampled in altitude. Their vertical resolution drops from about 100 m in the boundary layer to about 500 m in the lower stratosphere. In contrast, radiosondes have a vertical resolution of 10 m at all heights. Therefore, we selected four typical heights and vertically interpolated the radiosonde to the same height resolution as ERA5 for comparison. The four height intervals are 0.8–1.3, 2.2–3.2, 6–15, and 20–21 km a.g.l., as shown in Table 2. In these height intervals, the vertical resolution of ERA5 is about 100, 200, 300, and 400 m, respectively. Even at the same vertical resolution, ERA5 still seriously underestimates the value of OF($\mathit{\text{Ri}}<\mathrm{1}/\mathrm{4}$) (OF denotes occurrence frequency) at all heights and in all climate zones. These results indicate that the greatest difficulty in evaluating subcritical Ri with ERA5 is that its simulation of wind shears might be seriously underestimated compared with radiosonde. As illustrated in Table 3, even accounting for the fact that ERA5 has a comparable vertical resolution to the radiosonde, wind shears in ERA5 reanalysis are still underestimated by around 50.3 %, 48.7 %, 43.6 %, and 62.2 % at 0.8–1.3, 2.2–3.2, 6–15, and 20–21 km a.g.l., respectively. In order to obtain an occurrence frequency of subcritical Ri from ERA5 reanalysis that is comparable with radiosonde-based OF($\mathit{\text{Ri}}<\mathrm{1}/\mathrm{4}$), the Ri_t for ERA5 should be set to larger than $\mathrm{1}/\mathrm{4}$. For instance, at 0.8–1.3 and 2.2–3.2 km a.g.l., Ri_t equals 1 could be an appropriate choice for ERA5 reanalysis, rather than $\mathrm{1}/\mathrm{4}$ (Table 2). More generally, $\mathrm{0.5}<{\mathit{\text{Ri}}}_{\mathrm{t}}<\mathrm{1.5}$ could be more suitable for ERA5 reanalysis, compared to ${\mathit{\text{Ri}}}_{\mathrm{t}}=\mathrm{1}/\mathrm{4}$.

Table 2 The occurrence rate of low Ri at 0.8–1.3 (a), 2.2–3.2 (b), 6–15 (c), and 20–21 km a.g.l. (d). The critical Ri (Ri_t) is $\mathrm{1}/\mathrm{4}$ for radiosonde, and it increases from $\mathrm{1}/\mathrm{4}$ to 2 for ERA5 reanalysis. Note that HVRRS data were vertically resampled to 100, 200, 300, and 400 m at these four height intervals to match the ERA5 reanalysis. In addition, the moving-average number in Eq. (1) is 0. RS stands for radiosonde.

Table 3 Vertical wind shears at 0.8–1.3 (a), 2.2–3.2 (b), 6–15 (c), and 20–21 km a.g.l. (d). Note that HVRRS data were vertically resampled to 100, 200, 300, and 400 m at these four height intervals to match the ERA5 reanalysis. RS stands for radiosonde.

Figure 6 The altitude variation in the occurrence frequency of Ri below certain thresholds (0.25, 0.5, 1, 1.5, and 2) in ERA5 reanalysis in various climate zones. The ERA5-derived Ri is spatially and temporally collocated with that of the HVRRS. The occurrences of $\mathit{\text{Ri}}<\mathrm{1}/\mathrm{4}$ in the HVRRS are overlapped with red lines.

Due to the huge change in the vertical resolution of ERA5, it could be difficult to interpolate ERA5 into uniform data vertically with a relatively high resolution. Therefore, the question posed here is, what is the proper threshold value of Ri in predicting the occurrence of KHI when using the ERA5 reanalysis, compared to the HVRRS? The occurrence frequency of $\mathit{\text{Ri}}<\mathrm{1}/\mathrm{4}$ indicated by the HVRRS and the ERA5-determined occurrence frequencies produced by Ri<0.25, Ri<0.5, Ri<1, Ri<1.5, and Ri<2 at all heights up to 30 km a.g.l. are demonstrated in Fig. 6. It is notable that over all climate zones and in the free atmosphere, occurrence frequencies of Ri<0.25 and Ri<0.5 obtained from the ERA5 reanalysis are underestimated, but the frequencies of Ri<1.5 and Ri<2 are generally overestimated. The occurrence frequency of Ri<1 gives a close estimation for both magnitude and altitude variations compared to the HVRRS over all climate zones.

Figure 7 The correlation coefficients between the monthly averaged occurrence frequency of $\mathit{\text{Ri}}<\mathrm{1}/\mathrm{4}$ in the HVRRS and the monthly occurrence frequency of Ri below certain thresholds (0.25, 0.5, 1, 1.5, and 2) in ERA5 reanalysis. The ERA5-derived Ri is spatially and temporally collocated with that of the HVRRS. The coefficients in various climate zones are estimated in increments of 2 km.

Furthermore, the correlation coefficients between HVRRS-determined occurrence frequencies and the ERA5-determined frequencies indicated by different threshold values of Ri at height levels of 0 to 30 km are illustrated in Fig. 7. It is worth noting that, in the troposphere, the ERA5-based frequencies indicated by Ri<1, Ri<1.5, and Ri<2 are highly positively correlated with those from the HVRRS, with a correlation coefficient of around 0.6 over all climate zones. In the lower stratosphere, however, these coefficients rapidly decline to 0.1, which can be explained by the low occurrence frequency in this height regime.

Combining the findings in Figs. 6 and 7, in the free troposphere, we can conclude that the ERA5-determined occurrence frequency of Ri<1 is closest to the frequency of $\mathit{\text{Ri}}<\mathrm{1}/\mathrm{4}$ based on the HVRRS. In the free atmosphere, KHI is the dominant source for clear-air turbulence (CAT), which is a well-known hazard to aviation. Therefore, the global characterization of KHI occurrence frequency in the free atmosphere obtained from ERA5 reanalysis could be of importance for understanding the spatial–temporal variation in CAT. In the following sections, the occurrence frequency of subcritical Ri (hereinafter OF(Ri<Ri_t)) is based on Ri<1 in ERA5 reanalysis and $\mathit{\text{Ri}}<\mathrm{1}/\mathrm{4}$ in the HVRRS, unless otherwise noted.

Figure 8 The percentage of OF(Ri<0) relative to OF(Ri<Ri_t) in the HVRRS (red) and ERA5 reanalysis (blue).

Finally, it is noteworthy that OF(Ri<Ri_t) includes the component of Ri<0 that indicates potential for convective instability. However, both ERA5 and the HVRRS are difficult to totally avoid Ri<0 when calculating Ri. Therefore, we evaluated the proportion of Ri<0 of all Ri<Ri_t values across two datasets to evaluate the possible contribution from convections, as shown in Fig. 8. For the HVRRS, the proportion of OF(Ri<0) drops sharply from about 40 % in the low troposphere to about 18 % at 5–15 km a.g.l. Similarly, for ERA5, its proportion drops from about 40 % in the lowermost part of the atmosphere to about 2 % at 5–16 km a.g.l. These findings indicate that, in the free atmosphere, OF(Ri<Ri_t) is mainly composed of OF($\mathrm{0}<\mathit{\text{Ri}}<{\mathit{\text{Ri}}}_{\mathrm{t}}$).

3.3 The OF(Ri<Ri_t) climatology

For a first indication, the global distributions of OF(Ri<Ri_t) provided by the ERA5 reanalysis at 0–2, 5–10, 10–15, and 15–20 km a.g.l. are displayed in Fig. 9. OF(Ri<Ri_t) in the low troposphere is considerably spatially heterogeneous. In subtropical oceans, the intense OF(Ri<Ri_t) can be noticed and has a magnitude of around 50 % (Fig. 9a). In addition, in the Sahara OF(Ri<Ri_t) reaches as high as 65 %. Interestingly, the spatial variation in mean OF(Ri<Ri_t) during the years 2017–2022 is highly consistent with that of the planetary boundary layer height (PBLH) over oceans, such as the Pacific Ocean near Japan and the Atlantic Ocean near the USA, as shown in Fig. S4 in the Supplement. However, at 0–2 km a.g.l., the spatial variation in OF($\mathrm{0}<\mathit{\text{Ri}}<{\mathit{\text{Ri}}}_{\mathrm{t}}$) exhibits a large difference from that of OF(Ri<Ri_t) in terms of magnitude, as shown in Fig. S5 in the Supplement. It is around 40 % (20 %) lower than that of OF(Ri<Ri_t) in subtropical oceans (Australia and North Africa). At heights of 5–10 km a.g.l., intensive OF(Ri<Ri_t) can be viewed in the subtropical regions and has a value of around 10 % (Fig. 9b), which is likely attributed to upper-tropospheric jets. In the upper troposphere above the tropical oceans, OF(Ri<Ri_t) is as high as 30 % (Fig. 9c), possibly as a result of the maximal heating effect by mesoscale convective systems (e.g., Houze, 1982). In the lower stratosphere, OF(Ri<Ri_t) sharply decreases to around 0.1 % (Fig. 9d).

Figure 9 The spatial distribution of the mean OF(Ri<Ri_t) in ERA5 reanalysis at 0–2 (a), 5–10 (b), 10–15 (c), and 20–25 km a.g.l. (d). Note that Ri_t is set to 1.

Figure 10 The monthly averaged OF(Ri<Ri_t) in the HVRRS (a–g) and ERA5 reanalysis (h–n) in seven climate zones. NH: Northern Hemisphere; SH: Southern Hemisphere.

In comparison, the spatial–temporal variability in free-tropospheric OF(Ri<Ri_t) indicated by the HVRRS is fairly consistent with that of ERA5 reanalysis over all climate zones (Fig. 10). Seasonal cycles can be detected by both the HVRRS and the ERA5 reanalysis over all climate zones, especially in the subtropics and midlatitudes. However, the ERA5-based OF(Ri<Ri_t) can only reflect the large-scale structure of the cycles, and it is hard to quantify the detailed variation with ERA5 like the HVRRS does.

Figure 11 The seasonally averaged OF(Ri<Ri_t) in the HVRRS (a–g) and ERA5 reanalysis (h–m) in seven climate zones. MAM, March–April–May; JJA, June–July–August; SON, September–October–November; DJF, December–January–February. NH: Northern Hemisphere; SH: Southern Hemisphere.

Furthermore, the seasonal variation in OF(Ri<Ri_t) with ${\mathit{\text{Ri}}}_{\mathrm{t}}=\mathrm{1}/\mathrm{4}$ for the HVRRS and Ri_t=1 for ERA5 in all climate zones is provided in Fig. 11. In the midlatitudes and subtropics, OF(Ri<Ri_t) exhibits maximum values in the low troposphere, as well as a local minimum in the middle troposphere and a local maximum at altitudes around 9 km. In the lower stratosphere, the occurrence frequencies decrease to values of the order of 1 % (Fig. 11b, c, e, and f). In the tropics, a primary peak can be clearly noticed at around 13 km, with a maximum of 12 % for the HVRRS and 20 % for the ERA5 reanalysis (Fig. 11d and k). The seasonality in the tropical region could be related to some large-scale flow features like the summer Asian monsoon and the tropical easterly jet (Roja Raman et al., 2009; Sunilkumar et al., 2015; Kaluza et al., 2021). In polar regions, the tropospheric OF(Ri<Ri_t) is significantly lower than that of other climate zones, with values ranging from around 4 % at heights of 2–8 km to 1 % in the lower stratosphere (Fig. 11a and g).

In Table 4, the mean OF(Ri<Ri_t) magnitudes over seven climate zones and in three typical altitude regimes are listed. At 0–2 km a.g.l., the ERA5-based OF(Ri<Ri_t) is about 24 % larger than that of the HVRRS-based one. At 10–15 km a.g.l., the ERA5-based OF(Ri<Ri_t) is reasonably well estimated, except that it is overestimated by around 5.92 % in the tropics. In addition, ERA5 underestimates OF(Ri<Ri_t) by around 0.5 % in the lower stratosphere.

Table 4 Similar to Table 1 but for the occurrence frequency of Ri<Ri_t. Note that Ri_t is indicated by $\mathit{\text{Ri}}<\mathrm{1}/\mathrm{4}$ in radiosonde, but it is identified with 1 in ERA5 reanalysis.

Figure 12 The association of HVRRS-determined OF(Ri<Ri_t) with different standard deviations of orography (dimensionless). Panels (a), (b), and (c) are for height ranges of 1–2, 3–4, and 5–6 km a.g.l., respectively. The correlation coefficients between OF(Ri<Ri_t) and standard derivation of orography are marked in the top-right corner, where the star superscripts indicate that values are statistically significant (p<0.05).

According to Fig. 9a, it seems that low-level continental OF(Ri<Ri_t) is dependent on underlying terrains. We investigate the association of low-level HVRRS-determined OF(Ri<Ri_t) with the standard deviation of orography (SDOR). At heights of 1–2 km a.g.l., the underlying terrain with a large SDOR generally corresponds to a high OF(Ri<Ri_t), with a correlation coefficient between OF(Ri<Ri_t) and SDOR of 0.24. Then, the coefficient decreases to 0.15 at 3–4 km a.g.l. (Fig. 12b), and eventually, it equals 0.14 at 5–6 km a.g.l. (Fig. 12c). These findings indicate that complex terrain may locally enhance OF(Ri<Ri_t).

Moreover, it is quite evident from Figs. 9b and S5 that both OF(Ri<Ri_t) and OF($\mathrm{0}<\mathit{\text{Ri}}<{\mathit{\text{Ri}}}_{\mathrm{t}}$) are largely enhanced above the tropical ocean associated with the El Niño–Southern Oscillation (ENSO). Most of the enhanced OF(Ri<Ri_t) can be identified over the Niño-3 region (5^∘ N–5^∘ S, 150–90^∘ W), and the time–height cross section of OF(Ri<Ri_t) during the years 2000–2022 is illustrated in Fig. 13. OF(Ri<Ri_t) values in the height region of 6–13 km are evidently large, at around 35 %, which is about 15 % larger than the climatological mean value (Fig. 10k). More specifically, OF(Ri<Ri_t) during time periods of La Niña events is obviously stronger than that during the El Niño periods. The identification of ENSO events is based on Ren et al. (2018), Li et al. (2022), and Lv et al. (2022). It is also worth recalling here that the wind shear does not exhibit such an anomaly over the Niño-3 region (Fig. 2c), implying that the OF(Ri<Ri_t) anomaly could likely be attributed to the ENSO-related tropical convective heating in the upper troposphere, leading to a low Brunt–Väisälä frequency.

Figure 13 The monthly averaged OF(Ri<Ri_t) in ERA5 reanalysis over the Niño-3 region (5^∘ N–5^∘ S, 150–90^∘ W). The blue and red shading on the time axis indicates the time periods with El Niño and La Niña events, respectively.

3.4 The dynamical environment of OF(Ri<Ri_t) in the free troposphere

In the free troposphere the percentage of OF(Ri<0) relative to OF(Ri<Ri_t) is generally less than 20 % (Fig. 8); KHI is preferentially generated from strong wind shears, which may be closely associated with mean flows and wave activities.

Figure 14 Geographical distribution of mean tropospheric GW total energy obtained from the HVRRS (a). The latitudinal variation in mean energy in a grid cell of 5^∘ latitude (b). The joint distribution of OF(Ri<Ri_t) with GW energy and wind shear (c). OF(Ri<Ri_t) and wind shear are derived from individual HVRRS profiles and vertically averaged over the tropospheric segment that is used for GW study. The numerical number in (c) indicates the matched profile number in each grid, using a bin size of 2 J kg⁻¹ along the x axis and 2 $\mathrm{m}{\mathrm{s}}^{-\mathrm{1}}{\mathrm{km}}^{-\mathrm{1}}$ along the y axis.

The propagation of GW could raise strong wind shear and therefore generate KHI. Therefore, we investigate the joint distribution of OF(Ri<Ri_t) with tropospheric GW total energy and wind shear (Fig. 14). The latitudinal variation in GW total energy exhibits a double-peak structure, with two peaks at around 30^∘ in the Northern Hemisphere and Southern Hemisphere (Fig. 14a). The joint distribution of OF(Ri<Ri_t) with GW energy and wind shear indicates that large OF(Ri<Ri_t) (for instance, >10 %) generally corresponds to GW energy larger than 10 J kg⁻¹ or to wind shear exceeding 14 $\mathrm{m}{\mathrm{s}}^{-\mathrm{1}}{\mathrm{km}}^{-\mathrm{1}}$ (Fig. 14b). Also, OF($\mathrm{0}<\mathit{\text{Ri}}<{\mathit{\text{Ri}}}_{\mathrm{t}}$) exhibits a similar distribution (Fig. S6 in the Supplement). Overall, OF(Ri<Ri_t) obviously increases with GW total energy (Fig. S9a in the Supplement), possibly implying that the propagation of GWs could enhance wind shears and, therefore, the burst of KHI.

In addition, the interaction between low-level winds and mountain barriers could be a source of orographic GWs (Zhang et al., 2022a). We take orographic GW dissipation in ERA5 reanalysis, which is the accumulated conversion of kinetic energy in the mean flow into heat over the whole atmospheric column, as an indicator of the strength of orographic GWs. It is interesting to note that monthly averaged orographic GW dissipation and monthly ERA5-determined OF(Ri<Ri_t) at heights from near the ground up to 30 km demonstrate a close association (Fig. S7 in the Supplement). For instance, in the middle troposphere, they are positively associated over mountainous areas such as the Rocky Mountains and the European Alps, with correlation coefficients of around 0.5. These findings also suggest that during months with strong unresolved orographic gravity wave activity, the flow and stability parameters of the resolved flow are then modified, leading to a low Richardson number. Nevertheless, it is hard to quantify the effect of resolved orographic GWs on Ri here.

Figure 15 Joint distribution of HVRRS-derived wind speed, wind shear, and OF(Ri<Ri_t), with a bin size of 5 m s⁻¹ along the x axis and 5 $\mathrm{m}{\mathrm{s}}^{-\mathrm{1}}{\mathrm{km}}^{-\mathrm{1}}$ along the y axis. Note that all the relationships are based on the mean result of individual profiles at heights of 10–15 km a.g.l. Numbers indicate the matched profile number in each grid.

At jet heights (10–15 km a.g.l.), large shears can be easily induced by strong wind speeds. Figure 15 demonstrates the joint distribution of OF(Ri<Ri_t) with wind speed and wind shear. Generally, OF(Ri<Ri_t) larger than 10 % can be easily found when the wind shear exceeds 20 $\mathrm{m}{\mathrm{s}}^{-\mathrm{1}}{\mathrm{km}}^{-\mathrm{1}}$. In addition, OF($\mathrm{0}<\mathit{\text{Ri}}<{\mathit{\text{Ri}}}_{\mathrm{t}}$) gives a similar conclusion (Fig. S8 in the Supplement). In the middle and upper troposphere, OF(Ri<Ri_t) almost linearly increases with wind speed (Fig. S9b).

In a short conclusion, in the free troposphere, the occurrence of KHI would favor the dynamical environment with intensive orographic or non-orographic GW activities and large mean flows.

4 Conclusion and remarks

The occurrence of KHI is potentially crucial in many ways, such as for aircraft safety and mass transfer, but it is very hard to be globally understood due to its fine structure. The subcritical Richardson number is commonly used as an indictor for KHI. This study uses ERA5 as the latest reanalysis product from ECMWF as well as a comprehensive dataset of high-resolution radiosonde to globally characterize the distribution of low Richardson numbers as a proxy for the occurrence of KHI, for the years 2017 to 2022.

Vertical wind shears are considerably underestimated at almost all heights and over all climate zones by the ERA5 reanalysis compared to the HVRRS. It is noteworthy that vertical wind shear in the ERA5 reanalysis at heights of 10–15 km a.g.l. is dramatically underestimated by around 7.65 $\mathrm{m}{\mathrm{s}}^{-\mathrm{1}}{\mathrm{km}}^{-\mathrm{1}}$, especially at middle latitudes. However, the spatial distribution of the ERA5 shear exhibits a statistically significant positive correlation with the HVRRS shear. As a result, the ERA5-determined occurrence frequency of $\mathit{\text{Ri}}<\mathrm{1}/\mathrm{4}$ is significantly underestimated. In addition, it is weakly correlated with HVRRS-determined ones at most heights and over most climate zones.

However, the vertical resolution of ERA5 reanalysis sharply decreases with altitude, which is not comparable with the HVRRS. Thus, to match ERA5 reanalysis at specified height intervals, the HVRRS was vertically interpolated with resolutions spanning 100 to 400 m. Even at a comparable resolution, ERA5-derived shear is underestimated by around 50 %, leading to a considerable underestimation of OF($\mathit{\text{Ri}}<\mathrm{1}/\mathrm{4}$) compared to that derived from radiosondes.

Interestingly, the ERA5-determined occurrence frequency of Ri<1 is generally consistent with the frequency of $\mathit{\text{Ri}}<\mathrm{1}/\mathrm{4}$ obtained from the HVRRS, in terms of magnitude and temporal variation. Rather than $\mathit{\text{Ri}}<\mathrm{1}/\mathrm{4}$, we argue that the threshold value of Ri<1 could be more appropriate when using ERA5 reanalysis for KHI study, especially in the middle and upper troposphere over midlatitude and subtropical regions in the Northern Hemisphere and Southern Hemisphere, where a high consistency between the HVRRS and ERA5 has been found in terms of OF(Ri<Ri_t) magnitude. In other words, for a similar occurrence frequency, the identification of vertical segments with Ri<1 in ERA5 is equitable with the identification of vertical segments with $\mathit{\text{Ri}}<\mathrm{1}/\mathrm{4}$ using the HVRRS. It is worth highlighting that the HVRRS experiences a 200 m vertical moving-average procedure to inhabit measurement noises and turbulence fluctuations. Without this procedure, the threshold Ri for the ERA5 reanalysis would be even higher than 1.

The climatology of OF(Ri<Ri_t) exhibits significant seasonal cycles over all latitudes. A poleward decrease can be clearly identified in the middle and upper troposphere. In addition, over mountainous areas, complex terrains may locally enhance low-level OF(Ri<Ri_t). Moreover, it is immediately obvious that both OF(Ri<Ri_t) and OF($\mathrm{0}<\mathit{\text{Ri}}<{\mathit{\text{Ri}}}_{\mathrm{t}}$) in the middle and upper troposphere of the Niño-3 region are considerably enhanced, probably by the tropical convective heating.

Moreover, both OF(Ri<Ri_t) and OF($\mathrm{0}<\mathit{\text{Ri}}<{\mathit{\text{Ri}}}_{\mathrm{t}}$) exhibit a close relationship with GW activities and background mean flows. Large OF(Ri<Ri_t) favors intensive GW activities and strong mean flow. Over complex terrains, the orographic GW breaking could locally enhance OF(Ri<Ri_t).

Those findings are valuable for pointing out the performance of the ERA5 reanalysis in terms of resolving low Richardson numbers as a proxy for KHI in comparison with a near-global high-resolution radiosonde measurement. In addition, the spatial–temporal variability in OF(Ri<Ri_t) over different climate zones from near the ground up to 30 km is quantitatively characterized by ERA5 and the HVRRS, which could provide new insights that increase our understanding of the fine structure of upper air.