Effects of Transient Stellar Emissions on Planetary Climates of Tidally Locked Exo-Earths

Howard Chen; Paolo De Luca; Assaf Hochman; Thaddeus D. Komacek

doi:10.3847/1538-3881/add33e

1. Introduction

Space weather events originating from a star—including X-ray and extreme UV radiation, coronal mass ejections, and stellar energetic particle precipitation—will likely impede or disturb the atmospheres of rocky exoplanets (see, e.g., C. Dong et al. 2017; K. France et al. 2020). Analyses of the Kepler data show that the average flare energy detected is ∼10³⁵ erg (J. R. A. Davenport 2016), and the average flare amplitude does not change substantially between G-K-M stars. In contrast, flare frequency strongly depends on spectral type (T. Van Doorsselaere et al. 2017). From the Transiting Exoplanet Satellite Survey (TESS) catalog, it has been shown that 97% of the most intense flares released energies of >10³³ erg (K. Vida et al. 2021; W. S. Howard & M. A. MacGregor 2022; B. Stelzer et al. 2022), with a typical TESS star having daily flares with optical energies between 10³² and 10³³ erg (J. N. Ealy et al. 2024). Meanwhile, the level of magnetic activity in terms of the energy and frequency of super-flare events with energies between 10³³ and 10³⁵ erg drops down within 600 Myr for G-type stars and 1 Gyr for K dwarfs and M dwarfs (J. R. A. Davenport et al. 2019; E. Ilin et al. 2021). The latter hosts include the well-characterized TRAPPIST-1 system, which, due to their deep convective zones and longer spin-down timescales, are known to retain a high level of magnetic and chromospheric activity for Gyr (e.g., I. Baraffe et al. 2015). This can affect atmospheric environments of close-in exoplanets on long timescales, inducing water loss via photolysis and hydrogen escape (e.g., A. J. Louca et al. 2023; E. F. Fromont et al. 2024; L. N. R. do Amaral et al. 2025). However, the detailed flare characteristics of the youngest of the samples (<300 Myr) are still being investigated (see, e.g., A. D. Feinstein et al. 2024; W. S. Howard et al. 2025).

Stellar eruptive events have been shown to dramatically influence the abundance of major and minor chemical species in the atmospheres of habitable zone planets, particularly for those around M dwarfs. The most significant effects on Earth-like, or strongly oxygenated, planets are those on the ozone column abundances (H. Chen et al. 2021; R. J. Ridgway et al. 2023). For a single large AD Leonis-like flare, the impact of ozone erosion is minimal as the ozone layer is recovered over just a few Earth years (A. Segura et al. 2010). However, the ozone layer is rapidly eroded even with moderate flare energy and proton fluence assumptions (M. A. Tilley et al. 2019) for repeated events with frequencies typical for Kepler stars. . The effects on other photochemical constituents, for example, nitric acid or HNO₃ (e.g., K. Herbst et al. 2019), depending on the assumed energetics of the associated primary protons precipitated in the atmosphere of the planet, typically between 16 MeV and 0.5 TeV (M. Scheucher et al. 2020).

Others, using a combination of observational reconstruction and numerical modeling, have found that secular injection of energetic particles can lead to enhancements of greenhouse gases such as water vapor and nitrogen oxides, i.e., NO₂ and N₂O, during large stellar proton events (V. Airapetian et al. 2016; H. Chen et al. 2021). Such enhancement may have important implications for surface climate, for example, by causing abrupt surface warming that could alter previous conclusions regarding surface habitability. As implied by other 1D models, sudden and periodic enhancements in the abundance of other greenhouse gases or UV absorbers by stellar flares can be substantial (up to 4–5 orders of magnitude change in their mixing ratios), but further efforts to quantify these processes using modeling tools at the global (three-dimensional; 3D) level are needed. In CO₂-rich conditions, an increase in upper atmospheric CO₂ caused by repeated flaring in oxidizing atmospheres can dramatically sculpt the temperatures of upper atmospheres around M dwarfs through CO₂ radiative cooling (R. Wordsworth & R. Pierrehumbert 2014). In these environments, repeated flaring may initiate an entirely different chemical reaction chain that may have substantial implications for surface climate. Ammonia, for example, which has been seen as a potential greenhouse gas to warm the surfaces of early Earth and potentially Earth-like planets (C. Sagan & G. Mullen 1972), would have severely limited abundances due to the absence of strong UV absorbers between 200 and 300 nm that would prevent photodestruction of ammonia in the atmosphere.

These 1D investigations show, at least on a theoretical level,⁸ that photochemical effects as a result of time-dependent stellar flares cannot be neglected in efforts to more realistically characterize the composition of astrobiological significant trace gases in rocky exoplanet atmospheres. It is, therefore, essential to further examine their impacts on habitability and exoplanet temporal variability and climate. Planetary climate dynamics can be an important factor in determining the observable properties of exoplanets because they set the three-dimensional atmospheric temperature structure, mixing of chemical species and aerosols, and resulting cloud distributions. An understanding of planetary climate dynamics is needed to interpret spectroscopic observations of exoplanet atmospheres. The precise characteristics of an exoplanet's atmosphere are linked to its orbital and rotation period around an M-dwarf star. Due to the dependence of the rotation period on the stellar spectral type, different rotational states lead to diverse climate evolution pathways of their surface temperatures (e.g., T. M. Merlis & T. Schneider 2010; J. Yang et al. 2013; R. Kopparapu et al. 2016, 2017; M. Turbet et al. 2016; M. J. Way et al. 2016; L. Carone et al. 2018; T. D. Komacek & D. S. Abbot 2019; M. J. Way & A. D. Del Genio 2020; F. He et al. 2022; H. Chen et al. 2023; T. Hammond et al. 2025), as well as divergent chemical distributions (e.g., H. Chen et al. 2018; M. Cohen et al. 2024; P. De Luca et al. 2024; M. Braam et al. 2025).

