Simulated Long-term Evolution of the Thermosphere during the Holocene: 1. Neutral Density and Temperature

. In the previous work of Yue et al. (2022), the ionospheric evolution during the Holocene (9455 BC to 2015 AD) 10 was comprehensively and carefully investigated for the first time using the Global Coupled Ionosphere-Thermosphere-Electrodynamics Model developed at the Institute of Geology and Geophysics, Chinese Academy of Sciences (GCITEM-IGGCAS), driven by realistic geomagnetic fields, CO 2 levels, and solar activity derived from the ancient media records and modern measurements. In this study, we further quantify the effects of the three drivers on thermospheric neutral density and temperature variations during the Holocene. We find that the oscillations of solar activity contribute more than 80% of the 15 thermospheric variability, while either CO 2 or the geomagnetic field contributes less than 10%. The effect of CO 2 on the global mean neutral density and temperature is comparable to that of the geomagnetic field throughout the Holocene but is more significant after 1800 AD. In addition, thermospheric density and temperature show approximately linear variations with the dipole moment of the geomagnetic field, CO 2 , and F10.7, with only the linear growth rate associated with the geomagnetic field varying significantly in universal time and latitude. The increasing dipole moment and CO 2 cool and contract the 20 thermosphere, while solar activity has the opposite effect. The higher the altitude, the greater the influence of the three factors


Introduction 25
Global glaciers have been melting in the recent century due to climate warming, and this melting has been accelerating in the last 20 years, leading to rising sea levels and elevating natural disasters (Hugonnet et al., 2021;Zemp et al., 2019). The main cause of climate warming is the use of fossil fuels as a source of anthropogenic greenhouse gases (Tollefson, 2021). The simulations of Roble and Dickinson (1989) show that the increase in greenhouse gases warms the troposphere, but cools the thermosphere. The long-term trend studies in subsequent decades largely support this consensus (Laštovička, 2009;Laštovička 30 et al., 2006;Laštovička et al., 2008). Tropospheric warming seriously affects human life, while changes in the thermosphere https://doi.org/10.5194/egusphere-2023-233 Preprint. Discussion started: 14 February 2023 c Author(s) 2023. CC BY 4.0 License. affect various human-launched satellites, space stations, and spacecraft, such as the SpaceX Starlink satellites destruction event on 4 February 2022 (Dang et al., 2022;Lin et al., 2022), therefore is also relevant to human life. The current understanding of the thermosphere is based on modern satellite observations over the last ~70 years. The energy sources driving the variability of the thermosphere include mainly solar irradiance and geomagnetic activity generated by the interaction between the solar 35 wind and the Earth's magnetic field (Knipp et al., 2004). These two energy sources are responsible for thermospheric temporal variations on time scales ranging from minutes to decades. On longer time scales, however, the effects of greenhouse gases and the geomagnetic field must be taken into account (Laštovička et al., 2006). Several review papers have summarized the knowledge of thermospheric variations and their driving mechanisms (Laštovička, 2017;Laštovička et al., 2012;Qian et al., 2011;. 40 Long-term trends of neutral density at different altitudes between 200 and 600 km have been extensively investigated using satellite orbit measurements since the 1960s (Emmert, 2015;Emmert et al., 2004;Emmert et al., 2008;Keating et al., 2000;Marcos et al., 2005;Saunders et al., 2011). These studies suggested that the trend is mainly attributed to a dramatic increase in anthropogenic greenhouse gases and becomes stronger with increasing altitude. A summary of the trends derived from satellite orbit data can be found in Emmert (2015) and Solomon et al. (2018). Overall, the observed long-term trend of the 45 thermospheric neutral density at 400 km ranges from -2% to -5% per decade. These observed characteristics are qualitatively consistent with the model-predicted effects of increasing CO 2 concentrations. (Qian et al., 2006;Roble & Dickinson, 1989;Solomon et al., 2015). The effect of geomagnetic field strength and configuration on the thermosphere has also been paid enormous attention from observations and simulations Cnossen, 2014;Cnossen, 2022;Cnossen & Maute, 2020;Cnossen et al., 2012;Cnossen et al., 2011;Förster & Cnossen, 2013). However, none of these effects are as significant as the 50 effect of solar activity on the thermosphere. The amplitude of the solar-driven variation increases with height, by a factor of two for temperature and an order of magnitude for density in the upper thermosphere Solomon et al., 2019). Overall, these observations and simulations of the thermosphere are limited to the recent 100 years, while the "ancient" thermosphere has never been investigated on longer time scales such as during the Holocene. Therefore, we propose to reconstruct the paleo-thermosphere since the Holocene using the first principle numerical model driven by these indices 55 including solar activity derived from tree ring (Solanki et al., 2004), paleo geomagnetic field model (Korte et al., 2011), and greenhouse gas concentration derived from polar ice core (Lüthi et al., 2008). Based on the reconstructed simulations, we can understand in detail the evolution of the thermosphere over 10,000 years and the effects of changes in the geomagnetic field, CO 2 , and solar activity on the thermosphere, which can provide a foundation for the future effect of climate change on human life. 60 The rest of the paper is organized as follows. Section 2 will briefly describe the numerical model and driving parameters used in this study. Section 3 will show the simulation results and discuss them. Finally, we draw conclusions in Section 4.

