Trends in the Airglow Temperatures in the MLT Region—Part 1: Model Simulations

: Airglow intensity-weighted temperature variations induced by the CO 2 increase, solar cycle variation (F10.7 as a proxy) and geomagnetic activity (Ap index as a proxy) in the Mesosphere and Lower Thermosphere (MLT) region were simulated to quantitatively assess their inﬂuences on airglow temperatures. Two airglow models, MACD-00 and OHCD-00, were used to simulate the O( 1 S) greenline, O 2 (0,1) atmospheric band, and OH(8,3) airglow temperature variations induced by these inﬂuences to deduce the trends. Our results show that all three airglow temperatures display a linear trend of ~ − 0.5 K / decade, in response to the increase of CO 2 gas concentration. The airglow temperatures were found to be highly correlated with Ap index, and moderately correlated with F10.7, with the OH temperature showing an anti-correlation. The F10.7 and Ap index trends were found to be ~ − 0.7 ± 0.28 K / 100SFU and ~ − 0.1 ± 0.02 K / nT in the OH temperature, 4.1 ± 0.7 K / 100SFU and ~0.6 ± 0.03 K / nT in the O 2 temperature and ~2.0 ± 0.6 K / 100SFU and ~0.4 ± 0.03 K / nT in the O1S temperature. These results indicate that geomagnetic activity can have a rather signiﬁcant e ﬀ ect on the temperatures that had not been looked at previously. nT in the O1S temperature. This indicates that the increase of CO 2 gas concentration tends to amplify airglow temperature variation, yielding larger F10.7 and Ap index trends in the O 2 and O1S temperatures, and smaller trends in the OH temperature. The O 2 temperature shows that it has the largest response to the inﬂuences of F10.7, the Ap index variations, and CO 2 increase. Our simulation results show that the airglow temperatures would decrease linearly in response to the increase of anthropogenic gas emissions. Further, geomagnetic activity can a ﬀ ect airglow temperatures and kinetic temperature rather signiﬁcantly, in addition to what is already known about how the solar cycle variations can a ﬀ ect the atmospheric temperature. More work needs to be done in the area of observations, modeling and data analysis to learn more about the e ﬀ ects on the atmosphere caused by the increase of anthropogenic gas emission, solar cycle variation and geomagnetic activity.


Introduction
Global warming has been recognized as being responsible for extreme weather events in recent years. The cause of global warming has been attributed to be of anthropogenic origin. The pioneering work by [1], who did a numerical study to investigate the impact on the atmospheric structure caused by the CO 2 and CH 4 concentration change in the lower atmosphere, has revealed that gas concentration and temperature would change in response to the change of anthropogenic gas emissions. In particular, the study predicted a 10 K (50 K) cooling in the mesosphere (thermosphere) and a 30-50% decrease in the number densities of major gas species and atomic oxygen at 300 km, for a doubling of CO 2 and CH 4 gas concentrations.
Several studies have attempted to identify global warming footprint from various sources. For example, satellite drag measurements have been used to find if there exists any decreasing trend in neutral density [2][3][4][5]. See [6] for a comprehensive comparison of these studies. One big challenge for the trend analysis was that solar cycle variation would cause neutral density variation as well. Needless to say, in order to quantify global warming of anthropogenic origin, we need to separate the effects by natural variabilities from those by anthropogenic gas emissions.
As shown in [7], airglow is very sensitive to global change, which is not surprising, because airglow is known to be sensitive to atmospheric conditions. However, as pointed it out by [8], "using airglow measurements to deduce the change by the increase of anthropogenic gas emissions presents some challenges . . . " The challenges come from the fact that there are variations in the measurements by other influences (be it dynamical or natural variability) that need to be distinguished from the variation caused by the increase of anthropogenic gas emissions. The study by [7] also concluded that solar variability and dynamics are two main factors that affect the airglow variations in the measurements. Therefore, the current work does not include the dynamics intentionally.
That we need to quantitatively assess the variations of airglow intensities and temperatures by the known influences was the motivation for the current study and previous numerical studies by [8][9][10]. As pointed out in [8], the advantage of a numerical approach is that we can do controlled simulations to better assess the individual impact of the influences. The author of [9] did a simulation study to investigate how airglow intensities and Volume Emission Rates (VERs) would change in response to hypothetical scenarios of CO 2 doubling, CO 2 decreased by half and CO 2 increased by 10%. The work by [8] simulated the time series of airglow intensities, VERs, and airglow peak heights due to the influences of CO 2 concentration increase, using realistic CO 2 measurements and solar cycle variation for over a solar cycle time period. The recent study by [10] simulated the time series of the aforementioned airglow quantities due to CO 2 increase, solar cycle variation and geomagnetic activity over a much longer time period of 55 years.
In the current study, we use OHCD-00 and MACD-00 models, which employ the NRLMSISE-00 model [11] as a reference model [12] for the simulations of airglow temperatures. This work is a follow-up of [10], focusing on airglow temperatures but with two major differences-the current work uses the NRLMSISE-00 model, whereas [10] used MSISE-90 in their OHCD-90 and MACD-90 models, and the current work expands the years of interest from 1960 to 2019. The airglow intensity-weighted temperatures are obtained from the NRLMSISE-00 kinetic temperatures, weighted by the airglow VERs in the airglow models. The variations induced by the influences (CO 2 gas concentration increase, F10.7, and Ap index variation) are simulated with OHCD-00 and MACD-00 models, and then the trends in the airglow temperatures by these influences are deduced.
The current work is focused on the model simulations of airglow temperatures under the influences of CO 2 increase, F10.7 and Ap index variations. We will defer the comparisons to observation to a follow-up paper ("Trends in the Airglow Temperatures in the MLT Region -Part 2: Observations"), in which we will compare our model simulations and the reanalysis of some other observations to the Sounding of the Atmosphere using Broadband Emission Radiometer (SABER) observations. The paper is organized as follows. The models and data sources (CO 2 , F10.7, and Ap index) are described in Section 2. Results are presented in Section 3. Discussion is in Section 4 and Conclusions are in Section 5.