There has yet to be a systematic study of the effect of stellar flaring events on the planetary climate and atmospheric dynamics of planets orbiting M-dwarf stars, as the majority of previous work has focused on understanding effects on stratospheric chemistry using 1D photochemistry models. Exoplanet climate forcings due to external agents such as astrophysical events, a relatively unexplored area of research, could have major ramifications for both habitability and observations. Studies on the relationship between Earth and the Sun during strong geomagnetic storms found a ∼0.5 K cooling in the middle stratosphere over the tropics and up to 2 K over southern high latitudes (E. Rozanov et al. 2005). The most recent investigations into the effects of extreme solar energetic events underscored the key role of the geomagnetic field strength in moderating instantaneous and long-term effects on the atmosphere (P. Arsenović et al. 2024). Specifically, during historical periods of weakened planetary magnetic fields, their global climate models have shown that odd nitrogen species can remain in the atmosphere on the order of several years in the lower latitudes, leading to ozone losses on a decadal timescale and hence markedly increased ground-level UV radiation.

Another study on the effects of the AD 775 great solar event on Earth’s climate has shown ozone depletion followed by a decrease in near-surface air temperature by up to 4 K (T. Sukhodolov et al. 2017). This temperature anomaly then affects the Polar Night Jet Oscillation and accelerates the tropospheric mean flow soon after the event in December. Their model also found that a more pronounced meridional circulation leads to warming in the Earth’s polar regions and the northern hemisphere. The most recent investigations into the effects of extreme solar energetic events on Earth’s particle and atmospheric environment underscored the key role of the geomagnetic field strength in moderating instantaneous and long-term effects on the atmosphere (P. Arsenović et al. 2024), echoing earlier results on strongly to weakly magnetized exoplanets (H. Chen et al. 2021). These terrestrial studies indicate that climate effects may also be important for extrasolar environments, in particular for activities from cool stars extending up to higher energy domains compared to solar measurements.

Here, we simulate the putative climate dynamics of TRAPPIST-1e analogs, assuming an Earth-like atmosphere with different levels of flaring activity based on observational measurements of M dwarfs. We analyze the flare-driven variability in the thermal structure and dynamics of the atmosphere and their correlation with the variability in the abundance and distribution of key chemical constituents. This manuscript is organized as follows. We describe our GCM modeling framework and statistical analyses in Section 2. We present results statistically quantifying the impact of varying levels of flaring activity on climate and atmospheric chemistry in Section 3. Finally, Section 4 discusses our findings and places this work in the broader context of climate and atmospheric evolution models.

2. Methods and Numerical Model

2.1. 3D Global Chemistry-climate System Modeling

We use the National Center for Atmospheric Research's Community Earth System Model (CESM) version 1.2.2 (R. B. Neale et al. 2010; D. R. Marsh et al. 2013) to simulate the atmospheres of template exo-Earths under the influence of stellar eruptive events, including both the effects of stellar UV flares and stellar energetic particle precipitation. The Whole Atmosphere Community Climate Model (WACCM), a supercomponent set of CESM, is a 3D global climate model that simulates atmospheric chemistry, radiation, thermodynamics, and dynamics. In particular, the model includes chemical speciation, network reaction, and transport processes necessary to fully account for the effects of time-dependent stellar events that are of interest here. We follow previous uses of WACCM in exoplanet studies (e.g., H. Chen et al. 2019; G. P. Afentakis et al. 2023; Y. Luo et al. 2023) and specify the Community Atmosphere Model v4 (CAM4) as the atmosphere component of WACCM. CAM4 employs the native Community Atmospheric Model Radiative Transfer radiation scheme (J. Kiehl & V. Ramanathan 1983), the Hack scheme for shallow convection (J. J. Hack 1994), the Zhang-McFarlane scheme for deep convection (G. J. Zhang & N. A. McFarlane 1995), and the Rasch-Kristjansson scheme for condensation, evaporation, and precipitation (M. Zhang et al. 2003). Note that the specific atmospheric component and the convection parameterization scheme have been shown to substantially influence the simulated climates of planets around M dwarfs (J. Bin et al. 2018; D. E. Sergeev et al. 2020).

We conduct four independent simulations initialized from the same steady state solution. These simulations have TRAPPIST-1e's planetary mass and radius (E. Agol et al. 2021), Earth-similar compositions, Earth’s present-day insolation, but with slightly different transient stellar emissions; that is, time-varying UV spectra and input proton fluence. We also assumed modern Earth’s surface albedo, surface topography, and landmass distribution to facilitate comparison with previous work studying ocean-covered (exo)planets. Previous exoplanet applications of WACCM simulations have also assumed Earth’s surface properties in part due to the numerical challenges of replacing ocean bathymetries of aquaplanets with the existing model land masks (see, e.g., H. Chen et al. 2019; Y. Luo et al. 2023).

We use an open-source Python flare toolkit developed by R. O. P. Loyd et al. (2018) to compute flare spectra and flare light curves. This stochastic flare model is based on the MUSCLES Treasury Survey V, which we use to generate UV light curves. The MUSCLES survey (HST observing program 13650) characterized the low-mass stellar radiation environment, including X-ray, XUV, FUV, and NUV fluxes (K. France et al. 2016). Model flares are stochastically generated, depending on the range of user-specified energies and frequencies, from observed data drawn from two stellar populations: the MUSCLE M-dwarf sample and four active stars Proxima Centauri, AD Leo, EV Lac, and AU Mic. The UV flare spectra are obtained by multiplying the quiescent UV spectra of each stellar SED with the active-to-quiescent flux ratio given by the flare generator. However, we note that energetic particles, or stellar proton events, are found to be responsible for ∼99% of the photochemical effects (A. Segura et al. 2010; J. L. Grenfell et al. 2012; M. A. Tilley et al. 2019). As the estimated color temperatures of M-dwarf flares range between 7700 and 14,000 K (A. F. Kowalski et al. 2013), here all of our modeled flaring events assume blackbody temperatures of 9000 K. M-dwarf flares at much hotter temperatures have been observed (e.g., R. A. Osten et al. 2016).

