© N. Jakowski & M. Mainul Hoque, Published by EDP Sciences 2025

1 Introduction

Climate change due to enhanced greenhouse gas concentration concerns the entire “Earth System”, including the hydrosphere, atmosphere, thermosphere, ionosphere, and magnetosphere. Being aware of a permanently changing environment of the Earth, the question arises how the Magnetosphere-Ionosphere-Thermosphere (MIT) system may change over long time series. Considering the dominating impact of solar irradiation on the structure and dynamics of the MIT system, separating a potential anthropogenic impact on ionospheric observables from their natural variability is an extremely difficult task.

More than 30 years ago, initial studies began investigating how surface climate change due to greenhouse gases could affect the upper atmosphere and thermosphere, including related coupling processes within the MIT system. Early studies by Roble & Dickinson (1989), Rishbeth (1990), Rishbeth & Roble (1992) were crucial in addressing the broader impacts of climate change on our planet’s atmosphere. Up to now, the objective is to find out how strong the potential impact of climate change might be in highly interacting geospheres. The analysis of long-term ionospheric databases is one way to figure out long-term trends in the ionospheric structure and plasma density that may be linked to climate change on the ground. Thus, many attempts have been made to separate long-term trends of ionospheric parameters such as the F2 layer plasma critical frequency (foF2) at the peak density height (hmF2) or E-layer peak plasma frequency (foE) from the dominating impact of solar radiation (Bremer, 1992; Rishbeth, 1997; Mikhailov & Marin, 2001; Laštovička et al., 2006; Cnossen & Richmond, 2008; Bremer et al., 2012; Mielich & Bremer, 2013; Cnossen & Franzke, 2014; Sivakandan et al., 2023, Tan Jun Rios et al., 2025).

Long-term effects of the total electron content (TEC) of the ionosphere can be studied separately for most data sets starting in 1995 when the Global Positioning System GPS became operational (e.g., Lean et al., 2016; Emmert et al., 2017; Laštovička et al., 2017; Jakowski et al., 2024). This period of available data, however, is too short to study the relationship of TEC with solar radiation. As will be shown in Section 4.2, a substantial increase in ground temperature occurred already between approximately 1965 and 1975.

In addition to the ionospheric database, solar radiation characteristics must be considered when examining how ionospheric parameters relate to solar radiation proxies in long-term data analysis. Several studies (e.g., Dudok de Wit & Bruinsma, 2017; Mielich & Bremer, 2013; Zossi et al., 2025a) highlight differences among solar activity indices. As stated by Zossi et al. (2025a), these discussions highlight the need for ongoing evaluation of solar indices to ensure accurate long-term predictions. Additionally, trends in indices reported by Clette (2021), Mursula et al. (2022, 2024), and Zossi et al. (2025a, 2025b) affect relationships between atmospheric and solar parameters. Drawing on recent research (Laštovicka & Burešová, 2023; Laštovička, 2024; Elias et al., 2024), this study employs the solar radio flux index F30 at a 30 cm wavelength to assess long-term links between solar radiation and ionospheric parameters, replacing the previously used F10.7 index. Dudok de Wit & Bruinsma (2017) also recommend F30 for investigating its connection with neutral gas density modeling. Data preparation and analysis methods are critical for obtaining accurate long-term trend results. Several studies (e.g., Ulich et al., 2003; Laštovička et al., 2012, Laštovicka & Jelínek, 2019; Cai et al., 2019) have reported weak or inconsistent long-term trends due to varying analytical approaches. Laštovička et al. (2006) compared different methods for studying foF2 plasma frequency trends, and Danilov (2017) noted that methodological uncertainties may also stem from variations in data selection or preparation.

To extract climate change signals from ionospheric data by filtering out solar radiation effects using solar activity indices, it is important to account for MIT’s complexity related to greenhouse gas emissions. It is evident that the existence of several physical impact factors, which may vary over time, can complicate explanations of long-term trend estimates. Substantial impact factors include changes in geomagnetic field configuration (Cnossen & Richmond, 2008; Elias et al., 2010), geomagnetic activity (Mikhailov & Marin, 2001), solar irradiation (Tapping & Valdés, 2011; Clette, 2021; Mursula et al., 2022, 2024; Zossi et al., 2025b), dynamic coupling processes such as gravity waves (Oliver et al., 2013), and alterations in thermospheric density and composition from increasing greenhouse gas concentrations (Qian et al., 2006, 2009; Solomon et al., 2015; Parker et al., 2025).

Roble & Dickinson (1989) found that global cooling reduces the F2 layer peak height (hmF2) by 15–20 km, with little effect on electron density. However, Qian et al. (2009), using a model with double CO₂, showed that both hmF2 and peak electron density (NmF2) can decrease by up to 40%, depending on factors such as location, time, season, and solar activity. This result supports the need for further investigations into the long-term behavior of ionospheric observation data. In this paper, we investigate the long-term behavior of NmF2 in more detail by including ionosonde measurements starting before 1960. In Section 4, we compare the long-term behavior of NmF2 with the temperature increase on the ground.