Models and Data Sources
In this study, airglow intensity-weighted temperatures of OH (8,3), O 2 (0,1) atmospheric band, and O( 1 S) greenline emissions were simulated for a location at latitude 18 N and longitude 290 E, the same as in [10]. Two airglow models, OHCD-00 and MACD-00 models, were used with the values of CO 2 concentration level, F10.7 and Ap index input to the models for a time period from year 1960 to year 2019.
The OH Chemistry Dynamics (OHCD) model considers essential chemical reactions for OH airglow [13,14], whereas the Multiple-Airglow Chemistry Dynamics (MACD) model considers chemical reactions for O 2 atmospheric bands and O( 1 S) greenline [15,16]. The extension of the airglow models indicates which reference model is used. For example, the OHCD-90 and MACD-90 models use the MSISE-90 model [17] as a reference model, and they were used in [10] to simulate airglow intensities and VERs. Detailed information about the different versions of the OHCD and MACD models and the references therein can be found in [12]. Suffice to say that the OHCD-00 and MACD-00 models use NRLMSISE-00 as the reference model to obtain the vertical profiles of major gas species' number density, atomic oxygen, atomic hydrogen and atmospheric temperature. Vertical profiles of minor gas species' number densities of O 3 , OH(v = 0), OH(v = 8), HO 2 , O 2 ( 1 C), O 2 (0,1) and O( 1 S) are obtained from the OHCD-00 and MACD-00 models, based on the chemical equilibrium requirement.
Three scenarios have been simulated. The first scenario, S1, simulates the time series of airglow temperatures under the influence of F10.7 and Ap index variations, when the CO 2 gas concentration level is kept fixed at its value in year 1960. The second scenario, S2, simulates the time series of airglow temperatures under the influence of the CO 2 gas concentration increase, when F10.7 and Ap index values are kept fixed at their values in year 1960. The third scenario, S3, simulates the time series of airglow temperatures under the influence of F10.7, Ap index and CO 2 changes.
The CO 2 gas concentration (in ppm) was obtained from the NOAA's website and the F10.7 index values (in Solar Flux Unit, SFU) were taken from NASA's website. The F10.7 index was used as a proxy for the 11-year solar cycle variation. The Ap index (in nT) values were obtained from the World Data Center for Geomagnetism, Kyoto website. These data from 1960 to 2019 were further averaged to produce the annual means.
The solar cycle and geomagnetic activity variations of the atmospheric species' number densities and temperature are from the F10.7 and the Ap index dependence of these quantities output from the NRLMSISE-00 model. Note that the F10.7 and Ap variations have not been carried below 110 and 90 km, respectively, where their coefficient became insignificant in the MSISE-90 or NRLMSISE-00 models. We should also point it out that the MSISE-90 and NRLMSISE-00 are empirical models, meaning that the incorporated observations inherently contained the signatures of solar cycle variation and geomagnetic activity. So even when the F10.7 is not explicitly "turned on" below 110 km in the model, the solar cycle influence is in the data implicitly. Also, the advantage of using NRLMSISE-00 over MSISE-90 is that it contains more and recent data with more coverage in locations and data sampling. The change of the atmospheric gas concentrations and temperature due to the increase of the CO 2 gas concentration were found by using the analytic expressions listed in Table 1 of [8], following the procedure described in [8]. Figure 1 shows the OH(8,3) airglow intensity-weighted temperature (OH temperature hereafter) as a function of year for the three scenarios. Figure 1a shows the temperatures (axis on left) and F10.7 values (axis on right), and Figure 1b shows the temperatures (axis on left) and Ap index (axis on right). Our simulation results of S1 (red line with dots) show a strong anti-correlation between the OH temperature and Ap index. The correlation between the OH temperature and F10.7 is not as obvious. In some years, it showed a strong positive correlation during the solar maximum in solar cycle 20 (1968-1970), solar cycle 21 (1978-1982) and solar cycle 23 (1999-2002), and in some years, especially during the solar minima, it showed an anti-correlation. The temperature varies between 193.9 and 197.9 K, a 4.0-K variation over the 60 years. The S2 scenario (blue line with upright triangles) shows that the OH temperature decreases linearly with the increase of CO 2 concentration. The temperature change is about 3.2 K in the 60 years, indicating a −0.53 K/decade trend. For the S3 scenario (green line with squares), the temperature varies between 192.2 and 196.8 K, a 4.6-K variation. In looking at the S1 and S3 results, we can see that they share the same features, but the S3 results are shifted downward due to the CO 2 increase, and the separation between S1 and S3 temperatures gets larger in time, when the influence of CO 2 increase gets larger.