To compute the associated energetic particle precipitation in the ambient atmosphere, we follow previous studies using solar scaling of an event’s proton fluence based on near-Earth satellite data, assuming a linear relationship between the proton flux and production rate of ion pairs that should be valid for the most Earth-like atmospheres (C. Jackman et al. 2008). We assume that all of the particles are protons. We compute the expected peak proton fluences from the SiIV energy of stellar flares using the prescription of A. Youngblood et al. (2017). The proton energy spectra for each event are assumed to be fixed and consistent with the 2003 Halloween solar proton event (B. Funke et al. 2011). Although inappropriate for the most extreme events and M-dwarf activity (K. Herbst et al. 2019; J. Hu et al. 2022), this assumption should not deviate from the typical events expected from our chosen flare energy range. Across the four numerical experiments conducted, we assume three different levels of flaring activity of TRAPPIST-1. The template “quiescent” simulation uses a TRAPPIST-1 spectra with its UV regime (λ < 125 nm) replaced with an inactive PHOENIX BT-Settl UV model. The “moderate” simulation uses the out-of-the-box flare frequency for sample MUSCLES stars (α = 0.76; K. France et al. 2016; R. O. P. Loyd et al. 2018). The “extreme” simulations assume a cumulative flare index α = 0.58, indicative of more magnetically active stars (K. France et al. 2016; R. O. P. Loyd et al. 2018). Both these simulations have flare rates μ = 6 day⁻¹. Lastly, the “active” used a power-law index of α = 0.52 and a reduced flare rate of μ = 4 day⁻¹. This choice in flare statistic is in closer agreement with the sample of actively flaring stars (Table 6 in R. O. P. Loyd et al. 2018). The latter simulation set aims to test the effect of a slightly higher frequency in sampling more energetic flares, but a lower flare rate overall across all energies. For the three nonqueiscent simulations, we restrict their input flare energies to between 10³⁰ and 10^34.5 erg, and all model data are outputted at daily cadences equal to one Earth day. Table 1 summarizes the stellar and planetary parameters in this study.

Table 1. Ranges of Initial Stellar and Planetary Properties Explored in the Numerical Simulations

Stellar or Planetary Parameter	Value/Range (units)
Planetary orbital period	4.62 (days)
Stellar T_eff	2600 (K)
Incident flux	900 (W m⁻²)
Planet radius	5733 (km)
Proton energy spectra	2003 Oct–Nov SPE¹
Cumulative flare index	0.52, 0.58, 0.76²
Flare rate	4–6 day⁻¹
Flare energy	10³⁰–10^34.5 erg
Flare blackbody temperature	9000 K
Blackbody SiIV energy ratio	160
SiIV quiescent flux	0.1 erg s⁻¹ cm⁻²

Note. (2) The “moderate” experiment refers to those with the cumulative flare index set to α = 0.76; the quiescent refers to those with the variability of the stellar spectral energy distribution turned off, and the “extreme” is with a lower flare index α = 0.58. Lastly, the “active” case uses the lowest cumulative flare index of 0.52, with a slightly lower flare rate of flare rate of 4 day⁻¹. These estimates of their flare frequency and flare indices are based on the HST MUSCLES survey: (K. France et al. 2016; R. P. Loyd et al. 2018). (1) M. López-Puertas et al. (2005), C. Jackman et al. (2008).

Download table as: ASCII Typeset image

J. A. Jackman et al. (2023) recently showed that the assumed blackbody flare temperature value (i.e., 9000 K⁹ ) underestimates the flare energies of field-age M dwarfs by a factor of 6.5 in the NUV and 30 in the FUV range. Therefore caution must be placed when interpreting our numerical results, as higher flare temperatures would impact photodissociation rates and enhance/mute the impacts of stellar protons. For example, the NUV regime is largely responsible for ozone photolysis, whereas the FUV contains key bands for ozone production. The precise spectral shape of each flare determines the production and loss of cooling/warming agents such as NO and N₂O. How different flare temperatures might change spectral windows and how UV irradiation intersects with the effects of ionizing particles remains to be seen. Exploring alternative blackbody temperatures and an expanded parameter space covering flare frequency distributions of earlier type stars, is an important next step for both 1D photochemistry and 3D climate models.

2.2. Statistical Analyses

We test the statistical significance of the field medians by applying a two-tailed Wilcoxon rank-sum test (H. B. Mann & D. R. Whitney 1947). The test does not assume data to be normally distributed, and it checks the null hypothesis that the field median of the extreme and moderate simulation is equal to the field median of the quiescent simulation at the 1% significance level (p value < 0.01).

For the above four sets of climate model simulations, we calculate boxplots of the anomalies for regions as follows: equatorial substellar (−15^∘–15^∘N, 165^∘–195^∘E), equatorial antistellar (−15^∘–15^∘N, 345^∘–15^∘E), midlatitude substellar (30^∘–60^∘N, 165^∘–195^∘E), and midlatitude antistellar (30^∘–60^∘N, 345^∘–15^∘E). We select these subregions because they divide the dayside and nightside hemispheres, including the terminator regions, better to identify the effects of tidal locking on flare activity. The anomalies were computed concerning the quiescent simulation. Specifically, we first take the field medians of the extreme, moderate, and active simulations, and then subtract from them the field and temporal medians of the quiescent simulation. In this case, we apply a Kruskal–Wallis H test (or one-way ANOVA on ranks; W. H. Kruskal & W. A. Wallis 1952), which is an extension of the Wilcoxon rank-sum test to three groups, to check whether at least one population anomaly median of one group is different from the population median of at least one other group at the 1% significance level.

Along with field medians and boxplots, we also compute monthly 8 yr median time series for all four simulated cases over the entire planet, dayside (−90^∘–90^∘N, 90^∘–270^∘E) and nightside (−90^∘–90^∘N, 270^∘–90^∘E). Lastly, we calculate the vertical difference profile of wind speed (m s⁻¹) over the region −15^∘–15^∘N, 0^∘–50^∘E. This region was chosen as it showed one of the largest anomalies in the field median analyses. As for the other cases, we apply a two-tailed Wilcoxon rank-sum test to infer whether the temporal and field medians are statistically different at the 1% significance level.

