Changes in Global Aviation Turbulence in the Remote Sensing Era (1979–2018)

: Atmospheric turbulence primarily originates from abrupt density variations in a vertically stratified atmosphere. Based on the prognostic equation of turbulent kinetic energy (TKE), we here chose three indicators corresponding to the forcing terms of the TKE generation. By utilizing ERA5 reanalysis data, we investigate first the maximum achievable daily thickness of the planetary boundary layer (PBL). The gradient Richardson number ( R i ) is used to represent turbulence arising from shear instability and the daily maximum convective available potential energy ( CAPE ) is examined to understand the turbulence linked with convective instability. Our analysis encompasses global turbulence trends. As a case study, we focus on the North Atlantic Corridor (NAC) to reveal notable insights. Specifically, the mean annual number of hours featuring shear instability ( R i < 0.25) surged by more than 300 h in consecutive 20-year periods: 1979–1998 and 1999–2018. Moreover, a substantial subset within the NAC region exhibited a notable rise of over 10% in the number of hours characterized as severe shear instability. Contrarily, turbulence associated with convective instability ( CAPE > 2 kJ/kg), which can necessitate rerouting and pose significant aviation safety challenges, displays a decline. Remote sensing of clouds confirms these assertions. This decline contains a component of inherent natural variability. Our findings suggest that, as air viscosity increases and hence a thickened PBL due to a warming climate, the global inflight turbulence is poised to intensify.


