© L. Spogli et al., Published by EDP Sciences 2025

1 Introduction

The long-term variations in electron concentration of the Earth’s ionosphere are very small and have minimal practical importance in terms of impact on technological systems. However, they may suggest links between the lower (troposphere, stratosphere) and the upper layers (mesosphere and thermosphere/ionosphere) of the atmosphere. Although the long-term variations of the ionospheric parameters are closely related and reflect corresponding variations in thermospheric parameters, the trends in ionised and neutral components of the upper atmosphere are not identical and must be treated separately.

Currently, more than thirty years have passed since the first publications by Roble & Dickinson (1989) and Rishbeth (1990) devoted to studying, using model simulations, the possible impact of greenhouse gases (mainly CO₂) on Earth’s upper atmosphere. According to these studies, the increase of CO₂ abundance should result in cooling and subsiding of the thermosphere, i.e., a decrease of neutral gas density at fixed heights and a consequent shrinking.

Negative trends in the ionosphere have been found in electron concentration at the height of F2-layer (Alfonsi et al., 2001, 2002; Bremer, 1998, 2001, 2008; Danilov, 2006; Danilov & Konstantinova, 2013; Jarvis et al., 1998; Laštovička, 2013; Laštovička, 2017; Laštovička et al., 2012; Mielich & Bremer, 2013; Roininen et al., 2015; Sharma et al., 1999; Ulich & Turunen, 1997). Similar negative trends are also observed in thermospheric temperature and neutral gas density, as indicated by satellite drag and Incoherent Scatter Radar (ISR) observations (Elias et al., 2021; Emmert, 2015a, 2015b; Laštovička, 2023; Oliver et al., 2014; Ogawa et al., 2014; Zhang et al., 2011; Zhang and Holt, 2013a, and the references therein). The literature extensively discusses the contradiction between long-term trends in thermospheric temperature inferred from satellite drag measurements and those derived from ISR observations (Emmert, 2015a; Oliver et al., 2014; Zhang et al., 2011, and references therein).

