Future Global Convective Environments in CMIP6 Models

The response of severe convective storms to a warming climate is poorly understood outside of a few well studied regions. Here, projections from seven global climate models from the CMIP6 archive, for both historical and future scenarios, are used to explore the global response in variables that describe favorability of conditions for the development of severe storms. The variables include convective available potential energy (CAPE), convection inhibition (CIN), 0–6 km vertical wind shear (S06), storm relative helicity (SRH), and covariate indices (i.e., severe weather proxies) that combine them. To better quantify uncertainty, understand variable sensitivity to increasing temperature, and present results independent from a specific scenario, we consider changes in convective variables as a function of global average temperature increase across each ensemble member. Increases to favorable convective environments show an overall frequency increases on the order of 5%–20% per °C of global temperature increase, but are not regionally uniform, with higher latitudes, particularly in the Northern Hemisphere, showing much larger relative changes. The driving mechanism of these changes is a strong increase in CAPE that is not offset by factors that either resist convection (CIN), or modify the likelihood of storm organization (S06, SRH). Severe weather proxies are not the same as severe weather events. Hence, their projected increases will not necessarily translate to severe weather occurrences, but they allow us to quantify how increases in global temperature will affect the occurrence of conditions favorable to severe weather.

This degree of certainty is not common to all impactful climate extremes (National Academies of Sciences, Engineering, & Medicine, 2016). Severe thunderstorms (ones producing hail, tornadoes, or damaging straight-line winds) hold greater uncertainty, with confidence low in both observed and projected changes. This is despite a long-term awareness of the need to understand how convective extremes are changing since the earliest climate change assessments (Griffiths et al., 1993). The large annual impacts of these events motivate much of this need. Every year hail causes widespread damage to crops and property, while tornadoes in particular can produce intense localized destruction, and are more frequent and violent in the U.S. than in any other country (LLoyd's, 2013). In the last decade, Contiguous United States (CONUS) average annual losses from severe thunderstorms have trended upward and due to their regularity have exceeded those due to hurricanes (Sander et al., 2013). While evidence of changes is growing for the historical period, the challenge for observed trends is that these phenomena can be difficult to assess directly through time as historical data are inconsistent or incomplete (Doswell, 2007). This has lead to the use of environmental proxies (Brooks et al., 2003) and convective allowing simulations as a way to explore current and future trends. In fact, the observed trends have been accompanied by evidence of increases to severe thunderstorm hazards through proxies for the frequency of tornadoes and tornado outbreaks (Gensini & Brooks, 2018;, hailstorms (Mohr & Kunz, 2013;Rädler et al., 2018;Tang et al., 2019), and the overall frequency of severe thunderstorm events (Koch et al., 2021;. The main challenge for understanding the impact of climate change on severe thunderstorms is that they are not resolved in traditional climate modeling frameworks (Bindoff et al., 2013;Kunkel et al., 2013). Early studies on the effects of climate change on severe convection were limited by horizontal model resolution and limited computational resources (McMaster, 1999;Niall & Walsh, 2005). Increasing spatial and temporal resolution since the CMIP3 generation led to several studies that explored projected changes in convective environments in response to increasing CO 2 concentration (Allen et al., 2014b;Brimelow et al., 2017;Diffenbaugh et al., 2013;Glazer et al., 2020;Hoogewind et al., 2017;Mohr et al., 2015;Púčik et al., 2017;Rädler et al., 2019;Seeley & Romps, 2015;Singh et al., 2017;Trapp et al., 2007Trapp et al., , 2009, focusing primarily on the CONUS and more recently Europe. These convective environmental proxies are functions of quantities such as convective available potential energy (CAPE), convection inhibition (CIN), and wind shear (S, in particular for 0-6 km AGL atmospheric depth, S06, and 0-1 km AGL depth, S01, or in the form of storm relative helicity, SRH) and co-variate products for severe weather days. Generally, these studies have highlighted increasingly available instability due to warming of the lower troposphere, with the largest increases resulting from increasing mixing ratios proximal to warm tropical oceans (Allen et al., 2014b;Diffenbaugh et al., 2013;Glazer et al., 2020;Rädler et al., 2019;Singh et al., 2017). This increase has been shown to play a role in extending the duration of the year favorable to severe convective storms, and perhaps reducing the frequency during summer when CIN is maximized (Diffenbaugh et al., 2013;Hoogewind et al., 2017). While early studies highlighted the potential for S06 to decrease (Trapp et al., 2007), these have been tempered regionally by evidence that these decreases occur on days that have low CAPE and vary between continents, while appearing to be outweighed by increases during the spring and fall (Allen et al., 2014b;Diffenbaugh et al., 2013;Púčik et al., 2017). However, while such assessments of change can provide useful context, not all favorable environments result in convective storms.
To address this issue, studies using dynamical downscaling , 2015Hoogewind et al., 2017;Prein et al., 2015;Rasmussen et al., 2020;Trapp & Hoogewind, 2016;Trapp et al., 2019) have highlighted that, over North America, environment-based assessments indicate higher rates of occurrence in both present and projected frequency as compared to the realized frequency using convection allowing models, as these methods explicitly include the frequency of storm initiation. The general consensus of these projections has been toward an increase to severe convective frequency in the transition seasons, and a shift toward more extreme precipitation events, echoing results from environment-based studies. The approach also allows consideration of hazard specific projections, as expected changes to hail have a greater number of competing factors (Trapp et al., 2019), while tornadoes also rely on low-level instability (Trapp & Hoogewind, 2016). These simulations have suggested that while environments may increase in instability and hence would have higher potential updraft velocity (predicted through max = √ 2 ⋅ ), the resulting updrafts in a warmed climate appear to yield lower updraft velocity increases than would be expected due to processes excluded from theoretical parcels, such as entrainment or condensate loading (Trapp & Hoogewind, 2016). While explicitly modeling convection does allow for increased insight into the non-conditional occurrence of these storms, the method is computationally expensive as it requires resolutions of 4 km or higher, and thus is generally considered impractical beyond continental scale projections, or confined to sub-samples of the future environment (Prein et al., 2015). Recent reviews of 10.1029/2021EF002277 3 of 21 our understanding of severe convective storm climatology and future projections have highlighted the paucity of our global understanding of the expected changes of severe thunderstorm events to climate change (Allen, 2018;Raupach et al., 2021), suggesting that greater insight on a global scale is necessary. In lieu of appropriate global simulations at the convective scale, the best available approach is to consider changes to convective environments in climate model output, while understanding the limitations of this approach.
The main focus of this manuscript is to explore the projected changes in convective environments in response to increasing atmospheric CO 2 concentration at the global scale, and to explore the competing physical processes that are known to govern severe weather occurrence. Individual convective parameters (e.g., CAPE, SRH, and S06) are computed at a 6-hourly temporal resolution for seven CMIP6 global climate models which, at the time of writing, provide the necessary output for the calculation. They are first explored individually, to infer a greater understanding of the contributions to changes to organized storms. We then include a suite of broadly accepted severe weather indices in an effort to gain a clearer picture of the interaction of convective environments and to quantify the effect that our changing climate will have on the occurrence of severe weather at the global scale.
In using convective proxies, we acknowledge that they represent a useful tool in severe weather modeling, but we do realize they have several limitations. The first is that despite a "perfect" severe thunderstorm environment, there is no guarantee that an environment will actually produce a storm, which leads to proxies being an overestimate of frequency both in the present and future (Hoogewind et al., 2017;Trapp et al., 2009). Second, by using a fixed proxy for occurrence, this assumes that the relationships between the proxies and the events remain stationary throughout time, which may not necessarily be the case given the limitations of parcel theory (Raupach et al., 2021). Other caveats can be due to inherited properties of the global climate models themselves, related to the type of convective precipitation schemes, mistimed diurnal peaks, and coarse spatial resolution (Allen et al., 2014a;Diffenbaugh et al., 2013;Seeley & Romps, 2015).
Despite these limitations, it is still possible to draw significant and robust conclusions with respect to the future trajectories of severe weather activity in response to climate change. In fact, to address the uncertainties listed above, we consider a range of appropriate steps, novel as compared to prior studies. First, we focus on quantifying percent changes as a function of global temperature change, in addition to time, as each model projection and scenario has a different response to forcing and a different underlying climatology. Second, variables are calculated at 6-hourly resolution, but properties are examined on a seasonal basis, and where appropriate aggregated at regional scale. The advantage of such an approach is that it smooths the "weather-scale" variations that are not the goal of climate modeling and focuses on the overall projected changes. Third, we consider a suite of proxies, rather than a single approach to estimate changes in frequency, which use different set of convective parameters and are formulated in different ways to encompass a wide set of possibilities, allowing for a multi-index based view of the projected changes. These calculations, based on huge volume of CMIP6 data, are computationally complex and costly, and to accomplish them, we also demonstrate a novel approach to data processing using Pangeo cloud computing.
The study is structured as follows; in Section 2, we present the convective environments we computed and the severe weather proxies included in the study. In Section 3, we present the data included in our analysis, and describe the relevant parameters, experiments, models, and temporal and spatial resolution, together with the methods and tools used in the computation and analysis of the data. Section 4 presents a global view of the calculated convective parameters and their expected changes in a future scenario. Section 5 discusses the approach and benefits in evaluating the results as a function of temperature trajectories, and Sections 6-10 summarizes the findings.