Results
The O 2 (0,1) airglow intensity-weighted temperature (O 2 temperature) for the three scenarios is shown in Figure 2. Figure 2a shows a moderate correlation of the O 2 temperature with F10.7. In general, the O 2 temperature follows the overall F10.7 variation, and a close correlation is observed during the solar maximum in solar cycles 20 (1967-1970), 22 (1988-1991) and 23 (2000-2002). Figure 2b shows that the O 2 temperature is highly correlated with Ap index. The temperature varies between 193.4 and 204.4 K, showing a 11-K variation for S1. The O 2 temperature also shows a decreasing trend, −0.55 K/decade, in response to the increase of CO 2 concentration for S2. For the S3 scenario, the temperature varies from 191.1 to 203.8 K, showing a 12.7-K variation. The presence of CO 2 increase amplifies the O 2 temperature variation with time.  Figure 3 shows a similar plot for O( 1 S) airglow intensity-weighted temperature (O1S temperature). The O1S temperature displays a similar response to the influences of F10.7, Ap index, and CO2 increase as the O2 temperature for the three scenarios. The O1S temperature overall follows the F10.7 variation, with an exception that occurred during the solar maximum in solar cycle 21 (1979-1981). This is also seen in the O2 temperature. The O1S temperature also shows a high correlation with Ap index. The O1S temperature in S1 changes within a 6.9-K range. The increase of CO2 concentration leads to a decreasing temperature with a trend of −0.52 K/decade. For the S3 scenario, our results show that the temperature fluctuates in a 9.0-K range. The results shown in Figures 1-3 indicate that the airglow temperatures will decrease in response to the increase of CO2 gas concentration with a magnitude of ~0.5 K/decade. Our results also indicate that the O2 temperature shows a larger response to the F10.7 and Ap index variations than the O1S temperature and both temperatures show a positive correlation with F10.7 and Ap index. The OH temperature has the smallest response and is found to be anticorrelated with F10.7 and Ap index. In addition, the increase of CO2 gas concentration tends to amplify airglow temperature variation.  The O1S temperature displays a similar response to the influences of F10.7, Ap index, and CO 2 increase as the O 2 temperature for the three scenarios. The O1S temperature overall follows the F10.7 variation, with an exception that occurred during the solar maximum in solar cycle 21 (1979-1981). This is also seen in the O 2 temperature. The O1S temperature also shows a high correlation with Ap index. The O1S temperature in S1 changes within a 6.9-K range. The increase of CO 2 concentration leads to a decreasing temperature with a trend of −0.52 K/decade. For the S3 scenario, our results show that the temperature fluctuates in a 9.0-K range. The results shown in Figures 1-3 indicate that the airglow temperatures will decrease in response to the increase of CO 2 gas concentration with a magnitude of~0.5 K/decade. Our results also indicate that the O 2 temperature shows a larger response to the F10.7 and Ap index variations than the O1S temperature and both temperatures show a positive correlation with F10.7 and Ap index. The OH temperature has the smallest response and is found to be anticorrelated with F10.7 and Ap index. In addition, the increase of CO 2 gas concentration tends to amplify airglow temperature variation.
To see how well the airglow temperatures are correlated with either F10.7 or Ap index, we did a linear regression fit of the airglow temperatures as a function of F10.7 and as a function of Ap index using Equations (1) and (2), respectively.
(1)  To see how well the airglow temperatures are correlated with either F10.7 or Ap index, we did a linear regression fit of the airglow temperatures as a function of F10.7 and as a function of Ap index using Equations (1) and (2), respectively.  Figure 4 shows the simulated airglow temperatures as a function of Ap index, along with the fits for S1. They all display a very high correlation with the Ap index, with the correlation coefficients 0.97 and the fitting lines fit the data very well. The OH temperature shows an anticorrelation, whereas the O 2 and O1S temperatures show a positive correlation with Ap index. The study by [10] reveals that major gas species' number density and kinetic temperature ratios are anticorrelated with the Ap index in the OH airglow region, e.g., Figure 2 of [10]. In the O 2 and O1S airglow region, the major gas species' number density ratios remain anticorrelated with the Ap index, but the temperature ratio is found to be positively correlated with the Ap index. Our airglow temperatures, which are NRLMSISE-00 kinetic temperature weighted by the airglow VERs, show that the Ap index trends are similar to the trends in kinetic temperature, in that the sign of the trend is altitude-dependent. The simulated airglow temperatures, however, do not show such high correlation with F10.7, as shown in Figure 5. The data are scattered along the fitting lines. That being said, the airglow temperatures do show a moderate correlation with F10.7, as indicated by the correlation coefficients varying between 0.49 and 0.69. The OH temperature decreases with increasing F10.7, whereas the O 2 and O1S temperatures increase with increasing F10.7.
The results of the linear regression analysis for each of the airglow temperatures for the three scenarios are listed in Table 1. For S1 and S3, the airglow temperatures were fitted as a function of F10.7 and as a function of the Ap index. For S2, the temperatures were fitted as a function of CO 2 gas concentration. The correlation coefficient, R, is also listed to show the goodness of the fit for each fitting. Based on the correlation coefficients, our results show that the airglow temperatures are well correlated with the Ap index, and moderately correlated with F10.7 for S1. The presence of CO 2 in S3 decreases the fidelity of the fits with F10.7 or the Ap index. The O1S and O 2 temperatures for S3 still show a very high correlation coefficient for the Ap index trend, but the OH temperature is affected the most, with the R value changing from~0.996 for S1 to~0.545 for S3. This implies that it is difficult to identify an Ap index trend in the OH temperature if the CO 2 trend has not been removed first. As for the O 2 and O1S temperatures, the Ap index trends are found to be a little bit larger when the CO 2 gas concentration increases. For S2, the airglow temperatures all show a CO 2 trend with a slope of −0.03, with the O 2 temperature having a slightly larger slope (−0.035). Using the trend equations for S2, it should be fairly easy to infer the airglow temperature change in response to the change of CO 2 gas concentration.  To see how well the airglow temperatures are correlated with either F10.7 or Ap index, we did a linear regression fit of the airglow temperatures as a function of F10.7 and as a function of Ap index using Equations (1) and (2), respectively. Atmosphere 2020, 10, x 6 of 12 Figure 4 shows the simulated airglow temperatures as a function of Ap index, along with the fits for S1. They all display a very high correlation with the Ap index, with the correlation coefficients ~0.97 and the fitting lines fit the data very well. The OH temperature shows an anticorrelation, whereas the O2 and O1S temperatures show a positive correlation with Ap index. The study by [10] reveals that major gas species' number density and kinetic temperature ratios are anticorrelated with the Ap index in the OH airglow region, e.g., Figure 2 of [10]. In the O2 and O1S airglow region, the major gas species' number density ratios remain anticorrelated with the Ap index, but the temperature ratio is found to be positively correlated with the Ap index. Our airglow temperatures, which are NRLMSISE-00 kinetic temperature weighted by the airglow VERs, show that the Ap index trends are similar to the trends in kinetic temperature, in that the sign of the trend is altitude-dependent. The simulated airglow temperatures, however, do not show such high correlation with F10.7, as shown in Figure 5. The data are scattered along the fitting lines. That being said, the airglow temperatures do show a moderate correlation with F10.7, as indicated by the correlation coefficients varying between 0.49 and 0.69. The OH temperature decreases with increasing F10.7, whereas the O2 and O1S temperatures increase with increasing F10.7.   The results of the linear regression analysis for each of the airglow temperatures for the three scenarios are listed in Table 1. For S1 and S3, the airglow temperatures were fitted as a function of F10.7 and as a function of the Ap index. For S2, the temperatures were fitted as a function of CO2 gas concentration. The correlation coefficient, R, is also listed to show the goodness of the fit for each fitting. Based on the correlation coefficients, our results show that the airglow temperatures are well correlated with the Ap index, and moderately correlated with F10.7 for S1. The presence of CO2 in S3 decreases the fidelity of the fits with F10.7 or the Ap index. The O1S and O2 temperatures for S3 still show a very high correlation coefficient for the Ap index trend, but the OH temperature is affected the most, with the R value changing from ~0.996 for S1 to ~0.545 for S3. This implies that it is difficult to identify an Ap index trend in the OH temperature if the CO2 trend has not been removed first. As for the O2 and O1S temperatures, the Ap index trends are found to be a little bit larger when the CO2 gas concentration increases. For S2, the airglow temperatures all show a CO2 trend with a slope of ~ −0.03, with the O2 temperature having a slightly larger slope (−0.035). Using the trend equations for S2, it should be fairly easy to infer the airglow temperature change in response to the change of CO2 gas concentration.