Table 2. Summary of Simulated Results for the Quiescent, Moderate, and Extreme Stellar States

Value	Atm. Level	Atm. Variable	Case Name
	(mbar)
654.88 K	1.223 × 10⁻⁵	Temperature	Moderate
204.87 K	0.008	Temperature	Moderate
62.39 m s⁻¹	1.223 × 10⁻⁵	Wind Speed	Moderate
28.35 m s⁻¹	0.008	Wind Speed	Moderate
4.10 × 10⁻¹⁰	1.223 × 10⁻⁵	OH vmr	Moderate
1.18 × 10⁻⁸	0.008	OH vmr	Moderate
1.83 × 10⁻⁶	1.223 × 10⁻⁵	H₂O vmr	Moderate
1.64 × 10⁻⁵	0.008	H₂O vmr	Moderate
3.04 × 10⁻⁷	1.223 × 10⁻⁵	N₂O vmr	Moderate
2.55 × 10⁻⁶	0.008	N₂O vmr	Moderate
1.09 × 10⁻¹³	1.223 × 10⁻⁵	NO₂ vmr	Moderate
1.90 × 10⁻⁷	0.008	NO₂ vmr	Moderate
2.1 × 10⁻⁴	1.223 × 10⁻⁵	NO vmr	Moderate
3.18 × 10⁻⁶	0.008	NO vmr	Moderate
5.14 × 10⁻¹⁰	1.223 × 10⁻⁵	O₃ vmr	Moderate
2.59 × 10⁻⁸	0.008	O₃ vmr	Moderate
761.76 K	1.223 × 10⁻⁵	Temperature	Quiescent
210.25 K	0.008	Temperature	Quiescent
44.39 m s⁻¹	1.223 × 10⁻⁵	Wind Speed	Quiescent
24.97 m s⁻¹	0.008	Wind Speed	Quiescent
2.25 × 10⁻¹⁰	1.223 × 10⁻⁵	OH vmr	Quiescent
2.10 × 10⁻⁹	0.008	OH vmr	Quiescent
1.95 × 10⁻⁶	1.223 × 10⁻⁵	H₂O vmr	Quiescent
6.85 × 10⁻⁸	0.008	H₂O vmr	Quiescent
1.55 × 10⁻⁹	1.223 × 10⁻⁵	N₂O vmr	Quiescent
1.81 × 10⁻⁰⁸	0.008	N₂O vmr	Quiescent
5.34 × 10⁻¹⁴	1.223 × 10⁻⁵	NO₂ vmr	Quiescent
1.34 × 10⁻¹⁰	0.008	NO₂ vmr	Quiescent
1.36 × 10⁻⁴	1.223 × 10⁻⁵	NO vmr	Quiescent
1.29 × 10⁻⁷	0.008	NO vmr	Quiescent
4.28 × 10⁻¹⁰	1.223 × 10⁻⁵	O₃ vmr	Quiescent
3.09 × 10⁻⁷	0.008	O₃ vmr	Quiescent
603.99 K	1.223 × 10⁻⁵	Temperature	Extreme
240.58 K	0.008	Temperature	Extreme
62.61 m s⁻¹	1.223 × 10⁻⁵	Wind Speed	Extreme
39.92 m s⁻¹	0.008	Wind Speed	Extreme
1.43 × 10⁻⁹	1.223 × 10⁻⁵	OH vmr	Extreme
1.34 × 10⁻⁷	0.008	OH vmr	Extreme
1.87 × 10⁻⁶	1.223 × 10⁻⁵	H₂O vmr	Extreme
7.9 × 10⁻⁴	0.008	H₂O vmr	Extreme
5.71 × 10⁻⁶	1.223 × 10⁻⁵	N₂O vmr	Extreme
2.03 × 10⁻⁵	0.008	N₂O vmr	Extreme
2.02 × 10⁻¹³	1.223 × 10⁻⁵	NO₂ vmr	Extreme
4.34 × 10⁻⁶	0.008	NO₂ vmr	Extreme
3.81 × 10⁻⁴	1.223 × 10⁻⁵	NO vmr	Extreme
2.95 × 10⁻⁵	0.008	NO vmr	Extreme
4.44 × 10⁻¹⁰	1.223 × 10⁻⁵	O₃ vmr	Extreme
3.77 × 10⁻⁹	0.008	O₃ vmr	Extreme

Note. Values recorded at pressures 1.223 × 10⁻⁵ mbar correspond to the lower thermosphere, while those recorded at 0.008 mbar correspond to the mesosphere. On Earth-like planets and Sun-like stars, the latter region is where energetic particle precipitation is expected to be particularly important (C. Jackman et al. 2008). The labeling of each input flare characteristic is described in Table 1. vmr = volume mixing ratio.

Download table as: ASCII Typeset image

3. Results

First, we present the temporal medians of our four simulated cases. We then examine regional anomalies between different hemispheres and vertical domains before delving into the temporal evolution of atmospheric temperature and relevant chemical species. We conclude the results by presenting possible regimes of flare-induced dynamic variability in the atmosphere. Table 2 shows selected results, including atmospheric temperature, wind speed, and mixing ratios of gas phase constituents (OH, H₂O, N₂O, NO₂, and O₃). We focus on the variables most relevant to climate dynamics in the the main figures, whereas the Appendix provides information on the photochemically important species.

3.1. Temporal Medians

