© L. Perrone et al., Published by EDP Sciences 2017

1 Introduction

Long-term variations of ionospheric and thermospheric parameters are widely discussed in the literature especially in relation with the last deep and prolonged solar minimum in 2008–2009. Although ionospheric long-term trends are very small and have no practical importance, they are closely related to the upper atmosphere parameter variations and may serve as an indicator of the thermosphere long-term changes. The latter is very interesting and important as we live on the Earth surrounded by the neutral atmosphere. However, there are no direct observations of the thermospheric parameters compared on their duration to the ionospheric ones which are available for the period of 5–6 solar cycles and even longer at some ionosonde stations. The interest to long-term changes in the ionospheric and thermospheric parameters has been initiated by Roble & Dickinson (1989), Rishbeth (1990), and Rishbeth & Roble (1992) who predicted the ionospheric effects of the atmosphere greenhouse gas concentrations (mainly CO₂) increase. They have shown that even under the double CO₂ increase scenario (which we are very far from) the predicted ionospheric effects should be small. But their results have stimulated researchers to relate the observed ionospheric long-term trends to the thermosphere greenhouse cooling (Ulich & Turunen, 1997; Sharma et al., 1999; Alfonsi et al., 2001, 2002; Laštovička et al., 2008; Qian et al., 2008, 2009; Laštovička et al., 2012; Danilov & Konstantinova, 2013; Mielich & Bremer, 2013; Konstantinova & Danilov, 2015; Roininen et al., 2015). The mechanisms of the thermospheric and ionospheric trends may be different and serious contradictions with the CO₂ hypothesis confirm this (Perrone & Mikhailov, 2016). It is well-known that the ionospheric F-layer is strongly controlled by geomagnetic activity and nobody has denied yet the geomagnetic control concept of ionospheric long-term trends (Mikhailov, 2002). The analysis by Perrone & Mikhailov (2016) using all available (including recent ones) f_oF₁ and f_oF₂ observations on Slough/Chilton and Juliusruh ionosonde stations has confirmed that the geomagnetic control of the f_oF₂ and f_oF₁ long-term variations was still valid in the 21st century. Moreover the dependence on geomagnetic activity has become more pronounced and explicit after 1990.

Along with pure morphological analyses of the observed f_oF₂, f_oF₁, and f_oE long-term variations we have proposed a so called “self-consistent approach” to the analysis of the thermospheric and ionospheric long-term trends (Mikhailov & Perrone, 2016a). The idea of this approach is in using the observed f_oF₁ long-term variations to retrieve a consistent set of the main aeronomic parameters responsible for these f_oF₁ variations. Keeping in mind the same scheme of photo-chemical processes and common neutral composition in the daytime mid-latitude F₁ and F₂ regions it is possible to perform a simultaneous analysis of long-term variations in the two ionospheric region.

The main link in our approach is the method to extract thermospheric parameters and it needs a thorough testing. We have used a comparison with the excellent CHAMP/STAR neutral gas density (ρ) observations for this testing (Mikhailov & Perrone, 2016a). Neutral gas density is the integral characteristic which includes three retrieved (O, O₂, N₂) neutral concentrations as well as neutral temperature which is used to reduce ρ from heights of F₁-layer to the height of the CHAMP satellite. A comparison with CHAMP/STAR neutral gas density observations is the only opportunity to test the method using direct observational data. Another possibility is a comparison with the empirical models like MSISE00 (Picone et al., 2002) and JB2008 (Bowman et al., 2008) (both are used in this paper), but empirical models are climatologic ones describing average values and the priority should be given to the comparison with CHAMP/STAR neutral gas density observations.

With the new method to retrieve thermospheric neutral composition (O, O₂, N₂) and temperature Tex from routine f_oF₁ ionosonde observations the mechanism of f_oF₁ and f_oF₂ long-term variations (daytime, mid-latitudes) can be specified. Such analysis conducted with Slough/Chilton and Juliusruh observations (Mikhailov & Perrone, 2016a) has shown that f_oF₁ and f_oF₂ long-term variations are controlled via two channels: [O] and [O]/[N₂] variations. Both channels, in their turn, are controlled by solar and geomagnetic activity long-term variations.

Keeping in mind that our approach to long-term trends analyses is a new one and the results obtained with this method are not in the mainstream additional tests are needed using new observations.

Rome, with manually scaled ionosonde observations for the period of ∼5 solar cycles, was chosen for such testing. On one hand, Rome is a lower latitude station compared to Slough and Juliusruh with a different response to geomagnetic activity variations – only strong geomagnetic storms result in negative F₂-layer disturbances at Rome. On the other hand, long-term h_mF₂ variations have not been analyzed at Slough and Juliusruh using the retrieved thermospheric parameters and it would be interesting to check the geomagnetic control in h_mF₂ long-term variations. According to theory (e.g. Ivanov-Kholodny & Mikhailov, 1986), N_mF₂ and h_mF₂ are closely related via the F₂-layer formation mechanism and the geomagnetic control should be seen in the h_mF₂ long-term variations as well. Therefore, the aims of the present paper may be formulated as follows:

• to reveal long-term variations in monthly median f_oF₂, h_mF₂, f_oF₁, including the recent observations at Rome, and to check the existence of the geomagnetic control in their variations;

• to retrieve neutral composition (O, O₂, N₂) and temperature Tex from daytime monthly median f_oF₁ observations and to analyze the long-term variations in the retrieved thermospheric parameters in the light of the geomagnetic control concept;

