Future Decreases in Thermospheric Neutral Density in Low Earth Orbit due to Carbon Dioxide Emissions

Increasing carbon dioxide causes cooling in the upper atmosphere and a secular decrease in atmospheric density over time. With the use of the Whole Atmospheric Community Climate Model with thermosphere and ionosphere extension (WACCM‐X), neutral thermospheric densities up to 500 km have been modeled under increasing carbon dioxide concentrations. Only carbon dioxide and carbon monoxide concentrations are changed between simulations, and solar activity is held low at F10.7 = 70 throughout. Neutral density decreases through to the year 2100 have been modeled using four carbon dioxide emission scenarios produced by the Intergovernmental Panel on Climate Change (IPCC). The years 1975 and 2005 have also been simulated, which indicated a historic trend of −5.8% change in neutral density per decade. Decreases in the neutral density relative to the year 2000 have been given for increasing ground‐level carbon dioxide concentrations. WACCM‐X shows there has already been a 17% decrease in neutral densities at 400 km relative to the density in the year 2000. This becomes a 30% reduction at the 50:50 probability threshold of limiting warming to 1.5°C, as set out in the Paris Agreement. A simple orbital propagator has been used to show the impact the decrease in density has on the orbital lifetime of objects traveling through the thermosphere. If the 1.5°C target is met, objects in Low Earth Orbit (LEO) will have orbital lifetimes around 30% longer than comparable objects from the year 2000.