Temporal medians of temperature for the upper and middle atmosphere in the “extreme,” “active”, “moderate,” and “quiescent” simulations are shown in Figure 1. With the inclusion of the effects of UV flares and energetic particles, the moderate simulation leads to 100 K of cooling in the upper atmosphere over eight years (Figure 1(a)). This effect is maximized with the extreme flare cases. In the middle atmosphere regions (Figure 1(b)); however, the temperature response is reversed as we find warming (30–50 K) in the extreme case. In contrast, very minimal temperature changes are found between the quiescent and the moderate. The peak atmospheric heating occurs at the substellar hemisphere due to the stellar UV and proton injected energies and subsequent chemical reactions that affect the distributions of nitrogen oxides and water vapor (Figures A1 and A2). We find that middle-to-lower atmospheric N₂O and H₂O, which contribute to near-surface warming, increase in the moderate and maximize in the extreme case. In the upper atmosphere, dayside enhancement of these species is suppressed by UV photodissociation of H₂O. Thus, an increase in the nightside is slightly more evident for H₂O. In the middle atmosphere, we find increased moisture across the entire planet in both extreme and moderate. We also see localized elevated wind speeds (Figure A2).

Nitrogen species’ mixing ratios (NO₂, N₂O) all show a steady increase as the stellar flaring activity is dialed up (Figure A3), where O₃ behavior is reversed and its mixing ratio scales nonmonotonically with increasing flare activity (Figure A4). Stratospheric ozone, which influences the degree of the temperature inversion above the tropopause, reduces in its mixing ratios with increasing levels of stellar flaring activity. These latter photochemical responses generally echo those of previous 1D modeling results. In the middle atmosphere, midlatitude enhancement of wind speeds and H₂O mixing ratios are found, with the latter particularly pronounced outside of the equatorial regions. Both types of perturbations, to regional atmospheric dynamics and water vapor abundance, may be affected by transport processes and three-dimensional atmospheric thermal structure.

Figure 1. Refer to the following caption and surrounding text. — **Figure 1.** Temporal medians of “extreme,” “moderate,” and “quiescent” simulations for Temperature (T in K). Medians are taken from 8 yr and for (a) 1.223 × 10⁻⁵ mbar and (b) 0.008 mbar atmospheric levels. Extreme and moderate medians are only taken for days when flares occur. Stippling uses a two-tailed Wilcoxon rank-sum test to represent “extreme” and “moderate” areas that are not statistically different from the “Qquiescent” one at p < 0.01.
Download figure:
Standard image High-resolution image

3.2. Regional Anomalies

Regional statistical analyses from the box and whiskers plots reveal negative temperature anomalies in the upper atmosphere (Figure 2(a)). The most substantial cooling is seen in the antistellar equatorial region, where we find the greatest temperature anomalies between the moderate and extreme cases (Figure 2(a)). In the middle atmosphere, the temperature anomalies are much greater in the extreme than in the moderate cases (Figure 2(b)). All four regions of the moderate experience minor cooling, and those of the extreme cases experience similar magnitudes of warming. Flare-induced cooling is highest for the extreme cases (up to 100–150 K). In the middle atmosphere (Figure 2(b)), the highest (but positive) temperature anomalies are also seen in the extreme, suggesting mesospheric heating. Simulations that experience frequent high energy flares but lowest overall number of flares (i.e., the active case) have intermediate responses in antistellar hemosphere (Figure 2(a)), but the smallest thermal responses in the substellar region(s).

Figure 2. Refer to the following caption and surrounding text. — **Figure 2.** Boxplots of extreme, active, and moderate anomalies for temperature (T in K), nitric oxide (NO, in mol mol⁻¹), and hydroxide (OH, in mol mol⁻¹). Anomalies are computed by subtracting the quiescent temporal field median from the daily medians of the three sets of simulated cases: extreme, moderate, and active boxplots in (a), (c), and (e) are computed for the 1.223 × 10⁻⁵ mbar middle atmospheric level, whereas the ones in (b), (d), and (f) represent the 0.008 mbar lower level. Each panel also represents six key regions: equatorial substellar (Eq Sub), equatorial antistellar (Eq Ant), midlatitude substellar (Mid Sub), midlatitude antistellar (Mid Ant) In all panels, at least one of the extreme, moderate, and active anomaly medians significantly differ from one another at the 1% level.
Download figure:
Standard image High-resolution image

Chemically, except for NO mixing ratios at the middle atmosphere (Figure A5(b)), NO and OH have greater anomalies for the extreme cases (Figures A5(a), (b), (e)) due to enhanced dayside production rates of NO and OH driven by the proton fluxes associated with the flare events. Similar nightside NO in the extreme and moderate (Figure A5(a)) suggests horizontal transport, whereas OH experiences little transport due to its much shorter chemical lifetime (Figure A5(e)).

The importance of flares and stellar proton events on atmospheric flow is also indicated in Figures 2(c) and (d). The positive wind velocity anomalies across all three sets of simulations show that flare-driven dynamics, particularly in the equatorial antistellar and equatorial west regions for the upper and middle atmosphere, can be considerable (Figures 2(c), (d)).

Enhancements in oxidized greenhouse gases such as N₂O (Figures 2(e)–(f)) in upper and middle atmosphere are most evident in the extreme case. This suggests that greenhouse warming is particularly relevant for climates around the most active stars (reflected in Figure 2(b)). Despite having a few outliers due to sporadically large events in the moderate and active, on average the effects on N₂O are small. A similar trend is also found for effects on water vapor and hydroxyl radical (Figure A5).

Ozone is well-studied constituent that can alter atmospheric temperature by scattering and absorbing UV radiation through the Chapman cycle. Our results show that O₃ abundances are sensitive to stellar flare inputs (see also A. Segura et al. 2010; J. L. Grenfell et al. 2012; M. Scheucher et al. 2020). The moderate scenario shows an overall positive O₃ mixing ratio anomaly, while the extreme shows a negative anomaly in the upper atmosphere (Figure A6(c)). This nonlinearity between input UV and the resultant ozone abundance, hinted at in the time-averaged maps in Figure A4, is absent when examining the middle and stratospheric atmospheric abundances. For the middle atmospheric O₃ in particular, only negative anomalies are found in both extreme and moderate at the substellar hemisphere, while some areas of the planet experience almost no net losses (Figure A6(d)).