Methodology
This study will use the Global Coupled Ionosphere-Thermosphere-Electrodynamics Model developed at the Institute of Geology and Geophysics, Chinese Academy of Sciences (GCITEM-IGGCAS) (Ren et al., 2009), which is the same as Yue et 65 al. (2022). This model self-consistently solves the energy, momentum, and continuity equations of neutrals and ions in altitude coordinate rather than pressure level coordinate between 90 and 600 km, and solves the electrodynamic equations using magnetic apex coordinate (Richmond, 1995) based on provided spherical harmonic coefficients of any dipole-dominated geomagnetic field, such as the International Geomagnetic Reference Field (IGRF) model (Alken et al., 2021). This model has a good performance that has been confirmed by several ionospheric and thermospheric weather and climate simulations (Ren 70 et al., 2011;Ren et al., 2020;Ren et al., 2010;Yue et al., 2022;Zhou et al., 2022). To remove the effects of initial conditions, each simulation is run for an interval of 15 days, and the final day results are used for further analysis.
Three drivers associated with geomagnetic field, CO 2 level, and solar activity during the Holocene, 9455 BC to 2015 AD, are used to drive the GCITEM-IGGCAS. These drivers have been summarized in Figure 1 of Yue et al. (2022). The geomagnetic field is CALS10k.2 (Constable et al., 2016) for the period 9455 BC to 1900 AD and IGRF after 1990. The geomagnetic field 75 has undergone complex and nonlinear changes during the Holocene, with a dipole moment variation of ~40%, much larger than the ~7% variation since 1900. The CO 2 concentration evolution that is derived from Antarctica Vostok and EPICA Dome C ice cores (Lüthi et al., 2008), Antarctica Law Dome ice cores (Macfarling Meure et al., 2006), and direct atmospheric measurement at Mauna Loa Observatory, Hawaii (Keeling et al., 1995). The CO 2 concentration increases roughly linearly from 250 ppm (parts per million) around 10,000 BC to 402 ppm around 2015 AD, with a major increase occurring after 1800. 80 The F10.7 index evolution converted from the tree ring derived sunspot number (SSN) (Solanki et al., 2004) and the group SSN (Hoyt & Schatten, 1998), and modern instrument measurement (Tapping, 2013), which reveals the long-term oscillations in solar activity with relatively high solar activity around 9000 BC and 1960 AD. A detailed description of the three drivers can be found in Yue et al. (2022).
In this simulation, four control runs (CR1-CR4) have been implemented, which is the same as Yue et al. (2022) (summarized 85 in their Table 1). CR1 is used to identify the effect of geomagnetic field variation on thermospheric evolution. CR2 is used to reveal the effect of CO 2 . CR3 is used to diagnose the combined effect of geomagnetic field and CO 2 . CR4 is used to determine the combined effects of geomagnetic field variation, solar activity, and CO 2 . In addition, the combined analysis of the simulations of CR3 and CR4 allows for discerning the effect of solar activity.  Figure 1 gives the corresponding F10.7 index, CO 2 level, and dipole moment of the geomagnetic field. It is clear that the global mean neutral density and temperature are mainly controlled by solar activity as they show essentially the same temporal evolution and oscillation as solar activity. The higher the solar activity, the larger the neutral density and temperature in general, with a relatively larger value around 9000 BC and 1960 AD. The global distribution of neutral density and temperature shows a remarkable feature during the dayside with two peaks on either side of the geomagnetic inclination equator, so-called 100 equatorial mass anomaly (EMA), like the equatorial ionization anomaly (EIA) in the ionosphere (Appleton, 1946;Balan et al., 2018), indicating a strong coupling between the thermosphere and ionosphere modulated by the geomagnetic field at low and middle latitudes (Hedin & Mayr, 1973;H. Liu et al., 2005;Huixin Liu et al., 2007;Raghavarao et al., 1993;Raghavarao et al., 1991). Comparing the simulations for the considered four years can basically reveal the influence of CO 2 , geomagnetic field, and solar activity on the evolution of the thermosphere. The F10.7 index and CO 2 level used to drive GCITEM-IGGCAS are 105 similar in 9005 BC and 7635 BC, except that the dipole moment is about 20% larger in 7635 BC than in 9005 BC. It is clear that the simulated neutral density and temperature for these two years are not significantly different in global distribution pattern and magnitude, which indicates that the ~20% change in dipole moment has a weak effect on the thermosphere at 400 km altitude. Comparing the simulations of 9005 BC and 3015 BC reveals that a ~25% reduction of the F10.7 index leads to a ~50% decrease in neutral density and a ~170 K decrease in neutral temperature. Furthermore, the cooling effect of CO 2 on the 110 thermosphere can be found by comparing the simulations of 9005 BC and 2005 AD. An increase of 110 ppm (~40%) of CO 2 concentration causes a temperature reduction of ~40 K and a density reduction of ~21%. In addition, we also checked the simulations at other UTs and during the June solstice as well. In summary, the thermosphere is primarily controlled by solar activity, with secondary controlled factors being CO 2 and geomagnetic field. Figure 2 shows the global mean for all UTs and grids of the neutral density and temperature as a function of altitude and year 115 (left column) and the zonal mean neutral density and temperature versus latitude and year at 400 km (right column) in March equinox during the Holocene. According to (Afraimovich et al., 2008), a latitude-dependent area weighting factor was used in the calculation of the global mean to make it more representative. As shown in Figure 2, both the global mean and zonal mean display significant oscillations throughout the whole Holocene, which is consistent with the oscillations of the solar activity whose relative change was marked by the red lines in the left column. When F10.7 reaches its relatively higher value before 120 8000 BC and in the recent century, a larger value of neutral density and temperature also appear. This feature is not clear in the altitude profile of the global mean of neutral density (top left panel) due to the span of about 10 orders of magnitude.