• to analyze the role of solar and geomagnetic activity in the revealed ionospheric and thermospheric parameter long-term variations and to make a conclusion on their nature.

2 A method to extract ionospheric parameter long-term variations

A standard simple method using a regression of monthly median ionospheric parameter with an index F of solar activity $$f_{\text{reg}}=a_0+a_{1}F+a_2{F}^2\text{,}$$(1) is used to find monthly relative deviations $$\delta f=(f_{\text{obs}}-f_{\text{reg}})\text{/}{f}_{\text{obs}}\text{.}$$

Depending on the analyzed parameter indices of solar activity, F used in the regression may be different and this should be checked each time. Table 1 gives an example of such analysis applied to June noontime monthly median f_oF₂ and f_oF₁ values observed at Rome in 1957–2015. Indices of solar activity: monthly F_10.7, 3-month F_10.7, 12-month running mean F_10.7, 11-month running mean weighted F_10.7 (Mikhailov & Perrone, 2016a), and 12-month running mean sunspot number R₁₂ have been compared to find the best correlation coefficient.

Table 1 shows that all indices provide a good correlation but F_3mon and F_{11mon_w} are the best and they may be used in the further analysis.

It was also checked (the result earlier stressed by many researches) that an addition of Ap indices: either monthly or annually or smoothed Ap values to (1) does not improve the regression accuracy. Although F-layer parameters depend on geomagnetic activity, this dependence cannot be removed by a regression of this type and it was stressed repeatedly in our earlier publications (e.g. Mikhailov & Marin, 2000; Mikhailov, 2006). Furthermore, differently from our earlier approach (Mikhailov, 2002), here we use monthly median f_oF₂, f_oF₁, h_mF₂ values for individual months instead of annual mean ones. For further analysis monthly relative deviations δf, monthly Ap and F_10.7 indices should be smoothed using running mean weighted smoothing with an 11-year gate (only June values are used) $$P={\displaystyle \sum _{k=-5}^0}{p}_k\times (6+k)+{\displaystyle \sum _{k=1}^5}{p}_k\times (6-k)\text{.}$$(2)

The selection of summer months (June) is due to the following reasons. On one hand, due to a seasonal peculiarity of the thermospheric circulation, the geomagnetic control is the best seen in summer. On the other hand, the method to retrieve the thermospheric parameters (used for physical interpretation) can be applied only to summer conditions when F₁-layer is reliably observed by ground-based ionosondes.

3 Thermospheric parameter retrieval

We start with the thermospheric parameters as they will be further used for the h_mF₂ trend analysis. This step is a very important link in our approach. On one hand, the thermospheric parameters allow us to understand the mechanism of the ionospheric parameter long-term variations on the other hand, their long-term variations are interesting by themselves as they manifest thermospheric long-term trends widely discussed in the literature in relation with the thermosphere cooling due to CO₂ abundance increase in the Earth's atmosphere.

A new method to retrieve thermospheric parameters (Tex, O, O₂, N₂) and solar EUV flux with λ < 1050 Å from routine f_oF₁ ionosonde observations was proposed by Mikhailov & Perrone (2016a). The method is applicable only to summer months and around noon hours, when f_oF₁ is regularly and reliably observed, but even with these limitations the method has turned out to be useful for trend analyses (Mikhailov & Perrone, 2016b).

Observed f_oF₁ is the input information to the method, therefore its quality is crucial for the final results. Unfortunately, the quality of f_oF₁ measurements is different at different stations especially after the introduction of the automatic scaling of ionograms. Rome ionosonde observations have a long history and experience in ionogram scaling and such manually scaled ionospheric parameters can be used for our analyses.