Human activities, principally through emissions of greenhouse gases, have unequivocally caused global warming, with global surface temperature reaching 1.1 °C above 1850–1900 in 2011–2020. Global greenhouse gas emissions have continued to increase, with unequal historical and ongoing contributions (IPCC, 2023). These changes have far-reaching consequences for ecosystems on the ground but also for processes in the atmosphere up to the thermosphere and related interactions with the ionosphere. Future climate is partly determined by the magnitude of current emissions of greenhouse gases, aerosols, and other natural and man-made forcings. This also concerns future structure and processes in the thermosphere and ionosphere to a certain extent, as addressed in this paper.

2 Database and data processing method

The peak electron density NmF2 used in this study has been derived from ground-based vertical sounding observations providing the plasma frequency foF2 at the peak electron density height. We established a database of hourly foF2 vertical sounding data from 10 ionosonde stations located in mid-latitudes at quite different longitude sectors, in order to obtain continuous long time series starting whenever possible from 1952. All stations used in this paper are listed in Table 1.

Juliusruh data were obtained directly from the Juliusruh observatory, but can also be obtained from the ionospheric section of the WDC at the Australian Space Weather Forecasting Centre (ASWFC). From ASWFC, we also downloaded hourly foF2 data for the ionosonde stations Slough, Brisbane, and Port Stanley. Manually scaled foF2 data from the Japanese station Kokubunji (TO535/TO536) and Wakkanai (WK546) are obtained from the National Institute of Information and Communications Technology. American foF2 data from the US-station Boulder (BC840) were obtained via the World Data Center C1 (1958/07-2002/12, manually scaled) and automatically scaled data from the Global Ionospheric Radio Observatory (Reinisch & Galkin, 2011) (2004/03-2022/12).

The foF2 data quality is different, ranging from 0.1 MHz accuracy for manually scaled data to about 0.5 MHz for automatically scaled ionograms data (e.g., Jakowski et al., 2024). The electron density is derived from the foF2 data by the well-known relation:

$$\mathrm{NmF}2\mathrm{ }=1.24\times 10^{-2}(\mathrm{foF}2{)}^2,$$(1)

where foF2 is given in MHz and the F2-peak plasma density NmF2 is provided in 10¹² m⁻³.

To reduce the impact of diurnal plasma dynamics, we averaged noontime NmF2 data in the period 12–14 LT. This period is characterized by nearly equilibrium conditions in the ionosphere around local noon. Thus, if NmF2 is mentioned in the text, it means the noontime average value of NmF2 unless there is another meaning indicated.

To reject data outliers, commonly appearing in observation data, e.g., caused by solar storms and related effects in the MIT system, we have extracted monthly medians of NmF2 for further data analysis. No additional filtering was done to exclude geomagnetically disturbed days or other perturbations. Monthly medians of NmF2 have been related to corresponding arithmetic mean values of the solar radio flux indices F10.7 and F30, given in solar flux units (1 sfu = 10⁻²² W Hz⁻¹ m⁻²). To further reduce the strong seasonal variation of NmF2, the arithmetic means of their monthly medians were calculated to get yearly averages (cf. panels a and c in Fig. 1). The yearly NmF2 data were then used to mitigate the solar cycle effect by computing 11 years of running arithmetic means for all station data sets.

Figure 1 Long-term behavior of NmF2 over the German ionosonde station Juliusruh in comparison with solar radio flux indices F10.7 (a, b) and F30 (c, d). Panels a and c: Yearly means of NmF2 monthly medians and solar indices F10.7 (a) and F30 (b). Panels b and d: Sliding solar cycle averages (previous 11 years) of NmF2 estimates (NmF2*) and linear models of NmF2* (until 1984) run by sliding solar cycle averages of F10.7 (F10.7*) and F30 (F30*). The differences between observations and linear NmF2* models, ΔNmF2, are plotted in the lower part. scaled at the right side. The dashed lines indicate the RMS ± error (RMSE) range of the linear models taken as reference.

Jakowski et al. (2024) have shown that 11-year sliding averages of NmF2 (NmF2*) can be modeled by the associated F10.7 radio flux index (F10.7*) with around 1.0 – 1.3 × 10¹⁰ m⁻³ accuracy prior to 1985. We used a similar approach, replacing F10.7 with the more solar-sensitive F30 index (Laštovicka & Burešová, 2023; Laštovička, 2024), and analyzed NmF2 data from 1952 to 1984, including pre-1990 data from Port Stanley and Brisbane. In analogy to former computations, here the F30* based model of NmF2* until 1984 is taken as a reference for comparing subsequent observation values of the selected ionosonde station. Consequently, the long-term behavior of NmF2* was then computed for each ionosonde station by subtracting model values from corresponding observations, providing ΔNmF2.

To get statistically representative results for the Northern and Southern hemispheres, the long-term trends from different ionosondes are averaged by computing arithmetic means of related ΔNmF2 values. The results are described by linear regression lines for absolute and percentage values. Finally, these linear ΔNmF2 fits for both hemispheres were adapted to corresponding temperature anomaly data from the ground.

3 Relationship of vertical sounding data and solar radiation