Results and Discussions
However, this does not affect the conclusion, after all, the fixed height of the zonal mean is more revealing of this feature.
Furthermore, the thermospheric neutral density and temperature also show significant long-term decreases from 6000 BC to 3500 BC and the most famous grand solar minimum (Usoskin et al., 2007), the Maunder Minimum between 1645 and 1715 125 (Eddy, 1976). Only a weak latitudinal variation can be seen in the zonal mean neutral temperature in the March equinox.

The effects of the three drivers on the evolution of thermospheric neutral density and temperature
In this section, the changes in thermospheric neutral density and temperature caused by the geomagnetic field, CO 2 , and solar activity will be diagnosed by subtracting the beginning year (9455 BC) of the simulations, as shown in Figures 3-5 for the results of the global mean and zonal mean. 130 Figure 3 shows the global mean neutral density profile variations in percentage caused by the three drivers in the left column.
The black, magenta, and red lines represent the relative changes of the diploe moment, CO 2 , and F10.7 index, respectively, during the Holocene. In general, the higher altitude, the larger the effect of the three drivers on the neutral density, because the neutral density decreases exponentially with altitude causing all effects to be amplified at higher altitudes. For the effect of the geomagnetic field, its nonlinear variation causes a nonlinear change in the neutral density, and a decrease in its intensity 135 represented by the weakening of the dipole moment (around 5500 BC) generally leads to an increase in the neutral density. In addition, an increase in the dipole moment would make the neutral density increase weaker, which might be related to the decrease in Joule heating due to the strong dipole moment (Cnossen et al., 2012;Cnossen et al., 2011;Glassmeier et al., 2004;Wang et al., 2017). When the dipole moment increases beyond ~3×10 22 Am 2 (in ~7500 BC and from ~1500 BC to ~1000 AD), the density will decrease in turn between 150 and 250 km. For the effect of CO 2 , the neutral density decreases during the 140 increase phase of CO 2 level (before ~8000 BC and after ~4000 BC). This is because greenhouse gases can cool and contract the thermosphere (Qian et al., 2011), as shown in Figure 4 for temperature reduction. In turn, when CO 2 decreases between ~8000 BC and ~4000 BC, the thermosphere neutral density increases. It is worth noting that the effect of CO 2 has been more significant since 1800 AD due to the much larger growth rate of CO 2 . For the effect of solar activity, the overall change in neutral density due to solar activity is more than 10 times larger than that of CO 2 and the geomagnetic field, so it is the dominant 145 factor in neutral density change. The neutral density has increased by more than 100% in the recent century and around 9000 BC, which corresponds to relatively greater solar activity. The right column of Figure 3 is the zonal mean results at 400 km.
The grey lines in panel (d) mark the latitude of the north and south magnetic poles for the corresponding year. The effects of the three drivers on the temporal evolution of the neutral density at all latitudes are similar to those characterized in the left column, with no significant latitude variations in the effects of CO 2 and solar activity, and a weak latitude variation in the 150 effect of the geomagnetic field, which is stronger at high latitudes. The region with greater geomagnetic field effects corresponds exactly to the magnetic pole locations, such as the south magnetic pole at ~64° around 5500 BC, implying the importance of the magnetic pole locations and further supporting the contribution of Joule heating in the polar region. Figure 4 shows a similar pattern as Figure 3 except that the neutral temperature variations are shown. Overall, the effects of CO 2 and solar activity on temperature are essentially the same as those on neutral density, only the temperature change is more 155 significant. The dramatic increase in CO 2 level in the past century has led to a decrease in global mean neutral temperature of more than 20 K, which is well in line with previous understanding (Cnossen, 2014;Roble & Dickinson, 1989). In addition, the global mean neutral temperature increases within 5 K due to CO 2 reduction between 8000 and 4000 BC. Solar activity remains the dominant factor in the neutral temperature variability, which leads to neutral temperature changes in the range of ±200 K with its own oscillations. Furthermore, the effect of the geomagnetic field on neutral temperature differs significantly 160 from the effect on neutral density, but can still be explained by the geomagnetic field structure, dipole moment strength, and Joule heating. The zonal mean temperature changes contributed by the geomagnetic field show a clear latitude variation. The neutral temperature increases ~24 K at south latitude ~65° around 5500 BC, caused by a reduction in the dipole moment resulting in stronger Joule heating around the south magnetic pole. Conversely, between 1500 BC and 1000 AD, the neutral temperature in the polar regions dropped by up to 22 K due to the weakening of Joule heating caused by the increase of the 165 dipole moment. In addition, the effect of the geomagnetic field is significantly weaker at north latitude ~10°, perhaps owing to tides in the lower atmosphere. Although the neutral temperature changes in the polar regions are large, the change in global mean neutral temperature due to the geomagnetic field is essentially within ±10 K at all altitudes during Holocene. The June simulations have also been carefully analyzed (not shown here), and the effects of CO 2 and solar activity are similar to those of March, and the geomagnetic field effects differ greatly from those of March, which lead to a weakening of both the neutral 170 density and temperature except for an increase near 5500 BC. The contributions of geomagnetic field structure, dipole moment, and Joule heating are still evident in the June results.