• Thermospheric neutral density at 500 km altitude lowers by over 80% with a high ground-level carbon dioxide concentration of 890 ppm • Meeting the 1.5°C Paris Agreement target limits the reduction in neutral density at 400 km since the year 2000 to around 28% • Objects in Low Earth Orbit (LEO) will have orbital lifetimes around 30% longer at the 1.5°C target than comparable objects from the year 2000 with increased extreme ultraviolet emission which heats the thermosphere, causing oscillatory variation in neutral densities of an order of magnitude. These peaks also vary in magnitude from cycle to cycle, with the most recent (solar cycle 24) maxima having one of the smallest peak sunspot numbers observed since records began in 1749. There is a wide range of literature trying to predict the peak of the next solar cycle, with these predictions ranging from lower than the last, to being one of the largest recorded (McIntosh et al., 2020). The length of each cycle also varies, with an average period of around 11 years and the majority falling between 10 and 12 years. Part of this increased period can be attributed to the length of time taken for the Sun to start recovering from low solar activity to high, with the last 2 solar minima being particularly prolonged (Lockwood et al., 2012). Both of these minima were also signifcantly lower than previous historic minima, resulting in historically low neutral densities.
The secular decrease in thermospheric density caused by an increase in CO 2 concentration is of particular concern for all objects orbiting in the Low Earth Orbit (LEO) region. With decreasing atmospheric density, an orbiting object experiences less atmospheric drag. This becomes particularly relevant at altitudes below 500 km where drag is a dominant perturbing force and the effect of a secular density trend is substantial. The smaller drag force leads to a reduced rate of semi-major axis decrease, and hence a longer orbital lifetime. Observations of changes in the semi-major axis within the tracking data of historical LEO objects have been used to measure density trends (Emmert, 2015;Emmert et al., 2008;Keating et al., 2000;Saunders et al., 2011). Numerical atmospheric models have also been used to simulate historical time periods and estimate the magnitude of density trends associated with an increase in CO 2 concentration (Qian et al., 2006;Solomon et al., 2015Solomon et al., , 2018. Results from both groups of studies, at an altitude of 400 km, found under low solar activities are summarized in Table 1. All of these studies have found the historical secular trend to be negative, with values ranging from −2.5% to −7.2% per decade. While historical observations relied on tracked satellites, more recently satellites such as CHAMP, GRACE, and GOCE have been launched with accelerometers, allowing for higher resolution thermospheric density data (Doornbos, 2011;. Further improvements can be made to the data by bettering the drag models of these satellites, whether that be gas-surface interactions or the geometry and aerodynamic models, as detailed in . In the long term, this will lead to better estimations of the density trend over multiple decades. The magnitude of the trend is inversely proportional to the level of solar activity (Emmert, 2015;Solomon et al., 2019), so the largest secular changes in neutral density are seen during low solar activity.
Space debris in LEO also have their orbital lifetimes increased due to the declining thermospheric density. Made up of all the discarded components and collision fragments left over from human activity in orbit, space debris poses a substantial risk to operational spacecraft (ESA Space Debris Office, 2020). The number of trackable objects (greater than 10 cm in size) intersecting the LEO region continues to increase, reaching nearly 17,500 in 2020, of which around 2300 are active spacecraft. As the amount of debris increases, BROWN ET AL.  (2015) changes depending on the period it is calculated over. b k q , CO 2 -O collisional deactivation rate, of ∼1.5 × 10 −12 cm 3 s −1 or 3.0 × 10 −12 cm 3 s −1 .

Table 1
Observations and Models of the Historic Density Trend at 400 km and for Low Solar Activities Only so does the probability of a collision which would create further debris. This poses further risk to active spacecraft and disruption to operations as spacecraft operators have to respond to possible conjunctions.
Space debris models have been created to investigate the evolution of the debris environment. However, the decreasing thermospheric density trend has not been included in the majority of these models. Simulations can run many decades, even centuries into the future, over which secular density trends will have an important cumulative effect. The models are also used to assess possible ways to reduce the amount of debris and the risk it poses. One method for operators to reduce the risk posed by space debris in the future is to follow the voluntary debris mitigation guidelines introduced by the United Nation's Committee on the Peaceful Use of Outer Space (UN COPUOS), which includes the limiting of orbital lifetime to 25 years once a satellite's mission is complete. Another method known as Active Debris Removal (ADR) is the launching of satellites with the aim of removing the pieces of space debris which pose the greatest risk to the environment.
The Intergovernmental Panel on Climate Change (IPCC)'s Fifth Assessment Report (AR5) published four Representative Concentration Pathways (RCPs), with each giving a CO 2 trajectory through to the year 2100. These are shown in Figure 1 and summarized in the Synthesis Report from the IPCC (2014). These are entitled RCP2.6, RCP4.5, RCP6.0, and RCP8.5, where the number refers to the radiative forcing in W/m 2 in the year 2100 for each scenario. While these are not meant to be predictions of the future, they provide a limited number of baseline scenarios and CO 2 concentration trajectories from which modeling can be performed and results across studies compared. The RCPs are used to simulate future density reductions and trends within this work. For more information on how four independent groups arrived at the CO 2 concentration profiles used within the RCPs, see van Vuuren et al. (2011).

Model Simulations
The Community Earth System Model (CESM) from the National Center for Atmospheric Research (NCAR) allows for simulation of the whole, coupled Earth climate, using separate modules for each major system (Hurrell et al., 2013). In this work, the Whole Atmosphere Community Climate Model with thermosphere and ionosphere extension (WACCM-X) module is used (Liu et al., 2010). This models the atmosphere numerically from ground level through to 4 × 10 −10 hPa, which corresponds to an altitude of 350 km at CO 2 levels from the year 2000 and low solar activity. The model has a resolution of 1.9° in latitude and 2.5° in longitude (using a 96 by 144 grid), with 81 vertical pressure levels and a resolution of one-quarter scale height above 1 hPa. Solar and geomagnetic activity are parametrized in the model using established indices of activity such as the solar 10.7 cm radio flux, F 10.7 , and the planetary Kp index respectively. For a full description of the model, including the chemistry and radiative transfer within WACCM-X, see Liu et al. (2010) and Liu et al. (2018). For this work, CESM 1.2.2 was run on the University of Southampton computing cluster, Iridis 4.
There are only small variations in carbon dioxide concentration from ground level up to around 70 km where the atmosphere is well mixed. Above this altitude molecular diffusion takes over and the concentration reduces towards 0 ppm, as seen in Figure 2c. New initial conditions for WACCM-X under varying CO 2 concentrations were created by scaling the CO 2 values at each pressure level in the original initial history files (here chosen as the WACCM-X default files for the year 2000) by the relative increase in ground-level CO 2 concentration. At the higher altitudes, CO 2 and carbon monoxide (CO) exist in chemical equilibrium. To account for this, CO values are scaled similarly to CO 2 . This accounts for over 99.7% of the carbon in the thermosphere and minor constituents containing carbon such as methane (CH 4 ) were not scaled. However, BROWN ET AL.  the effect of changes in methane on thermospheric density is expected to be much smaller than the effect of the increase in CO 2 concentration (Roble & Dickinson, 1989).
Nitric oxide (NO) also plays an important role in thermospheric cooling during solar maximum (Mlynczak et al., 2016), and could be increased by the anthropogenic emission of the greenhouse gas nitrous oxide (N 2 O). However, the large amount of Nitrogen (N 2 ) in the lower atmosphere acts as a reservoir, keeping NO in the thermosphere at a relatively stable level. To remove the effects of solar activity (and therefore also NO) on cooling and the secular density trend, the F 10.7 and Kp indices were held fixed at 70 sfu and 0.33 respectively in all simulations. While Kp values can be used as integers between 0 and 9, WACCM-X uses a Kp decimal value which provides a finer scale for geomagnetic activity. A Kp value of 0.33 is similar to the value of 0.3 used in Solomon et al. (2015), allowing for easier comparison between the historical trends. Other increases in gas concentrations such as O 3 , CH 4 , and H 2 O vapor and CFCs which would impact the stratosphere and mesosphere regions were not considered.
The Earth's magnetic field also changes over time, affecting the ionosphere and in turn the thermosphere (Cnossen, 2014). It was found that for the period 1908 to 2008, the historic changes in the magnetic field do contribute to cooling at 300 km, however, the increasing CO 2 concentration dominates the thermospheric cooling. A recent study by Cnossen and Maute (2020) confirmed the predicted changes in the Earth's magnetic field up to the year 2065 will result in at most a 1% to 2% increase in neutral density from 2015 to 2065. This is equivalent to an increase of 0.2% to 0.4% per decade over the next 5 decades, which is an order of magnitude smaller than the decreasing historical trends of −2.5% to −7.2% per decade summarized in Table 1. The impact of the changing magnetic field will therefore not be considered in the simulations for BROWN ET AL.  this study. The magnetic field within WACCM-X, namely the International Geomagnetic Reference Field (IGRF-12) described in Thébault et al. (2015), is held fixed at the year 2000 level.
The CO 2 -O collisional deactivation rate (quenching rate), k q , has a major impact on the cooling within atmospheric models (Akmaev, 2003). This can be seen by the change in density trend of Solomon et al. (2015) from −4.9% per decade at k q = ∼1.5 × 10 −12 cm 3 s −1 to −6.8% at k q = 3.0 × 10 −12 cm 3 s −1 . Recent measurements of the value vary between 1.5 × 10 −12 and 8.0 × 10 −12 cm 3 s −1 and are summarized in Feofilov et al. (2012). The value of k q used for this study has been left as the default in WACCM-X, namely 3.0 × 10 −12 cm 3 s −1 as this was the value WACCM-X was calibrated with during development.
WACCM-X performed simulations of 10 different ground-level CO 2 concentrations, each for 16 model months total. These are the CO 2 concentrations of RCP8.5 at 10 years intervals starting from 2005 and ending on 2095, along with the year 1975 (for validation against previous studies). All simulations were performed as free runs, and hence the lower atmosphere (< ∼50 km) was not constrained by data. Lower boundary conditions were set following the Coupled Model Intercomparison Project (CMIP5) implementation of RCP8.5 (Taylor et al., 2012). Ground-level CO 2 concentrations were allowed to vary with season and propagate through the atmosphere. Therefore, annual mean CO 2 concentrations are stated through this work. The seasonal sinusoidal variation had a maximum amplitude of 1% of the annual mean. Numerical models (like WACCM-X) require time for the user-set, initial conditions of the modeled atmosphere to stabilize into a physical equilibrium. While other studies such as Solomon et al. (2018) allowed one year for minor chemical constituents to equilibrate, we found temperature and density appeared to equilibrate after one month. As a precaution, the first 4 months of each data set were ignored, leaving 12 model months for each CO 2 concentration modeled. Key values from the simulations such as temperature and chemical concentrations were stored with a temporal resolution of 3 hours.
Along with the 10 simulations stated earlier, an additional 64 month long simulation of the year 2000 has also been completed. This provided a more robust reference point against which the later simulations could be compared. 64 months was chosen as 12 months are simulated 5 times, with the end of one simulated year being used as the initial conditions of the next 12 months cycle. Again 4 months are added to the start of this for equilibration. Figure 3 plots the differences for each of five cycles from the overall globally averaged temperature and globally averaged density altitude profiles for the 64 months simulation of the year 2000 (with the first 4 months ignored). This shows the impact the initial conditions have on the simulation of each year. These longer simulations could not be performed for the other concentrations due to limited computing time.

Processing
A geopotential height, h, is output at each of the grid points of WACCM-X. This was converted into a geometric altitude, z, via where r E is the average radius of the Earth, 6,371 km. The geometric altitude is equivalent to a satellite's orbiting altitude. All interpolation between the 81 altitude levels of WACCM-X for each latitude, longitude pair was performed via 1-D monotonic cubic interpolation.
WACCM-X models the thermosphere down to a pressure level of 4 × 10 −10 hPa. This minimal pressure level varies in height, with the maximum altitude decreasing as the thermosphere cools. At high ground-level CO 2 concentrations the maximum model altitude is only 280 km. As the motivation of the work was to understand the reduction in density at commonly used LEO altitudes, the WACCM-X model data was extrapolated to an altitude of 500 km. It was found that the function that best fit the densities at points above 175 km for all cases was calculated by where ρ is atmospheric density and a, b, c, and d are coefficients fit to the modeled density with non-linear least squares. This allowed a neutral density profile as a function of altitude to be obtained for each latitude-longitude combination. It should be noted that from 175 km to the maximum modeled altitude, WACCM-X is dominated by atomic oxygen. Helium becomes the dominant species within the thermosphere above around 700 km, which would not be accounted for by this extrapolation method. Applying the extrapolation method to the NRLMSISE model (Picone et al., 2002) which does account for helium showed an increasing deviation with altitude, reaching an overestimate of 3% at 500 km. As the relative differences in densities was of primary concern in this study, and the overestimate was consistently applied, this extrapolation method proved acceptable when used up to 500 km. Temperatures in the thermosphere tend towards the exospheric temperature with increasing altitude. It is assumed the temperature at the highest modeled altitude is the exospheric temperature, and therefore this temperature is used for all higher altitudes up to 500 km.
To obtain global mean annual mean values of temperature and density, data was transformed from the model's pressure levels to geometric altitudes using the above methods, and then averaged temporally to calculate annual means at each altitude, latitude, and longitude. The global mean annual mean is then obtained by averaging over latitude (with cosine latitude weighting) and longitude.
BROWN ET AL.

Figures 4a and 4b
show the WACCM-X global mean temperature plotted against ground-level CO 2 concentration at the fixed pressure level of 2.9 × 10 −7 mbar (around 200 km altitude) during the March equinox and June solstice. This was done to compare the WACCM-X data at large CO 2 concentrations against three implementations of the Coupled Middle Atmosphere Thermosphere version 2 (CMAT2) model (Cnossen et al., 2009), where each implementation has different gravity wave parameterization schemes. While WACCM-X is cooler than all three implementations of the CMAT2 model, WACCM-X's temperature decreases with increasing CO 2 as expected, and at a similar rate to that of CMAT2.
Using the methods described in this study, the trend of decreasing neutral density due to increasing carbon dioxide for the period 1975-2005 has been calculated to be −5.8% per decade at 400 km altitude. This places it in the middle of the range predicted by other models and observations which were summarized in Table 1.
Stating trends for the reduction in neutral thermospheric density is useful for studies of historical periods, as the change in carbon dioxide concentration is fixed. However, for future reductions in density, CO 2 concentration will be variable and any future density trend will be dependent upon it. Figure Figure 6 shows the density change relative to the year 2000 for each RCP scenario. These data were generated by combining the data shown in Figures 1 and 5. Under the high emissions RCP8.5 scenario, the neutral density relative to the year 2000 reduced by 75% at 400 km by the year 2095. Under the RCP2.6 scenario, neutral densities at 400 km reach their largest change between the years 2040 and 2060 as CO 2 concentrations peak and begin to decline. This peak sees global ground-level temperature rise limited to 2°C or below and a 25% reduction in neutral density compared to the year 2000. Neutral densities recover to a 20% reduction by 2100 as the CO 2 concentration declines (see Figure 1). Regardless of scenario, the WACCM-X model shows that neutral densities have dropped by 17% at 400 km in 2020 relative to the year 2000.

Discussion
The targets and predictions set out within the Paris Agreement can provide some further context to this work's results. The Emissions Gap Report from the United Nations Environment Programme (2019) states that for a 50% probability of limiting global warming to 1.5°C, the carbon budget from 2018 onward is 580 GtCO 2 . Taking the 2017 globally averaged CO 2 concentration of 405.0 ppm (Le Quéré et al., 2018), this sets a target of 480 ppm above which there is higher than 1.5°C warming at ground level. This target is independent of time and is reached in different years in the RCP scenarios. By looking at the density reductions in Figure 5, it can be seen that the 1.5°C warming target is equivalent to a 32% drop in neutral density at 400 km relative to the year 2000. For a 66% probability of limiting warming to 1.5°C, there is a remaining budget from 2018 of 420 GtCO 2 . This is equivalent to a concentration of 459 ppm and a drop of 27% in neutral density at 400 km. The Emissions Gap Report also gives a prediction that under the current uncon- The empirical atmospheric model NRLMSISE-00 (Picone et al., 2002) is one of the thermospheric models used in space debris models and is calibrated with observational data up to the year 1997. The relative density drop seen in Figure 5 can therefore be used as a scaling factor for the neutral density output from NRLMSISE-00 to account for the density trends from increases in CO 2 . This has been done for a few objects in LEO in Table 2. A simple numerical propagator was used, with drag as the only perturbing force and NRLMSISE-00 as the density model. The changes presented here are to provide initial implications of the work for the debris environment. The lifetimes at 468 ppm was chosen as the closest modeled point using WACCM-X to the 1.5°C target CO 2 of 480 ppm. The orbital lifetimes of the chosen objects increase by between 26.9% and 32.6% relative to the year 2000. Under the highest modeled CO 2 concentration of 890 ppm, the orbital lifetime approximately triples when compared to the year 2000.
If these increases in orbital lifetime are experienced by all objects, they will have a substantial impact upon the space debris environment in at least two ways. First the number of orbiting objects will increase. The number of objects can only increase by adding new objects (most commonly via launches), or having orbiting objects collide and create collision fragments. Longer lifetimes lead to a greater cumulative chance of collison, resulting in further fragmentation events. Second, operators who meet the 25-year guideline by using atmospheric drag to passively remove their satellite from orbit will have to take the reduction in density into account by lowering the perigee of their final orbit or by using re-entry technologies such as drag augmentation devices. In both cases, further satellite mass has to be dedicated to the change.
These results have been computed under low solar activity only. Therefore, it is only appropriate to apply these density reductions at solar minima. There will still be a secular density reduction when looking at solar maxima, albeit smaller (see Emmert (2015)). It is unknown to what scale this will be and this will be the subject of future work. The density reduction is expected to scale with the solar activity between the minima and maxima, albeit the exact relationship is unknown. It is acknowledged that a number of assumptions have had to be made to allow for the modeling of future densities. These include the value of k q , extrapolation to higher altitudes and ignoring the density change due to the Earth's changing magnetic field. The initial state of each object was given by the two line element set on November 2, 2020 with use of the North American Aerospace Defense (NORAD) satellite catalog number. Coefficient of drag for all objects was set as Cd = 2.2. a Reference lifetime given is for the year 2000, equivalent to 379 ppm.

Conclusion
Increasing ground-level CO 2 concentrations result in decreasing thermospheric densities. These have been historically observed and modeled as summarized in Table 1. The impact on neutral densities under further increasing CO 2 concentrations have been modeled and presented within Figure 5. These neutral density reductions were also estimated for the four IPCC RCP scenarios in Figure 6. Density reductions reported throughout the paper have been given relative to the density in the year 2000 (under low solar activity) to provide a scaling ratio which can be applied to fast, empirical atmospheric models such as NRLMSISE-00. It has also been shown that limiting warming to the 1.5°C target of the Paris Agreement will still result in a 27% to 30% reduction in density relative to the year 2000 at 400 km, with orbital lifetimes increasing by around 30% as a result. Even under the low CO 2 emissions of the RCP2.6 scenarios, the drop in neutral density will have a substantial impact on LEO object orbital lifetimes, as well as the debris environment within this region. Any CO 2 concentrations higher than this will have an even more substantial impact on thermospheric densities and the space industry as a whole, particularly on constellations relying on atmospheric drag to passively remove their satellites from orbit.