For our long-term trends study, we selected ionosonde observations covering at least six solar cycles with a small number of data gaps. Ionosonde data available for shorter periods from the nineteen fifties until the mid-eighties are also utilized to strengthen the model fit. In order to get a global view on long-term trends of NmF2, we analyzed ionosonde data obtained at European, Asian, and American longitude sectors. Furthermore, to avoid severe geomagnetic perturbation impact at high latitudes and strong solar impact combined with equatorial crest-related plasma dynamics at low latitudes, we limited our analysis to data from ionosondes located at mid-latitudes at both hemispheres.

To provide an impression of how relationships of NmF2 with solar indices F10.7 and F30 differ, we demonstrate computations using both indices for Juliusruh data in Figure 1. It should be mentioned that in order to utilize data as much as possible, we used foF2 data from Juliusruh and other stations principally without additional filtering. We rely on the fact that extreme values, e.g., due to enhanced geomagnetic activity, are filtered out in monthly medians. Yearly arithmetic means of monthly medians of NmF2, and related data of solar indices were computed to mitigate the strong seasonal variation, as can be seen in the panels a and c for both indices. To mitigate the remaining strong solar cycle dynamics, running arithmetic means of yearly averages were computed over a full solar cycle length of 11 years.

The NmF2 model results are shown in panels b and d of Figure 1 for both F10.7* and F30*. It is obvious that the subsequent development of the deviation of NmF2* from the initial linear approach ΔNmF2 is clearer and stronger for F30* than for F10.7* indices over the time until 2023. Laštovička (2024) has recommended F30 as the best solar proxy for studying long-term trends of foF2 at middle latitudes. Dudok De Wit & Bruinsma (2017) found that F30 improves the response of the thermospheric density to solar forcing compared with F10.7. Consequently, F30 was chosen as the solar activity proxy for the following calculations and trend assessments.

In principal agreement with former results obtained by Jakowski et al. (2024), the fit of 11 years running-averaged radio flux data to corresponding NmF2* values is very good with an RMS deviation of less than 1 × 10¹⁰ m⁻³ that corresponds with about 0.05 MHz for foF2 between 1969 and 1984. This is better than the accuracy of about 0.1 MHz for manually scaled ionosonde data of Juliusruh (Jakowski et al., 2024). The linear model (LM) of the running mean of NmF2 derived for the period 1969–1984 is given by LM(NmF2*) = F30*/0.8 – 20.4 in 10¹⁰ m⁻³ (RMSE: 0.9 × 10¹⁰ m⁻³). As panel d shows, the deviation ΔNmF2 = NmF2* − LM(NmF2*) grows permanently after about 1985, reaching −14.2 × 10¹⁰ m⁻³ in 2023, which corresponds with a percentage decrease of 20.8%.

Besides Juliusruh, also data from the ionosonde stations Slough (SL051) and Chilton (RL052) provide a continuous data set in Europe from 1952 to 2021 (cf. Table 1). Since the location of both stations is very close (Slough: 51.5° N; −0.6° E, Chilton: 51.1° N; −0.6° E), a combined data set has been arranged to fully cover the long period 1952–2021. Whereas the measurements from Slough cover the period 1952–1993, the Chilton ionosonde station started operation in 1993. To get a continuous time series of NmF2, the monthly medians of Slough and Chilton data have been merged at the end of 1993. As the panel of Figure 2 shows, the yearly averages of monthly medians of NmF2 show no irregular feature during the transition of data sources from Slough to Chilton. The assumed linear model LM(NmF2*) of the 11 years running NmF2 data derived for the period 1963–1984 is given by LM(NmF2*) = F30*/0.9 – 9.4 in 10¹⁰ m⁻³ (RMSE: 1 × 10¹⁰ m⁻³). The deviation ΔNmF2 grows up to −14.8 × 10¹⁰ m⁻³ in 2021, which corresponds with a reduction of 22.5%.

Figure 2 Long-term behavior of NmF2* over the ionosonde stations Slough/ Chilton from 1952 to 2021 (a, b) and Boulder from 1958 to 2022 (c, d). Panels a and c: yearly means of NmF2 and solar index F30. Panels b and d: Sliding solar cycle averages (previous 11 years) of NmF2 estimates (NmF2*) and linear models of NmF2* (until 1984) run by sliding solar cycle averages of F30 (F30*). The differences between observations and corresponding linear NmF2* models, ΔNmF2 are plotted in the lower part, scaled at the right side. Dashed lines indicate the RMS ± error range.

The data coverage for the American ionosonde station Boulder is slightly reduced in the early phase. Later, a data gap includes the full year 2003. To perform the analysis, the gap was filled with a linearly interpolated value for NmF2. The assumed linear approach of NmF2* for the period 1969–1984 is given by LM(NmF2*) = F30*/0.77 − 9.6 in 10¹⁰ m⁻³ (RMSE:1.0 × 10¹⁰ m⁻³). The deviation ΔNmF2 of observations from the linear model starts around 1980–1985 and reaches a value of −17.6 × 10¹⁰ m⁻³ in 2022 that corresponds with a reduction of 21.7%.