The long-term trends generated by the variations of the geomagnetic field, CO2, and solar activity
From Figures 3-5, it can be found that the effects of the three factors vary approximately linearly. Therefore, we calculated linear growth rates for the effects of the three factors on neutral density and temperature, as shown in the text of each panel in Figure 6, which reveals the long-term trends of the thermosphere generated by the variations of the geomagnetic field (dipole moment and the colatitude of the north magnetic pole), CO 2 , and solar activity, respectively. The grey dots in Figure 6  of the global mean neutral density at 400 km over the entire simulation time interval is about 2.26×10 -12 kg m -3 . Based on this value and the linear growth rate shown in Figure 6, the effects due to changes in the geomagnetic dipole moment, the colatitude of north magnetic pole, CO 2 , and solar activity during the entire Holocene can be calculated to be about -1.4%, -1.6%, -31% and +250%, respectively. While the effects are about -7 K, -4 K, -48 K, and +557 K for the neutral temperature, respectively. 195 Solomon et al. (2019) pointed out that the temperature change from solar minimum to maximum increases by about 500 K at 400 km based on the simulation of the Whole Atmosphere Community Climate Model-eXtended (WACCM-X), which is generally consistent with our results. Since the effect of magnetic pole position is not as large as that of dipole moment, only the dipole moment is used later to quantify the effect of the magnetic field. March and June essentially increases with altitude, but it is close to constant above 300 km. The effect of the geomagnetic field on the neutral temperature is slightly greater in June than in March. For every 10 22 Am 2 increase in dipole moment, the neutral temperature in June decreases by ~2.7 K, compared to ~1.7 K in March. In contrast, CO 2 and solar activity have the 205 opposite effect on neutral temperature, with a greater effect in March. For every 10 ppm increase in CO 2 , the global mean neutral temperature above 200 km decreases by ~3.5 K in March and ~2.3 K in June. Akmaev and Fomichev (1998) suggested a trend of about -3.1 K per 10 ppm at 200 km in the thermosphere in April due to increasing CO 2 , while Cnossen (2014) and Solomon et al. (2018) reported trends of about -1 K and -1.8 K per 10 ppm above 200 km, respectively. Therefore, it is reasonable that our results are 0.8-3.5 K per 10 ppm between 150 and 600 km. For each 1 sfu increase in F10.7, the neutral 210 temperature increase is ~6.2 K and ~4.4 K in March and June, respectively. The effect of the three factors on neutral density is similar to that on neutral temperature. Since neutral density decreases exponentially with altitude, the effect of the three factors on neutral density (absolute value change) also decreases with an altitude above 150 km, as shown in the left column of Figure 7. To present more clearly the effect of the three factors on neutral density, the solid lines display the corresponding percentage changes using 2005 simulations of CR4 as a reference. This reveals that the effects of solar activity and dipole 215 moment on neutral density increase significantly with increasing altitude in March and June with a stronger effect in June, while the effect of CO 2 is basically unchanged with altitude in both March and June with a stronger effect in March. In addition, the increasing dipole moment and CO 2 decrease the neutral density, while the rising solar activity increases the neutral density.
As mentioned in Section 3.2, only the effect of the geomagnetic field on the thermosphere displays the latitude and UT variations. Figure 8 shows the linear growth rate of neutral density (top row) and temperature (bottom row) at 400 km attributed 220 to the geomagnetic field versus latitude and UT (left column) or longitude (right column). In general, an increase in the dipole moment attenuates the thermospheric neutral density and temperature at all latitudes, longitudes, and UTs. The linear growth rate is greater at high latitudes in the 0-8 and 18-24 UT and can reach -3.6×10 -14 kg m -3 /10 22 Am 2 or -9.5 K/10 22 Am 2 , while it is about -1×10 -14 kg m -3 /10 22 Am 2 or -2 K/10 22 Am 2 at other latitudes and UTs. In addition, a larger linear growth rate of up to -4.4×10 -14 kg m -3 /10 22 Am 2 or -8.7 K/10 22 Am 2 is seen in all longitudes above ±60° latitude, and it is also about -1×10 -225 14 kg m -3 /10 22 Am 2 or -2 K/10 22 Am 2 in other regions. Overall, the thermosphere in the polar regions of the Southern Hemisphere is more influenced by the geomagnetic field than that in the Northern Hemisphere, while the influence of the magnetic field is weaker and of about a similar magnitude in the middle and low latitudes.