As highlighted by Mikhailov et al. (2021), there is no direct experimental confirmation that the 29% increase in CO₂ in the Earth’s atmosphere observed at that time with respect to the 1960 level (317 ppm, source: https://www.co2.earth) can account for the observed ionospheric trends. Simulations using the Thermosphere-Ionosphere-Electrodynamics General Circulation Model (TIE-GCM) by Cnossen (2014) have clearly demonstrated “how little influence the increase in CO₂ concentration has had on foF2 and the daily magnetic perturbations associated with the Sq current system compared to changes in the Earth’s magnetic field”, aligning with theoretical and previous model predictions by Rishbeth (1990) and Rishbeth & Roble (1992).

Cnossen (2020) found that spatial patterns of trends in hmF2, NmF2, and total electron content indicate a superposition of CO₂ and geomagnetic field effects, with the latter dominating trends in the region of ∼50–20°N, ∼60°W to 20°E, such as the South Atlantic Anomaly. Model simulations consistently support the CO₂ concept (e.g. Cnossen, 2014; McInerney et al., 2024; Qian et al., 2011, 2008; Solomon et al., 2018). Utilising the Whole Atmosphere Community Climate Model-eXtended (WACCM-X), Solomon et al. (2018) demonstrates a cooling rate in ion temperature from the 1970s to 2000s of 2.7 K/decade, while McInerney et al. (2024) indicate a decrease of 3.7 K/decade for the same time range near the same altitude. The disparity arises from the use of different versions of the model. It is noteworthy that satellite drag observations do not substantiate the theoretically predicted magnitude of CO₂ thermospheric cooling, as indicated, e.g., by Emmert (2015a).

The neutral gas density at thermospheric heights is influenced by both neutral composition (primarily atomic oxygen) and exospheric temperature, T_ex. Emmerts’s (2015a) analysis lacked data on atomic oxygen, and the observed decrease in neutral density was ascribed to a decrease in T_ex. However, the study performed by Perrone & Mikhailov (2019) did not reveal any statistically significant decrease in the abundance of atomic oxygen over some decades. Therefore, the decrease in neutral density observed by Emmert (2015a) was appropriately attributed to a decrease in T_ex. The deducted decline in neutral temperature can be compared to the predicted cooling of the thermosphere resulting from the increase in CO₂ concentration.

Some papers on ionospheric and thermospheric long-term trends have shown the fundamental concept of the “geomagnetic control”. According to this, ionospheric and thermospheric long-term variations have a natural (not anthropogenic) origin associated with long-term variations in solar and geomagnetic activity (see e.g., Mikhailov & Perrone, 2016b; Perrone & Mikhailov, 2016, 2017).

The aim of this study is to investigate the long-term trends in the ionosphere and in the thermosphere using data from an ionosonde located at the headquarter of the Istituto Nazionale di Geofisica e Vulcanologia (INGV) in Rome, Italy. The objective is to assess whether the obtained results align with the natural or anthropogenic hypothesis.

The paper is organised as follows: Section 2 introduces the utilised data, encompassing ionosonde, solar flux, and geomagnetic data, along with the methods employed, including THERION and Fast Iterative Filtering. The adopted regression methods are also introduced in this section. Section 3 focuses on the results, covering the analysis of scales involved in parameter variability, trend retrieval, identification of drivers for the long-term variability of the examined parameters, and role of CO₂ concentration. Lastly, conclusions are drawn in Section 4, focusing on strengths and limitations of the proposed analysis.

2 Data and methods

This section provides details on the data, models, and signal processing techniques employed in our analysis. Specifically, the considered data include ionospheric parameters obtained from the ionosonde situated in Rome (41.8°N, 12.5°E); solar flux at 10.7-centimeter wavelength (F10.7); the Planetary Ap Index, serving as the indicator of the Earth’s overall geomagnetic activity; and thermospheric parameters derived from the THERmospheric parameters from IONosonde observations (THERION) method (Perrone & Mikhailov, 2018). THERION utilises ionosonde data in conjunction with solar and geomagnetic activity indices. The data span of the period is from 1976 to 2020, encompassing approximately four complete solar cycles.

Additionally, the methods employed include the Fast Iterative Filtering (FIF) algorithm (Cicone & Zhou, 2021) and its associated spectrogram, referred to as the “IMFogram” (Cicone et al., 2024). Furthermore, regression techniques such as linear regression and second-degree polynomial regression are utilised. More detailed explanations of the aforementioned data and methods can be found in subsequent sections.

2.1 Ionosonde data

Considering the extensive historical background and expertise in manually scaling ionospheric parameters from Rome ionosonde observations, data from the Rome Ionospheric Observatory, spanning the period from 1976 to 2020, are included in this study. This observatory has functioned as an ionospheric research facility since 1956, with the availability of digitalized ionograms starting from 1976. Measurements of the critical frequency of the F1 layer (foF1) and the F2 layer (foF2) from 10:00 to 14:00 Local Time (LT) are derived through manual interpretation of ionograms obtained from various ionosondes across different periods: (a) Union Radio Mark II recorded type ionosonde (Somoye, 2009) until 1979, (b) Bibl 128P type ionosonde (Bibl & Reinisch, 1978) from 1980 to 1995, (c) Barry ionosonde (Zuccheretti, 1998) from 1996 to 2004, and (d) AIS-INGV (Zuccheretti et al., 2003) since 2005.

To facilitate the identification of long-term trends, monthly median values were utilised to mitigate transient deviations in the examined parameters. The manually scaled ionosonde data constitute part of the comprehensive dataset accessible through the electronic Space Weather Upper Atmosphere (www.eswua.ingv.it) data portal managed by INGV (Upper Atmosphere Physics and Radiopropagation Working Group et al., 2020).

2.2 Solar flux and geomagnetic data

The F10.7 serves as a widely used solar activity indicator, akin to the sunspot number (as mentioned by Tapping, 2013). An analysis by Perrone & Mikhailov (2016) and Perrone et al. (2017) demonstrated that F10.7 correlates more strongly with the monthly median values of foF2 than other solar indices, such as the 12-month running mean sunspot number (R₁₂).

Each F10.7 value represents the cumulative emission at a wavelength of 10.7 cm, corresponding to a frequency of 2800 MHz, from all sources on the solar disk, collected during a 1-hour period centred around the specified epoch. It is worth noting that this quantity is actually a flux density rather than a flux, although it is conventionally referred to as such. The selection of the 10 cm wavelength range for monitoring solar activity is based on its high sensitivity to conditions in the upper chromosphere and the solar corona’s base. In contrast to numerous other solar indices, the F10.7 radio flux can be consistently and conveniently measured on a daily basis from the Earth’s surface, regardless of prevailing weather conditions. It is employed to drive both statistical and first-principles models of the ionosphere and thermosphere, with diverse applications across various domains (see, e.g., Elvidge et al., 2023 and references therein). The F10.7 measurements considered here are on a daily basis and are available through the National Research Council Canada in partnership with Natural Resources Canada (https://www.spaceweather.gc.ca/forecast-prevision/solar-solaire/solarflux/sx-en.php).

The Ap index is a measure of the general level of geomagnetic activity over the globe for a given day in universal time (UT) (Matzka et al., 2021; Rostoker, 1972). It is derived from measurements made at a number of stations worldwide of the variation of the geomagnetic field due to currents flowing in the Earth’s ionosphere and, to a lesser extent, in the Earth’s magnetosphere. The daily Ap-value is obtained by averaging the eight 3-hour values of ap for each day. To get Ap values, 3-hourly Kp-values must be converted to ap-values. The Ap values, along with other related indices of geomagnetic activity like the three-hour Kp and ap indices, are computed by the GeoForschungsZentrum Helmholtz Centre in Potsdam, Germany (https://kp.gfz-potsdam.de/en/data). The solar flux F10.7 and the geomagnetic disturbance index Ap are expected, based on physical principles, to be the primary indicators of variability in thermospheric and ionospheric parameters over long time scales.

2.3 THERION method

In the absence of routine thermosphere observations co-located with ionosondes, long-term values of the thermospheric parameters can be retrieved from the ionospheric observations themselves. This is the purpose of the THERmospheric parameters from IONosonde observations (THERION) method, designed to extract aeronomic parameters from ground-based ionosonde observations (Mikhailov et al., 2012). The THERION method has two versions: a general one, which may be used for any daytime conditions when foF2 and five plasma frequencies at 180 km height (fp₁₈₀) at (10, 11, 12, 13, 14) LT are available (Perrone & Mikhailov, 2018), and a version when, instead of fp₁₈₀, five foF1 values are used (Mikhailov & Perrone, 2016a). Plasma frequencies fp₁₈₀ read from automatically scaled fp(h) vertical profiles are available only for the last two decades at the best, while historical observations are needed, when only foF1 were available as F1-layer is systematically present in the summer daytime ionograms. For such a reason, this second version of the method (Mikhailov & Perrone, 2016a) is used in our analysis to retrieve thermospheric parameters under sunlit conditions (i.e., 12:00 LT) from ionosonde observations for the month of June for each year included in the study. As detailed in Mikhailov et al. (2012), the input solar and geomagnetic activity indices encompass F10.7 for the current and previous days, as well as the 81-day running mean F10.7 value, a background F10.7 value used in the solar EUV model, and the daily Ap index used in the MSIS-86 model (Hedin, 1987). The background F10.7 is a smoothly varying function whose values change from approximately 65 during solar minimum to about 120 during solar maximum.

The conceptual workflow of THERION is reported in Figure 1. As reported, from ionosonde observations and solar and geomagnetic indices as input the method provides a consistent set of noontime aeronomic parameters such as neutral composition ([O], [O₂], and [N₂]), exospheric temperature T_ex, total solar extreme ultraviolet (EUV) flux, and vertical plasma drift W mainly related to thermospheric winds.

Figure 1 Conceptual workflow of the THERION method as used in this study.

In this study, the THERION method, utilising monthly median values of ionospheric parameters as inputs, has been applied to obtain an extensive series of thermospheric parameters suitable for long-term trend analyses. The retrieved parameters are utilised to assess whether their variations align with natural or anthropogenic hypotheses. Specifically, [O] and T_ex values are directly considered in this study, while [N₂] and [O₂] are employed to obtain [O]/[N₂] ratio and the total thermospheric density ρ = m₁[O] + m₂[O₂] + m₃[N₂], where m₁ is the atomic mass of the atomic oxygen, and m₂ and m₃ the molecular masses of O₂ and N₂, respectively.

2.4 CO₂ concentration data

To investigate the role of the atmospheric CO₂ concentration ([CO₂]), we leverage the deseasonalised monthly mean recorded by the Mauna Loa (Hawaii) observatory managed by the National Oceanic and Atmospheric Administration (NOAA), which has the longest continuous record of direct atmospheric CO₂ measurements, dating back to 1958 (Keeling et al., 1976). We preferred to use the observations at Mauna Loa rather than the globally averaged at the surface, because the latter is available only since 1980. Despite some small differences are present in the two datasets, as Mauna Loa is located at an altitude of 3400 m in the northern subtropics, for the purposes of our analysis, based on the trend retrieval, we do not expect meaningful differences, also because the CO₂ is a well-mixed gas in the atmosphere.

2.5 Fast iterative filtering and IMFogram

To examine the spectral properties of the considered time series and to extract trends, we employ the Fast Iterative Filtering (FIF) technique (Cicone & Zhou, 2021). FIF, along with its precursor, Adaptive Local Iterative Filtering (ALIF) (Cicone et al., 2016), or its multivariate version (Cicone & Pellegrino, 2022), has been successfully applied in the ionospheric field for performing time-frequency, spectral, and multi-scale analyses (see, e.g., Ghobadi et al., 2020; Piersanti et al., 2018; Spogli et al., 2021; Urbar et al., 2023). The FIF technique offers an efficient decomposition of non-linear, non-stationary signals into functions oscillating around zero, known as Intrinsic Mode Functions (IMFs) in mathematics. Thus, every signal s = s(t) of that kind can be expressed in the following way:

$$s\left(t\right)=\sum _{i=1}^{N_{\mathrm{IMF}}}{\mathrm{IMF}}_i\left(t\right)+s_{\mathrm{res}}\left(t\right),$$(1)in which N_IMF is the total number of the found IMFs and s_res is the residual of the technique, identified as the trend of the signal. Each IMFs is a quasi-stationary component. Regarding the retrieval of the trend, the FIF parameter “ExtPoints” defines the maximum number of extrema (relative maxima) allowed in the remainder of the signal after extracting the IMFs. This determines the identified trend at the end of the process. In applying FIF technique, we manually selected this parameter to obtain a meaningful trend. Our choice, ExtPoints = 6, results in the trend that we believe better captures the variability and provides a non-trivial, non-linear representation. Bearing this in mind, in the proposed approach we assume the trend reconstructed with FIF a reliable proxy for the trend of investigated quantities over the considered period and by using the specified time window (start and end points).

We highlight that the FIF method assumes periodicity at the boundaries of the signal under consideration. To overcome this limitation, we employ the Signal Extension Algorithm as described by Stallone et al. (2020). This algorithm, through an anti-symmetric (anti-reflective) extension at the signal’s boundaries before decomposition, significantly reduces boundary errors.

The FIF method has demonstrated several advantages compared to other decomposition techniques (Cicone & Pellegrino, 2022; Cicone & Zhou, 2021). Notable among these advantages are its low computational complexity, making it the fastest technique of its kind; the assurance of obtaining unique decompositions; a comprehensive mathematical framework; and full adaptivity to the signal under investigation, eliminating the need to predefine either the number of components to be extracted or the basis to be used in the process.

As an illustration of the FIF method’s output, Figure 2 depicts the decomposition of the F10.7 time series from 1976 to 2020 (in black) using FIF. The red curves represent the various IMFs ordered by their descending quasi-stationary frequency values. The bottom plot represents the trend component. To emphasise FIF’s capability in capturing quasi-stationary modes effectively, Figure 3 showcases selected IMFs extracted from the F10.7 time series (also displayed in black), covering a longer period from 1950 to 2020. The extended dataset provides a clearer demonstration of the expected periodicities in the F10.7 data, such as the solar cycle-related ~11-years periodicity and its multiples at ~22 and ~44 years, which are distinctly identified through FIF, as indicated by the red curves in Figure 3.

Figure 2 Decomposition with FIF of the F10.7 time series considered in this study (black). Red curves represent the various IMFs, ordered according to the descending value of the quasi-stationary frequency. The bottom plot represents the trend.

Figure 3 IMFs (in red) of the F10.7 extended time series (black) identifying the ~11-years variation and the harmonics at ~22 and ~44 years.

Furthermore, a characterization based on FIF for identifying frequencies embedded in the signal provides the finest time-frequency resolution compared to more commonly known and utilised techniques such as the Fast Fourier Transform (FFT) and Discrete Wavelet Transform (DWT) (Urbar et al., 2023).

To achieve finer time-frequency identification, the IMFogram can be derived from the FIF decomposition (Cicone et al., 2024). The IMFogram is an analogue of the spectrogram that can be computed on the IMF decomposition to simultaneously identify the local frequency and amplitude information of a signal. In other words, the IMFogram offers precise resolution for scale identification in both the time and frequency domains. The basic formulation of the IMFogram and the comparison against FFT and DWT spectrograms are recalled in Urbar et al. (2023), while the complete theory behind is described in Cicone et al., 2024. In this study, we utilise the IMFogram to generate periodograms that are related to it. For the sake of simplicity, we will refer to these periodograms as IMFograms, by generalising its meaning. An example of the IMFogram obtained from the F10.7 time series is provided in Figure 4. In this figure, the period at ~11 years due to solar cycle variations is found to be the most energetic one. As reported in the thorough review by Bhowmik et al. (2023) while citing the pioneering work by Waldmeier (1935), it is known since long-time that the length of the cycle is correlated with the cycle’s amplitude, meaning that more intense cycles tend to have a shorter duration, while weaker cycles tend to be longer. This is clearly depicted by the F10.7 IMFogram, as a scale shift to slightly larger timescales. This serves to demonstrate the power of the joint use of the FIF and IMFogram techniques.

Figure 4 Left plot: Time series of F10.7 from 1976 to 2020 (black) and of its trend retrieved by FIF (red). Right plot: Corresponding F10.7 IMFogram in terms of temporal scale in years.

In summary, we employed the FIF/IMFogram to achieve the following objectives:

1. Understand the temporal scales involved in the long-term variation of the investigated ionospheric and thermospheric parameters.

2. Reconstruct the trends of the investigated thermospheric and ionospheric parameters, as well as of F10.7, Ap, and of [CO₂].

2.6 Regression analysis

The trends identified by FIF of the ionospheric and thermospheric parameters have been modelled through regression in the following way:

$$y=F\left(x_1,x_2\right),$$(2)in which y are the ionospheric/thermospheric parameters trends, x₁ = F10.7_trend is the trend of the F10.7 time series, and x₂ = Ap_trend is the trend of the Ap time series.

The investigated models are the linear regression, having the form:

$$y=F\left(x_1,x_2\right)=A+B{x}_1+C{x}_2,$$(2a)and the second-degree polynomial, having the form:

$$y=F\left(x_1,x_2\right)=y_0+\alpha x_1+\beta x_2+\gamma \left(x_1{x}_2\right)+\delta x_1^2+\epsilon x_2^2.$$(2b)

The regression analysis aims at determining whether variabilities in Ap and F10.7, which are expected to nicely correlate with long-term variations in the studied ionospheric and thermospheric parameters, cover the full range of variability, and, to which extent, other sources, like [CO₂], can contribute. In other words, we investigate if the long-term trends in the upper atmosphere over Rome are exclusively attributable to solar flux and geomagnetic influences when a correlation through regression analysis is performed. To achieve this, we conducted an additional assessment of the coefficient of determination R² (i.e., the square of the Pearson’s correlation coefficient R) between the trends in ionospheric and thermospheric parameters and their respective regression models. An R² value approaching 1 would indicate that the variability is predominantly accounted for by long-term solar flux and geomagnetic variations. Moreover, to further investigate potential unmodeled minor effects, we performed a linear fit on the time profile of the ratio between the observed trends and the regressed models. Furthermore, we provide some insights about the application of our method whenever CO₂ concentration data are considered.

3 Results and discussion

The results are provided in terms of: (i) analysis of the involved scales and trend identification, (ii) regression analysis to investigate if the long-term trends have further drivers besides the solar flux and geomagnetic variability, (iii) role of the carbon dioxide concentration in the proposed approach.

3.1 Analysis of the involved scales and trend identification

Before applying FIF to foF1 and foF2 data, addressing data gaps in these datasets is essential. This step is crucial because FIF requires uninterrupted data series as input. To achieve this, we have opted to fill the data gaps using a smoothing spline. While this may introduce spurious IMFs, it predominantly affects components with high-frequency variations that fall outside the scope of our work, as our primary focus is on trend reconstruction. This aspect will be considered when discussing the analysis of the relevant scale. Additionally, we explored alternatives such as linear interpolation and various values of the smoothing parameter p for the spline fitting (0.25, 0.5, 0.75, 1). We remind the reader that the spline smoothing parameter p controls the trade-off between the fidelity to the data and the smoothness of the spline: p = 0 produces a least-squares straight-line fit to the data; p = 1 produces a cubic spline interpolant, i.e., the spline fits the data points exactly; intermediate values of p produce a spline that balances fitting the data points closely while maintaining smoothness. From our analysis, we concluded that different values of p do not yield significant changes in the resulting IMFogram and the identified trend (not shown) for both foF1 and foF2. Consequently, we have chosen a smoothing parameter value of 0.5, and the reported results are based on this choice. Figure 5 reports the time series of foF2 (green dots) and the related smoothing spline adopted to fill the data gaps (black), while the red curve represents the identified trend. In the case of the considered signals and the use of the FIF technique, it is important to note that the identification of the trend depends on the lowest identifiable frequency and the length of the dataset. Concerning the former, as mentioned in Section 2.5, for each considered dataset we set the FIF options on the maximum number of extrema allowed in the remainder of the signal (see Cicone & Zhou, 2021 for details) to get the first trend which is not the trivial, linear trend. Additionally, it is worth noticing that the choice of starting and ending points within the dataset may result in slightly different trends. The dependence of the trends on the dataset selection has been pointed out in the work by Mikhailov et al. (2021). In that work, both increasing and decreasing phases have been considered. Different authors report either positive or negative trends depending on the duration of the period they examine. Furthermore, the lack of a standardised method to remove the effects of geomagnetic activity from foF2 long-term variations leads to different foF2 trends inferred for the same period when using different methods (Laštovička et al., 2006). This is a feature of our analysis affecting all the datasets under consideration and that must be taken into account when interpreting the results.

Figure 5 Time series of foF2 considered in this study (green dots) and the related smoothing spline adopted to fill the data gaps (black). Red curve represents the identified trend.

Figure 6 shows the results of the trend identification, depicted as a dashed red line, for all considered parameters: Ap (panel a), F10.7 (panel b), foF1 (panel c), foF2 (panel d), [O]/[N₂] (panel e), [O] (panel f), ρ (panel g), T_ex (panel h). The number of reconstructed IMFs for each quantity is reported in Table 1, illustrating how the considered parameters exhibit different ranges of modes, as revealed by FIF, indicating the complexity of the relationship among the considered quantities. The corresponding right plots depict the associated IMFogram in terms of the time scale in years. The colour bar is omitted, as the intensity is just supposed to indicate the relative intensity of the various temporal scales. Panel b reports the same content as Figure 4, illustrating trend identification with FIF and related IMFogram.

Figure 6 Same as Figure 4, but for all the investigated quantities: Ap (panel a), F10.7 (panel b), foF1 (panel c), foF2 (panel d), [O]/[N2] (panel e), [O] (panel f), ρ (panel g), Tex (panel h).

There are many parameters that determine the number of IMFs that must be carefully tuned to obtain a set of meaningful IMFs, which are detailed in Cicone & Zhou (2021). As mentioned in Section 2.5, a stopping criterion for FIF in the retrieval of the trend concerns the “ExtPoints” parameter, that we manually selected equal to 6 for each of the considered quantities in order not to get a non-trivial, i.e., linear trend.

It is noteworthy that the temporal plots shown in panels from a to d have finer temporal resolution compared to the thermospheric parameters plots (panels from e to h), as the latter values were retrieved only for the month of June and not for all the months of each year. Despite all the considered datasets featuring the main ~11-years periodicity (except for [O]/[N₂]), the IMFogram reveals a more complex structure for all of them.

Specifically, F10.7 (Fig. 6b) exhibits the most energetic IMFs in the range between 10 and 20 years, with a clear peak around 10 years and a slight increase over the years. The wiggles in the IMFogram are attributed to a not ideal behaviour of FIF in disentangling the various periodicities below 10 years, which conversely does not affect the trend reconstruction and the main features of the expected long-term variability in the solar cycle.

For Ap (Fig. 6a), the energy of the IMFs in the timescale range around 10 years is similar to harmonics at larger periods, and the energy is more distributed across all IMFs. Concerning the ionospheric parameters foF1 and foF2 (Figs. 6c and 6d, respectively), the most energetic components are again located at about 10 years and above. We remind here that the low-periods part of the foF1 and foF2 related IMFograms is influenced by the application of the smoothing spline to remove the data gaps and, therefore, must be considered unreliable for interpreting the results. However, our analysis is not leveraging at all on this part of the IMFograms and it is not affected by the method we adopted to fill the data gaps.

The parameters retrieved using THERION further confirm that the observed variability within the examined time frame results from multiple modes exchanging their energy, leading to a more intricate structure rather than a straightforward response to solar flux conditions influenced by the approximately 11-year solar cycle and its related harmonics. Also, the [O]/[N₂] ratio, for which the ~11-year was not clearly depicted in the time-series, presents a set of IMFs with periods ranging around 10 years. It is also worth noting that the behaviour of the [O]/[N₂] trend ratio deviates from that of other parameters, especially due to its noticeable increasing trend, which warrants further investigations. However, as a first consideration, the [O]/[N₂] ratio does not exhibit a clear ~11-year periodicity likely because it is derived from the ratio of two parameters that individually manifests this periodicity.

3.2 Drivers of the long-term variability

In this section, we present the findings related to the reconstructed trends of the examined thermospheric and ionospheric parameters, as well as the F10.7 and Ap indices, along with the associated regression analysis. The purpose of the regression analysis is to investigate whether the trends in Ap and F10.7 are the exclusive factors influencing the long-term changes in the studied ionospheric and thermospheric parameters whenever addressed with a regression problem on parameter trends. Figure 7 shows the trends reconstructed with FIF for all considered quantities, forming the basis for the results reported in this section.

Figure 7 Trends reconstructed with FIF of all considered quantities.

As described in Section 2.5, we model the trends identified by FIF of the ionospheric and thermospheric parameters by applying a regression having F10.7 and Ap as the independent variables, as per equation (2). The investigated functional forms are the linear and the second-degree polynomial models. Perrone & Mikhailov (2016), Mikhailov & Perrone (2016b) and Mikhailov et al. (2021) already showed that a linear regression is sufficient to explain long-term variations in the European sector. We also adopt second-degree polynomials as they were found to be even better in modelling the identified variations, as reported below. In general, the mixed terms in the regression models allow including the possible role of the interactions among the various terms. Moreover, others studies used a dependence on the square of some solar activity index (e.g., Emmert, 2015a; Zhang & Holt, 2013) and the same result has been obtained by Cnossen (2020) who showed that a multilinear regression model including an F10.7a² term (where F10.7a is the 81-day average of the F10.7 index) gives the most reliable upper atmosphere trends based on 66 years of model data.

The left panel of Figure 8 shows the foF1 trend and associated models with F10.7 and Ap derived through linear and second-degree polynomial fitting. A scatter plot depicting foF1 and its associated models is reported in the middle panel, while the right panel reports the time profile illustrating the ratio between the foF1 trend and the regressed models. The dashed lines correspond to the linear fits assessed for both curves. From the left panel of Figure 8, it is possible to deduce with the naked eye how the second-degree polynomial outperforms the linear model in depicting the foF1 trend reported in black. This is quantitatively confirmed by the R² values reported in the legend of the middle panel, which pass from 0.840 for the linear model to 0.996 for the second-degree polynomial model. The latter also aligns with the ideal conditions represented by the bisector, indicated by a black dashed line in the middle panel. Furthermore, the ratio between foF1 and its regressed model values, as shown in the right panel, further substantiates that the second-degree polynomial model accounts for all the variations in the long-term foF1 trend. In fact, when examining the time-dependent linear fit of this ratio, it remains nearly constant and consistently remains around 1 for the second-degree polynomial model, while some residual unmodeled effects are apparent in the red curve, which represents the linear regression model.

Figure 8 Left panel: foF1 trend (black) and related models obtained by applying a linear fit (red curve) and second-degree polynomial fit (green). Middle panel: Scatter plot of foF1 trend and related models obtained by applying a linear fit (red curve) and second-degree polynomial fit (green). The R2 of the models is reported in the legend, while the black dashed line indicates the bisector. Right panel: Time profile of the ratio between foF1 trend and the regressed models, according to the same colour code of the previous panel. The related dashed lines indicate the linear fit evaluated for the two curves.

Figures 9–13 report the same plot of Figure 8 but for the remaining investigated ionospheric and thermospheric quantities, i.e., foF2, [O], [O]/[N₂], ρ, T_ex, respectively. The results confirm the findings observed for foF1, indicating that the second-degree polynomial model outperforms the linear model, accurately capturing the variability in the trends of the investigated ionospheric/thermospheric quantities. In fact, all R² values associated with the second-degree polynomial model, as summarised in Table 2, consistently approach unity and are always larger than the corresponding values of the linear fit model. The sole exception is the [O] trend, which is featured by a R² of 0.929 for the second-degree polynomial fit. This warrants further investigation, as discussed further in this paper.

Figure 9 Same as Figure 8, but for foF2.

Figure 10 Same as Figure 8, but for [O].

Figure 11 Same as Figure 8, but for [O]/[N2] ratio.

Figure 12 Same as Figure 8, but for thermospheric density ρ.

Figure 13 Same as Figure 8, but for exospheric temperature Tex.

According to the theory, the same scheme of photo-chemical processes and common neutral composition in the F2 and F1 regions results in foF1 manifesting similar long-term variations as foF2. This is confirmed by our previous analysis (Perrone et al., 2017), which found a correlation coefficient of 0.83 between (δfoF2)11y and (δfoF1)11y variations, where monthly relative deviations (δfoF2)11y and (δfoF1)11y are smoothed using running mean weighted smoothing with an 11-year gate (Perrone & Mikhailov, 2016).

This means that foF1 and foF2 long term trends are similar because they depend on the same thermospheric parameters (neutral composition and temperature) and solar EUV radiation. Therefore, all long-term variations in foF2 and foF1 (used in the THERION method) will be reflected in the retrieved thermospheric parameters regardless of the nature of these foF1 and foF2 long-term variations.

Table 3 reports a summary of the results, presenting the coefficients of the second-degree polynomial fit with F10.7 and Ap for the investigated quantities. These coefficients are characterised by larger R² values compared to the linear fit. We recall here that the second-degree polynomial fit is given in the form $y_0+\alpha x_1+\beta x_2+\gamma \left(x_1{x}_2\right)+\delta x_1^2+\epsilon x_2^2$ as per equation (2b), in which x₁ = F10.7_trend and x₂ = Ap_trend.

It is noteworthy that the correlation coefficients obtained through the FIF analysis are similar to those calculated in previous studies, providing confirmation, albeit using a different methodology, of the corresponding results, as indicated in Table 4. As the works cited in Table 4 did not use any detrending techniques, this also suggests that the assumption that FIF-retrieved trends capture the entire long-term variability is physically sound, rather than being influenced by any periodicity longer than the 11-year cycle.

The obtained results show that the unexplained variability in ρ is below 1%, which is significantly less than the expected 14.8% from model simulations (Solomon et al., 2018) considering the actual concentration of CO₂. Similarly, the same amount of variability is observed for T_ex, a finding close to the 3% variability found at the same station in a previous study (Perrone et al., 2017). This level of variability is also consistent for the ionospheric parameters foF1 and foF2. However, in the case of [O] at 300 km, the remaining 7% of variability may be attributed to other factors, or it could reflect the inaccuracy of the input data or limitations of the applied method.

Similar results have been reported in Perrone & Mikhailov (2019), in which the unexplained variability of columnar [O] was reported as 7%. Additionally, Oliver et al. (2014) observed a very small negative trend of 0.0 ± 1.5% per decade in atomic oxygen density at 400 km altitude over the period from 1976 to 2013 at Millstone Hill.

In Mikhailov et al. (2021), correlation coefficients were computed for the Rome station in the period 1964–2020, focusing on the T_ex and ρ parameters retrieved from THERION and those approximated with a linear regression (Eq. (2a)). On the other hand, in Mikhailov & Perrone (2016b), correlation coefficients were calculated between the retrieved parameters (T_ex, ρ and [O]) and the corresponding values derived from the MSIS-86 thermospheric model (Hedin, 1987) over the period 1957–2015.

In summary, the variations in [O] account for 93%, while those in ρ and T_ex range between 99% to 100%, all attributed to solar and geomagnetic activity variations. Remarkably, the same results were obtained considering MSIS-86. Given that MSIS-86 is driven by solar F10.7 and geomagnetic Ap indices, it can be concluded that the variations in the retrieved parameter are mainly controlled by solar and geomagnetic activity.

3.3 Concentration of the carbon dioxide

In this section, we report about further analyses conducted to investigate how the role of [CO₂] can be directly highlighted by the proposed analysis. Figure 14 shows the deseasonalised monthly mean [CO₂] data from Mauna Loa (black), its trend reconstructed with FIF (red) and fit of the [CO₂] with a linear (blue) and second-degree polynomial (green) for the period 1976–2021. The R² of the two fits are reported in the legend. As reported by the fits on the [CO₂] trend, a linear increase is statistically meaningful (R² = 0.98) to describe the rise in carbon dioxide concentration. The ideal fit is given by a second-degree polynomial, whose R² is compatible with the unity. This functional behaviour must be taken into account when understanding the cause-effect relationship addressed with the proposed approach. Our analysis relies on regression, which implicitly assumes that correlation implies causation. While the causation is well consolidated by first principles for what concerns the forcing of which F10.7 and Ap are proxies, the causation chain linked with [CO₂] variations is not completely understood and consolidated. This aspect is discussed in detail later in this section.

Figure 14 Deseasonalised monthly mean [CO2] data (black), its trend reconstructed with FIF (red) and fit of the [CO2] with a linear (blue) and second-degree polynomial (green) for the period 1976–2021. The R2 of the two fits are reported in the legend.

We investigate if the [CO₂] trend may account for the missing variability highlighted in Section 3.2. To this scope, we include the [CO₂] trend as a third independent variable in the regression analysis based on the second-degree polynomial reported in Figures 8–13. Therefore, we compare the models with and without the [CO₂] trend. Specifically, the expression of a second-degree polynomial with three independent variables is the following:

$$y=F\left(x_1,x_2,x_3\right)=y_0+\alpha x_1+\beta x_2+\gamma x_3+\delta \left(x_1{x}_2\right)+\epsilon \left(x_1{x}_3\right)+\zeta \left(x_2{x}_3\right)+\eta x_1^2+\theta x_2^2+\iota x_3^2,$$(3)in which the extra independent variable is x₃ = [CO₂] trend. Figure 15 reports the output of such an analysis on [O] trend, which was the quantity featured by the larger missing variability (about 7%) when only F10.7 and Ap trends are considered as drivers. For this specific case, the inclusion of [CO₂] leads to meaningful improvements: the blue curves in Figure 15 visibly represent the [O] trend more accurately and R² improves significantly, increasing from 0.929 to 0.998. This better performance of the model with [CO₂] trend may suggest that [CO₂] is exactly the missing part in our model formulation.

Figure 15 Left panel: [O] trend (black) and related models obtained by applying the second-degree polynomial fit with (blue) and without (green) considering [CO2] trend. Middle panel: Scatter plot of [O] trend and related models obtained by applying the second-degree polynomial fit with (blue) and without (green) considering [CO2] trend. Right panel: Time profile of the ratio between [O] trend and the regressed models, according to the same colour code of the previous panel. The related dashed lines indicate the linear fit evaluated for the two curves.

The performance of the model with and without [CO₂] trends for all the ionosphere/thermosphere parameters are summarised in Table 5, in terms of the relative behaviour of R² and Root Mean Square Error (RMSE) evaluated on the scatter plots, as per middle panel of Figure 15. Additionally, Figure 16 reports the percentage increase in R² of the model with [CO₂] trend with respect to the model without [CO₂] trend (left plot) and the corresponding percentage decrease in RMSE (right plot). The improvement reported for [O] trend is confirmed for all the considered quantities, and the model with [CO₂] trend seems to provide the better performances, suggesting again that [CO₂] is exactly the missing part in our model formulation. The increase in R² for [O] is at the level of 7.5%, while for the other quantities it is below 0.5%, which is seemingly insignificant for an already well-performing model.

Figure 16 Left panel: Percentage increase in R2 of the model with [CO2] trend with respect to the model without [CO2] trend. Right panel: Corresponding percentage decrease in RMSE.

Despite promising, the outcome of this analysis must be carefully interpreted. In fact, it is assumed that we know a-priori that [CO₂] is a driver of the variability to move from correlation to causation. This is consolidated for F10.7 and Ap, but not thoroughly assessed for [CO₂], especially for what concerns the amount of possible variability it brings into play.

Thus, we further tested the significance of the [CO₂] trend in the unmodeled part of each second-degree polynomial model. To the scope, we perform a Fisher’s regression test on the residual from the second-degree polynomial model with F10.7 and Ap as independent variables against the [CO₂] trend for each considered quantity. An example of the input to the Fisher’s test is provided in Figure 17, which reports the scatter plot of the residual between [O] trend and its modelled counterpart based on F10.7 and Ap trends against [CO₂] trend in the period 1976–2021. The red line indicates the related linear fit.

Figure 17 Scatter plot of the residual between [O] trend and its modelled counterpart based on F10.7 and Ap trends against [CO2] trend in the period 1976–2021. The red line indicates the related linear fit.

Table 6 reports the results of the Fisher’s regression test for all the considered quantities. In the table, R² values relate to the linear fits against the [CO₂] trend, as the one reported in the red curve of Figure 17 for [O], while N is the number of data points, F is given by:

$$F=\frac{R^2(N-2)}{1-R^2},$$(4)and P(X > x) indicates the probability that the null hypothesis is true or, in other words, the probability that there is not a linear relationship between the [CO₂] trend and the unmodeled part of the second-degree polynomial model. If P(X > x) < 0.05, it is usually assumed that the linear relationship does exist. Therefore, the P(X > x) values reported in Table 6, which are all well above the threshold of 0.05, indicate that, for all investigated quantities, the unmodelled part of each trend is likely not dependent on the [CO₂] trend.

The further analyses conducted to investigate the role of [CO₂] reveal that, despite being seemingly very promising in explaining the residual variability, some doubts remain. Indeed, although the regression analysis suggests the residual variability to be actually brought by [CO₂], it can be not statistically meaningful. It should be also noted that, as confirmed by further tests (not shown), the effect of [CO₂] on the regression is not distinguishable to any other linear trend one may consider in the same period for whatever quantity. This is a known feature, as time can be used as a proxy for [CO₂] in interannual scales. Additionally, our results do not substantiate the theoretically predicted magnitude of [CO₂] thermospheric cooling. These aspects limit the possibility of drawing conclusions on [CO₂] and further studies are needed in that regard.

4 Conclusions

In this work, we investigate the nature of the long-term changes in the upper atmosphere at mid-latitude by leveraging on ionosonde data and corresponding ionospheric and thermospheric retrieved parameters under sunlit conditions. Our analysis is based on a regression of the trends identified through the application of the FIF technique (Cicone & Zhou, 2021). Ionospheric parameters (foF1, foF2) are obtained from ionosonde data recorded digitally at the Rome observatory from 1976 to 2020, while thermospheric parameters ([O], [O]/[N₂], ρ, T_ex) are retrieved through THERION method (Perrone & Mikhailov, 2016) under summer daytime conditions. Both linear and second-degree polynomial models are here tested for the regression analysis. The results indicate that the variability in the trends of the investigated quantities is almost completely covered by the variability in F10.7 and Ap trends with differences at or below the percent level or even less, except for [O], which displays a variability at a few percent. The second-degree polynomial model outperforms the linear model, with achieving a coefficient of determination R² close to 1 within a percent margin except for [O], where R² is 0.93. Although statistically validated, the significant increase in [O]/[N₂] observed in this study warrants further in-depth investigation in future research.

This leads to the following conclusions:

• The long-term variations observed in the ionospheric and thermospheric parameters investigated in this study seem to be almost entirely governed by the solar and geomagnetic activity long-term variations, when the proposed methodology is applied.

• According to the FIF-based approach, the long-term [CO₂] increase suggests to explain the long-term variability of upper atmosphere characteristics, which is not explained by solar and geomagnetic activity. This concerns especially the [O] trend. However, the impact of [CO₂] apparent in our study is statistically not significant and requires further investigation.

• The analysis depends on the length of the considered datasets, and additional data are required to validate and confirm our findings. The issue of the length includes also the selection of the time window over which trends are evaluated, which may result in different functional forms of the found trends.

• The results found in this study with a novel FIF-based approach are consistent with previous analyses.

Acknowledgments

This work is part of the “A Multidisciplinary Analysis of Climate change indicators in the Mediterranean And Polar regions” (MACMAP) project, a Strategic Department Project funded by Istituto Nazionale di Geofisica e Vulcanologia. The project focuses on studying climate evolution in Mediterranean and polar regions by analysing and integrating new and existing data from models, observations, and historical qualitative information. The authors are grateful to Antonio Cicone (University of l’Aquila and INGV, Italy) for his valuable help with the FIF algorithm and useful discussions. The authors are also grateful to Annalisa Cherchi (CNR, Italy), Antonio Guarnieri (INGV, Italy) and Simona Simoncelli (INGV, Italy) for their support with the CO₂ concentration data. The editor thanks two anonymous reviewers for their assistance in evaluating this paper.

Data availability statement

Rome ionosonde data are available at the eSWua website: www.eswua.ingv.it.

F10.7 data are available at the following link: https://www.spaceweather.gc.ca/forecast-prevision/solar-solaire/solarflux/sx-en.php.

Fast Iterative Filtering and all related software is available at www.cicone.com.

The Ap values are available at https://kp.gfz-potsdam.de/en/data.

The CO₂ data are available at https://gml.noaa.gov/ccgg/trends/mlo.html.

References

All Tables

All Figures

	Figure 1 Conceptual workflow of the THERION method as used in this study.
In the text

	Figure 2 Decomposition with FIF of the F10.7 time series considered in this study (black). Red curves represent the various IMFs, ordered according to the descending value of the quasi-stationary frequency. The bottom plot represents the trend.
In the text

	Figure 3 IMFs (in red) of the F10.7 extended time series (black) identifying the ~11-years variation and the harmonics at ~22 and ~44 years.
In the text

	Figure 4 Left plot: Time series of F10.7 from 1976 to 2020 (black) and of its trend retrieved by FIF (red). Right plot: Corresponding F10.7 IMFogram in terms of temporal scale in years.
In the text

	Figure 5 Time series of foF2 considered in this study (green dots) and the related smoothing spline adopted to fill the data gaps (black). Red curve represents the identified trend.
In the text

	Figure 6 Same as Figure 4, but for all the investigated quantities: Ap (panel a), F10.7 (panel b), foF1 (panel c), foF2 (panel d), [O]/[N2] (panel e), [O] (panel f), ρ (panel g), Tex (panel h).
In the text

	Figure 7 Trends reconstructed with FIF of all considered quantities.
In the text

Figure 8 Left panel: foF1 trend (black) and related models obtained by applying a linear fit (red curve) and second-degree polynomial fit (green). Middle panel: Scatter plot of foF1 trend and related models obtained by applying a linear fit (red curve) and second-degree polynomial fit (green). The R2 of the models is reported in the legend, while the black dashed line indicates the bisector. Right panel: Time profile of the ratio between foF1 trend and the regressed models, according to the same colour code of the previous panel. The related dashed lines indicate the linear fit evaluated for the two curves.

In the text

	Figure 9 Same as Figure 8, but for foF2.
In the text

	Figure 10 Same as Figure 8, but for [O].
In the text

	Figure 11 Same as Figure 8, but for [O]/[N2] ratio.
In the text

	Figure 12 Same as Figure 8, but for thermospheric density ρ.
In the text

	Figure 13 Same as Figure 8, but for exospheric temperature Tex.
In the text

	Figure 14 Deseasonalised monthly mean [CO2] data (black), its trend reconstructed with FIF (red) and fit of the [CO2] with a linear (blue) and second-degree polynomial (green) for the period 1976–2021. The R2 of the two fits are reported in the legend.
In the text

Figure 15 Left panel: [O] trend (black) and related models obtained by applying the second-degree polynomial fit with (blue) and without (green) considering [CO2] trend. Middle panel: Scatter plot of [O] trend and related models obtained by applying the second-degree polynomial fit with (blue) and without (green) considering [CO2] trend. Right panel: Time profile of the ratio between [O] trend and the regressed models, according to the same colour code of the previous panel. The related dashed lines indicate the linear fit evaluated for the two curves.

In the text

	Figure 16 Left panel: Percentage increase in R2 of the model with [CO2] trend with respect to the model without [CO2] trend. Right panel: Corresponding percentage decrease in RMSE.
In the text

	Figure 17 Scatter plot of the residual between [O] trend and its modelled counterpart based on F10.7 and Ap trends against [CO2] trend in the period 1976–2021. The red line indicates the related linear fit.
In the text