The analysis described in the previous section is repeated in the same way for several stations in the Asian longitude sector. Although NmF2 data from the Japanese ionosonde station Kokubunji suffer from a data gap of 4 months in 2003, no interpolation has been applied. The linear model of the 11-year running NmF2 average, derived for the period 1963–1984, is given by LM(NmF2*) = F30*/0.73 + 3.4 in 10¹⁰ m⁻³, RMSE: 0.7 × 10¹⁰ m⁻³). The deviation ΔNmF2 starts around 1985–1990, reaching a strong decrease of −23.2 × 10¹⁰ electrons per m³ in 2022, which corresponds with a reduction of 24.0%.

As shown in Figure 3, NmF2 data from the Japanese ionosonde station Wakkanai suffer from a large data gap from 1988 to 2000, as indicated in panels c and d, making a closed analysis over the full time span impossible. However, due to a good data coverage from 1950 to 1987, it was possible to derive a linear model for NmF2* data until 1984 in the same way as for other ionosonde stations. Despite the large gap, it was possible to compute the reduction of observed NmF2 against the linear model as shown in Figure 3 from 2001 onward. The linear approach of the 11-year running mean of NmF2, derived for the period 1963-1984, is given by LM(NmF2*) = F30*/0.62 – 40.5 in 10¹⁰ m⁻³ (RMSE: 1.3 × 10¹⁰ m⁻³). The absolute deviation ΔNmF2 starts at around 1985–1990, reaching a decrease of −12.1 × 10¹⁰ m⁻³ in 2023 that corresponds with a reduction of 17.4%.

Figure 3 Long-term behavior of NmF2 over Kokubunji from 1958 to 2022 (a, b) and Wakkanai from 1950 to 2023 (c, d). Panels a and c: yearly means of NmF2 and F30. Panels b and d: NmF2* estimates and linear models of NmF2* (until 1984) as a function of F30*. The differences between observations and corresponding linear NmF2* models, ΔNmF2, are presented in separate plots scaled at the right side. Dashed lines indicate the RMS ± error range of the linear models taken as reference.

To include data from the Southern hemisphere, long-term datasets obtained for Canberra and Hobart were used. Additional data from ionosonde stations Port Stanley and Brisbane were used to estimate NmF2* models for observations ranging up to 1984.

For Canberra, the data coverage is essentially reduced mainly in the early phase, with some monthly data gaps. Nevertheless, no interpolation has been applied. The assumed linear approach of the 11 years running NmF2 means for the period 1963–1984 is given by LM(NmF2*) = F30*/1.15 + 11.5 in 10¹⁰ m⁻³, (RMSE: 1.0 × 10¹⁰ m⁻³). The difference ΔNmF2 between the observations and the model is about −15.0 × 10¹⁰ m⁻³ in 2022, which corresponds to a reduction of 21.1%.

Unfortunately, the data coverage for the Australian ionosonde Hobart is essentially reduced in the early phase in 1959 and 1960. To save the analysis, the gap for 1960 was filled with a linearly interpolated value. The assumed linear approach of NmF2* is then given by LM(NmF2*) = F30*/1.05 – 4.6 in 10¹⁰ m⁻³(RMSE: 2.3 × 10¹⁰ m⁻³). The less accurate fit up to 1984 is probably due to a less reliable data quality. The deviation of observations ΔNmF2 starts around 1980–1985 and reaches a level of −11.1 × 10¹⁰ m⁻³ in 2023 that corresponds with a reduction of 17.6%.

4 Discussion

4.1 Average deviations from references at Northern and Southern hemispheres

As can be seen in the Figures 1–4, there is a consistent behavior of the NmF2*-F30* relationship. In particular, there is a clear linear relationship between the 11 years running averaged noontime values of the peak electron density and a linear model approach until 1984, of course, different for each ionosonde station. The fitting results are summarized for each station in Table 2.

Figure 4 Long-term behavior of NmF2 over the Australian ionosonde stations Canberra from 1952 to 2022 (a, b) and Hobart from 1950 to 2022 (c, d). Panels a and c: yearly means of NmF2 and F30. Panels b and d: NmF2* estimates and linear models of NmF2* (until 1984) as a function of F30*. The differences between observations and corresponding linear NmF2* models, ΔNmF2, are presented in separate plots scaled at the right side. Dashed lines indicate the RMS ± error range of the linear models taken as reference.