Introduction
Aircraft have evolved into an essential and highly efficient mode of global transportation for both people and goods (see Figure 1 in Ref. [1]).This rapidly expanding sector of the industry has consistently shown sensitivity to shifts in the global environment [1][2][3][4][5][6].Currently, the aviation industry confronts the dual challenges of ensuring the safety and comfort of passengers and crew, while simultaneously striving to enhance efficiency and minimize environmental impacts.These efforts are particularly important due to the projected post-pandemic doubling of aviation traffic by ~2030 [5,7,8] from the traffic levels of 2010.
A longstanding issue within aviation has been turbulence-induced safety concerns (WMO [9] aviation and environmental workshop on 12 March 2021; Refs.[2,10]).Recently, endeavors to mitigate climate change impacts on trans-Atlantic flights are in progress [11] as passenger safety is always the foremost priority.Studies of this critical topic have grown rapidly both in volume and scope (Ref.[12] and references therein).Environmental factors that optimize payload [6] and fuel efficiency [1,13] inherently bear implications for aviation safety.For commercial airlines, advancement in predictive wind-shear technology and enhanced pilot training programs have substantially curbed fatal accidents stemming from this cause.This study focuses on more precise quantification of atmospheric turbulence elements that hold the greatest relevance for enhancing passenger comfort and safety.
In addition to providing enhanced visibility, aircraft cruising near the tropopause have the benefits of encountering the coldest air and consequently least viscous air as the air viscosity is proportional to the temperature [1].This advantageous condition leads to increased fuel efficiency and reduced environmental footprints.Because of thermal wind balance [14], the temperature gradient generates westerly winds that amplifies with altitude.Vertically, the wind speeds maximize near the tropopause, a region where the sign of the meridional temperature gradient reverses (Figure 1a).While the westerly atmospheric jet streams furnish ideal tailwind conditions, they concurrently represent a notable source of turbulence [12] due to shear and inertial instability [15].Under the influence of jet stream, vertical wind shear alters surfaces of constant angular momentum, potentially leading to the occurrence of conditional symmetric instability.Here, we focus on the explanation of changes in aviation turbulence, assessed through the classical Gradient Richardson number (R i ) to quantify shear instability, boundary layer thickness (H) for near surface turbulence, and the convective available potential energy (CAPE) to gauge convective instability.These three indices correspond respectively to the turbulence kinetic energy (TKE) generation terms [16].Aside from the static topography caused gravity wave break and wake vortices, the three covered the full physical ranges of aviation turbulence generation.Various operational aviation turbulence indices (e.g., Ref. [10]) are based on and are various combination of the three.However, without explicitly being guided by the TKE generation equation, previous studies are heavily biased on clear air turbulence (i.e., the R i based shear instability [17,18]).True energy sources for the increase in aviation turbulence under a warming climate thus needs further elucidating.This study intends to complement the previous ones by providing a physically complete picture of aviation turbulence and the sensitivity to warming.The net effect is equivalent to a superimposed circulation that counters the present climate mean thermal winds at the lower troposphere (1000-600 hPa) and enhances the present climate mean thermal winds at the upper troposphere (600-200 hPa).The commercial aviation operating domain is hatched in (b).The illustrations were produced using GrADS (available at http://cola.gmu.edu/grad.php;accessed on 16 June 2018).

Data (Sources and Processing) and Methods
The datasets utilized in the prior work by Ren et al. (2020) [1] have also been incorporated in this study.The following is a brief description of the finest resolution available, quality controlled hourly meteorological reanalysis provided by the European Centre for Medium-Range Weather Forecasts (version 5; ERA5 henceforth [21]).

Quality-Assured Hourly ERA5 Reanalysis
Airline flights typically span durations of 1 h for domestic routes and up to 17 h for international journeys.Given this temporal range, stability indices are analyzed at hourly or sub-hourly intervals and this analysis relies on the high-quality ERA5 data during the remote sensing era available from 1979 onwards.The needs for finer temporal resolution The net effect is equivalent to a superimposed circulation that counters the present climate mean thermal winds at the lower troposphere (1000-600 hPa) and enhances the present climate mean thermal winds at the upper troposphere (600-200 hPa).The commercial aviation operating domain is hatched in (b).The illustrations were produced using GrADS (available at http://cola.gmu.edu/grad.php;accessed on 16 June 2018).
Although these selected indicators may not comprehensively cover the entire spectrum of aviation turbulence, they establish connections between three aspects of aviation turbulence and the corresponding meteorological phenomena.Considering that some previous studies [19] used 20 indicators for clear-air turbulence alone, our approach is apparently succinct and avoided the overlap of indicators.In fact, previously used empirical indicators are intercorrelated and tend to distort the true source/origin of aviation turbulence.To assess trends in various forms of aviation turbulence, reanalysis data of the atmospheric parameters over a 40-year period (1979-2018) are utilized.Justification of the use of these indicators derived from low spatial-temporal resolution datasets is given in Section 2, supported by recent developments in turbulence research [12,20].Section 3 shows the result and Section 4 provides the conclusion.

Data (Sources and Processing) and Methods
The datasets utilized in the prior work by Ren et al. (2020) [1] have also been incorporated in this study.The following is a brief description of the finest resolution available, quality controlled hourly meteorological reanalysis provided by the European Centre for Medium-Range Weather Forecasts (version 5; ERA5 henceforth [21]).

Quality-Assured Hourly ERA5 Reanalysis
Airline flights typically span durations of 1 h for domestic routes and up to 17 h for international journeys.Given this temporal range, stability indices are analyzed at hourly or sub-hourly intervals and this analysis relies on the high-quality ERA5 data during the remote sensing era available from 1979 onwards.The needs for finer temporal resolution data arises from the dynamic and more variable nature of the turbulence, which surpass the variability of the convective systems or jet streams that give rise to the turbulence.As a result, even flights scheduled within hours of each other can encounter highly distinct meteorological conditions.
Hourly ERA5 data are available on a 0.25 degrees lat/lon (~30 km) horizontal grid with 137 vertical levels spanning from the surface to ~80 km above the sea level.The stability indices discussed below are generated using standard pressure level data which is interpolated from the original model levels.Single level products of the planetary boundary layer (PBL) depth H and CAPE also are examined.
Using reanalysis data to estimate aviation turbulence indicators is justified by Jaeger and Sprenger (2007) [22], who posit that mechanisms behind turbulence generation originate from large scale forcing (resolvable by ERA5), thereby contributing to atmospheric instability [23][24][25].Aviation turbulence indicators are typically defined as vertical derivatives of atmospheric parameters.With the fine vertical resolution provided by ERA5, critical values of turbulence indicators (e.g., 0.25 for critical R i ) can refer to those from archived PBL experiments.It is important to note that the scale of aviation turbulence is ~1000 m and winds at a horizontal resolution about 30 km may not fully capture the extreme values encountered by the aircraft.Hence, the critical values of the turbulence indicators should be regarded as qualitative, serving as conservative estimates.
In addition to ERA5, we also examine four other reanalysis products to minimize uncertainty in our analyses: (i) the European Centre for Medium-Range Weather Forecasts interim reanalysis (ERA-Interim reanalysis [26]), (ii) the NCEP/NCAR reanalysis [27]; (iii) the NASA's Modern Era Retrospective analysis for Research and Applications (MERRA [28]); and (iv) the JRA-55 reanalysis [29], at horizontal resolutions of 0.25 • , 0.75 • , 2.5 • , 0.625 • × 0.5 • and 1.25 • respectively.Their vertical resolution near cruising altitudes (~200 hPa) varies from 25 hPa (for ERA5, ERA-Interim and JRA-55) to 50 hPa (for NCEP/NCAR and MERRA).In locations where regular atmospheric soundings are available (https://weather.uwyo.edu/upperair/sounding.html; accessed on 15 March 2022), CAPE and PBL heights show good agreement between sounding analyses and the reanalyses.The root mean squared errors in PBL and CAPE in May are respectively 55 m and 4 J/kg for Norman, OK.It is routine in NWP models to calculate CAPE from atmospheric sounding profiles.Here we used plotskew.gsfrom https://moe.met.fsu.edu/~rhart/plotskew.gs (accessed on 2 February 2017).The method used for calculating PBL depth is outlined in the Appendix A (Figure A1).

Pathfinder Atmospheric Extended (PATMOS-x) Cloud Properties Products
By assimilating remotely sensed atmospheric parameters, especially over oceanicdominate southern hemisphere, the quality of reanalysis data after 1979 are reliable.Additionally, aviation turbulence has salient remote sensing fingerprints, in wind fields as well as cloud characteristics [30].In this study, cloud products already processed using the PAT-MOS algorithm developed by Heidinger [31] of the National Environmental Satellite, Data, and Information Service (NESDIS) group located in the Space Science and Engineering Center (SSEC) at the University of Wisconsin.
The PATMOS-x cloud properties climate data record (CDR) provides data for multiple cloud properties and brightness.In addition to being open to public, the dataset has advantage that it uses a consistent set of spectral bands over the whole period from 1979 to present.Further, new satellites when launched are overlapped in time with existing on-orbit satellites so that there is a cross-satellite radiometric calibration of the whole time series.Partly due to this unique quality control feature, the PATMOS series of products has been anointed by GEWEX as a high-quality global data set.In this study, unique convective anvil clouds (e.g., CAPE indicated aviation turbulence) as well as the cirrus clouds associated with polar jet streams (e.g., Ri indicated CAT) are examined to estimate trends in their occurrence frequency and areal extent.Both cirrus (Ci) and cumulonimbus (Cb) have low cloud top temperatures.Cloud thickness (i.e., the difference between cloud top and base pressures) and total vapor content are further considered in identifying cloud types.For instance, very dense Ci, that appears frequently in early spring season in the Eurasia continent, can be easily separated from other lower level stratiform clouds, which have similar visible signature for satellites.In case of multiple-layer clouds, environmental atmospheric conditions (from simultaneous ERA5 reanalyses) are also considered to avoid ambiguities.Due to limited spatial resolution, remote sensing cloud-type identification is not expected to identify the individual cloud cells but only overall optical properties of them.Automatic schemes [32] have evolved over the past thirty years and now become rather mature.For this study, a simplified joint histograms of cloud optical depth (τ) and cloud top pressure (p c ), as in Figure 2 in Ref. [33], suffices during the daytime (because the optical depth is for 0.65-micron visible channel).Since the IR ~1.375 micron channel is ideal for cirrus detection during nighttime, other datasets such as MODIS (Moderateresolution Imaging Spectroradiometer), AIRS (Atmospheric Infrared Sounder), CALIPSO (Cloud-Aerosol Lidar and Infrared Pathfinder Satellite Observations), and CloudSat were also examined for clear air turbulence.This strategy would be an important supplement to aviation turbulence analysis using reanalysis data, benefited by the advance in remote sensing technology.

Definition of Indicators of Aviation Turbulence
Atmospheric turbulence is a major aviation hazard [10], sometimes leading to aircraft damage and injuries among passengers and crew.The socio-economic implications of such occurrences were reviewed by Storer et al. [12].Their work affirms that wind shear (e.g., clear air turbulence, CAT), mountain lee waves (e.g., mountain wave turbulence, MWT), and convective instability (e.g., convectively induced turbulence, CIT, or in-cloud turbulence, ICT) collectively constitute the primary sources of aviation turbulence, either directly or indirectly through the influence of gravity waves [34,35].A more nuanced understanding of aviation turbulence may help stakeholders in aviation sector prepare for future incidences with more targeted maneuver codes and refining corresponding techniques to reduce the risk of injuries.
Atmospheric turbulence is produced at regions of intense winds fluctuation, usually has apparent remote sensing signature (e.g., unique cloud types and areal extends; concentrated flow patterns at specific altitudes).Although this study is diagnostic of air turbulence using atmospheric parameters, it would be beneficial to have a examine of the generic TKE equation: where PBL means the primary eddy vertical diffusion (i.e., TKE transport T) occurs in the planetary boundary layer-as energy provider for aviation turbulence, the SHS means horizontal shear eddies dominates the shear production term (S), the CAPE means convective induced turbulence (i.e., CIT and ICT) dominates the buoyancy production term (B), and the MWT represents topography caused wakes of vortices and lee-way, gravity wave breaking turbulence.Based on Mellor and Yamada's level 1.5 closure scheme, the eddy dissipation rate (ε) for aviation turbulence can be parameterized as proportional to TKE This study aims at understanding the changes of aviation turbulence as climate warms.Consequently, the MWT term, or the major contributor to ε, is not a focus of this study.Our philosophy is that a situation that favors TKE generation also favors aviation turbulence.In this sense, Equation (1) forms the basis of our indicator selection for aviation turbulence.This also provides a unified rubric for the twenty plus widely used empirical indicators [10,12,17].Among the three indices used here, the first two-the PBL depth (H) and the CAPE are primary outputs of the ERA5 hourly reanalysis.The third one, R i , is calculated using the hourly reanalysis wind and temperature at pressure levels (at increment of 25 hPa), serving as an indicator of turbulence linked to the large-scale gravity waves resolvable within the reanalysis data.Importantly, these three aviation turbulence indicators exhibit minimal overlap.The choice of these three specific indicators is rooted in the division of the atmospheric general circulation into two overarching classes, both originating from the uneven distribution of diabatic heating within the atmosphere.One class is driven by horizontal heating gradients within a stably stratified atmosphere and accounts for over 98% of the atmospheric kinetic energy (e.g., Ref. [14]; lecture notes of H. Bluestein 2002 on convective storms at the University of Oklahoma; and https://atmos.uw.edu/academics/ classes/2014Q1/545/545_Ch_1.pdfaccessed on 8 April 2024).Turbulence associated with the vertical shear of the horizontal wind field is typically linked to the R i .The other class is driven by convective instability and accounts for the remainder of the atmospheric kinetic energy.Convection is continuous process unfolding within distinct atmospheric regions due to the vertical gradient of diabatic heating.CAPE emerges as an appropriate descriptor of aviation turbulence associated with this type of air motion.The resultant motions span spatial scales ranging from ~30 km in the largest thunderstorms as minute as 1 mm in microscale movements within the surface layer.In this context, H becomes influential.
It is important to note that despite contributing relatively modestly to the total atmospheric kinetic energy, convectively driven motions play an important role in vertically transporting latent and sensible heat, hence linking to the daily maximum H. Furthermore, H holds relevance for aviation turbulence not only due to Earth's atmosphere obtains its energy primarily from its base (through interception of longwave radiation emitted from the Earth's surface), but also because, on the planetary scale, thermal gradients regulate wind fields, and H is a measure of the thermal condition of the lower atmosphere.
From fluid mechanics, R i is an indicator of Kelvin-Helmholtz Instability (KHI, e.g., clear air turbulence CAT in the vicinity of jets).The R i is readily calculable from ERA5 hourly data, using its standard definition: where → V is the horizontal wind vector, and z is the vertical coordinate.In Equation ( 2), the numerator, N 2 , is the Brunt-Väisälä frequency squared (N 2 ), which is a function of temperature, and the denominator is the square of the vertical wind shear, which is directly related to the horizontal temperature gradient.This approach to estimating atmospheric turbulence is widely used, including for aviation impact studies [12,36,37].Turbulence, originating from sources such as the PBL, mountain lee waves, or in the vicinity of jet streams, is generated by substantial density contrasts.In this context, R i serves as a broad-spectrum aviation turbulence indicator.The chosen critical value for R i , set to 0.25 (implying the needs for passengers to fasten seat belts), have been informed by discussions with X. Wang (2018).Meanwhile, CAPE is an indicator for Rayleigh-Taylor instability (e.g., buoyancy driven instability in the vicinity of storms).The selection of a threshold value for CAPE > 2 kJ/kg to denote severe convective turbulence (i.e., prompting flight rerouting) is to some extent arbitrary.But the critical value of 2 kJ/kg is also used in Ren et al. ( 2020) [1] in the investigation of turbulence effects on aviation fuel efficiency.
For a commercial aircraft cruising at ~0.85 M, the relevant aviation turbulence arises from eddies of sizes ranging between ~50-5000 m, corresponding to frequencies of 0.1-10 Hz (i.e., 10 3 -10 5 wave numbers for the planetary atmosphere).Both reanalysis data and climate model outputs lack the resolution necessary to directly capture the scales of these eddies.Traditionally, there is an absence of temporally coherent, high horizontal resolution measurements spanning a sufficiently extended timeframe for analytical purposes.Energy spectrum analyses drawn from along track measurements (e.g., GASP aircraft measurements [38]) reveal that even the canonical Kolmogorov −5/3 slope is no longer fully valid across the log-log energy spectrum diagram (kinetic energy density versus wave number).A slope of -3 extends comfortably into wave number regimes exceeding 10 3 , which corresponds to the aviation turbulence scale.As such, there exists strong selfresemblance across the eddy size spectrum [22].Recent statistical analysis of atmospheric turbulence using Thorpe analysis [20,39,40] align with this assertion.With the sufficient vertical resolution, a wide spectrum of turbulence can be adequately represented [34,35,41].
If the focus is not extensively on the turbulence encountered by an aircraft during its specific flight, but rather on the environmental conditions giving rise to aviation turbulence, stability indicators estimated from coarse resolution data remain adequate for analytical purpose.

Results and Discussion
Climate warming effects on aviation safety arise from the fluxes-coupling of the planetary boundary layer with the overlying free atmosphere (in communication with P. Williams, 2014).Upper atmospheric features (e.g., jet streams) are connected to the development of the PBL [42].PBL thicknesses have been collected over the past 40 years.From ERA5 hourly data, the maximum reachable PBL depths of each day are collected and averaged into monthly values (Figure 2).
Figure 3 elucidates alterations in March PBL depths between two consecutive 20-year periods 1979-1998 and 1999-2018.The most striking feature is the marked increase in PBL depth.Examination of PBL depth time series at selected locations reveals a prevailing trend of thickening PBL depths across the entire reanalysis period, punctuated by fluctuations stem from a blend of both nature and anthropogenic drivers.For example, Guo et al. [43], based on radiosonde data over China, found a spatially uniform increase in PBL depth from 1979-2003, followed by a spatially non-uniform decrease of PBL depth from 2003-2016.The decrease was ascribed to alterations in land-surface characteristics affecting the surface energy budget.
This overall increasing trend of PBL depth is explicable using fundamental physical principles.The PBL, being the layer most affected by friction, often has its depth parameterized as proportional to the square root of the eddy viscosity coefficient [44].This coefficient functions as a macroscopic representation of air's kinematic viscosity [45].This provides a pathway through which a warming climate could influence aviation turbulence, as discussed in the Appendix B. Multiple factors can govern PBL depth, including surface heating, surface energy partition, the free atmospheric lapse rate, among others.This overall increasing trend of PBL depth is explicable using fundamental physical principles.The PBL, being the layer most affected by friction, often has its depth parameterized as proportional to the square root of the eddy viscosity coefficient [44].This coefficient functions as a macroscopic representation of air's kinematic viscosity [45].This provides a pathway through which a warming climate could influence aviation turbulence, as discussed in the Appendix B. Multiple factors can govern PBL depth, including surface heating, surface energy partition, the free atmospheric lapse rate, among others.
Viscosity-related feedbacks persist as a discernible signal within a warming climate context.However, changes in ground surface heat fluxes stemming from changes in Bowen ratio [46] also acts as a regulatory factor.These internal (viscosity-related) and external (surface heat fluxes) factors, operating in synergy or opposition, collectively modulate the evolution of actual PBL depth with climate warming.PBL emerges as a notable source of diverse turbulence, which can pose challenges for aircraft flights.As shown in Figure 2, elevated plateaus can lead to the development of extraordinarily thick boundary layers.For example, the PBL depth over the Tibetan Plateau reaches ~6000 m.Considering that the elevated ground surface with a mean elevation ~4 km above sea level, the PBL reaches up to ~10 km above sea level.This proximity to the tropopause height facilitates an exchange between the stratosphere and the PBL, involving momentum as well as latent and sensible heat.This largely explains why commercial  (1979)(1980)(1981)(1982)(1983)(1984)(1985)(1986)(1987)(1988)(1989)(1990)(1991)(1992)(1993)(1994)(1995)(1996)(1997)(1998).Details of the calculations are described in Figure 2. Except for very limited areas in south Africa and Tibetan Plateau, all regions of increase PBL depth passed 95% confidence interval of a t-test of degree-of-freedom (dof) 40.Regions indicating a decrease of PBL depth (e.g., Western Australia) did not pass the same t-test.The map background and colour shades were produced using GrADS (available at http://cola.gmu.edu/grad.php;accessed on 1 October 2022).
Viscosity-related feedbacks persist as a discernible signal within a warming climate context.However, changes in ground surface heat fluxes stemming from changes in Bowen ratio [46] also acts as a regulatory factor.These internal (viscosity-related) and external (surface heat fluxes) factors, operating in synergy or opposition, collectively modulate the evolution of actual PBL depth with climate warming.
PBL emerges as a notable source of diverse turbulence, which can pose challenges for aircraft flights.As shown in Figure 2, elevated plateaus can lead to the development of extraordinarily thick boundary layers.For example, the PBL depth over the Tibetan Plateau reaches ~6000 m.Considering that the elevated ground surface with a mean elevation ~4 km above sea level, the PBL reaches up to ~10 km above sea level.This proximity to the tropopause height facilitates an exchange between the stratosphere and the PBL, involving momentum as well as latent and sensible heat.This largely explains why commercial flights do not typically traverse directly over the Tibetan Plateau.
While March serves as an example in this context, the same overall trend is exhibited in other months too, with winter seasons usually showcasing more pronounced increasing trends in PBL depths.This can also be seen across various months, including January (boreal winter), April (spring), July (summer), and October (fall), as depicted by the percentage changes in Figure 4. Particularly, the Siberian flight corridor has encountered the most substantial increase in PBL depth over the last 40 years.Concurrently, this region experienced the most salient warming.This alignment reinforces the connection between PBL depth and warming, with air viscosity serving as the intermediary.The consequences for turbulence are manifested in R i values and other pertinent indices throughout the air column.
As indicated by Ref. [42], atmospheric circulations extending up to 100 hPa within the tropical region (30 • S-30 • N) are sensitive to PBL thickness.The most intense turbulence can manifest either above or below aircraft cruising altitude (~200 hPa).However, it is the turbulence in proximity to the cruising altitude that is directly experienced by the aircraft.Hence, in the subsequent discussion, the indices are defined in relation to the cruising altitude unless explicitly stated otherwise.
To achieve optimal flight conditions including factors like enhanced visibility, strong tailwinds, and best thermal efficiency, commercial airliners often choose to cruise in the vicinity of the tropopause.The core of jet stream, the primary source of turbulence, is typically located just below the tropopause [47][48][49].Consequently, the tropopause altitude holds implications for virtually all aspects of civil aviation [1,6].Based on ERA5 data, Santer et al. [47] shows a continuous increase of tropopause height globally.
Disregarding minor changes in atmospheric composition, an elevated tropopause height is necessitated in response to an augmented PBL depth which is in line with the convective-radiative balance and influenced by planetary-scale features such as jet streams.The atmospheric thermal expansion resultant from global warming is manifested near the surface as a reinforced well-mixed boundary layer.To maintain energy balance, the tropopause must be elevated (in accordance with the inertial stability requirement [50].With exceptions in regions poleward of 60 • S over Southern Hemisphere oceans, which witness fewer flight routes (see Figure 1 in Ref. [1]), the increased tropopause height is apparent.Most notably, this increase is prominent over the mid-high latitudes and attains rates of up to 1 m/yr in certain aeras (e.g., over Europe and the North Pacific Ocean).
Apart from the relatively infrequent convective instability (Rayleigh-Taylor instability), the turbulence experienced by aircraft is primarily related to KHI.The onset of KHI instability is classically given by the value of the R i .Typically, a value of R i < 0.25 is required for layer instability.In subsequent sections, the focus first revolves around determining if there are changes in atmospheric turbulence activity using ERA5 reanalysis.Subsequently, investigation into global patterns of atmospheric turbulence relevant to aviation safety is undertaken.Further exploration is conducted using CAPE to quantify the turbulence associated with severe convective activity.To achieve optimal flight conditions including factors like enhanced visibility, strong tailwinds, and best thermal efficiency, commercial airliners often choose to cruise in the vicinity of the tropopause.The core of jet stream, the primary source of turbulence, is typically located just below the tropopause [47][48][49].Consequently, the tropopause altitude holds implications for virtually all aspects of civil aviation [1,6].Based on ERA5 data, Santer et al. [47] shows a continuous increase of tropopause height globally.
Disregarding minor changes in atmospheric composition, an elevated tropopause height is necessitated in response to an augmented PBL depth which is in line with the convective-radiative balance and influenced by planetary-scale features such as jet streams.The atmospheric thermal expansion resultant from global warming is manifested near the surface as a reinforced well-mixed boundary layer.To maintain energy balance, the tropopause must be elevated (in accordance with the inertial stability requirement

More Frequent Turbulent Conditions
Here, wind speed and potential temperature centered at 200 hPa have been employed to estimate grid point values of R i .These values have been used to construct an hourly time series for the period 1979-2018, covering the entire region of interest.Through an exploration of various critical values R i , a threshold of R i < 0.25 has been chosen to define severe turbulence (i.e., requiring seat belt fastening).At each grid point, the annual cumulative count of hours surpassing this threshold serves as an indicator of its impact on aviation.A Student t-test was applied to the full 40-year annual time series.
By comparing two consecutive 20-year period (namely 1979-1998 and 1999-2018), we found an increased occurrence of severe turbulence by greater than 300 h a year in the second period than that in the first period in a substantial area (about 1.66 × 10 8 km 2 ) within the global belt between 40 • S and 70 • N (Figure 5).For example, over the North Atlantic Corridor (NAC), there are two segments (at ~30 • N and at ~55 • N, corresponding to the polar and subtropical jets respectively), experiencing this dramatic increase in turbulent activity.Moreover, about 20.5% of the total area (~8.2 × 10 7 km 2 ) has experienced at least 400 additional hours (with ~7.4% experiencing >500 h) turbulent conditions over a typical year in the second period than that in the first period.A significant portion of the Southern and Pacific Oceans experienced a > 10% increase in total hours of severe turbulence (blue shading in Figure 5).Over the Pacific Ocean, the pattern closely resembles the changes in the Madden-Julian Oscillation (MJO) activity.
lative count of hours surpassing this threshold serves as an indicator of its impact on aviation.A Student t-test was applied to the full 40-year annual time series.
By comparing two consecutive 20-year period (namely 1979-1998 and1999-2018), we found an increased occurrence of severe turbulence by greater than 300 h a year in the second period than that in the first period in a substantial area (about 1.66 × 10 8 km 2 ) within the global belt between 40°S and 70°N (Figure 5).For example, over the North Atlantic Corridor (NAC), there are two segments (at ~30°N and at ~55°N, corresponding to the polar and subtropical jets respectively), experiencing this dramatic increase in turbulent activity.Moreover, about 20.5% of the total area (~8.2 × 10 7 km 2 ) has experienced at least 400 additional hours (with ~7.4% experiencing >500 h) turbulent conditions over a typical year in the second period than that in the first period.A significant portion of the Southern and Pacific Oceans experienced a > 10% increase in total hours of severe turbulence (blue shading in Figure 5).Over the Pacific Ocean, the pattern closely resembles the changes in the Madden-Julian Oscillation (MJO) activity.Figure 6 shows a representative NAC flight trajectory (marked by red dashed line), linking New York and London, with stippling denoting subregions with a change of severe turbulence days excessing 95% significance level based on a student t-test (p-value < 0.05).Regions with an increase in severe turbulence days for more than 10% (i.e., days with at least two hours exceeding the criterion) between two consecutive 20-year periods centered in 1993 and 2013 are highlighted with blue shading.The expansiveness of these areas, coupled with corroborating evidence from recent flight logs, underscores the growing concern surrounding the heightened occurrence of in-flight severe turbulence instances over the near-to-mid-term.
The distribution of wind kinetic energy is spatially heterogeneous and concentrated in proximity to the polar and subtropical jet streams (figures not shown for brevity).The Polar Jet Stream (PJS) is particularly impactful on commercial aircraft operations because of its contribution to strong headwinds or tailwinds along busy mid-latitude flight routes.However, it also gives rise to clear-air turbulence due to the pronounced vertical wind shear associated to PJS.Polar amplification (i.e., Arctic amplification, Appendix B) played an important role in the increased vertical shear around the PJS.This is partly because, a wavier PJS as a consequence of polar amplification (in communication with M. Wehner, 2023) expands the airspace influenced by PJS.0.05).Regions with an increase in severe turbulence days for more than 10% (i.e., days with at least two hours exceeding the criterion) between two consecutive 20-year periods centered in 1993 and 2013 are highlighted with blue shading.The expansiveness of these areas, coupled with corroborating evidence from recent flight logs, underscores the growing concern surrounding the heightened occurrence of in-flight severe turbulence instances over the near-to-mid-term.The distribution of wind kinetic energy is spatially heterogeneous and concentrated in proximity to the polar and subtropical jet streams (figures not shown for brevity).The Polar Jet Stream (PJS) is particularly impactful on commercial aircraft operations because of its contribution to strong headwinds or tailwinds along busy mid-latitude flight routes.However, it also gives rise to clear-air turbulence due to the pronounced vertical wind shear associated to PJS.Polar amplification (i.e., Arctic amplification, Appendix B) played an important role in the increased vertical shear around the PJS.This is partly because, a wavier PJS as a consequence of polar amplification (M.Wehner, personal communication 2023) expands the airspace influenced by PJS.In contrast, the increased turbulence in the vicinity of subtropical jets (Figure 7) is due primarily to the expansion of the Hadley cell.In Figure 7, a vertical transect along the 29 • N is taken as an example.The 'core' of the subtropical jet has been elevated consistently over the past 40 years.The trend, as indicated in the annual time series, satisfies the criteria of a t-test at the 95% confidence level.This trend is expected to persist as the climate continues to warm.Cconsequently, a far more pronounced increase in wind speed above 150 hPa and a reduction in wind speeds below 150 hPa (or a less significant increase around the 120 • W) signifies an increased vertical wind shear within the aviation domain (250-100 hPa).From Equation (2), this (increased denominator in the R i indicator) means a reduced R i value and increased likelihood of clear air turbulence.Corroborating this, Prosser et al. [18] systematically examined the Northern Hemispheric winter CAT and found CAT causing >0.5 g vertical acceleration of the aircraft could increase by 40-170% by 2050.
Trends in the reanalysis aligned with the interpretation that ongoing climate warming is progressively exerting an influence on aviation turbulence.Anthropogenic climate change is modifying the thermal structure of the stratified atmosphere (Figure 1b).Due to Arctic amplification, the temperature gradients reduce at the lower portion of the tropopause.In the upper troposphere, due to the sedimentation feature of CO 2 , the meridional temperature gradient increases.This equates to a superimposed zonal circulation that aids the upper-level westerlies but retards the lower-level westerlies.Consequently, the vertical shear is more concentrated at commercial aviation altitudes [51,52].In addition to the traditional polar amplification mechanism, the Polar Amplification of the Second kind (PAOSK), caused by air viscosity increasing with temperature, also contributes to the persistence of the increasing in aviation turbulence (Appendix B).With PAOSK, the reduced strength of PJS favors more frequent winter weather extremes.The associated aviation turbulence mostly is describable by a reduced R i .
29°N is taken as an example.The 'core' of the subtropical jet has been elevated consistently over the past 40 years.The trend, as indicated in the annual time series, satisfies the criteria of a t-test at the 95% confidence level.This trend is expected to persist as the climate continues to warm.Cconsequently, a far more pronounced increase in wind speed above 150 hPa and a reduction in wind speeds below 150 hPa (or a less significant increase around the 120°W) signifies an increased vertical wind shear within the aviation domain (250-100 hPa).From Equation (2), this (increased denominator in the Ri indicator) means a reduced Ri value and increased likelihood of clear air turbulence.Corroborating this, Prosser et al. [18] systematically examined the Northern Hemispheric winter CAT and found CAT causing >0.5 g vertical acceleration of the aircraft could increase by 40-170% by 2050.Trends in the reanalysis aligned with the interpretation that ongoing climate warming is progressively exerting an influence on aviation turbulence.Anthropogenic climate change is modifying the thermal structure of the stratified atmosphere (Figure 1b).Due to Arctic amplification, the temperature gradients reduce at the lower portion of the tropopause.In the upper troposphere, due to the sedimentation feature of CO2, the meridional temperature gradient increases.This equates to a superimposed zonal circulation that aids the upper-level westerlies but retards the lower-level westerlies.Consequently, the vertical shear is more concentrated at commercial aviation altitudes [51,52].In addition to the traditional polar amplification mechanism, the Polar Amplification of the Second kind (PAOSK), caused by air viscosity increasing with temperature, also contributes to the persistence of the increasing in aviation turbulence (Appendix B).With PAOSK, the reduced strength of PJS favors more frequent winter weather extremes.The associated aviation turbulence mostly is describable by a reduced Ri.

The Decline in Turbulence due to Severe Convective Instability
The sources of atmospheric turbulence impacting aviation are diverse and include regions characterized by strong vertical shear, severe convection and associated microbursts, and mountain lee waves [53] propagating both horizontally and vertically (depending on the Scorer parameter).Both severe convection and mountain waves can extend into the upper troposphere and lower stratosphere (UTLS), which coincides with the typical cruising altitude for commercial aircraft [12,36].As the spectrum of aviation turbulence is wide, it can be challenging to differentiate whether exceeding the R i threshold alone will result in passenger discomfort, route adjustments, or rerouting to evade severe storms.
To identify turbulent conditions associated with severe storms, we examine the CAPE value using ERA5 reanalysis data at a temporal-spatial resolution of hourly and quarter degree.CAPE values >2 kJ/kg is used as the threshold for convective instability.Around 80% of the globe experiences fewer than 300 h per year with CAPE values surpassing this threshold.However, certain regions, such as the tropical oceanic areas (e.g., the Bay of Bengal, the Gulf Mexico warm pool, and the Timor Sea) and portions of Central Africa (e.g., the Congo River Basin) can encounter 1200 to 2000 h annually with CAPE > 2 kJ/kg.
The change of the CAPE index is defined here as the annual mean total hours of CAPE > 2 kJ/kg for 1999-2018 minus those for 1979-1998, as shown in Figure 8a.Generally, there is a reduction in convective instability globally in the latter 20-year period than the earlier one.If 4 kJ/kg is used as the CAPE threshold, the areas with decreased convective activity then are much larger than those with increased convective activity, although the total global area satisfying the 4 kJ/kg criterion is much smaller.This trend also exists from the remote sensing products.For example, using the same cloud identification scheme of Chedzey et al. [30] and additionally setting cloud top temperature <−55 • C (and visibly thick) as the criteria for severe convective clouds (e.g., the cumulonimbus Cb clouds), the PATMOS-x series of radiance products provide ~6-hourly cloud cover (aerial coverage) maps.If Cb clouds are presence, it is assumed that they are present for the full 6 h.For a grid point, when there are more than five pixels detected having convective clouds, it is assumed to have convective clouds.The annual number of hours with convective clouds can be estimated for the 1979-2018 period, except for the limited polar hole and also tropical areas in the orbital gaps.For the tropical belt (30 • S-30 • N), the overall trend is a decrease.There are locations experiencing a ~17% reduction in cloud cover hours for the two time periods considered.Therefore, over the past four decades, the increase in KHI instability, indicated by low R i values, has been substantial, yet its impact is primarily confined to passenger discomfort.However, instances of severe turbulence, which could pose significant safety concerns, seem to be on a declining trend.This decline may be related to the global warming or due to decadal-to-multi-decadal climate natural variability.

Summary and Conclusions
Modern aircraft are engineered to endure severe turbulence and wing icing, so occurrences of these events are no longer being the primary triggers of fatal incidents.However, even mild turbulence can trigger discomfort among passengers and crew, and in certain instances, pose safety concerns leading to distress and injuries.Our investigation indicates that comfort and safety of air travel in relation to weather patterns are likely to gain increased significance in a warming climate.
Here, the canonical Richardson number (Ri) is calculated over the ERA5 reanalysis period and leveraged to identify regions experiencing increased turbulence.Beyond Ri as a comprehensive turbulence indicator, the convective available potential energy (CAPE) is examined as an index for turbulence induced by severe convective weather.By examining multiple indices associated with the atmospheric instability under present day conditions, we suggest that the likelihood of encountering increased turbulence during flights Figure 8a reveals that the pronounced reduction in CAPE has primarily occurred within the tropical belt (30 • S-30 • N).Studies show changes in the tropical basin interactions among the Atlantic, Indian and Pacific Oceans since the late 1990s [54,55].One of the significant changes on decadal timescale is the phase transition of the Interdecadal Pacific Oscillation (IPO [56]) from positive phase to a negative phase in late 20th century which led to an unprecedented acceleration of the Walker circulation.This in turn has amplified the wind-stress curl, resulting in enhanced ocean heat uptake [56].This enhanced Walker circulation also intensifies the Indonesian Throughflow (ITF [57]).The increased heat uptake in the Pacific was partly transported to the Indian Ocean through ITF.Because the temperature threshold for deep convection increases with the mean tropical temperature [58], the relative warming between different basins can be used to represent the different warming rate among the three tropical ocean basins.This inter-basin relative warming is defined as Trans-basin Variability (TBV) which effectively get the mean tropical warming removed.TBV is defined for both Indian-Pacific basins and Atlantic-Pacific basins (Figure 8b).Both Indian and Atlantic TBVs, operating through a Gill-type response in the atmosphere and an oceanic Bjerknes feedback, have contributed to a La Nina-like state in the Pacific [59,60].Through these mechanisms, deep convective activity over tropical Pacific has been suppressed over the past ~20 years.From Figure 8b, the IPO phase transition may occur again which may lead to opposite changes in Walker Circulation, wind stress curl and ITF to those during IPO negative phase.This may lead to enhanced CAPE index, thus raising the concerns of aviation safety in the coming decades.
Safety initiatives heavily depend on investigating past incidents and utilizing available data to detect less conspicuous or emerging patterns before accidents materialize.The insights from Figures 5 and 8 serve this purpose by shedding light on the evolving landscape of aviation turbulence.

Summary and Conclusions
Modern aircraft are engineered to endure severe turbulence and wing icing, so occurrences of these events are no longer being the primary triggers of fatal incidents.However, even mild turbulence can trigger discomfort among passengers and crew, and in certain instances, pose safety concerns leading to distress and injuries.Our investigation indicates that comfort and safety of air travel in relation to weather patterns are likely to gain increased significance in a warming climate.
Here, the canonical Richardson number (R i ) is calculated over the ERA5 reanalysis period and leveraged to identify regions experiencing increased turbulence.Beyond R i as a comprehensive turbulence indicator, the convective available potential energy (CAPE) is examined as an index for turbulence induced by severe convective weather.By examining multiple indices associated with the atmospheric instability under present day conditions, we suggest that the likelihood of encountering increased turbulence during flights would increase under future warmer climate.
Fortunately, the occurrence of severe turbulence episodes necessitating re-routing are becoming less common, partially due to the negative interdecadal Pacific oscillation (IPO) since late 1990s.An underlying rationale for the rise of clear air turbulence lies in the atmospheric dynamics under global warming.Specifically, as the atmosphere warms, the minimal Kolmogorov eddy size (d min ~ν3/4 , where ν is the kinematic viscosity of air) increases, while the maximum eddy size remains relatively unaffected by air viscosity.Thus, the distribution of turbulent kinetic energy across the eddy spectrum bin-levels intensifies with atmospheric warming.This phenomenon is particularly notable in the eddy size range pertinent to aviation turbulence, largely due to the proximity to d min , the lower limit of the eddy size spectrum.Our findings reinforce the ongoing influences of the warming climate on aviation.
Furthermore, aside from the rise of identified clear air turbulence, the reduced polar jets characterized by more concentrated shear and localized warming can increase PBL height.This holds turbulence-related consequences beyond the impacting maximum aircraft payload reduction [6].Thus, the interplay of maximum payload, fuel efficiency, environmental footprint, and safety is intertwined, both globally in relation to the vigor of polar and subtropical jet streams and locally with regard to PBL depth.form.Meanwhile, warmer air possesses a larger viscosity and hence affects the flow's shear structure.These two mechanisms are proven to significantly fortify the classical polar amplification but to different extents.A simple, force-restore type, box model sheds light on potential consequences of extreme weather events.
Polar amplification [64,65] sets up the background of climate warming impacts, especially for regional impact and mitigation and hence, formulates the basis of climate warming research.Here, using a proof-of-concept box model, we illustrate currently overlooked processes that can further reinforce positive feedback mechanisms.We consider the fluid medium of earth (plus the thermally active layer underneath) of two belt zones 32-45 • N and 45-80 • N.For the land surface, this would compose the atmosphere plus a soil layer of about 1 m depth.For oceanic regions, it would comprise the atmosphere plus the well mixed layer (a quarter of the thermocline ~50 m depth).earth-atmosphere system is vertically stratified, with temperature decreasing steadily upward within the troposphere.Above the tropopause, midlatitudes possess the lowest air temperature meridionally.Thus, for mid-to-high-latitudes, the tropopause is the location for both vertical and longitudinal reversals of air temperature.From the thermal wind relationship, this also is the location of jet streaks (i.e., jet core).Here we consider the domain of the troposphere (e.g., below 200 hPa) for two reasons: (1) the anthropogenic warming is vertically uniform within this regime; and (2) there is a steady meridional temperature gradient throughout this domain.
The second property makes it reasonable to assume a gross/equivalent temperature representing the entire troposphere in our idealized discussion of polar amplification ( [64][65][66]; reduction of the equator-to-pole temperature gradient by albedo and lapse rate feedbacks).The radiative-convective (RC) equilibrium temperatures of these two domains are assumed to be T e (for the extra-tropical region) and T 0 (for the polar domain).The drastic difference between T e and T 0 is due primarily to uneven solar insolation.Following existing research studying the Atlantic meridional overturning circulation (AMOC) [67], for both land surface processes [68,69], and for the El Niño-Southern Oscillation (ENSO) (e.g., Sun and Liu [70], SL96 henceforth), we propose the following coupled forcing-restoring system to illustrate the effects of poleward heat transportation by winds and oceanic flows.Consider, where T 1 and T 2 are respectively the temperatures of the extra-tropical and polar regions, α is the relaxation time scale for temperatures converging to their radiative-convective (RC) equilibrium values, and β is the dynamic coupling strength.Simpler than in SL96 is that the exchange between the two domains is internal, hence there is only a sign difference in Equations (A1) and (A2) for the second RHS term.The third term on the RHS of Equation (A2) is singled out to delineate the polar amplification effect on the solution.This profoundly changed the property of the feedback system, although all indicators of feedback strength (e.g., α * = γ(T e −T 0 ) α ) can be similarly defined as in SL96.While α * still measures the relative importance of thermodynamics and dynamics (feedback strength), arbitrary values of α * have unique solutions for stable temperatures T 1 and T 2 .That is, the exchange of heat between the two regimes, no matter the magnitude, stabilizes the Earth's climate from the impact of uneven solar insolation.
Transient events aside, the two temperatures have a tendency to being restored to T e and T 0 respectively (at the time scale of 1/α) but the interchange of thermal energy forces the two temperatures to stabilize at two values between T e and T 0 .The following discuss focuses on the stable solutions of T 1 and T 2 for the coupled system.

Appendix B.3. Results and Conclusions
The aim is to examine the response of the meridional temperature contrast at midlatitudes arising from a warming climate, using as simple a model as possible.The numbers here are used to illustrate mechanisms, the exact values should not be taken literally.In Figure A2, we have adjusted the relief at year 2000 so all experiments have the same starting values of stable temperatures T 1 and T 2 .The control experiments (CNTL results are the red lines in Figure A2) set the background of the canonical ice/snow albedo feedback induced polar amplification.As climate warms, T 2 increases faster (6.7 • C increment; from 288.89 K in 2000 to 295.59 K in 2100) than T 1 (11.3 • C increment; from 257.28 K in 2000 to 268.58 K in 2100) but their relative differences at the endpoints reduce with time (31.6 • C in 2000 to 27 • C in 2100).From comparison between the EXP01 and CNTL experiments, it is apparent that, the weakened PJS actually transfer heat more effectively, assuming a curvature appropriate for the reduced temperature contrast between the polar and extratropical domains.This further reduces the temperature contrast between the two domains.By 2100, the T 1 −T 2 is reduced to ~20.5 • C, about an additional 7 • C reduction from the CNTL.With the reduced T e and T 0 difference, the requirement of heat transport between the two regimes is reduced.However, the PJS system may resume a more efficient mechanism for transferring heat.In reality, it may involve an extended discharge of polar airmass to lower latitudes, such as via extreme winter storms.Increased air viscosity contributes in this same direction (yellow-brown lines in Figure A2), but is smaller than, but on the same order of magnitude as the curvature effect.The difference between the two stable temperatures is further reduced to ~16.7 • C by 2100.So, once the ice albedo-polar amplification effect initiates, the atmospheric circulation has a positive feedback mechanism to amplify, as warming continues.Here the PJS system also conceptualizes the wind driven ocean currents and their fair share in the poleward relay of heat transfer.As PJS becomes wavier, the heat exchanges between the two regimes may take more extreme forms.driven ocean currents and their fair share in the poleward relay of heat transfer.As PJS becomes wavier, the heat exchanges between the two regimes may take more extreme forms.For the fluidic sphere of the climate system, the equatorial (polar) region is the source (sink) of energy/heat and (angular) momentum.Advection by winds and oceanic flows is a fast response mechanism that reduces the meridional gradients, permi ing a milder climate for the overall ecosystem.As climate warms, polar amplification is changing the overall efficiency of this mechanism and may involve weather extremes as compensation.Aside from the traditional ice/snow albedo feedback, changes in internal physical properties of air and oceanic water also contribute.For the fluidic sphere of the climate system, the equatorial (polar) region is the source (sink) of energy/heat and (angular) momentum.Advection by winds and oceanic flows is a fast response mechanism that reduces the meridional gradients, permitting a milder climate for the overall ecosystem.As climate warms, polar amplification is changing the overall efficiency of this mechanism and may involve weather extremes as compensation.Aside from the traditional ice/snow albedo feedback, changes in internal physical properties of air and oceanic water also contribute.

21 Figure 1 .
Figure 1.Temperature profiles over polar and midlatitude stations (a) at present and (b) for a warming climate.In (a), soundings taken at 12Z, 25 November 2021 over the south pole and mid-to-high latitude stations: Amundsen-Scott (90°S; 0°E), Syowa (69°S; 39.58°E), SAVC (−45.78°S;67.5°W) and NZNV (46.41°S; 168.31°E).At the tropopause, temperature gradients reverse both vertically and meridionally, explaining why the jet core occurs near the tropopause.Panel (b) illustrates how anthropogenic climate change differentially affects the polar and mid-latitude temperature profiles.The net effect is equivalent to a superimposed circulation that counters the present climate mean thermal winds at the lower troposphere (1000-600 hPa) and enhances the present climate mean thermal winds at the upper troposphere (600-200 hPa).The commercial aviation operating domain is hatched in (b).The illustrations were produced using GrADS (available at http://cola.gmu.edu/grad.php;accessed on 16 June 2018).

Figure 1 .
Figure 1.Temperature profiles over polar and midlatitude stations (a) at present and (b) for a warming climate.In (a), soundings taken at 12Z, 25 November 2021 over the south pole and mid-tohigh latitude stations: Amundsen-Scott (90 • S; 0 • E), Syowa (69 • S; 39.58 • E), SAVC (−45.78 • S; 67.5 • W) and NZNV (46.41 • S; 168.31 • E).At the tropopause, temperature gradients reverse both vertically and meridionally, explaining why the jet core occurs near the tropopause.Panel (b) illustrates how anthropogenic climate change differentially affects the polar and mid-latitude temperature profiles.The net effect is equivalent to a superimposed circulation that counters the present climate mean thermal winds at the lower troposphere (1000-600 hPa) and enhances the present climate mean thermal winds at the upper troposphere (600-200 hPa).The commercial aviation operating domain is hatched in (b).The illustrations were produced using GrADS (available at http://cola.gmu.edu/grad.php;accessed on 16 June 2018).
Remote Sens. 2024, 16, x FOR PEER REVIEW 7 of 21 communication 2014).Upper atmospheric features (e.g., jet streams) are connected to the development of the PBL [42].PBL thicknesses have been collected over the past 40 years.From ERA5 hourly data, the maximum reachable PBL depths of each day are collected and averaged into monthly values (Figure2).

Figure 2 .
Figure 2. The maximum global PBL depths over the globe for March.Using hourly data, daily maximum PBL depths are extracted and used as an index.Maximum reachable March PBL depths (color shading is in meters) at each grid point are obtained by averaging the largest daily maximum PBL depths.Then the 20-year mean value is calculated and subtracted to obtain Figure 3.The background and color shades are produced using GRADS (available at http://cola.gmu.edu/grad.php;accessed on 15 May 2022).

Figure 3
Figure 3 elucidates alterations in March PBL depths between two consecutive 20-year periods 1979-1998 and 1999-2018.The most striking feature is the marked increase in PBL depth.Examination of PBL depth time series at selected locations reveals a prevailing trend of thickening PBL depths across the entire reanalysis period, punctuated by fluctuations stem from a blend of both nature and anthropogenic drivers.For example, Guo et al. [43], based on radiosonde data over China, found a spatially uniform increase in PBL depth from 1979-2003, followed by a spatially non-uniform decrease of PBL depth from 2003-2016.The decrease was ascribed to alterations in land-surface characteristics affecting the surface energy budget.This overall increasing trend of PBL depth is explicable using fundamental physical principles.The PBL, being the layer most affected by friction, often has its depth parameterized as proportional to the square root of the eddy viscosity coefficient[44].This coefficient functions as a macroscopic representation of air's kinematic viscosity[45].This provides a pathway through which a warming climate could influence aviation turbulence, as discussed in the Appendix B. Multiple factors can govern PBL depth, including surface heating, surface energy partition, the free atmospheric lapse rate, among others.Viscosity-related feedbacks persist as a discernible signal within a warming climate context.However, changes in ground surface heat fluxes stemming from changes in Bowen ratio[46] also acts as a regulatory factor.These internal (viscosity-related) and external (surface heat fluxes) factors, operating in synergy or opposition, collectively modulate the evolution of actual PBL depth with climate warming.

Figure 2 .
Figure 2. The maximum global PBL depths over the globe for March.Using hourly data, daily maximum PBL depths are extracted and used as an index.Maximum reachable March PBL depths (color shading is in meters) at each grid point are obtained by averaging the largest daily maximum PBL depths.Then the 20-year mean value is calculated and subtracted to obtain Figure 3.The background and color shades are produced using GRADS (available at http://cola.gmu.edu/grad.php; accessed on 15 May 2022).Remote Sens. 2024, 16, x FOR PEER REVIEW 8 of 21

Figure 4 .
Figure 4.The differences of maximum PBL depths (m) between two consecutive 20-year periods; (1999-2018) versus (1979-1998), for four months representing respectively winter (January), spring (April), summer (July), and fall (October).The lowest row is the climatology of the maximum PBL height over the 1979-2018 period (mean over this period).The central row is absolute changes (same as Figure 3, (1999-2018) minus (1979-1998)).The upper row is the percentage changes, or the central row divides the lowest row.

Figure 4 .
Figure 4.The differences of maximum PBL depths (m) between two consecutive 20-year periods; (1999-2018) versus (1979-1998), for four months representing respectively winter (January), spring (April), summer (July), and fall (October).The lowest row is the climatology of the maximum PBL height over the 1979-2018 period (mean over this period).The central row is absolute changes (same as Figure 3, (1999-2018) minus (1979-1998)).The upper row is the percentage changes, or the central row divides the lowest row.

Figure 5 .
Figure 5. Increased instability over the globe (a 40°S-70°N latitude belt is shown for clarity).The severe atmospheric turbulence locations were generated from the gradient Richardson Numbers (Ri) calculated from ERA5 hourly atmospheric parameters during the period 1979-2018.The critical Ri value was set as 0.25.The total annual number of hours with Ri below the critical value is plotted.Contour lines are increases in annual mean number of hours with Ri < 0.25, between1979-1998 and

Figure 5 .
Figure 5. Increased instability over the globe (a 40 • S-70 • N latitude belt is shown for clarity).The severe atmospheric turbulence locations were generated from the gradient Richardson Numbers (R i ) calculated from ERA5 hourly atmospheric parameters during the period 1979-2018.The critical R i value was set as 0.25.The total annual number of hours with R i below the critical value is plotted.Contour lines are increases in annual mean number of hours with R i < 0.25, between 1979-1998 and 1999-2018.The blue shaded areas identify regions with increases of at least 10% in annual-average number of hours of severe atmospheric turbulence between the two 20-year periods.The map background and color shading were produced using GrADS (available at http: //cola.gmu.edu/grad.php;accessed on 1 October 2022).

Figure 6 .
Figure 6.Similar to Figure 5, but for the North Atlantic Corridor (NAC).The number of hours below the critical Ri (0.25) in a given year is defined as the index.Contour lines show increases in the annual mean number of hours with Ri < 0.25, between 1979-1998 and 1999-2018.The blue shaded regions have annual-average numbers of hours with severe atmospheric turbulence increased by at least 10%, between the two consecutive 20-year periods 1979-1998 and 1999-2018.A t-test is performed for each grid of a 40-year annual time series of this index (i.e., the number of hours with Ri < 0.25 in a year).Stippled regions passed the t-test at the 95% confidence level.More than 60% of the traditionally defined NAC passed the t-test.The red dashed line indicates the JFK-LHR route.The map backgrounds and contours were produced using GRADS (available at http://cola.gmu.edu/grad.php;accessed on 1 October 2022).

Figure 6 .
Figure 6.Similar to Figure 5, but for the North Atlantic Corridor (NAC).The number of hours below the critical Ri (0.25) in a given year is defined as the index.Contour lines show increases in the annual mean number of hours with Ri < 0.25, between 1979-1998 and 1999-2018.The blue shaded regions have annual-average numbers of hours with severe atmospheric turbulence increased by at least 10%, between the two consecutive 20-year periods 1979-1998 and 1999-2018.A t-test is performed for each grid of a 40-year annual time series of this index (i.e., the number of hours with Ri < 0.25 in a year).Stippled regions passed the t-test at the 95% confidence level.More than 60% of the traditionally defined NAC passed the t-test.The red dashed line indicates the JFK-LHR route.The map backgrounds and contours were produced using GRADS (available at http://cola.gmu.edu/grad.php;accessed on 1 October 2022).

Figure 7 .
Figure 7. Annual mean trends (color shades) in wind magnitude in the subtropical jets as a function of latitude and atmospheric pressure over the period 1979-2018.Values were calculated as the temporal mean over the (1999-2018) period minus the temporal mean over the (1979-1998) period, based on ERA-5 reanalyses.Dashed and solid lines are, respectively, the contours of the magnitude of winds for the second (1999-2018) and first (1979-1998) periods.Wind speeds generally increased aloft (above 150 hPa), whereas they mainly decreased below 150 hPa.The map backgrounds and contours were produced using GRADS (available at h p://cola.gmu.edu/grad.php;accessed on 20 March 2022).

Figure 7 .
Figure 7. Annual mean trends (color shades) in wind magnitude in the subtropical jets as a function of latitude and atmospheric pressure over the period 1979-2018.Values were calculated as the temporal mean over the (1999-2018) period minus the temporal mean over the (1979-1998) period, based on ERA-5 reanalyses.Dashed and solid lines are, respectively, the contours of the magnitude of winds for the second (1999-2018) and first (1979-1998) periods.Wind speeds generally increased aloft (above 150 hPa), whereas they mainly decreased below 150 hPa.The map backgrounds and contours were produced using GRADS (available at http://cola.gmu.edu/grad.php;accessed on 20 March 2022).
Remote Sens. 2024,16,  x FOR PEER REVIEW 14 of 21Safety initiatives heavily depend on investigating past incidents and utilizing available data to detect less conspicuous or emerging patterns before accidents materialize.The insights from Figures5 and 8serve this purpose by shedding light on the evolving landscape of aviation turbulence.

Figure 8 .
Figure 8. Trends in convective instability.Panel (a) uses the same approach as for Figure 5 but identifies the annual number of hours with CAPE exceeding 2 kJ/kg as the index.Color shades show changes in convective turbulence over the past 40 years, defined as the mean annual number of hours during (1999-2018) minus the mean annual number of hours during (1979-1998).Note that, globally, regions with reductions in convective instability are dominant.Panel (b) is a sliding-window of 20-year trends of inter-basin SST differences (TBV: Trans-basin Variability) for Atlantic-minus-Pacific trans-basin SST gradient (ATL-PAC) and Indian-minus-Pacific TBV (IND-PAC).SST is from the Extended Reconstructed Sea Surface Temperature (ERSST).Indian (IND), Pacific (PAC), and Atlantic (ATL) basin SST are calculated between 20°S to 20°N, respectively for longitudinal ranges of 21°E-120°E, 121°E-90°W, and 70°W-20°E.

Figure 8 .
Figure 8. Trends in convective instability.Panel (a) uses the same approach as for Figure 5 but identifies the annual number of hours with CAPE exceeding 2 kJ/kg as the index.Color shades show changes in convective turbulence over the past 40 years, defined as the mean annual number of hours during (1999-2018) minus the mean annual number of hours during (1979-1998).Note that, globally, regions with reductions in convective instability are dominant.Panel (b) is a sliding-window of 20-year trends of inter-basin SST differences (TBV: Trans-basin Variability) for Atlantic-minus-Pacific trans-basin SST gradient (ATL-PAC) and Indian-minus-Pacific TBV (IND-PAC).SST is from the Extended Reconstructed Sea Surface Temperature (ERSST).Indian (IND), Pacific (PAC), and Atlantic (ATL) basin SST are calculated between 20 • S to 20 • N, respectively for longitudinal ranges of 21 • E-120 • E, 121 • E-90 • W, and 70 • W-20 • E.

Figure A2 .
Figure A2.Stable temperatures (°C) for the two regimes under different assumptions of the heat exchange coefficient: control (red lines), with curvature effects (blue lines) and with further consideration of sensitivity of viscosity to temperature (yellow brown lines).

Figure A2 .
Figure A2.Stable temperatures ( • C) for the two regimes under different assumptions of the heat exchange coefficient: control (red lines), with curvature effects (blue lines) and with further consideration of sensitivity of viscosity to temperature (yellow brown lines).
1.5/L/15.Here L is Blackadar turbulence length scale.A quick view of global ELLROD (a popular aviation turbulence indicator) values indicate that downstream (eastward) regions of the Rockies, Iranian Plateau, and the Tibetan Plateau are prone to MWT.