The only direct way to test the efficiency of the method is to compare the retrieved neutral gas density with CHAMP/STAR neutral gas density measurements (http://sisko.colorado.edu/sutton/data.html) which have been conducted for many years under various geophysical conditions. Neutral gas density, ρ is an integral characteristic which includes the retrieved neutral concentrations (O, O₂, N₂) and temperature Tex. The latter is used to reduce neutral concentrations retrieved at F₁-layer heights to the height of CHAMP for a comparison. This is a strict type of a comparison which gives an objective estimate of the method efficiency.

A comparison with the empirical thermospheric models like MSISE00 (Picone et al., 2002) and JB2008 (Bowman et al., 2008) may serve as an independent check of our method.

Summer (June–July) daytime CHAMP/STAR observations in 2003, 2006–2008 in the European sector were used for testing. June–July 2003 was the period of elevated solar activity (monthly F_10.7 ∼ 130), and magnetically it was a very disturbed period with monthly Ap = 20–24. About half of tested days belonged to 2003 and some of them were strongly disturbed, with daily Ap up to 40–60. Another half of the tested days present low with F_10.7 = 76–73 (in 2006–2007) and extremely low with F_10.7 = 66 (in 2008) solar activity. Geomagnetic activity was low or slightly elevated for the second half of the selected dates. Observed CHAMP/STAR neutral gas densities were reduced to the locations of Rome (41.9N; 12.5E) and 12 LT using MSISE00 (Picone et al., 2002) thermospheric model and the following expression: $${\rho }_{\text{station}}={\rho }_{\text{CHAMP}}\frac{\text{MSISE0}{0}_{\text{station}}}{\text{MSISE0}{0}_{\text{CHAMP}}}$$

The height of CHAMP orbit changed from ∼400 km in 2003 to ∼335 km in 2008. The reduction height should be close to the satellite height to minimize possible errors due to the MSISE00 imperfectness. Three successive observations close to the latitude of the ionosonde station (after the reduction) were averaged to give the neutral gas density for our comparison. Normally, the reduced values of ρ at three points are close so the average value is reliable.

The retrieved from f_oF₁ neutral gas density ρ = m₁[O] + m₂[O₂] + m₃[N₂] does not include the contribution of He and N, therefore the observed neutral gas densities were corrected using MSISE00. Normally this correction is small (≤2%) at the reduction height, but it was applied.

Overall 48 comparisons between the retrieved and observed ρ have been done. We calculated the distribution of the R = ρ_cal/ρ_obs ratio where ρ_cal is the neutral density retrieved from the observed f_oF₁ values and ρ_obs are the corresponding CHAMP/STAR measurements reduced to the ionosonde location and 12 LT. Left panels of Figure 1 give the histograms of R. R_ave gives the average shift of the calculated ρ with respect to the observed ones. Middle panel gives a comparison with the JB2008 model (Bowman et al., 2008), and bottom panel of Figure 1 shows a comparison with the MSISE00 model (Picone et al., 2002).

Along with the histograms, we provide some statistical metrics (mean relative deviation (MRD), standard deviation (SD), and the bias with respect to the observed values) for a comparison between the retrieved neutral gas densities and two thermospheric models. The proposed method gives MRD = 12.3%, SD = 0.548 × 10⁻¹⁵ g cm⁻³ and the bias = −0.058 × 10⁻¹⁵ g cm⁻³, the JB2008 gives MRD = 13.1%, SD = 0.558 × 10⁻¹⁵ g cm⁻³ and the bias = −0.063 × 10⁻¹⁵ g cm⁻³, while MSISE00 gives MRD = 14.8%, SD = 0.574 × 10⁻¹⁵ g cm⁻³ and the bias = 0.207 × 10⁻¹⁵ g cm⁻³. The testing results show that the proposed method provides better accuracy than the modern empirical models. MSISE00 demonstrate a large positive bias, while JB2008 is well-centered and manifests less relative and SDs. It should be also stressed that the uncertainty of the retrieved neutral gas density coincides with the announced absolute uncertainty ±(10–15%) of the neutral gas density observations with the CHAMP satellite (Bruinsma et al., 2004). For a quick visual inspection, the plots of the retrieved and model ρ versus observed neutral gas densities are given in Figure 1 (right column). MRD and SD values along with the bias are given for a comparison. This graphical representation and the statistical results show that the retrieved densities are more centered with respect to the observations. Although the JB2008 model manifests a good distribution of R and is well-centered, both models have tails with large R values.

The undertaken testing shows that Rome f_oF₁ observations provide acceptable results in a comparison with CHAMP/STAR measurements and empirical models, therefore such f_oF₁ observations and the retrieved thermospheric parameters can be used for further long-term trend analyses.

Fig. 1 Left panels – distributions of R = ρcal/ρobs ratio for the retrieved neutral gas densities and those based on the JB2008 and MSISE00 models. Average Rave and the number of analyzed cases are given. Right panels – retrieved and model neutral gas densities versus the observed values. MRD and SD deviations along with the bias are given for a comparison.

4 F-layer parameter long-term variations

Long-term (δf_oF₂)_11y, (δf_oF₁)_11y, and (δh_mF₂)_11y variations for June 12 LT calculated with our method are given in Figure 2. Usual monthly hourly median f_oF₁ and f_oF₂ were used in our calculations, but the method of getting (δh_mF₂)_11y needs explanations. Only M(3000)F₂ routinely observed values are available for the whole pre-digisonde historical period of ionospheric observations.

Traditionally, these M(3000)F₂ are converted to h_mF₂ using the Shimazaki (1955) formula or more sophisticated expressions (e.g. Dudeney, 1974). Anyway, such h_mF₂ are not directly observed F₂-layer maximum heights but their approximation. An analysis by Ulich (2000) has shown that the overall inaccuracy of such conversion is about 20 ± 10 km depending on geophysical conditions, however such expressions are widely used in trend analyses (e.g. Bremer, 1998, 2001; Jarvis et al., 1998; Cnossen & Franzke, 2014; Roininen et al., 2015). It should be stressed that various improvements of the initial Shimazaki (1955) expression include the f_oF₂/f_oE ratio which by itself manifests long-term variations therefore the usage of such expressions for h_mF₂ trend analyses is questionable.

Upper panel of Figure 3 gives h_mF₂ long-term variations calculated from the observed monthly median M(3000)F₂ using the Shimazaki (1955) expression. F₂-layer maximum heights are seen to be unrealistically large, especially under solar minimum conditions, and corresponding long-term h_mF₂ trend calculated with such h_mF₂ also looks unreal which does not correspond to f_oF₂ long-term variations. Both parameters are related by the unique F₂-layer formation mechanism and their variations should agree with this mechanism. Therefore the thermospheric neutral composition (O, O₂, N₂) and temperature Tex, retrieved from June noontime monthly median f_oF₁, were used in the analytical expression for h_mF₂ obtained from a solution of the continuity equation for the electron concentration in the stationary daytime mid-latitude F₂-layer (Ivanov-Kholodny & Mikhailov, 1986, p. 43) $${\text{h}}_{\text{m}}{\text{F}}_2=\frac{H}{3}[\ln ({\beta }_1{[\text{O}]}_1)+\ln (H^2\text{/}0.54d)]+h_1$$(3) where H = kT/mg – scale height and [O]₁ concentration of neutral atomic oxygen at a fixed height h₁ (300 km in our case), β = γ₁[N₂] + γ₂[O₂] – linear loss coefficient at h₁,

d = 1.38 × 10¹⁹√ T/1000. Results of such analytical calculation of h_mF₂ are shown in the upper panel of Figure 3.

Theoretical h_mF₂ are seen to demonstrate quite different long-term variations with reasonable values both under solar minimum and maximum. The difference with the Shimazaki (1955) h_mF₂ values reaches 50–70 km rather than 20 ± 10 km (Ulich, 2000). A new global monthly median h_mF₂ empirical model by Shubin (2015) was used as a reference to compare with the theoretically calculated h_mF₂ variation (Fig. 3). This model is based on COSMIC radio-occultation observations and digisonde h_mF₂ data. Figure 3 (upper panel) shows that the empirical model by Shubin (2015) perfectly coincides with our theoretical h_mF₂ variations even in details (cf. the period 1971–1972 or 2012–2015). For this reason we used theoretical h_mF₂ (Eq. (3) + retrieved thermospheric parameters) for our long-term trend analysis (Fig. 2). For further discussion theoretical h_mF₂ variations are also compared to calculated ones (Fig. 3, bottom panel) when model MSIS-86 thermospheric parameters are used in equation (3).

According to theory due to the same scheme of photo-chemical processes and common neutral composition in the F₂ and F₁ regions, f_oF₁ manifests similar f_oF₂ long-term variations, while (δh_mF₂)_11y should demonstrate anti-phase variations with (δf_oF₂)_11y and (δf_oF₁)_11y as is seen in Figure 2. The correlation coefficient between (δf_oF₂)_11y and (δf_oF₁)_11y variations is 0.826, with the 99% confidence level according to Fisher F-criterion. The similarity in f_oF₁ and f_oF₂ long-term variation was stressed earlier (Mikhailov, 2008), but that time we did not have the required thermospheric parameters to explain this correlation.

Fig. 2 June noontime 11-year running mean weighted (δfoF2)11y, (δfoF1)11y, and (δhmF2)11y long-term variations at Rome.

Fig. 3 Long-term hmF2 variations at Rome, June 12LT: squares – calculated from observed M(3000)F2 values using the Shimazaki (1955) formula, asterisks – theory (expression 3), triangles – global empirical model by Shubin (2015). Bottom panel – the same theoretical hmF2 values but in a comparison with hmF2 when model MSIS-86 thermospheric parameters were used in equation (3).

5 Thermospheric parameter long-term variations

Monthly median f_oF₁, usable for our analysis, are available for ∼5 solar cycles (1957–2015). Thermospheric neutral composition (O, O₂, N₂), retrieved at heights of F₁-layer was reduced to 300 km altitude using the MSIS-86 neutral temperature T_n(h) profile with the retrieved Tex value. The same procedure was used in Section 3 when the retrieved neutral gas density was compared to CHAMP/STAR observations. The retrieved neutral composition and temperature are compared to the MSIS-86 thermospheric model (Hedin, 1987). On one hand, this is done for an additional control of the performance of our method. On the other hand, this is done to show that the retrieved and modelled thermospheric parameters manifest similar long-term variations indicating the origin of these variations.

The retrieved exospheric temperature Tex, neutral gas ρ and atomic oxygen [O] densities at 300 km altitude versus empirical MSIS-86 model values are given in Figure 4. To provide a correct comparison the model monthly Tex, ρ, and [O] medians were calculated for each June of all years using the observed 3-hour Ap and daily F_10.7 indices for each day of June and 12 LT. Along with the plots we provide some statistical metrics: the MRD, the bias with respect to MSIS-86 model and correlation coefficients between the retrieved and model values. The correlation coefficients are seen to be large for all parameters but there are some systematic shifts between the retrieved and model values: 6.2% for Tex, 16.5% for ρ, and 22.3% for [O], i.e. MSIS-86 gives larger values with respect to the retrieved ones.

To estimate the residual trends solar and geomagnetic activity effects should be removed from the retrieved parameter variations. The retrieved parameters manifest a good correlation with 3-month mean F_10.7 (Fig. 5, left panel), therefore it is possible to remove these solar activity effects and to check the residual variations. If they bear the geomagnetic activity effects they should be also removed. However, an addition of any Ap indices (monthly, annually or 11-year smoothed) to the regression practically does not affect the results. The obtained variations of δ were smoothed using 11-year running mean weighted smoothing (Fig. 5, right panel).

The residual 11-year running mean weighted smoothed δ manifest well-pronounced long-term variations (Fig. 5, right panels), which may be related to long-term variations in geomagnetic activity (see later).The magnitude of the revealed variations is small: ±1.5% for Tex, ±6% for ρ, and ±5% for [O]. They manifest both positive and negative phases and depending on the selected time interval the estimated trends will demonstrate different signs and magnitudes. Linear trends estimated over all available years (1962–2010) are very small (<0.5% per decade) for ρ and [O] and even less for Tex being statistically insignificant according to Fisher criterion. This means that practically all variations in the retrieved Tex, ρ, and [O] are due to solar activity variations.

Summarizing the results of undertaken analysis one may conclude that the retrieved Tex, ρ₃₀₀ and [O]₃₀₀ do not manifest any significant long-term trends estimated over a 57 year time period. However, it should be stressed that this conclusion has been obtained for June noontime mid-latitude conditions.

Fig. 4 Retrieved exospheric temperature, neutral gas density, and atomic oxygen concentration at 300 km versus MSIS-86 model values. See text for statistical metrics in the plots.

Fig. 5 Retrieved Tex, ρ and [O] at 300 km versus 3-month mean F10.7. Correlation coefficients are given. Solid lines – polynomial approximation (left panel). Right panels give 11-year running mean weighted smoothed δ obtained after the regression of retrieved Tex, ρ300 and [O]300 with 3-month mean F10.7. Straight lines – linear trends estimated over all years. Trends can be quantified using expressions given in the plots. The error bars indicate the ±1σ uncertainties.

6 Mechanism of f_oF₁, f_oF₂, h_mF₂ long-term variations

Using the retrieved thermospheric parameters and the analytical expressions for the daytime mid-latitude F₂-layer maximum parameters, it is possible to understand the mechanism of N_mF₁, N_mF₂ and h_mF₂ long-term variations. Equation (3) for h_mF₂ and the following expression (4) for N_mF₂, the latter taken from (Mikhailov et al., 1995, Appendix B) $${\text{N}}_{\text{m}}{\text{F}}_2\propto \frac{J_{\text{o}}{{[\text{O}]}_1}^{4/3}{T}_{\text{n}}^{-5/6}}{{{[{\text{N}}_2]}_1}^{2/3}}\propto {{[\text{O}]}_1}^{2/3}{\left(\frac{{[\text{O}]}_1}{{[{\text{N}}_2]}_1}\right)}^{2/3}{T}_{\text{n}}^{-5/6}\text{,}$$(4) were used for this analysis. The formation mechanism of the mid-latitude F₁-layer considered by Mikhailov & Schlegel (2003) shows that N_mF₁ is mainly controlled by the q(O⁺)/β ratio which is proportional to [O]/[N₂]. Ionospheric observations are taken directly from Rome ionosonde database (http://www.eswua.ingv.it/), and from the Lowell DIDBase via GIRO (Reinisch et al., 2004). Figure 6 gives 11-year running mean weighted Ap and F_10.7 indices, calculated (δf_oF₂)_11y, (δf_oF₁)_11y, (δh_mF₂)_11y, as well as the retrieved (δTex)_11y, and [O]_11y, ([O]/[N₂])_11y long-term variations at 200 km. Here we use 200 km height, which is closer to F₁-layer maximum, while expression (4) is invariant with respect to h₁ selection in the isothermal atmosphere.

Equation (4) shows that N_mF₂ depends not only on the O/N₂ ratio but on the absolute [O] concentration as well. For this reason, atomic oxygen turns out to be the main thermospheric parameter controlling daytime N_mF₂ variations. Molecular nitrogen is a passive species which follows the T_n variations, but it determines the recombination rate both in the F₁ and F₂ regions. For this reason along with [O]_11y we show [O/N₂]_11y variations in Figure 6. This ratio is usually used as an indicator of the thermosphere perturbation in the F₂-layer storm mechanism (e.g. Prölss, 1995, 2004 and references therein). It is seen that ups and downs in [O/N₂]_11y variations mainly coincide with the corresponding ups and downs in (δf_oF₁)_11y variations, as F₁-region is totally controlled by photo-chemical processes. This coincidence does not always take place for (δf_oF₂)_11y. The difference is due on one hand, to the dynamical processes that are important in the F₂-region and, on the other hand, to the fact that N_mF₂ depends also on the [O] absolute concentration (Eq. (4)). The [O/N₂]_11y ratio is seen to vary anti-phase with Ap_11y. Further, the O/N₂ ratio with some delay manifests anti-phase variations with atomic oxygen which, in turn, varies in-phase with (F_10.7)_11y. Geomagnetic activity is known to lag behind solar activity in solar cycles. Therefore, the (O/N₂)_11y ratio, which mainly follows geomagnetic activity variations, demonstrates some time lag with respect to ([O])_11y variations, the latter may be related to solar activity represented by direct solar indices like F_10.7. Ups and downs in ([O])_11y variations coincide with ups and downs in (F_10.7)_11y, the correlation coefficient is 0.974 with the 99% confidence level. All ups in ([O])_11y variations correspond to solar maxima and all downs to solar minima in solar cycles.

The (δh_mF₂)_11y long-term variations given in Figure 6 (right panel) may be explained using equation (3). This expression indicates that h_mF₂ is linearly related to neutral temperature while the dependence on [O] and β is weaker (via logarithm). For this reason, the h_mF₂ long-term variations should reflect the corresponding variations in the retrieved Tex. To demonstrate this, solar activity variations should be removed from both h_mF₂ and Tex variations using the 3-month regression with F_10.7. The residual variations should be smoothed using 11-year running mean weighted smoothing. The two plots in Figure 6 (right panels) manifest the similarity in the two variations. The correlation coefficient between (δh_mF₂)_11y and (δTex)_11y is 0.897, with the 99% confidence level.

Solar and geomagnetic activities are two channels, which provide the control of f_oF₁, f_oF₂ and h_mF₂ long-term variations, but via different aeronomic parameters. The atomic oxygen, [O] and the [O]/[N₂] ratio control f_oF₁ and f_oF₂. The neutral temperature, Tex controls h_mF₂ long-term variations. The rising phase (1965–1985, Fig. 6) in the [O] long-term variation corresponds to positive (δf_oF₂)_11y and (δf_oF₁)_11y deviations, while the falling phase in (1985–2008) results in negative (δf_oF₂)_11y and (δf_oF₁)_11y deviations. During the (1965–1985) period, [O] and [O]/[N₂] mainly work in one direction (both are increasing), thus f_oF₂ and f_oF₁ are also increasing (see Eq. (4)) and (δf_oF₂)_11y and (δf_oF₁)_11y are positive. After ∼1985 the decrease in [O] (dashes in Fig. 6, left bottom panel) becomes dominating and (δf_oF₂)_11y and (δf_oF₁)_11y become negative. In the end we have a well-pronounced negative f_oF₁ and f_oF₂ trends (see linear trend in Fig. 6) over the whole (1962–2010) period commonly discussed in the literature.

It is interesting to estimate the f_oF₂ trend, as it is without removing the geomagnetic activity effects and compare to other estimates. Figure 6 gives 0.05 for the total change in (δf_oF₂)_11y over 48 years. Accepting average f_oF₂ ∼ 7 MHz for June 12 LT at Rome we find a linear f_oF₂ trend ∼−0.007 MHz/year. This is much less than the f_oF₂ trend = −(0.020 − 0.015) MHz/year found by Laštovička et al. (2006) but it is closer to a recent estimate ∼−0.003 MHz/year (Mielich & Bremer, 2013). The analysis of (δf_oF₁)_11y long-term variations shown in Figure 6 gives a linear f_oF₁ trend ∼−0.001 MHz/year. This is, also, much less than the trend of 0.019 ± 0.011 MHz/year (Laštovička et al., 2012), moreover the sign of the trend is different. The absolute value of the estimated f_oF₁ trend is close to 0.0019 MHz/year (Bremer, 2008), but the sign is opposite. It should be stressed that the (1962–2010) period includes both periods of positive and negative (δf_oF₂)_11y and (δf_oF1)_11y trends (Fig. 6, left middle panel). This suggests that the final linear trend depends on the selected time interval (e.g. Mikhailov & Marin, 2001; Konstantinova & Danilov, 2015).

Long-term (δh_mF₂)_11y variations (as it was mentioned earlier) are controlled by neutral temperature long-term variations, i.e. they reflect (after the removal of solar activity effects) the variations in geomagnetic activity (Fig. 6, right column). As long as neutral temperature does not manifest any significant trend (Fig. 5), daytime h_mF₂ also does not demonstrate any significant long-term trend.

Fig. 6 11-year running mean weighted Ap and F10.7 indices along with (δfoF2)11y, (δfoF1)11y, (δhmF2)11y, and retrieved (δTex)11y, ([O]200)11y, ([O]/[N2]200)11y long-term variations at Rome. Curves – polynomial approximations. Numbers are given to identify the ups and downs with the corresponding downs/ups in the Ap index long-term variations. Straight line is the linear trend in (δfoF2)11y estimated over the (1962–2010) period.

7 Discussion

Ground-based ionosonde observations provide valuable data to analyze long-term trends not only in the main ionospheric parameters such as f_oF₂, h_mF₂, f_oF₁ but also in the thermospheric ones retrieved with our recently developed method.

This approach is based on the theory of the ionosphere formation which relates the ionospheric parameters to the thermospheric ones suggesting, for instance, that trends in N_mF₂ and h_mF₂ cannot be arbitrary, being related by the unique F₂-layer formation mechanism. The same can be said about f_oF₂ and f_oF₁ long-term trends. Following the pioneer paper by Rishbeth and Roble (1992) that reads: “The largest density changes occur in the F₁-layer near 180 km, with 50% increase at mid-latitudes”, some researchers find positive trends in the F₁-layer (Bremer, 2008; Qian et al., 2008; Laštovička et al., 2008). But according to theory daytime mid-latitude F₂ and F₁-layers are mainly controlled by common neutral composition and photo-chemical processes and they should manifest similar long-term variations. There are many possible reasons for these contradictions: poor quality of data at some stations, wrong data selection for trend analyses, poor methods of data development.

Figure 6 shows that f_oF₂ and f_oF₁ long-term variations manifest a negative trend estimated over the whole (1962–2010) period without the removal of geomagnetic activity effects (as the majority of trend researchers do). However this trend is not related to the increase in the CO₂ concentration in the Earth's atmosphere (as it is commonly accepted), but reflects the variations of neutral composition (mainly atomic oxygen) which in its turn are due to solar and geomagnetic activity variations.

According to theory the variations of mid-latitude daytime N_mF₂ are controlled by atomic oxygen [O] and [O]/[N₂] ratio variations (Eq. (4) and Fig. 6 left bottom panel); the effects of solar EUV are presumably removed via the 3-month mean F_10.7 regression. N_mF₂ also depends on neutral temperature (mainly via the temperature dependence for the reaction O⁺ + N₂ rate constant), but the T_n contribution to N_mF₂ long-term variations is small due to a small trend in Tex (Fig. 5).

According to Rishbeth (1990) the CO₂ global cooling of the upper atmosphere was expected to have a negligible small effect on N_mF₂: “…the “global cooling” is unlikely to have any significant effect on daytime values of N_mF₂, or critical frequency f_oF₂”. The same result follows from TIE-GCM model simulations by Cnossen (2014) which show “very clearly how little influence the increase in CO₂ concentration has had on f_oF₂”. Therefore, negative f_oF₂ and f_oF₁ trends (Fig. 6 left panel) estimated over the whole (1962–2010) period should be attributed to atomic oxygen decrease after ∼1990 (dashes in Fig. 6), which has overpowered a general [O]/[N₂] increase over the same period (solid line in Fig. 6). It should be mentioned that Danilov and Konstantinova (2014) were the first who have proposed to relate negative f_oF₂ trends with the atomic oxygen decrease in the upper atmosphere, but they prescribe this decrease to the intensification of the eddy diffusion. It was suggested by Mikhailov and Perrone (2016) that atomic oxygen long-term variation in (1961–2010) was due to the corresponding long-term variations in solar activity.

The analysis of h_mF₂ long-term variations have shown the following. Firstly, h_mF₂ values found from the M(3000)F₂ parameter using the Shimazaki (1955) formula cannot be used for trend analyses, at least under daytime summer conditions, especially in solar minimum because this formula strongly overestimates h_mF₂. For instance, the Shimazaki formula under such conditions gives h_mF₂ = 275–325 km (Fig. 3), while the Millstone Hill ISR observations http://madrigal.haystack.mit.edu/madrigal/ give h_mF₂ = 220–230 km. However, this formula directly or with some corrections is widely used for h_mF₂ trend analyses (Bremer, 1998; Jarvis et al., 1998; Roininen et al., 2015). In Rome the usage of this formula gives a strong negative trend which does not correspond to f_oF₂ long-term variations. For this reason we used a theoretical expression (3) with the retrieved values of thermospheric parameters. Although this analytical expression was obtained for a daytime mid-latitude stationary F₂-layer, without including the thermospheric wind effects, it gives reasonable h_mF₂ values. This is confirmed by a comparison with a modern empirical monthly median h_mF₂ model by Shubin (2015) included to the last version of IRI (Bilitza et al., 2017). To check this result we have used the MSIS-86 monthly median values of the thermospheric parameters in equation (3). Figure 3 (bottom panel) illustrates a good closeness of h_mF₂ calculated with the two sets of thermospheric parameters. A comparison of the two h_mF₂ variations gives the following statistical metrics: SD = 14 km, MRD = 3.8%, and the bias = −5.4 km. This result tells us that the retrieved thermospheric parameters are close to MSIS-86 model ones (see also Fig. 4) and also that they are controlled only by solar and geomagnetic activity represented by the solar F_10.7 and geomagnetic Ap indices which drive the MSIS-86 model. Therefore, after a proper removal of these two dependences one should not expect any significant residual trends either in the thermospheric parameters or in h_mF_2.

On the other hand the CO₂ concentration increase is a reality and some related effects should be seen, if not in N_mF₂, then in h_mF₂ long-term variations as the main effect of the CO₂ increase is the lowering of neutral temperature (see Eq. (3)). Under CO₂ doubled increase scenario Tex is predicted to decrease by 50 K (Roble & Dickinson, 1989). Under a 20% CO₂ increase in the Earth's atmosphere (Houghton et al., 2001) the expected Tex decrease is ∼10 K assuming a linear dependence. Taking 1000 K as an average solar cycle estimate for Tex, the expected Tex decrease is ∼1.0%. An average daytime h_mF₂ is ∼300 km then an h_mF₂ decrease is ∼3 km. This is close to the result of model (TIE-GCM) simulations by Cnossen (2014) for a ∼28% increase in CO₂ concentration which gave a fairly uniform decrease in h_mF₂ of about 5 km. Assuming that CO₂ increase has started 30–40 years ago one may expect a trend in h_mF₂ of ≤1 km/decade. Such trend hardly can be reliably detected keeping in mind the inaccuracy of (14–17) km of h_mF₂ determination with modern digisondes, obtained in a comparison with ISR observations (Chen et al., 1994). Anyway a trend of ≤1 km/decade strongly contradicts a 30 km h_mF₂ decrease obtained at Sodankylä over the (1957–2014) period and attributed to the CO₂ cooling (Roininen et al., 2015).

In the end some words concerning the results on the thermospheric parameters long-term trends are worth mentioning. As it was said earlier, under a 20% CO₂ increase in the Earth's atmosphere one may expect a ∼10 K Tex decrease. This gives cooling rate of ∼(3–4) K/decade analyzing a period of 30–40 years. This is close to an exospheric temperature trend of −1 to −2 K/decade estimated from satellite drag observations (Emmert, 2015) and much smaller than Tex trends inferred from ground-based ISR measurements: −18 K/decade for noontime exospheric temperature at Millstone Hill (Oliver et al., 2014), −60 K/decade at 350 km for daytime hours at Saint Santin/Nancay (Donaldson et al., 2010), −10 to −15 K/decade at F₂-layer heights for day-time hours at Tromso (Ogawa et al., 2014), and −20 K/decade at 350 km for daytime hours at Millstone Hill (Zhang & Holt, 2013). A recent analysis by Zhang et al. (2016) of Sondrestrom and Chatanika/Poker Flat ISR observations has shown that the high latitude long-term trend results are compared to those from the Millstone Hill mid-latitude dataset.

The retrieved thermospheric parameter long-term trends at 300 km, estimated over a 57 year time interval in the present analysis were shown to be small (<0.5% per decade) and statistically insignificant. Large negative trends in neutral temperature (supposing T_n = T_i) obtained with ISRs look absolutely unreal and probably are due to the IS method as it was discussed by Perrone & Mikhailov (2017).

Summarizing the results of our analysis it is possible to conclude that long-term variations of the thermospheric parameters retrieved from monthly median f_oF₁ observations have their origin in the Sun, i.e. they are of natural (not anthropogenic) origin and are controlled by long-term variations in solar and geomagnetic activity. After removing these dependencies the residual trends are very small and statistically insignificant. However, it should be stressed that these results on the thermospheric trends have been obtained on the limited observations – European region, June daytime hours. In future similar analyses should be conducted in other regions and for other seasons, providing reliable f_oF₁ observations are available.

8 Conclusions

The results of our analysis may be formulated as follows:

• Due to the same scheme of photo-chemical processes f_oF₁ manifests similar f_oF₂ long-term variations. The correlation coefficient between (δf_oF₂)_11y and (δf_oF₁)_11y variations is 0.826 under the 99% confidence level. In accordance with the F₂-layer formation mechanism the mid-latitude daytime f_oF₂ and h_mF₂ manifest anti-phase long-term variations as a reaction to geomagnetic activity.

• A comparison of neutral gas density retrieved from Rome f_oF₁ routine observations to CHAMP/STAR measurements has demonstrated satisfactory results: the proposed method provides better accuracy than the modern empirical models MSISE00, JB-2008 and the uncertainty of the retrieved neutral gas density coincides with the announced absolute uncertainty ±(10–15%) of the neutral gas density observations with the CHAMP satellite. That was an additional test of the developed method using new observations.

• There are periods of positive and negative long-term trends in exospheric temperature, neutral gas density, and atomic oxygen concentration retrieved from f_oF₁ observations for the period of ∼5 solar cycles (1957–2015), which are related to the corresponding periods in solar activity. After the removal of solar activity effects the residual trends estimated over the period of ∼5 solar cycles (1957–2015) are very small (<0.5% per decade) and statistically insignificant.

• Solar and geomagnetic activities are two channels which provide the control of f_oF₁, f_oF₂ and h_mF₂ long-term variations but via different aeronomic parameters. Atomic oxygen, [O] and [O]/[N₂] ratio control f_oF₁ and f_oF₂, while the neutral temperature Tex controls the h_mF₂ long-term variations. A linear trend in (δh_mF₂)_11y estimated over the (1962–2010) period is very small and insignificant reflecting the absence of any significant trend in neutral temperature.

• The f_oF₂ and f_oF₁ long-term variations obtained without removing the geomagnetic activity effects (as the majority of trend researchers do) demonstrate a negative trend estimated over the (1962–2010) period. However this trend is not related to the CO₂ concentration increase in the Earth's atmosphere but should be attributed to atomic oxygen decrease after ∼1990, which has overpowered a general [O]/[N₂] increase over the same period.

Acknowledgements

The authors thank for: the CHAMP density measurements available online at http://sisko.colorado.edu/sutton/data.html; the Juliusruh data are kindly provided by Leibniz institute of Atmospheric Physics – Field Station Juliusruh, Germany. SPIDR, IPS, and the Lowell DIDBase through GIRO to provide ionospheric data.

The editor thanks two anonymous referees for their assistance in evaluating this paper.

References

All Tables

All Figures

	Fig. 1 Left panels – distributions of R = ρcal/ρobs ratio for the retrieved neutral gas densities and those based on the JB2008 and MSISE00 models. Average Rave and the number of analyzed cases are given. Right panels – retrieved and model neutral gas densities versus the observed values. MRD and SD deviations along with the bias are given for a comparison.
In the text

	Fig. 2 June noontime 11-year running mean weighted (δfoF2)11y, (δfoF1)11y, and (δhmF2)11y long-term variations at Rome.
In the text

	Fig. 3 Long-term hmF2 variations at Rome, June 12LT: squares – calculated from observed M(3000)F2 values using the Shimazaki (1955) formula, asterisks – theory (expression 3), triangles – global empirical model by Shubin (2015). Bottom panel – the same theoretical hmF2 values but in a comparison with hmF2 when model MSIS-86 thermospheric parameters were used in equation (3).
In the text

	Fig. 4 Retrieved exospheric temperature, neutral gas density, and atomic oxygen concentration at 300 km versus MSIS-86 model values. See text for statistical metrics in the plots.
In the text

Fig. 5 Retrieved Tex, ρ and [O] at 300 km versus 3-month mean F10.7. Correlation coefficients are given. Solid lines – polynomial approximation (left panel). Right panels give 11-year running mean weighted smoothed δ obtained after the regression of retrieved Tex, ρ300 and [O]300 with 3-month mean F10.7. Straight lines – linear trends estimated over all years. Trends can be quantified using expressions given in the plots. The error bars indicate the ±1σ uncertainties.

In the text

Fig. 6 11-year running mean weighted Ap and F10.7 indices along with (δfoF2)11y, (δfoF1)11y, (δhmF2)11y, and retrieved (δTex)11y, ([O]200)11y, ([O]/[N2]200)11y long-term variations at Rome. Curves – polynomial approximations. Numbers are given to identify the ups and downs with the corresponding downs/ups in the Ap index long-term variations. Straight line is the linear trend in (δfoF2)11y estimated over the (1962–2010) period.

In the text