Future projections
As the number of human space missions increases explosively, more and more spacecraft will operate in the thermosphere, so 230 projecting the future state of the thermosphere is also important to human life. Based on the IPCC projections of greenhouse gas emissions under different scenarios (IPCC, 2014), we can simply and reasonably assume that CO 2 concentrations will rise by 400 ppm over the next century. Therefore, according to the calculations shown in Table 2, the global mean neutral density will decrease by ~70% and ~50% in March and June, respectively, due to a 400 ppm increase in CO 2 . This is generally consistent with the trend of about -6.1±0.8% per decade throughout the 21st century projected by Cnossen (2022). Also, from 235 Figure 7, it can be concluded that the neutral temperatures in March and June will decrease by ~114 K and ~84 K, respectively. This is larger than the projections (~60 K) of Cnossen (2022).
The dipole moment decreases by about 3.5% over the next 50 years based on the prediction by Aubert (2015), causing an increase in global mean neutral density of up to 1% above 500 km according to the simulations of Cnossen and Maute (2020).
However, the increase of the global mean neutral density is projected to be ~0.08% and ~0.25% in March and June, respectively, 240 based on our calculated linear growth rate. In addition, we can also project that the temperature increase due to the decrease of dipole moment in March and June is about 0.5 K and 0.7 K, respectively. As shown in Figure 8, the effect of the geomagnetic field is strongly dependent on UT and geographic location, so the global mean state projection is provided for reference only.
Furthermore, the effect of geomagnetic field variation is negligible compared to the effect of rising CO 2 over the next 100 years. 245