3.3. Time Series

Much of the climate dynamic feedback suggests quantifying the temporal change of flare-affected atmosphereres on annual timescales may provide further insights. Figure 3 shows the time series of stratospheric temperature, wind speed, water vapor mixing ratios, and ozone mixing ratios, partly driven by direct UV heating and indirectly through chemical reactions. Both NO and OH are strongly perturbed (Figure A7), with the a more spatially varied abundance seen for NO and a higher dayside concentration for OH. The chosen altitude shows that a reversal inday–night temperature differences in the upper atmosphere can be driven by the external energy input provided by stellar events. In general, compared to the quiescent simulations, we see that the inclusion of stellar events pushes the atmospheric state of the planet beyond that suggested by its internal climate variability alone. Peak temperatures are found to be comparable between the extreme and the active (Figure 3(a)) due to the injection of similar levels of highly energetic events. On average, however, the extreme sees the highest temporal-mean temperature change across the 8 yr period. Except for one or two temperature spikes, muted thermal oscillations are found in both the moderate and active cases.

Figure 3. Refer to the following caption and surrounding text. — **Figure 3.** Global-median temperature (K) and horizontal wind speeds (m s⁻¹) at 0.008 mbar, or the middle atmosphere, showing the results for a planet around actively flaring stars (gold and red), around a moderately active star (blue), around a quiescent/inactive star (black), and a more energetic star with a lower overall flare rate (orange). Time series' are computed from the 8 yr simulation.
Download figure:
Standard image High-resolution image

The NO and OH abundance fluctuations often serve as indicators of initial effects of ion pair production due to incident charged particles and UV photons accompanied by stellar events. Differences between our results and those of previous 1D models, particularly for NO and other transport-prone constituents, are due to the 3D nature of the CCM employed. Nitrogen oxides are first produced on the dayside hemisphere and then transported to the terminators, where nightside sustenance can be maintained (H. Chen et al. 2021; R. J. Ridgway et al. 2023). The chemical residence time of NO is around 1–2 weeks in Earth’s middle atmosphere, and based on satellite measurements, the NO decay rate toward its preperturbed equilibrium state is crucial to determining the atmosphere’s mean temperature change during a geomagnetic storm. Our results show that NO lifetimes are comparable to that on present Earth. Still, the relative enhancements of day versus night sides depend on the ambient condition of the atmosphere at the time of the flaring event. The cascading influence upon key molecular species, including N₂O, NO₂, H₂O, and O₃ can be seen in Figures 3(c), (d), and A7.

For H₂O in the moderate, only the largest events appear to disrupt their time-evolving mixing ratios substantially. While the effects on temperature are relatively damped for the active case, its H₂O abundance change is much more considerable (Figure 3(c)), highlighting feedback processes even in the absence in atmospheric temperature changes. For the extreme case, the H₂O levels are at times enhanced into the moist greenhouse regime, or >10⁻³ mol mol⁻¹ (J. F. Kasting 1988; E. T. Wolf & O. B. Toon 2015). This enhancement is key to an exoplanet habitability and interpretation (see Section 4). For the quiescent run, the observed O₃ oscillation is tied to intrinsic variations in atmospheric temperatures and the longitudinally asymmetric stratospheric wind oscillation (LASO; M. Cohen et al. 2022). In the moderate run, such oscillatory behavior is suppressed by flares due to the mismatch of flare injection rates and the period of the LASO. In the extreme run, traces of periodic variations in ozone abundance are largely removed due to the rapid ozone erosion through catalytic reaction cycles initiated by NO and OH.

N₂O, NO₂, and H₂O mixing ratio variabilities in the quiescent scenario are muted or absent due to the dryness of the upper stratosphere/mesosphere and a lack of a significant N₂O/NO₂ source in the upper atmosphere of temperate planets (Figure 3(c) and Figure A7). We find that peak production of N₂O and NO₂ in the extreme and moderate simulations can be several orders of magnitude greater than planets around quiescent stars (which lack a distinct nitrogen chemistry source in the middle atmosphere). In the moderate run, NO₂ levels are not sustained by flares, while N₂O levels are quasisustained (Figure A7(a)). In the extreme run, N₂O levels are perturbed to a new regime where its mixing ratios fluctuate between ∼2 × 10⁻⁵ and ∼3 × 10⁻⁵ mol mol⁻¹, whereas NO₂ is rapidly returned to its preflare mixing ratios even after the largest of the flaring events.

The effects of stellar flares are particularly pronounced in the middle atmosphere. At 0.008 mbar (in the mesosphere), virtually all the observed atmospheric temperature oscillations are caused by the quasi-random stellar inputs (Figure 3(a)). At month one but before the flaring event at month three, the atmospheric temperatures of the extreme and quiescent cases are similar, with the moderate already showing some degree of cooling due to the 2 yr cutoff of our full 10 yr model data (see Section 2.1). The moderate shows much cooler temperatures throughout the period examined in this regime, with occasional temperature peaks at approximately months 40 and 73. Due to the greater assumed flare frequency and flare energies, the extreme scenario is consistently heated by 30–60 K above the temperature of the quiescent. In the final months, a complete departure from the preflare state is seen. In contrast to the widening gap in atmospheric temperatures here, the temperature difference between the moderate and quiescent has been reduced throughout the simulated duration. Flaring could instigate successive cooling rather than warming episodes at other altitudes, as the earlier results indicate (see, e.g., Figure 1(a)). These higher levels are extremely low in their molecular densities; however, the warming effects at the stratospheric and mesospheric levels will likely have a greater impact on near-term observations than those in the upper atmosphere.

Stellar flares affect planetary wind speeds at various levels and generate interesting variability patterns not seen in climates without the inclusion of time-dependent stellar inputs. In Figure 3(b) once again in the middle atmosphere at 0.008 mbar, we find that for the 50–60 months, the flare-affected global-median wind velocities across all experiments are comparable to those generated by the natural climate variability (quiescent curve in Figure 3(b)). Seen together with the temperature variations, the extreme case, after the ≈63 yr mark, has constantly higher middle atmosphere temperatures and wind speeds, implying a new dynamical equilibrium is being established even when flares have been removed. On the other hand, the wind velocity variability for the moderate is comparable with the quiescent across the entire 8 yr duration.