Whereas the relationship between NmF2* and F30* is linear until about 1984, the NmF2* values start to deviate more and more from the linear model level towards lower values afterwards. To statistically underline this statement, we have combined the differences ΔNmF2 between the running means of NmF2 and corresponding linear models for different ionosonde stations by arithmetic means. Data sets of ionosondes from the Northern and Southern hemispheres have been computed separately. The results are shown in Figure 5. To improve the statistics for the Southern hemisphere, we have added data from Port Stanley (PSJ5J)/South America and Brisbane (BR52P)/Australia, whose data sets cover the time period from 1952 until 1986, thus supporting the check of the linear relationship in the early phase of the time series. The high accuracy of the linear approach can be seen in the left panels a and b of Figure 5, providing a zoomed picture of the difference between the 11 years running means and corresponding linear models. Hereafter, the residuals are very small, varying between about ±1.5 × 10¹⁰ m⁻³ or ±1.5% until 1982. It is worth mentioning that the small residuals in both the Northern and Southern hemispheres display similar patterns, indicating a systematic effect not eliminated by the 11-year running mean of NmF2. Such a pattern might occur because the length of solar cycles changes for unclear reasons (e.g., Ptitsyna & Demina, 2023). ΔNmF2 shows periodic patterns, especially for the Southern hemisphere, as seen in panel d of Figure 5. Frequency modulation of the 11-year solar cycle by Gleisberg periods (about 60, 88, and 115 years) yields periods between 7 and 15 years, matching observed solar cycle variability (Ptitsyna & Demina, 2023). The periodic pattern of NmF2* – 7.1 years before 1985 and 9.8 years after 1980 – corresponds to these shorter modulation periods. Further research is needed to understand this finding.

Figure 5 Mean differences ΔNmF2 between estimates of NmF2* at ionosonde stations in the Northern (a and c) and Southern (b and d) hemispheres and corresponding models LM(NmF2*) as a function of F30*. Panels a and b zoom in on the early phase of measurements until 1982. Panels c and d show differences between the linear NmF2* models and observations over the full period 1963–2022, including linear regression lines starting in 1980.

Panels c and d of Figure 5 show the difference between NmF2* and corresponding linear models until 2022. Starting shortly after 1980, this difference is permanently growing, indicating a linear decrease of 4.0 × 10¹⁰ m⁻³ or 5% per decade for the Northern hemisphere. Overall, the peak electron density decreases by about 1.6 × 10¹¹ m⁻³ or about 20.6% until 2022 in the Northern hemisphere. In the Southern hemisphere, the decrease is about 2.7 × 10¹⁰ m⁻³ or 4% per decade, which amounts to up to a reduction of 1.2 × 10¹¹ m⁻³ or about 18.2% until 2022. It has to be mentioned that the choice of the break point has some impact on the fitting results. Here, we always selected 1980 as the starting point for the fitting computation. The ratio of Northern to the Southern hemisphere slopes for absolute and percentage values of ΔNmF2 is 1.7 and 1.4, respectively.

4.2 Comparison with the temperature anomaly at the Earth’s surface

In this section, we investigate how the strong decrease of the ionospheric peak electron density NmF2 observed after 1980 is correlated with temperature increases at the Earth’s surface. Roble & Dickinson (1989) and Rishbeth (1990) posed the question first of how the greenhouse effect might impact the thermosphere structure. After analyzing orbit data of five satellites during solar minimum years 1976, 1986, and 1996, Keating et al. (2000) derived thermospheric density changes of about 5% at about 350 km height. Emmert et al. (2004) report a density reduction in the order of 2–4% per decade in the height range of 200–700 km that includes the F2 layer peak density height. As pointed out by Akmaev et al. (2006), the cooling of the stratosphere and mesosphere due to radiative forcing of CO₂, O_3, and H₂O extends very well to thermospheric heights. Their simulations, made for medium solar activity, indicate a density reduction of 4% near 200 km height. In a more recent study, Emmert (2015) estimated density trends of about 2% per decade at 400 km height from satellite orbit data in the period 1967–2005. Solomon et al. (2015) have performed a 3D-simulation with the model NCAR TIME-GCM over 12 years, thus covering a full solar cycle. They found a density decrease of 5% at solar minimum and a decrease of 2% per decade at solar maximum conditions at 400 km height. It should be underlined that the simulations used a realistic increase of CO₂ in the period 1996–2008. Simultaneous estimations of NmF2 and hmF2 revealed a decrease varying between 0.7–1.5% and 0.8 and 1.9% per decade, respectively. The authors conclude that the study of long-term trends of the neutral atmosphere is much more straightforward for understanding the anthropogenic impact than long-term trend investigations of the ionosphere. Indeed, as considered in numerous publications, the level of ionospheric ionization depends on many drivers, coupling processes, and geophysical conditions. The main driver is the solar radiation in the extreme ultraviolet (EUV) wavelength range. The intensity of the ionizing spectrum of the solar radiation is usually approximated by solar indices such as F10.7, F30, or Mg II (Elias et al., 2024; Laštovička & Burešová, 2023). Different indices may behave in a different way, as considered in Section 1. Other factors, impacting the ionospheric plasma density such as thermospheric density and composition, changes of the geomagnetic field configuration and enhanced activity, varying electromagnetic forces during solar storms or dynamic coupling processes from below play a certain role as pointed out earlier (e.g., Mikhailov & Marin, 2001; Cnossen & Richmond, 2008; Elias et al., 2010, Oliver et al., 2013). However, considering these impact factors, there is no indication of a long-term impact reducing the peak electron density by about 20% on average since the early 1980s as reported here. On the other hand, modeling studies simulating an enhanced carbon dioxide level by a factor of 2, report the possibility of a large depression of the F2 layer ionization of up to 40% (Qian et al., 2009). Considering the possibility of a sensitive response of the ionosphere to greenhouse gases and considering the fact that global warming at the Earth’s surface has principally the same origin, a correlation is expected between the measured global surface temperature increase and the derived long-term decrease of the peak electron density. The question is only how strong such a correlation appears in the observations. To answer this question, we have compared two data sets, characterizing the permanent decrease of the observed noontime peak electron density illustrated in Figure 5, with corresponding data sets describing the temperature increase over land and oceans at the Earth’s surface. The used temperature anomaly (TA) data are provided by the Japan Meteorology Agency (JMA). The TA values are referenced to long-term averages over the 1991–2020 baseline for the Northern as well as for the Southern hemisphere. For defining a baseline for NmF2 data, one has to consider the strong relationship of NmF2 with the permanently varying solar activity, i.e., the reference is not a constant value but must follow the solar irradiance variation. Thus, in analogy to the definition of the temperature anomaly at the Earth’s surface, we define here the electron density anomaly (EDA) in the ionosphere by the deviation of NmF2* from a linear model approach as a function of F30* for the time period 1963–1984 (two solar cycles) considered as reference. The linear relationship between 11 years averaged NmF2 and F30 has been shown to be very stable in this reference period for all ionosondes with small deviations not exceeding 1.5% in average. It is interesting to note that Donaldson et al. (2010), estimating ion temperatures at 350 km height from Incoherent Scatter Radar (ISR) measurements at Saint Santin, found a rather good correlation with solar irradiance until 1980 but the period 1976–2001 showed a complete break in that relationship. Analyzing long-term trends of ion temperatures based on ISR measurements at Millstone Hill, Zhang et al. (2011) had the impression that the cooling emerged “only” since the early 1980s. A different behavior of foF2 between the periods 1958–1979 and 1998–2010 was also reported by Danilov & Konstantinova (2013). Danilov (2017) considered the period 1958–1980 as the “etalon” period, characterized by much lower temporal changes of foF2 than during later observation periods. Considering the fact that the deviations from the reference model shown in Figures 1–4 are based on arithmetic means of yearly NmF2 values of the previous 11 years, corresponding processes in the MIT system should be dated back by about 4–6 years. Considering this time shift and the results shown in Figure 5, the physical “break point” of ionospheric behavior should lie in the second half of the 1970s. This would agree in particular with observations made by Donaldson et al. (2010). It is interesting to see in Figure 6 that the temperature anomaly in the Northern hemisphere also indicates a kind of a break point at around 1975 (panels a and c). To determine such a break point in the Southern hemisphere (panels b and d) is more uncertain due to measurement noise and a lower increase of the TA.