Conclusions
In this study, the evolution of the thermosphere during the Holocene (from 9455 BC to 2015 AD) was simulated using the independently developed global ionosphere-thermosphere theoretical model GCITEM-IGGCAS, driven by the realistic geomagnetic field model, CO 2 level, and solar activity derived from modern measurements and ancient natural media.
Furthermore, through a series of control simulations, we quantify the thermospheric temperature and density changes due to 250 variations in the geomagnetic field, CO 2 levels, and solar activity. The main conclusions are presented below.
1. The climatological morphology of the global thermosphere during the Holocene is reconstructed for the first time.
Thermospheric neutral density is mainly controlled by solar activity and modulated by CO 2 and geomagnetic field.
Typically, the geomagnetic field configuration directly affects the morphology of the equatorial mass anomaly structure of the thermosphere, while CO 2 mainly affects the magnitude of the neutral density and temperature. In general, the 255 frequent oscillations of solar activity contribute more than 80% of the thermospheric variability, while the contributions of CO 2 and geomagnetic field are both less than 10%. The effect of CO 2 is comparable to that of the geomagnetic field throughout the Holocene for the global mean neutral density and temperature but becomes more significant after 1800 AD. Only the effect of the geomagnetic field is strongly dependent on the universal time and geographical location, and the weakening of the dipole moment leading to an increase in Joule heating in the polar region thus make the thermosphere 260 change more than the effect of CO 2 . Overall, the higher altitude, the larger the effect of the three drivers on the neutral density and temperature.
2. Both the thermospheric neutral density and temperature vary approximately linearly with the dipole moment of the geomagnetic field, CO 2 , and the F10.7 index of solar activity. The global mean variability of the neutral density at 400 km during the March equinox due to changes in the geomagnetic dipole moment, CO 2 , and solar activity during the entire 265 Holocene can be about -1.4%, -31%, and +250%, respectively. While the effects are about -7 K, -48 K, and +557 K for the neutral temperature, respectively. In addition, there is a clear altitude and seasonal variation in the thermosphere change due to an increase in per unit of dipole moment, CO 2 , and solar activity. Different factors produce different seasonal variations in thermosphere changes. The increasing dipole moment and CO 2 decrease the neutral density and temperature, while the rising solar activity increases them. 270 3. We project that a 400 ppm increase in CO 2 will result in a 50-70% reduction in global mean thermospheric neutral density depending on the season, while neutral temperatures will decrease by 84-114 K. This is enough to change the orbit and lifetime of spacecraft and space debris, which deserves the attention of future space missions. The effect of decreasing dipole moments of the geomagnetic field over the next 100 years on the global mean thermospheric neutral density and temperature is negligible, but the effect of changing magnetic field configurations (e.g., magnetic pole positions) on the 275 thermosphere should be considered, especially in the polar regions.

Competing interests
The contact author has declared that neither of the authors has any competing interests. 290