3.4. A Glimpse Into Atmospheric Dynamical Perturbations

Cumulative and time-averaged effects are most relevant to habitability and observations. Here, we analyze the temporal-median anomalies between the extreme and quiescent scenarios. In the results above, simulated median wind speeds throughout the 8 yr suggest that regional enhancements are pronounced (Figures A2 and 3(d)). Wind velocity anomalies also show that flare-affected atmospheres have modulated horizontal flows that are more pronounced in specific regions of the planet. Constructing a vertical profile over longitude 0° to 50° and latitude −15° to 15°, or the eastern nightside region, we find that the median horizontal WS are 25–50 m s⁻¹ faster between 21 and 40 km in altitude in the extreme case (Figure 4(a)). Conversely, a predominant decrease in WS across its vertical profile is seen for the active scenario (Figure 4(b)), reflecting the interplay between photochemical-, and photolytic-induced temperature changes. Our results suggest that sudden increases in stellar UV and proton fluence can, at least in some regions of the planet, substantially energize the atmosphere and potentially lead to long-term asymmetries in the atmospheric dynamics on tidally locked Earth-like exoplanets.

Figure 4. Refer to the following caption and surrounding text. — **Figure 4.** Vertical profile of wind speed (m s⁻¹) temporal and field median difference between two sets of simulations over −15^∘–15^∘N, 0^∘–50^∘E. Circles represent differences that are not statistically significant at the 1% level following a two-tailed Wilcoxon rank-sum test.
Download figure:
Standard image High-resolution image

4. Discussion and Conclusions

Observational surveys suggest that time-dependent solar and stellar activity manifested through UV flaring and Type-II bursts are common (J. D. Alvarado-Gómez et al. 2020; A. M. Veronig et al. 2021; D. Ó. Fionnagáin et al. 2022; P. Mas-Buitrago et al. 2025). Previous theoretical studies investigating these effects on the temperature and dynamics of attendant atmospheres are restricted to Earth–Sun relationships. Here, using a general circulation model with interactive chemistry, we study stellar flares’ coupled climate and chemical effects on synchronously rotating planets with Earth-like atmospheric compositions around the host star TRAPPIST-1. In general, our photochemical results echo those of previous efforts using lower-dimensional photochemical models without including atmospheric dynamics and chemical advection. However, the 3D nature of our modeling and analysis techniques reveal new insights on how planets around active and flaring stars could instigate enhanced climate anomalies as implied by their thermospheric temperature, stratospheric humidity, and wind velocity perturbations on regional and global scales. Our results suggest that, in addition to initiating key photochemical reaction pathways and inducing photochemical disequilibrium, these types of stellar events could affect atmospheric dynamics and entirely alter the planets’ circulation regime in the most extreme scenarios. For stellar environments that do not have extreme excess in X-ray or EUV irradiation, transient stellar emissions may be the dominant channel through which the dynamics of substellar atmospheres are driven. In the long term, these effects could determine the thermal evolution of high-mean molecular weight atmospheres on potentially terrestrial exoplanets. The results of our study will have implications for the observations of direct-imaging missions (e.g., HWO and LIFE) that may be able to probe the astrophysically-influenced weather systems on habitable zone planets.

It is worth noting a few potential sources of uncertainties that may affect our simulated—among these uncertainties, two are related to the time-varying stellar inputs. The first is the potentially inaccurate UV spectra, as we replace the observed UV spectra from D. J. Wilson et al. (2021) with synthetic flare spectra that are not (necessarily) representative of real TRAPPIST-1 flares or compatible with the observed optical and infrared spectra. Second, we use the same proton energy spectra for every flaring event. The 2003 halloween geomagnetic storm is a well-documented event that may have stemmed from a strong solar flare accompanied by a coronal mass ejection. A more self-consistent approach would be using proton spectra dependent on flare characteristics (e.g., K. Herbst et al. 2019). However, how the emitted proton energy of cooler, low-mass stars would scale the total flare energy for each event is still being determined. But that would be beyond the chief goals of this pilot study. The other two sources of uncertainty are related to the climate model itself. First, we use daily flare cadences, whereas real flares occur on timescales between minutes and hours (M. Engebretson et al. 2008). For larger flares, our approximation should be valid as the duration of energetic particle precipitation could persist for tens of hours to days. Conversely, higher temporal resolution would modify the effects of smaller but more frequent stellar emissions (e.g., J. Pettit et al. 2018). Additionally, the cutoff of the atmosphere column of our climate model occurs at the lower thermosphere (or at ∼10⁻⁶ mba, which is ∼150 km). Resolving the important flare-associated photoionization and photoexcitation processes would require at least 300–400 km model-top. Such work has already seen development for the upper atmospheres of sub-Neptune-sized planets (D. Kubyshkina et al. 2024) and for Earth-like atmospheres with 1D simplified hydrodynamic thermosphere models (F. Tian et al. 2005). Given the large uncertainties associated with measuring exoplanet upper atmospheres, it remains to be seen whether a full-fledged 3D treatment is necessary for these problems (e.g., J.-L. Baudino et al. 2017).

Due to the computational expense of the 3D model, we have limited the numerical experiments to Earth-similar N₂–O₂-dominated atmospheres, following the focus of previous works on strongly oxygenated conditions (A. Segura et al. 2010; J. L. Grenfell et al. 2012; M. A. Tilley et al. 2019). Alternative but commonly assumed exoplanet archetypes include CO₂- and steam-dominated atmospheres, which will require the removal of much of the molecular oxygen that is innate to the chemical scheme (Mozart; D. E. Kinnison et al. 2007) of WACCM, including the exploration of Protozoic Earth-like atmospheres (e.g., A. V. Young et al. 2024). At the other end of the spectrum of climate archetypes are those akin to Archean Earth, or atmosphere composition in reduced redox conditions. In these conditions, one could use fully coupled climate models to understand how flaring events of the young main-sequence star can influence photochemistry and aerosol production. These effects will be even more dominant due to the increased stellar activity at early times compared to older systems. At present, very few fully coupled general circulation models used to study rocky exoplanets are capable of simulating true Archean climates with the effects of photochemical haze (J. K. Eager-Nash et al. 2023), while the other well-versed models for early Earth are in 1D (see applications of such models by G. Arney et al. 2016; R. Seeburger et al. 2023).