Figure 6 Temperature anomaly (TA) at the Earth’s surface at both hemispheres (red dashed line). Linear regressions of the TA are fitted by the EDA anomaly averages (absolute at panels a and b and percentage at panels c and d) computed from the combined data of ionosonde stations Boulder, Juliusruh, Slough/Chilton Kokubunji and Wakkanai in the Northern hemisphere (panels a and c) and by the ionosonde stations Brisbane, Canberra, Hobart and Port Stanley from the Southern hemisphere (panels b and d). To get a coincidence of linear fits for TA and EDA, absolute EDA averages have been scaled by a divisor as indicated and offset by several years indicated by TD in the plots.

As panels c and d of Figure 5 show, the growing absolute and percentage deviations from the reference level are almost linear for the Northern as well as for the Southern hemisphere. Regression lines have been computed both for the EDA and for the TA data sets. In a subsequent step, the EDA regression data are matched to the TA time series shown in Figure 6. Good agreement was achieved by fitting slope factors and time delays for EDA regression lines representing the different data sets. The good match of the Earth’s temperature anomaly and the increasing reduction of the noontime peak electron density in the ionosphere cannot necessarily be considered as an approval of a common physical mechanism. However, considering the good match, there are several aspects that favor the assumption of a common causality:

1. The temperature anomaly at the Earth’s surface and thermospheric cooling associated with a permanent decrease of the peak electron density have a common physical mechanism – the increase of greenhouse gas concentration. The question is only how strong the correlation appears in the data. The answer given here indicates a strong correlation.

2. TA differs noticeably in the Northern and Southern hemispheres. The slopes of TA increase are stronger in the Northern hemisphere, i.e., by a factor of 1.23 for absolute and 1.46 for percentage values. As mentioned in the previous section, the corresponding ratios of EDA are about 1.67 and 1.42 for absolute and percentage slopes, respectively. Although the ratio of slopes for absolute values differs, the slopes of percentage deviations of EPA and TA behave similarly. We think that this analogy of TA and EDA behavior substantially supports the assumption of a common origin.