Discussion
Several studies have performed a trend analysis on temperature measurements to determine if there existed a linear trend in the temperatures and/or any solar cycle dependence of the temperatures. For instance, [18] found a trend of −0.89 ± 0.55 K/decade and of 4.2 ± 0.9 K/100SFU in the OH airglow temperature from 1988 to 2015, but noted that the multiple linear regression analysis that included linear and solar cycle terms was not enough to capture the observed long-term dynamics, and that there seemed to be a trend break in the data. The work by [19] found a cooling trend of 2.2 ± 0.9 K/decade in the OH airglow temperature from 2000-2012. The OH temperature in Svalbard from 1983 to 2013 showed a temperature response of 3.6 ± 4 K/100SFU to solar influence, and a linear trend of −0.2K ± 0.5 K/decade [20]. A cooling trend of 4 ± 2 K/decade at 90 km in the meteor radar temperature observations from 2001-2012 [21]. The OH temperature at South Pole

Discussion
Several studies have performed a trend analysis on temperature measurements to determine if there existed a linear trend in the temperatures and/or any solar cycle dependence of the temperatures. For instance, [18] found a trend of −0.89 ± 0.55 K/decade and of 4.2 ± 0.9 K/100SFU in the OH airglow temperature from 1988 to 2015, but noted that the multiple linear regression analysis that included linear and solar cycle terms was not enough to capture the observed long-term dynamics, and that there seemed to be a trend break in the data. The work by [19] found a cooling trend of 2.2 ± 0.9 K/decade in the OH airglow temperature from 2000-2012. The OH temperature in Svalbard from 1983 to 2013 showed a temperature response of 3.6 ± 4 K/100SFU to solar influence, and a linear trend of −0.2K ± 0.5 K/decade [20]. A cooling trend of 4 ± 2 K/decade at 90 km in the meteor radar temperature observations from 2001-2012 [21]. The OH temperature at South Pole from 1994 to 2004 showed a temperature response of 4 ± 1 K/100SFU, and a statistically insignificant linear trend of 1 ± 2 K/decade [22]. Most of the linear trends in these studies were found to show a cooling trend, but a large variation with values ranging from 0.2 K to 4 K exist among these studies. These studies all found an F10.7 trend with a magnitude of~4 K/100SFU. The trends from these studies are listed in Table 2. We also note that none of these studies examined the temperature response to geomagnetic activity, so no information is provided. Our simulated airglow temperatures show a linear trend of~−0.5 K/decade, which is due to the increase of CO 2 gas concentration (see S2 in Table 1). Our cooling trend in the airglow temperatures is smaller than the trends found by [19,21], but is close to the values obtained by [18,20]. Although our cooling trends might be different from some of the studies, it should be noted that our linear trends from the simulation results are due to the increase of CO 2 gas concentration alone. In reality, temperature measurements might contain more information than what has been identified in the data. Our simulated airglow temperatures for S3 shows a trend of −0.69 ± 0.28 K/100SFU for the OH temperature, 4.1 ± 0.7 K/100SFU for the O 2 temperature and 2.0 ± 0.5 K/100SFU for the O1S temperature. The F10.7 trend in the O 2 temperature is very close to the values in the aforementioned studies, and the trend in the O1S temperature is smaller. Our OH temperature shows a F10.7 trend that is much smaller in magnitude, compared to the studies in Table 2, and is negative. This might be due to the fact that the NRLMSISE-00 model that outputs the kinetic temperature does not explicitly include the F10.7 effect below 110 km in their model, and that there might be regions like in the OH airglow region that were under sampled. Our OH temperature also shows a weak negative Ap index trend −0.1 ± 0.02 K/nT when compared to the trend in the O 2 temperature (~0.6 ± 0.03 K/nT) or O1S temperature (~0.4 ± 0.03K/nT). Since the Ap index effect is included explicitly above 90 km in the NRLMSISE-00 model, our results show that the airglow temperatures, especially the O 2 and O1S temperatures, can respond rather significantly to geomagnetic activity.
To minimize the contamination of wave dynamics, we did not use the wave feature in the models, and the inputs to the models were the annual mean values, so that variations of dynamic nature on a shorter time scale would be averaged out. The airglow models (OHCD-and MACD) are sensitive to the inputs like the number densities and temperature from an atmospheric reference model. The simulation results can be different if these parameters show differences in the atmospheric reference models. Previous work by [8,10] used the OHCD-90 and MACD-90, with the MSISE-90 as the atmospheric reference model, whereas this work used NRLMSISE-00 in the OHCD-00 and MACD-00 models. It would be interesting to see how big a difference there would be when we compare the results from these simulations. An update on [10] is planned in the near future.
As is well known, the Earth's thermal structure is determined by the absorption of the incoming and re-radiated solar radiation. In the thermosphere from 80 to 600 km, it is predominately heated by solar EUV radiation. The variations of EUV directly cause the thermosphere to vary on similar time scales and magnitudes. Since F10.7 is a proxy for solar EUV radiation, it is expected that the solar variability will directly affect the temperature in the region of interest.
Reference [23] used four decades of F10.7 and Ap index data to study the correlation between solar cycle variation and geomagnetic activity. Their study shows that F10.7 and the Ap index tracked each other well, but there were times when the Ap index lagged F10.7 from 1 to 4 years. That there is a correlation between F10.7 and the Ap index suggests that the influences of F10.7 and the Ap index on the atmosphere could be correlated as well. As such, we cannot do a linear regression of temperature with both F10.7 and the Ap index simultaneously, since it is inherently nonlinear and nonorthogonal. Instead, we found the F10.7 and Ap index trends separately, so that we could gain a better understanding of how the atmosphere would respond to these influences individually.
That F10.7 and the Ap index tracked each other is interesting and bears a closer look at how they are related to each other. As described in [24], a 10.7 cm solar radio flux measurement comprises a mix of a slowly varying component and the quiet sun background, and sometimes radio bursts. These components may each vary independently with time. The slowly varying component originates primarily in active regions of solar magnetic activity and has a broad spectral peak at about a 10 cm wavelength. So F10.7 values indicate the strength of solar magnetic activity, which in turns indicates the solar wind variability. On the other hand, the impact of solar wind variability on the atmosphere can appear through geomagnetic storms, magnetospheric substorms and changes of geomagnetic activity [25]. That the F10.7 values are closely related to the solar magnetic activity, and that solar wind variability can cause geomagnetic storms, which indicates that there should exist a close correlation between F10.7 and the Ap index, as supported by the study of [23]. Given that the F10.7 measurements comprise more than just the solar magnetic activity component, this might be the reason why solar cycle variation (F10.7) and geomagnetic activity (Ap index) are found to track each other well in some periods, but not always. Regarding how the temperature is affected by geomagnetic activity, it was suggested that the energy deposition from geomagnetic storm-related precipitating particles can lead to temperature variation in the MLT region [25].
The authors of [22] did a correlation analysis between MSIS-90 temperature and F10.7 variation at the OH airglow peak height at 87 km, and they found that the MSIS-90 model does not reproduce the solar cycle dependence of the observed OH airglow temperature at the South Pole Station. This is not surprising because "the F10.7 and Ap variations have not been carried below 110 and 90 km, respectively, where their coefficient became insignificant" in the MSISE-90 model [17]. In this study, we used the NRLMSISE-00 model, and found that the simulated airglow temperatures in general do track F10.7 variation loosely and track the Ap index closely. We highly recommend that the F10.7 and Ap index variations be carried below the current altitudes in the NRLMSISE-00 model, to explicitly include their variations in the model output quantities.
Recent studies using meteor radar and Microwave Limb Sounder (MLS) on Aura satellite have shown that the mesospheric neutral density would change in response to geomagnetic storms [26,27]. The work of [10] also shows that atmospheric neutral density tracks geomagnetic storms closely. On the other hand, the MLS mesospheric temperature does not seem to show a clear response to geomagnetic storms [27], but our airglow temperatures do display a high correlation with geomagnetic activity. Yet another study using SABER zonal mean temperature in the lower thermosphere, in the altitude range of 100-120 km, shows that it seemed to strongly correlate with the recurrent geomagnetic activity at periods of 9 days and 13.5 days, more so than with the solar EUV variability [28]. These conflicting results reveal that the actual mechanism of the Mesosphere Ionosphere Thermosphere (MIT) coupling during geomagnetic activity is still unknown. More work needs to be done in the area of observations and modeling, in order to better understand the mechanisms in the MIT coupling during geomagnetic activity.