Convective Environments and Severe Weather Proxies
Convective available potential energy (CAPE; J kg −1 ) is a measure of the potential for instability of an atmospheric parcel (Emanuel, 1994). It corresponds to the integrated virtual temperature difference between a theoretical parcel and the environmental profile while the parcel is warmer than the environment, representing the positively buoyant portion of the parcel profile; convection inhibition (CIN; J kg −1 ) corresponds, instead, to the negatively buoyant portion. CAPE and CIN were calculated for two parcels (i.e., surface and mixed-layer approximately to about lowest 50mb, ML50mb), using a 700 Pa pressure increment in the integration procedure to improve processing speed while maintaining accuracy (Bryan, 2008). Moist adiabats were calculated using the pseudoadiabatic liquid only approach. Surface and ML50mb CAPE products were very similar for our application, so we limited our attention to the ML50mb product. Both CAPE and CIN are evaluated only for time steps with surface temperature greater than 0°C; we note that CIN is here expressed as a positive quantity, rather than negative.
SRH is a measure of the potential for updraft rotation in supercells, and is often applied to anticipate rotating and tornadic storms (Markowski & Richardson, 2011). It is a quantity proportional to streamwise vorticity and storm-relative winds, defined as: where V is the ground relative wind vector m s −1 , C is the storm motion m s −1 , ω is the horizontal vorticity vector s −1 , and the resulting parameter SRH has units of m 2 s −2 . SRH is typically calculated by an integration of some depth depending on the application, most commonly between the surface and 1 km AGL, and surface to 3 km AGL, both included in this study (and referred to as SRH03 and SRH01). To estimate the storm motion C, we follow the Bunkers et al. (2000) approach, which is based on separate empirical storm motion estimates for cyclonic and anticyclonically rotating storms. Linear interpolation with respect to the height field is used to estimate U and V at 13 vertical levels at 500 m increments from the surface to 6 km, and then the non-pressure weighted mean wind is calculated for this layer. From this mean wind, right moving (RM) and left moving (LM) storm motions are then determined using the respective relationships (see Bunkers et al., 2000 for details). RM and LM SRH values are then combined by keeping the highest value in absolute magnitude for each time step and grid point.
Wind shear between 0 and 6 km (S06) was calculated instead between the surface meridional and zonal wind components (vas, uas) and the average 450-550 hPa, which is able to approximate the 6 km winds for the majority of the global domain, except for few areas with highly elevated terrain (u 6 km , v 6 km ).
To investigate the expected changes in severe weather, we considered a suite of commonly used proxies, which combine two to four of the convective parameters. These proxies represent the occurrence of conditions that are favorable for severe weather. If we were to use these indices to quantify the expected number of severe weather occurrence across the globe, this would necessitate their recalibration for each region considered. This is however beyond the scope of this work, as we instead use them as a simple tool to quantify the relative occurrence of these types of environments in future scenarios with respect to the historical period.
We start with the product CAPE ⋅ S06 (SEV) conditional to S06 being greater than 5 m s −1 , and CAPE being greater than 100 J kg −1 (Diffenbaugh et al., 2013). In some formulations, this index has been coupled with the occurrence of convective precipitation, as a way to include thunderstorm initiation (Trapp et al., 2009). We then considered a modified SEV, for which the weight of S06 is increased such as, CAPE ⋅ S06 1.67 (SEV γ ) (Allen et al., 2011;Brooks et al., 2003), which has been applied for several continents. The two indices are then normalized by 10,000 and 25,000 respectively, and set to 1 when equal or greater than 1, and 0 otherwise. In fact, transforming all the indices in dichotomous variables (0s and 1s) allows us to focus on changes in their frequency, rather than their magnitude. These indices have been shown to be able to discriminate between environments conducive to the formation of severe and significant severe thunderstorms, or indicate the higher probability of occurrence for a significant severe thunderstorm (producing hail exceeding 5 cm, an F2+ tornado, or winds exceeding 120 km/hr). The normalizing factors for SEV and SEV γ are overall lower than the ones used in some literature, which was estimated based on proximity soundings, this results in a somewhat less strict definition of severe weather conditions.
We then include two helicity indices, a classic approach in the product of CAPE and SRH03 (EHI), and one that gives more weight to the shear component, equal to the product of CAPE, SRH03, and SRH01 (EHI2). Both approaches use the absolute value of SRH, and hence are valid for both right and left moving storms across the hemispheric divide. EHI is normalized by 1.60 × 10 5 and EHI2 by 3.60 × 10 6 , and once again they are set to 1 when equal or greater than 1, and to zero otherwise. EHI by combining helicity and instability is used within the CONUS domain to assess the potential for supercells capable of producing tornadoes (Rasmussen, 2003).
A version of EHI2 was also introduced to evaluate the number of tornadoes within an outbreak and assess the occurrence of tornadoes outbreaks ).
All the above indices, SEV, SEV γ , EHI, and EHI2, have also been evaluated conditional to CIN being less than 75 J kg −1 (not shown) and 125 J kg −1 , to include the potential role of inhibition on moderating the likelihood of occurrence of severe weather.
We also included a modified supercell composite parameter (SCP), a multiple component index that is intended to indicate conditions favoring supercell thunderstorms (Thompson et al., 2002). The index include CAPE, SRH03, S06, and CIN, and can take different forms in literature. For this application, we defined SCP as follows: conditioned to CIN being less than 125 J kg -1 . Once again, by using the absolute value of SRH 0-3 km, SCP is valid for both right and left moving storms, and is set to 1 when equal or greater than 1, and to 0 otherwise. It is worth noting than the original formulation of SCP includes the most unstable CAPE, rather than MLCAPE.
Finally, we include two indices that are formulated quite differently from these previous approaches: using monthly averages rather than subdaily values, and intended to model the number of monthly occurrences of tornado and hail events. The tornado environment index, TEI (Tippett et al., 2012(Tippett et al., , 2014, and the hail environment index, HEI (Allen et al., 2015), have been demonstrated to perform well in representing the monthly climatological frequency of hail and tornado occurrence and their interannual variability over the CONUS, for both reanalysis and model data. These indices have been carefully calibrated for the CONUS domain to relate the climatological monthly number of tornado and hail events over the US to climatological monthly averages of atmospheric parameters from reanalysis data using a Poisson regression.

Data and Methods
The CMIP6 ensemble is a collection of GCMs from 41 modeling groups from around the world. Each modeling group contributed to the CMIP6 archive with one or more models and a variety of experiments, as defined by the World Climate Research Program . In this work, we analyze two sets of experiments: historical and ssp585. The historical scenario is an all-forcing simulation of the recent past and ends in 2014. The ssp585 (SSP5-8.5) scenario, represents an update of the CMIP5 RCP8.5 scenario, based on shared socioeconomic pathway 5 (SSP5) (O'Neill et al., 2016), and it represents a high-end scenario in terms of emissions and global temperature increase.
For accurately computing convective environments, we require data at a sub-daily temporal resolution and with a relatively fine vertical resolution. The variables required are temperature, specific humidity, pressure, zonal and meridional wind components, for both surface and all vertical levels; the corresponding "variable id" are tas, huss, ps, uas, vas, for surface variables, and ta, hus, ua, va for the 3D profiles. The resolutions we included, "table id" within the CMIP6 database, are 6hrLev for the 3D data, with a temporal 6-hourly resolution and a vertical resolution of the original model sigma levels, and 3hr, with a temporal 3-hourly resolution for the surface variables. Both type of table ids report instantaneous measurements of the variables, therefore from the 3hr fields we extracted time stamps matching the 6hrLev time coordinate. Two variables, uas, vas, for model CNRM-ESM2-1 and scenario ssp585 at the time of the analysis were available for a table id of E3hr, which represents an average value for the 3-hr interval; given the overall goal of the project, and the nature of the variables, we still treated these data as instantaneous, and matched them to the 6-hourly data.
The choice of which models to use, together with these two scenarios, was based primarily on data availability at the start of the project, but it also allowed for our analysis to include the widest range of temperature increase to better quantify uncertainty (see Section 5). In particular, we computed the convective environments for the years 1980-2014 for the historical experiment, and 2015-2100 for the SSP5-8.5 experiment. Table 1 summarizes the models included in the analysis, together with details about the version, member, and the horizontal and vertical resolution of each model. For comparison and cross-validation of the historical climatological performance of these models, we computed the same variables from the ERA-Interim (Dee et al., 2011) reanalysis (see Supporting Information S1 for more details).
The data described above were uploaded to the Google Cloud CMIP6 Public Data set, as part of a broader collaboration with Google Cloud (https://cloud.google.com/blog/products/data-analytics/new-climate-model-data-now-google-public-datasets). In the process, they were also converted from NetCDF to the cloud-optimized Zarr format in order to facilitate more efficient processing. All computations and analysis were performed in python and on the Pangeo NSF-funded cloud based computing platform (Abernathey et al., 2021). This allowed for ease of navigation and management of hundreds of TB of data. The workflow used for this project leverages xarray (Hoyer & Hamman, 2017), dask (Rocklin, 2015), xcape, xgcm (Abernathey et al., 2020), xhistogram, xarrayutils, and xesmf (Zhuang et al., 2021).
In particular, the computation of CAPE, together with storm relative helicity (SRH), was facilitated by the xcape python package (https://github.com/xgcm/xcape). xcape wraps an existing Fortran CAPE routine (Bryan, 2008) in python. Combining this efficient low-level numerical implementation with Dask's ability to scale out over many nodes enabled us to effectively process hundreds of TB of data in parallel.
Overall, the workflow was as follows: first we calculated CAPE, CIN, SRHs (through xcape), and S06, from the input variables; then we computed the severe weather proxies and aggregated the data to monthly totals (for the indices) and monthly averages (for the environments), and to regional values, and wrote these results back to cloud storage. Finally, we computed statistics and histograms.
In terms of statistics, the significance of the difference in means for the spatial maps is tested through a two-tailed t-test based on the interannual seasonal variability, at a significance level of 0.05. Agreement across models is mentioned when all models (100%) agree in the sign of the changes, unless indicated otherwise (as in Figures 5  and 7).
All time series of the regional averages and 2D histograms are area weighted, and spatial maps are obtained after a regridding to a common 1 × 1°grid, to correctly compare models with different resolutions.

Convective Environments in CMIP6 Models
We  Climatologically, seasonal ensemble mean CAPE (Figure 1, first column) reflects the cycle of CAPE coinciding with the transition seasons and hemispheric summers/winters. Seasonally CAPE increases in the relevant convectively active regions including those associated with the Monsoons, consistent with reanalysis climatology ( Figure S1, Supporting Information S1). Changes by midcentury suggest strong agreement for increases to CAPE across the ensemble (hatching) and statistical significance at the 0.05 two-tailed level (stippling) for the vast majority of the domain (Figure 1 for CAPE, second and third columns). Relative to the historical period, increases are on the order of 0.5σ over the tropical regions, while increases in the extratropics range between 1 and 1.5σ. Over all continents with the exception of desert regions, there are robust increases with strong ensemble agreement, though relative increases are larger over the NH in comparison to the SH by midcentury. Over the high latitudes of the Arctic in the NH, there are robust increases already by mid-century for CAPE. This is consistent with recent research considering convection in these regions (Poujol et al., 2021). However, the magnitude of these changes reflects climatologically small values of mean CAPE, and thus, while significant, reflect comparatively small increases associated with the surface temperatures warming to above freezing and corresponding increased moisture capacity. Regions with little to no change are primary in the subtropical bands, and mostly confined to oceanic regions upstream of the continents. These changes are small and insignificant for the majority of seasons, with the exception of significant decreases over the eastern Atlantic during the boreal spring, and over the SH during winter. End century changes reflect a continuance of the signals seen for midcentury, with magnitudes increasing to 2-3σ for much of the NH, and 1-2σ for the SH. Increases are once again larger toward the poleward regions. Negative signals, mostly located on the ocean, also increase in robustness, with decreases reaching 1σ. These projections suggest that over all regions that currently experience thunderstorms and severe thunderstorms on a regular basis, CAPE will significantly increase by the end of century.
In contrast to CAPE, CIN is restricted to regions where warming of the midlevel column occurs. This occurs either by adiabatic descent during compensating subsidence in the tropical circulation, or downstream of high terrain where airmasses with elevated mixed-layer characteristics are generated by diabatic mixing and advected by the midlatitude flow; its spatial distribution is also consistent with reanalysis ( Figure S1, Supporting Information S1). Changes to CIN by midcentury are comparatively small compared to those of CAPE, with increases generally no more than 0.5-1σ across all seasons, though still with robust increases across the ensemble. The largest increases to CIN are seen in the respective transition seasons, in particular spring, and are more substantive over the NH. The main exception to this change is over the higher latitudes in both hemispheres, where CIN does not show robust increases. End century projections contrast those of the midcentury, with widespread robust increases to CIN across the ensemble. The magnitude of the change dramatically increases relative to midcentury, with changes exceeding 2σ for the majority of convective basins, and extending fully poleward in both hemispheres. This reflects a shift toward more resistance to convection across the domain, which could have significant implications for convective frequency (Taszarek, Allen, Marchio, et al., 2021).
Parameters characterizing shear and helicity are historically very consistent with reanalysis ( Figure S1, Supporting Information S1), but their changes, in comparison to the thermodynamic variables, are small. Through its definition SRH03 is predominantly positive for the NH and negative for the SH, and thus regions that are blue in the southern hemisphere reflect relative increases. Changes in SRH03 are small by midcentury, with the majority of signal across seasons lacking ensemble consensus despite an overall indication of change. The exceptions to this are in the SH during the winter and spring, where latitudinal alternating bands of increases and decreases with increasing latitude suggest an enhancement of SRH over the severe convective regions of South America.
While not significant, indications of similar patterns are found over both southern Africa and Australia. By late century, changes to SRH03 are more clear, with increases over the southern United States in both MAM and to a lesser extent, DJF. Similar increases are seen for Europe, particularly in the west, coinciding with the current peak of severe convective season. While these changes are significant, model consensus suggests that this trend is not present for all models. In the Southern Hemisphere, the midcentury pattern is amplified, with a weakening of the subtropical band, and robust increases to helicity poleward of the tropic of Capricorn, extending to ∼40°S. These suggest greater favorability during the winter months. Similar patterns are present in SRH01, suggesting that these changes are also manifest at lower levels of the atmosphere, which may reflect that the low level wind field is the driving cause.
Changes to S06 reflect the expected changes to storm track and the midtropospheric flow, and are similar to helicity with latitudinal bands of alternated increase and decrease. The significance is already present by mid century (stippling), as is the agreement across models (hatching). The decrease over the CONUS is confirmed for three out of four seasons, except for the MAM months, which presents a more localized increase and decrease variability which is significant by the end of the century, but this signal is not robust across models. It is worth to remind the reader that the results on very elevated terrain (i.e., Himalayan regions) can be affected by the approximation used in the computation of S06 (see Section 2).
In summary, these results suggest that increasing thermodynamic favorability on a global basis will coincide with greater convection inhibition, while regional increases to low-level shear will be offset by decreased deep-layer shear. This suggests that consideration of covariability of convective environments is important to assess the potential changes to severe weather in response to a warming climate.

Temperature Scaling of Convective Environments
The variability of the rate of temperature increase across CMIP6 models can be large (Meehl et al., 2020), and differences in temperature response to increasing CO 2 emissions is expected to propagate in dependent variables such as, for example, specific humidity and CAPE (Singh et al., 2017). With respect to the models included in this analysis, Figure 2 shows in the left panel the significant differences across models in global surface temperature trajectories compared to each model's pre-industrial, 1850-1900, state, particularly for CanESM5. Moreover, a similar inter-model variability is seen in the other scenarios (not shown), such as SSP3 7.0 and SSP2 4.5, which are expected to reach specific temperature range increases at different points in the future (Figure 2, right panel).
To address these significant differences across models and scenarios, we investigate the possibility of presenting our results as a function of global temperature increase (ΔT) relative to preindustrial period (1850-1900), rather than time.
There are two potential approaches to this re-scaling as a function of temperature. We could use the global ensemble mean temperature increase (i.e., black line in Figure 2), or use for each model their respective ΔT values (colored lines in Figure 2). In the first approach, all models x-axis are stretched in the same way, therefore they remain synchronized in time, but the temperature differences are not eliminated; in the second case, instead, the alignment across years is lost but the differences due to difference in temperature trajectories are reduced (see Supporting Information S1 and Figure S3 for details).
Scaling the results of each model based on its temperature increase, rather than the ensemble mean, not only addresses models differences in temperature response to CO 2 across models, but it also reduces the uncertainty of the projections by an average 25% (see Supporting Information S1, and Figure S4 for details). For this reason, all subsequent results are therefore presented as a function of ΔT, with models scaled based on their own global temperature trajectory.
Regional averages are also easily re-scaled as a function of ΔT (as shown in Figure S3). For the spatial analysis, the historical component of the analysis (used to calculate the difference in mean and the t-test analysis) remains the same as in the previous case, equal to the average of the 1980-2010 time period. For the future scenarios, instead, we extract the year where ΔTs 1.5°C, 2°C, and 3°C are reached within each model (Table 2) and use these as the center of 25-year averaging time windows, to maintain a significant sample size but limit the range of temperature increase. Moreover, for the sake of consistency we do not combine experiments (historical with ssp585); which means that for some combinations of model and ΔT, the sample size is not exactly 25 years (MRI-ESM2-0 has a 23-year sample for ΔT = 1.5°C, CanESM5 a 21-year sample for ΔT = 2°C) or leads to their exclusion (CanESM5 reaches ΔT = 1.5°C increase in 2011, within the historical experiment, so it is not included in the ΔT = 1.5°C analysis).

Convective Environments in CMIP6 Models as a Function of Global ΔT
Global changes to environments as a function of temperature suggests similar results to those obtained by temporal changes (Figure 3). These suggest that increases to CAPE are more confined to continental regions or the tropical Pacific for increases of 1.5°C and 2°C respectively, while large increases over NH landmasses are prevalent in all ΔT scenarios. While changes are generally statistically significant by the 1.5°C scenario, robustness across the ensemble is generally not evident until a change of at least 2°C is reached within global average temperature. Increases of 3°C, see significant responses that are latitudinally dependent, with the largest relative increases occurring in more poleward climates, particularly through Canada, Eurasia throughout the year. Wholesale changes in CAPE at 3°C generally reflect temperature changes within the model ensemble by around 2060-2070. While CAPE changes are relatively consistent between the two approaches, changes in CIN highlight the role of the temperature differences in its contribution signal, particularly in the mid-century projections. Changes to CIN at the 1.5°C ΔT point suggest that while increases are still prevalent, the statistical confidence and agreement across the ensemble is lower. Increases in CIN are smaller across the convective regions with lower confidence of a robust change across the ensemble. At the 2°C ΔT, increases show greater robustness and consistency in the ensemble, while by the 3°C ΔT point, changes are of a similar magnitude and signature to projections by end century. These projections suggest that these significant increases are likely to be impactful to thunderstorm occurrence particularly as the increase in global mean temperature approaches 2°C above historical levels.
Changes in both SRH01 and SRH03 are equally temperature dependent, with many of the regions of the largest relative changes only becoming evident as ΔT approaches 2°C-3°C, with the majority of the largest changes confined to period following a ΔT of 3°C. This suggests that while significant changes are likely to SRH03 following significant warming, these changes may be manifest closer toward the mid-end century period. A similar pattern of reduced changes until larger ΔT changes are reached is found for S06. These suggest that the rate of change across different convectively relevant parameters are not consistent, and that the temperature relative approach gives a different insight into projected changes that avoids individual ensemble member biases.

Regional Trajectories of CMIP6 Convective Environments
The results in Sections 4 and 6 show extended and significant increases in CAPE and CIN for the majority of the domain, whereas shear related variables have a more regional and time/temperature dependent variability.
To explore how these changes manifest regionally, in Figure 4 we consider the regional area-weighted seasonal averages of those quantities for 11 regions (see Figure S2 for details) for land-only. In particular, differently from the previous figures, each model time series is normalized by its own 1980-2010 average, to express the results in terms of percentage change, to reduce underlying parameter bias in any given simulation. This allows for an easy comparison across seasons and regions, which can have very different underlying climatologies. For SRH01 and 0-3 km, we show the regional trend for their absolute values. During the remapping of the data onto ΔT coordinates, we bin and average the values on a 0.5°C degree regular grid.
We note that percent changes are calculated with respect to the historical averages, that is, the 1980-2010 climatology, whereas to compute ΔT we use the temperature increase relative to pre-industrial times . During the historical scenario, global temperatures have already increased faster than in previous decades and the ΔT corresponding to the models included here and the years covered by the historical scenario is about 0.8°C. CAPE and CIN ratios display a linear increase as a function of ΔT, with slopes varying regionally. The steepest rate of change is found at higher latitudes, for the Canadian provinces (Can region) and across Europe and Asia (EUA). This stronger increase is in line with the well-established stronger increase in temperature for higher latitudes (Chen et al., 2021;Singh et al., 2017) and is consistent with trends seen in recent observational shifts (Poujol et al., 2021). Over North America (US, CAN), seasonal CAPE trends are stronger during SON, MAM, and DJF, relative to JJA. For the US this reflects 25% relative increases for 1-2°C change, rising to a 50%-75% increase following a 4°C increase, although uncertainties are large at greater temperature changes. Over the EUA region, increases are the strongest, reflecting increases where near-zero or low climatological values lead to larger magnitudes of relative change. For Africa, South America, Europe, India, Southeast Asia and Australia increases for the highest temperature changes are smaller, generally between 25%-40%. Confidence intervals, except for few cases (e.g., Can and EUA) are very narrow, and they depict a clear and significant increase by a ΔT = 2°C warming. The regions in the southern hemisphere (SA, SA_1, Af_1, AU) show smaller slopes in general, yet  Figure 1 but for a temperature dependent framework. Standard deviation is still computed on the historical period , but changes refer to a 25-year sample centered on 1.5°C, 2°C, and 3°C global ΔT increase relative to preindustrial levels . Stippling indicates significance at the 0.05 level, hatching 100% sign agreement in the response across all models. still significant increases. Accompanying the positive trends in CAPE are increases in CIN, which are similar in relative magnitude to those found in CAPE. This suggests that while overall regional instability appears to be consistently increasing around the world, this is also accompanied by significant growth in the resistance to convection, consistent with the argument that while more instability will likely be available, frequency will not necessarily increase.
In contrast, all shear quantities over different atmospheric depths, S06, SRH01, and SRH03, display no trend, except for SON and DJF seasons in Af, SEA, In, and EUA regions, as seen in the spatial results. These competing factors support the use of severe weather indices in which the co-occurring CIN is included in the analysis together with instability (CAPE) and shear (S06 and/or SRH) parameters to assess whether instability or resistance to convection more influence changes in the convective environment.

Conditional Distribution of CMIP6 Convective Environments
Severe weather indices are based on the product of multiple convective parameters. To more clearly explore how changes to bivariate environments will respond, here, we explore their conditional distributions at 6-hourly resolution. The results, as expected, vary across seasons and regions, especially in terms of the range of values and seasonal amplitude, but their changes as a function of temperature are very similar. In particular, Figure 5 shows the difference of the multi-model average probability distribution function between a 25-year period centered around 2° ΔT and the historical  period for the CONUS and Australia. Fall and Winter seasons (SON and DJF for US, MAM, and JJA for AU) are very similar in the shape of the expected increases (red areas), and similarly for the Spring and Summer months where increases are largest. There are, as expected, differences in the overall range of values (in particular the AU winter, JJA, season has lower CAPE values compared to US winter), but the overall relationship across environments is not overly different.
The increase in CIN values is across a range of CAPE values, however, peak changes display an inverse relationship with CAPE, with the greatest density of CAPE increases corresponds to lower CIN. This suggests that although CIN will modulate CAPE-dependent proxies, its effect may be smaller than might be expected given the magnitude of increases to this parameter. In the second column, the effects of increasing CAPE are much larger . Seasonal and regional area-weighted response of convective available potential energy (CAPE), CIN, S06, SRH03, and SRH01 as a function of global temperature increase, ΔT, as a percentage of historical period . Vertical lines indicate 1.5°C, 2°C, and 3°C degree. Shaded areas represent the intra-model 95% confidence intervals.
than any decrease seen in S06 values. The largest S06 decrease (blue area) is mostly centered on the vertical axis, at around 25 m s −1 but generally low CAPE values. Increases to CAPE however are considerably more likely for lower values of S06 (less than 25 m s −1 ), suggesting that although the overall frequency of favorable environments may increase, the likelihood that this corresponds to an increasing frequency of highly organized storms is less clear. This may imply a shift toward more frequent multicellular storm modes as opposed to supercells, with implications for the potential severe hazards. CAPE and SRH01 and SRH03 phase spaces display similar characteristics. The last two columns explore relationships for to EHI 2 and SEV, and hence represent a filtered subset of the second and third/fourth column. In the last two columns the largest signal is once again the increase in CAPE, visible especially for the CAPE 2 , SRH 2 2 case with a clear vertical red band. For the CAPE SEV , S06 SEV case, instead, we see a narrow red band propagating back toward low values of CAPE and high shear, although not robust across models. What is interesting however is the decrease in high-shear low-CAPE environments that are often found over the southeastern US and other thermodynamically constrained regions. These findings suggest that while increases to instability will lead to enhanced likelihood of severe convection, this may be at the expense of high-shear low CAPE environments that provide the majority of severe weather in some parts of the world. These results also indicate that the small decreases in shear and helicity do not play a large role in curbing the occurrence of severe weather proxies.
Overall, this type of analysis shows how the marginal increases and decreases seen in the previous figures are not able to tell the whole story with respect to changes in severe convective storm environments. The co-occurrence of two or more environmental parameters, and their magnitude, is in fact crucial in defining not only the possible occurrence of a severe weather event, but also its type and intensity . In future work, it will be interesting to quantify more details of the subtle shifts seen in Figure 5, and explore how these can affect the type and/or intensity of severe weather events.

Projected Changes in Severe Weather Indices
Severe weather proxies have various limitations. However, they represent the best available approximation for modeling severe weather activity in models that lack the capacity to resolve convective processes. To increase the robustness of our findings, we consider: (a) a variety of severe weather proxies that use different combinations of variables filtered at varying thresholds to reflect the lack of a universal proxy that works in all parts of the world, (b) results are aggregated to a coarser (regional) spatial and temporal (seasonal) resolution, and (c) to avoid climatological bias changes are quantified as percent change relative to the historical period. To encompass the wide scope of environments, we consider indices that allow the exploration of broad definitions for severe weather that have been widely applied in literature (e.g., SEV; Brooks et al., 2003;Allen et al., 2011;Blamey et al., 2017;Taszarek et al., 2018) and extend to strict criteria that are regionally specific (SCP). This is particularly important in areas where the underlying link between the presence of convective environments and severe weather activity is less clear. By aggregating regionally, we acknowledge that GCMs have limitations in that they represent only a possible realization of future climate and should not be looked at a temporal scale relevant to weather but rather on climatological mean scales. Similarly, finer spatial features can increase uncertainty in our findings compared to a regional aggregated analysis. Finally, we also acknowledge that the majority of these indices have been calibrated only for specific regions of the global domain (e.g., U.S., Europe, and Australia), and thus exploring their absolute magnitude or specified thresholds would be problematic, but quantifying the projected percentage changes compared to the historical baseline to some extent eliminate these potential underlying biases (Allen et al., 2014a). These three steps together with our temperature dependent framework, allow us to draw robust and significant conclusions related to the changes of severe weather indices. Figure 6 shows the multimodel (seven) and multiindex (seven) seasonal and regional variability of the projected changes for severe weather frequency as a function of global temperature increase from pre-industrial times , with the percent changes calculated with respect to the historical  period. Each trajectory is the average over 49 trajectories (seven models and seven indices), and the shaded area includes both intra-model and intra-index variability, with the former being the largest source of uncertainty. Because of the overall linear behavior of these trends, we fit a linear regression to the 0.8°C-3°C interval (the average ΔT value for the historical period  is ∼0.8°C for the seven models included in the analysis, see Figure 2), and compute the percentage change as a function of global temperature increase for 1.5°C, 2°C, and 3°C ΔT. To provide this insight in a more accessible form, the slope of the regression as percentage increase per unit °C increase is also calculated.
Examining the temperature dependent results, we find that increases are essentially linear (with few exceptions) for values of ΔT between 0.8°C and 3°C. Beyond 3°C the majority of trends tend to plateau or show more erratic responses in response to greater uncertainty at higher ΔTs as data points decrease (see Figure 2). It is also plausible that part of the change could be in response to process-based shifts that only occur at higher temperature increases. This potential response is intriguing warranting further analysis, but lies outside the scope of this paper. There are also strong regional differences in the slopes and increasing uncertainty with increasing temperature, together with significantly larger changes with increasing latitude, with the majority of the uncertainty due to intra-model variability.
Increases during the summer months, for the majority of domains, present the smallest overall changes to relative frequency, while increases are more pronounced during MAM and SON in both hemispheres, generally on the order of 10%-20%. As expected, closer to the equator changes are more uniform across the seasons. The regions with steepest increases, but also greater uncertainty, are located closer to the Polar regions, and particularly in the northern hemisphere. This reflects the largest temperature driven increases to instability (CAPE) particularly over Can and EUA (Figure 4). The magnitude of the projected trends in these two regions is very large, with up to a 50% increase in frequency per degree of warming. It is worth to note, however, that in these two regions, and particularly for DJF and MAM, the uncertainty is very large because of the very low or no occurrence of these proxies in the current climate (Etkin & Brun, 1999), and this drives exclusion of these regions from consideration for DJF. For the same areas, during JJA and SON, when climatological frequency is higher, increases are still high and have increased confidence, similar to other regions.
Another region with anomalous behavior is the EU. This partly stems from the smaller domain size relative to other regions, together with features that are commonly problematic at the coarse scale of GCMs, such as quickly varying orography and extensive land-ocean boundaries (Púčik et al., 2017). For example, the ocean masks of some models barely identify some regions of activity of severe convection as occurring over land (i.e., Italy). For EU, the seasons with robust linear increases in percentage with temperature are the cooler seasons, whereas spring and summer present essentially no significant linear trend, for MAM, or only slightly greater than zero, for JJA.
For the remaining regions, we find similar slopes and uncertainty in SA, SA 1 , Af, Af 1 , and AU, with slightly lower values for AF 1 region. The proximal warm oceans see In and SEA display similar behavior, especially given the comparable range of latitudes. To better contextualize the results, regional changes and their 95% confidence intervals are summarized in Table 3.
Finally, similarly to Figure 6, we combine the 49 model and index combinations to extract a spatial distribution of the projected percentage change as a function of global temperature increase (Figure 7). Each (model, index) combination is normalized by its own historical distribution, as in Figures 4 and 6, and then the multi-model and multiindex averages and standard deviations are calculated.
The spatial distribution, as expected, displays large, robust, and significant increases in areas identified as strongly increasing in the regional aggregates, but also known to be associated with active severe weather seasons in the current climate (Brooks et al., 2003). The highest values are present in the northern hemisphere and the higher latitudes. However, as indicated by the limited stippling, this is a result of the very rare rate of occurrence in the current climate. A similar influence is seen over the Sahel desert area, where the relative percentage increase is very large, but the current climatology is very low. In general, all the areas with very large percentage increase Figure 6. Regional and seasonal multimodel and multiindex response of severe weather activity to global temperature increase, ΔT. Vertical lines indicate 1°C, 1.5°C, 2°C, and 3°C degree, text located on the lines indicate the % change projected for that specific temperature increase. Shaded areas represent the intramodel and intraindex 95% confidence intervals. Linear fit (solid straight line) and the 95% confidence interval (dashed lines) are fit to a 0.8°C-3°C interval. Bold text on the top left of each panel represents the percentage change per unit degree increase in global temperature based on the linear fit.
(≫200%) correspond to areas with comparatively low occurrence in the current climate and wide confidence interval with greater uncertainty in the regional results. Overall, the spatial results shown in Figure 7 suggest that further studies are needed to understand the source of disagreement across models with respect to the processes relevant to severe convective storm environments, which results in very few hatched areas (Seeley & Romps, 2015).
There are significant and robust increases over North America for all temperature levels (1.5°C-3°C). Of particular note is the increase over the southern US for DJF, where a secondary peak in severe weather activity is typical for the cooler months (Lu et al., 2015). For the other seasons, the significance follows the known climatological spatial distribution (MAM for Southern US, and moving north for JJA to the Canadian provinces for SON). The signal becomes more significant, stronger, and more spatially robust as the temperature increases (1.5°C-3°C). Moving to the European continent, the only significant increase, and robust at higher temperature changes, are found in March-November. This differs somewhat from Figure 6, where the DJF months showed significant slopes. However, the large confidence intervals in that figure are reflected in the low confidence in the spatial results relative to this region. For Africa, as seen in the regional results, the significant results are widespread in the northern region (Af) for March-November, reflecting increases both along the Mediterranean coast and through central Africa. Over the southern parts of Africa, there are also significant increases, including over South Africa, where severe weather is commonly observed (Dyson et al., 2021).
Over South America, significant and robust increases appear evident in the higher temperature ranges, particularly over the La Plata basin and northern Argentina proximal to the Andes (2°C-3°C). These increases coincide with the September-March seasonal peak in frequent severe weather seen in this region (Beal et al., 2020;Martins et al., 2017). For the Asian continent, the months with widespread increments are found through the spring and summer across the temperature range, particularly over India, Bangladesh and eastern China. For Australia, the significant increases are limited to higher temperature ranges (2°C-3°C), in line with the regional results, that display more of a convex dependence on temperature, rather than linear, suggesting a delay in increases for this region compared to others as a function of global temperature increase.

Conclusions
Seven models from the CMIP6 archive are used to explore the global response of severe convective storm environments to a warming climate. Since models do not resolve severe convective storms, we computed a set of convective environments (CAPE, CIN, SRH, S06) at a 6-hourly temporal resolution. We first explored them individually, and then combined them in severe weather proxies to explore the competing physical processes that are known to be associated with severe weather occurrence.
The historical simulations are in line with current observational studies ( Figure S6) where changes in convective environments are small compared to variability and not always detected, except in northern latitudes where changes appear more extreme and easier to observe. As projections move into warmer scenarios, however, trends in convective environments emerge above the variability, and there is an increasing likelihood of changes to severe convection.
As the climate warms, projections show increases in severe convection that are clearly latitudinally and seasonally dependent, with increases generally tied to the respective seasonal peaks. Several regions show little to no detectable change until temperature changes exceed 2°C. The largest increases expressed as relative percentages are found over the higher latitudes of the Northern Hemisphere, where current climatological frequency is low. Nonetheless, except for the higher latitudes, projections suggest global increases on the order of 5%-20% per 1°C of global temperature increase. This rising frequency is driven predominantly by the substantial increase in instability (CAPE), which even when moderated by increasing inhibition (CIN) still drives the projected increase in severe weather proxies. The robustness of the results suggests that there is relatively high certainty that atmospheric instability will be more favorable to severe convection as a result of increasing temperature due to CO 2 emissions. The robustness and significance of these projections vary across regions and seasons, but no matter the continent or region, there are increases to intensity and frequency of favorable environments and indices respectively.
While these results suggest that severe convection is likely to be more frequent in a future climate in line with existing research (Chen et al., 2021;Diffenbaugh et al., 2013;Glazer et al., 2020;Singh et al., 2017;Trapp & Hoogewind, 2016;Trapp et al., 2019), we mentioned throughout the manuscript how the approach used herein has some limitations: severe weather proxies are not the same as severe weather events (Hoogewind et al., 2017). The occurrence of a perfect combination of convective environments, in fact, does not always translate to the occurrence of a severe weather event, and the "false alarm ratio" of these proxies can vary across the globe, with further work needed to explore the realized frequency. However, these limitations were ameliorated by multiple steps, such as normalizing each (model, index) combination by its own historical state, by exploring a multi-index multimodel combination, and by aggregating the results at the regional and seasonal scale. Finally, these indices are the most used tools to explore and model severe weather, from real time to seasonal forecasting, and given the current computational limits, they essentially represent the only possibility for exploring the effect of a warming climate on severe weather in climate models, which do not resolve convective processes.  There are many sources of uncertainty, with a variety of origins. Some are specific to the applications here, such as the choice of severe weather proxies or the extent to which the proxies capture the characteristics of actual storms. Others are common across CMIP applications, such as limitations in modeling schemes or/and temporal and spatial resolution. A common source of uncertainty across CMIP applications is that the climate sensitivity still remains somewhat uncertain, even for a specified scenario, and which scenario will match reality is still unknown. Convective environments, and the derived indices, essentially inherit this sensitivity to uncertainty in temperature changes, displaying at times large variability across models. An interesting possibility to reduce uncertainties in climate change related studies can be achieved by partitioning uncertainty and creating more credible conditioned results (Shepherd et al., 2018). We, in fact, applied this paradigm shift to our results, and by scaling changes in convective parameters and indices by temperature, we provide change estimates that are conditional on global mean temperature changes (Tebaldi & Arblaster, 2014). This approach not only reduces the uncertainty due to differences in climate sensitivity across models, but it also allows for conclusions that are not tied to a specific scenario.
While the insights gleaned using this approach are essential toward understanding the global response of convective environments to a warming climate, future work will need to focus on both regional and temporal variability of these events, to better understand these changes beyond the overall projected trajectory. Regional approaches that statistically and dynamically downscale the projections from the CMIP6 ensemble will also provide critical insights as to the realized fraction of change in frequency (Allen, 2018;Hoogewind et al., 2017;Raupach et al., 2021). It is worth to mention, in fact, that across the global domain similar percentages of increase in severe storm convective environments can translate to very different local absolute changes in occurrences, depending on the current regional and seasonal climatology of such events.
Moving forward, an ongoing challenge is the volume of the CMIP6 data set, and the best practices to effectively interrogate these big data. The authors note that none of the work presented here would have been possible without the Pangeo software ecosystem and developments in cloud computing that allowed the analysis to be carried out data-proximate. Further advances in these analysis techniques will be necessary to analyze data in the volume and scales needed to provide insight to changes in severe convective frequency.