3. It is expected that the propagation of climate changes recorded at the Earth’s surface requires some time to reach thermosphere/ionosphere altitudes. As the attempts to match TA and EDA data sets show, EDA effects are delayed, thus also supporting the assumption of a common physical origin. Considering that the match of EDA and TA data requires a time shift T_D of 5 and 7 years for absolute and percentage estimates, respectively, greenhouse gases would need up to 3 years to induce a measurable ionospheric reaction in the Northern hemisphere. For the Southern hemisphere, the fitting reveals much stronger delays of more than 10 years, which are probably not realistic. The reason for this result is probably a soft definition of a kind of break point for the Southern hemisphere. As Figure 6 indicates, the temperature increase in the Southern Hemisphere starts about 10 years earlier than in the Northern Hemisphere. This can explain the big T_D values over 12 to 15 years (Fig. 6, panels b and d). As Figure 5 indicates, the density decrease starts about 1–2 years earlier in the Southern hemisphere than in the Northern hemisphere. This qualitative analogy also supports the assumption of a common origin of global warming and electron density decrease since the 1980s. Since the estimation of time delays is imprecise in this analysis, it might be determined more precisely in future refined trend studies. A comparison with the results of atmospheric studies addressing the upward propagation of greenhouse gases up to thermospheric heights could improve our understanding of their impact on the entire atmosphere.

Considering the close relation of EDA and TA, we believe that greenhouse gases like CO₂ do not only cause global warming but also lead to a systematic and measurable decrease of the ionospheric electron density related to thermospheric cooling in the order of 10–15 K/decade (e.g. Ogawa et al., 2014) and contraction indicated by a decrease of the equivalent slab thickness (e.g. Jakowski et al., 2024). This solar cycle averaged trend is mostly independent of geomagnetic activity, local time, latitude, and season. If the decrease of the ionospheric ionization is projected to continue for next decades under similar solar-terrestrial relationships as described in this paper, we expect a 50% reduction of the ionospheric plasma density in mid-latitudes at the end of this century. The ionospheric effect at high and low latitudes also needs to be studied. The long-term reduction of the ionization level is an indication of a chain of complex coupling processes in the MIT system. Ionospheric range errors in space-based global navigation satellite systems (GNSS), radar systems for remote sensing, and radio astronomy observations will decrease. On the other hand, current ionospheric models utilizing coefficients that have been derived many years ago, like the coefficients used in the International Reference Ionosphere (Bilitza & Reinisch, 2008) or the Galíleo correction model NeQuick (Nava et al., 2008), should consider the long-term trend in updated versions. Other impact factors on technical systems will also change and need to be studied in detail. Thus, to give an example, thermospheric cooling will reduce satellite drag in low Earth orbit (LEO) constellations. Considering the rapid expansion in the number of satellites in low Earth orbit, their extended lifetime will lead to enhanced risks of collision hazards and cascading generation of more debris (Parker et al., 2025).

5 Summary and conclusions

To examine long-term trends of ionospheric ionization resulting from thermospheric cooling, composition changes, and contraction due to growing greenhouse gas concentrations, we have analyzed foF2 data collected from 10 ionosonde stations operating for more than 70 years. The long time series enabled separating peak electron density (NmF2) data into two periods with different characteristics. For the first-time interval 1952–1984 covering two solar cycles, we find a close linear relationship between the 11-year sliding means of NmF2 (NmF2*) and the corresponding 11-year sliding means of solar radio flux F30 (F30*). This strong linear relationship – representing a time period just before the significant intensification of global warming – serves as the baseline for subsequent changes in the ionospheric electron density. By defining this period as a reference, it is possible to quantify and track the changes in NmF2* in subsequent years. This approach allows for a clear comparison between the established baseline and the observed deviations in ionospheric electron density, thus facilitating the assessment of long-term trends and potential correlations with surface temperature anomalies. The use of this linear model as a reference seems to be crucial for distinguishing natural variability from anthropogenic influences in the ionosphere over the following decades.

NmF2* starts to deviate from this linear model approach at all ionosonde stations considered in this paper after about 1980. Since the ionosonde stations are located at both hemispheres, the obtained results refer to a global effect. The deviation of the 11-year sliding means of NmF2 from the corresponding linear approach reference may reach up to 1.6 × 10¹¹m⁻³ or about 20.6% until 2022 in the Northern hemisphere and up to 1.2 × 10¹¹ m⁻³ or about 18.2% until 2022 in the Southern hemisphere. In other words, the estimated rate of decrease is about 5% per decade in the Northern hemisphere and about 4% per decade in the Southern hemisphere. These rates of electron density reduction fit quite well with estimations of thermospheric density decreases in the order of about 5% per decade (e.g., Katsuda et al., 2023). The observed hemispherical difference in the ionospheric peak electron density mirrors the pattern seen in the temperature anomaly (TA) at the Earth’s surface. Specifically, the variations between the Northern and Southern hemispheres in electron density anomaly (EDA) align closely with the corresponding changes in surface temperature anomaly. Furthermore, a comparison between temporal EDA and TA-signatures suggests a time delay of several years that seems to be reasonable. Atmospheric research addressing upward transport times of greenhouse gases could evaluate this finding.

We believe the analogy between temperature anomaly and electron density anomaly EDA supports the assumption of a shared physical mechanism. This parallel is particularly significant in the context of long-term changes observed in the ionosphere.