Conclusions
Airglow intensity-weighted temperature variations in response to the influences of CO 2 gas concentration increase, solar cycle variation (F10.7 as a proxy) and geomagnetic activity (Ap index as a proxy) were simulated for a period of 60 years from 1960 to 2019. Two airglow chemistry dynamics models, OHCD-00 and MACD-00, were used in the simulations of OH(8,3) airglow, O 2 (0,1) atmospheric band and O( 1 S) greenline temperatures for three scenarios (S1: CO 2 gas concentration is fixed and F10.7 and the Ap index are allowed to change; S2: F10.7 and the Ap index are fixed and CO 2 gas concentration is allowed to change; S3: F10.7, the Ap index, and CO 2 gas concentration are allowed to change). Trends were deduced from the linear regression analysis of the simulated airglow temperatures.
Our S1 results show that the airglow temperatures are highly correlated with the Ap index and moderately correlated with F10.7. The OH temperature shows an anti-correlation with F10.7 and the Ap index. The trends are~−1 ± 0.24 K/100SFU and~−0.2 ± 0.02 K/nT in the OH temperature, 3.8 ± 0.5 K/100SFU and~0.51 ± 0.02 K/nT in the O 2 temperature and~1.7 ± 0.4 K/100SFU and 0.3 ± 0.01K/nT in the O1S temperature. Our S2 results show that all three airglow temperatures display a linear trend of~−0.5 K/decade in response to the increase of CO 2 gas concentration. Our S3 results show that the trends become~−0.7 ± 0.28 K/100SFU and~−0.1 ± 0.002 K/nT in the OH temperature, 4.1 ± 0.7 K/100SFU and~0.6 ± 0.03 K/nT in the O 2 temperature, and~2.0 ± 0.6 K/100SFU and~0.4 ± 0.03K/nT in the O1S temperature. This indicates that the increase of CO 2 gas concentration tends to amplify airglow temperature variation, yielding larger F10.7 and Ap index trends in the O 2 and O1S temperatures, and smaller trends in the OH temperature. The O 2 temperature shows that it has the largest response to the influences of F10.7, the Ap index variations, and CO 2 increase.
Our simulation results show that the airglow temperatures would decrease linearly in response to the increase of anthropogenic gas emissions. Further, geomagnetic activity can affect airglow temperatures and kinetic temperature rather significantly, in addition to what is already known about how the solar cycle variations can affect the atmospheric temperature. More work needs to be done in the area of observations, modeling and data analysis to learn more about the effects on the atmosphere caused by the increase of anthropogenic gas emission, solar cycle variation and geomagnetic activity.