Flare activity can affect the surface climate evolution pathways of M-dwarf planets mainly by inducing abundance changes in key photon absorbers and scatterers. Previous work studying the early atmospheric chemistry of early Earth under a young Sun found that possible enhancement of greenhouse gases such as N₂O by a younger, more active Sun (V. Airapetian et al. 2016). Our results suggest that the troposphere may experience persistent warming caused by secular flaring through a similar catalytic reaction initiated by NO production. Stellar proton events could also affect CH₄ abundance, another strong greenhouse gas. On the Earth, coupled climate models in conjunction with satellite observations (e.g., SABER/TIMED and MIPAS) indicate that solar storms and geomagnetic activities can drive regional changes in temperature and wind fields (T. Von Clarmann et al. 2013; S.-M. Päivärinta et al. 2013). In particular, E. Becker & S. L. Vadas (2018) have shown that temporal heating of the middle atmosphere (90 km) can reach up to 18 K during large solar proton events, compared to our average middle atmosphere warming of 25 K (Figure 3). Another channel through which flares could affect the thermal structure is seen in the relationship between solar (geomagnetic) storms and the so-called sudden stratospheric warming. The former can generate heat ejected into the upper atmosphere, which Rossby waves could then transfer into the polar stratosphere. Other models found that the polar lower mesospheric cooling could induce an anomalous eastward flow following the gradient wind balance. Based on these results from Earth–Sun studies, there is a need for a more detailed analysis of stellar flare-induced wave propagation with a broader range of flare behavior and flare-proton fluence relationships.

The modeled middle atmosphere effects of time-dependent stellar flaring may influence transmission or emission spectra via time-varying atmospheric circulation regimes and thermal structures. Sufficient energetic events may even penetrate down to the lower atmosphere and affect direct-imaging observations of otherwise unperturbed climates. These eruptive events may lead to vastly divergent characteristics depending on the specific nature of the event and stellar spectral type. Such suggestions have already been provided by reconstruction of the activities of the young Sun and extrapolations to other main-sequence stars (J. Hu et al. 2022). Therefore, apart from the use of a flare generator that provides light curves and EM spectra, a better estimate of energetic particle precipitation requires improved coronal mass ejection and magnetospheric transport models (F. Fraschetti et al. 2022; K. Scherer et al 2025). More generally, galactic-wide astrophysical events, beyond sources within the stellar/planetary system itself (e.g., H. Chen et al. 2018; A. Ambrifi et al. 2022; K. I. Sippy et al. 2025), will be substantially different in terms of their output (photon and proton) spectra and how they interact with the astrospheres of the planetary medium, relative to stellar type events.

Other coupled chemistry-climate studies using WACCM6 have found the photochemical variability driven by stellar UV emissions is probably within the detection limits of future flagship direct-imaging missions (G. Cooke et al. 2023a). However, due to stellar UV uncertainties, degenerate interpretations of key spectroscopic windows would make it even more challenging for near-term missions such as the JWST (G. Cooke et al. 2023b). Another caveat is that different subversions of WACCM (i.e., between versions 4–6) are expected to produce disparate results due to the different physical parameterizations used, and thus it may be necessary to start with comparing/contrasting the pure radiative effects of key photochemical species using a GCM (R. Deitrick et al. 2025), without the full suite of complex chemical reaction network For example, J. Bin et al. (2018) found that distinct microphysical and aerosol schemes within CAM4 and CAM5 can lead to different predictions in the climate states (and thus habitable zone limits) of synchronously rotating planets in the habitable zone. Given the diverse behavior of stellar emissions and the varied forms of photochemical and photolytic effects, careful model intercomparison efforts of this sort are crucial.

The suppression, or even the complete erosion, of a substantial ozone layer is found in the extreme simulated cases at the end of the 8 yr runtime, a behavior reminiscent of earlier 1D model predictions by M. A. Tilley et al. (2019). However, their study was focused on extremely long (Gyr) timescales. As shown in our simulations, climate variability on decadal/annual timescales could be an important factor for near-term observations that would integrate over multiple (10–20) transits (T. J. Fauchez et al. 2019; J. Lustig-Yaeger et al. 2019). Our results suggest that changes in temperature and other climate factors such as wind speed are almost entirely photochemically driven with relatively little contribution from direct heating imparted by UV photons. This is in contrast to XUV-irradiated planets around young main-sequence stars, wherein the absorbed X-ray and EUV energy is transferred to thermal and mechanical energy (C. P. Johnstone et al. 2021; P.-G. Gu & H. Chen 2023; G. W. King et al. 2024; J. Krissansen-Totton et al. 2024; J. G. Rogers et al. 2024). Initially, climate oscillations are the most pronounced for those around active stars (e.g., Figure 3). Still, after a few years, much of their oscillatory behaviors are dampened by reduced photochemical seasonality due to the near depletion of key species such as ozone. Thus, planet atmospheres around stars with sporadic rather than consecutive eruptive events (which permits more time between each event for chemical species such as ozone to recover photochemically) would likely experience the greatest degrees of variability. These candidates will likely be situated around moderately active stars.

Acknowledgments

T.D.K. acknowledges support from NASA Habitable Worlds Program grant No. 80NSSC24K0215.

Data Availability

The model data behind the figures and the scrips used to produce them can be found on https://zenodo.org/records/15289304, and is also accessible via https://doi.org/10.5281/zenodo.15289304.