The marked decrease in peak ionospheric electron density is both scientifically significant and practically relevant. It indicates that existing empirical ionosphere models, using coefficients from outdated observations and solar activity indices like sunspot numbers or F10.7, may be inaccurate if long-term trends are ignored. Therefore, these models should be reviewed and updated.

Acknowledgments

The authors are grateful to the operators of numerous ionosonde stations for maintaining high-quality vertical sounding over many years. We also thank numerous colleagues for making ionosonde, solar indices and global surface temperature data available via national and international data centers. The editor thanks Stephan Buchert and an anonymous reviewer for their assistance in evaluating this paper.

Funding

This research did not receive any specific funding.

Data availability statement

Ionosonde data for stations Slough, Brisbane and Port Stanley from https://downloads.sws.bom.gov.au/wdc/iondata/medians/Damboldt/. Manually scaled foF2 data from the Japanese station Kokubunji (TO535/TO536) and Wakkanai (WK546) were obtained from the National Institute of Information and Communications Technology https://wdc.nict.go.jp/IONO/HP2009/ISDJ/index-E.html. foF2 data from the US-station Boulder (BC840) were obtained via the World Data Center C1 https://www.ukssdc.ac.uk/cgi-bin/wdcc1/secure/iono_data.pl (1958/07-2002/12, manually scaled) and from the Global Ionospheric Radio Observatory (https://giro.uml.edu/didbportal/didb-web-portal/ (2004/03-2022/12, automatically scaled). F30 solar radio flux index data from CLS via https://lasp.colorado.edu/lisird/data/cls_radio_flux_f30.

Temperature Anomaly data via the Japan Meteorology Agency (JMA) web site https://ds.data.jma.go.jp/tcc/tcc/products/gwp/temp/list/year_wld.html.

References

All Tables

All Figures

Figure 1 Long-term behavior of NmF2 over the German ionosonde station Juliusruh in comparison with solar radio flux indices F10.7 (a, b) and F30 (c, d). Panels a and c: Yearly means of NmF2 monthly medians and solar indices F10.7 (a) and F30 (b). Panels b and d: Sliding solar cycle averages (previous 11 years) of NmF2 estimates (NmF2*) and linear models of NmF2* (until 1984) run by sliding solar cycle averages of F10.7 (F10.7*) and F30 (F30*). The differences between observations and linear NmF2* models, ΔNmF2, are plotted in the lower part. scaled at the right side. The dashed lines indicate the RMS ± error (RMSE) range of the linear models taken as reference.

In the text

Figure 2 Long-term behavior of NmF2* over the ionosonde stations Slough/ Chilton from 1952 to 2021 (a, b) and Boulder from 1958 to 2022 (c, d). Panels a and c: yearly means of NmF2 and solar index F30. Panels b and d: Sliding solar cycle averages (previous 11 years) of NmF2 estimates (NmF2*) and linear models of NmF2* (until 1984) run by sliding solar cycle averages of F30 (F30*). The differences between observations and corresponding linear NmF2* models, ΔNmF2 are plotted in the lower part, scaled at the right side. Dashed lines indicate the RMS ± error range.

In the text

Figure 3 Long-term behavior of NmF2 over Kokubunji from 1958 to 2022 (a, b) and Wakkanai from 1950 to 2023 (c, d). Panels a and c: yearly means of NmF2 and F30. Panels b and d: NmF2* estimates and linear models of NmF2* (until 1984) as a function of F30*. The differences between observations and corresponding linear NmF2* models, ΔNmF2, are presented in separate plots scaled at the right side. Dashed lines indicate the RMS ± error range of the linear models taken as reference.

In the text

Figure 4 Long-term behavior of NmF2 over the Australian ionosonde stations Canberra from 1952 to 2022 (a, b) and Hobart from 1950 to 2022 (c, d). Panels a and c: yearly means of NmF2 and F30. Panels b and d: NmF2* estimates and linear models of NmF2* (until 1984) as a function of F30*. The differences between observations and corresponding linear NmF2* models, ΔNmF2, are presented in separate plots scaled at the right side. Dashed lines indicate the RMS ± error range of the linear models taken as reference.

In the text

Figure 5 Mean differences ΔNmF2 between estimates of NmF2* at ionosonde stations in the Northern (a and c) and Southern (b and d) hemispheres and corresponding models LM(NmF2*) as a function of F30*. Panels a and b zoom in on the early phase of measurements until 1982. Panels c and d show differences between the linear NmF2* models and observations over the full period 1963–2022, including linear regression lines starting in 1980.

In the text

Figure 6 Temperature anomaly (TA) at the Earth’s surface at both hemispheres (red dashed line). Linear regressions of the TA are fitted by the EDA anomaly averages (absolute at panels a and b and percentage at panels c and d) computed from the combined data of ionosonde stations Boulder, Juliusruh, Slough/Chilton Kokubunji and Wakkanai in the Northern hemisphere (panels a and c) and by the ionosonde stations Brisbane, Canberra, Hobart and Port Stanley from the Southern hemisphere (panels b and d). To get a coincidence of linear fits for TA and EDA, absolute EDA averages have been scaled by a divisor as indicated and offset by several years indicated by TD in the plots.

In the text