Analysis and Attribution of Climate Change in the Upper Atmosphere From 1950 to 2015 Simulated by WACCM‐X

Monitoring climatic changes in the thermosphere and ionosphere and understanding their causes is important for practical purposes. To support this effort and facilitate comparisons between observations and model results, a long transient simulation with the Whole Atmosphere Community Climate Model eXtension (WACCM‐X) from 1950 to 2015 was conducted. This simulation used realistic variations in solar and geomagnetic activity, main magnetic field changes, and trace gas emissions, including CO2, thereby including all known drivers of upper atmosphere climate change. Analysis of the full 1950–2015 interval with a standard multilinear regression approach demonstrated difficulties in removing solar cycle effects sufficiently to obtain reliable trends. Results improved when an (F10.7a)2 was included in the regression model, in addition to terms for F10.7a, KP, and the trend itself. Comparisons with previous studies and analysis of spatial variations in trend estimates confirmed that the increase in CO2 concentration is the main driver of trends in thermosphere temperature and density, but at high (magnetic) latitudes effects of main magnetic field changes play a role as well, especially in the Northern Hemisphere. Spatial patterns of trends in hmF2, NmF2, and total electron content indicate a superposition of CO2 and geomagnetic field effects, with the latter dominating trends in the region of ∼50–20°N, ∼60°W to 20°E. Additional model experiments to investigate the indirect dynamical effects of climate change in the lower atmosphere (<50 km) on the upper atmosphere (>100 km) suggested that these effects are small and insignificant. However, current model limitations could mean that these effects are underestimated.


Introduction
Climatic changes have been observed to take place throughout the atmosphere system. While the troposphere shows a global warming trend, the middle and upper atmosphere (stratosphere, mesosphere, and thermosphere) have been cooling, with more complicated trend patterns occurring in the ionosphere (e.g., Cnossen, 2012;Laštovička et al., 2012;Qian et al., 2011). It is important to monitor long-term changes in upper atmosphere climate and understand their causes, so that predictions of the future state of this environment can be made. These long-term predictions are necessary for various practical applications, such as planning for new satellite missions, managing the risks of space debris (e.g., Lewis et al., 2011), the interpretation of long-term climate monitoring with satellite-based data (e.g., Scharroo & Smith, 2010), and assessing future space weather risks to assets and infrastructure in near-Earth space.
Known drivers of middle and upper atmosphere climate change include the increase in greenhouse gas concentrations, changes in ozone concentration, long-term solar and geomagnetic activity variations, and 10.1029/2020JA028623 the secular variation of the Earth's main magnetic field. Upper atmosphere trend analyses usually attempt to remove, or at least suppress, any signal of solar and geomagnetic activity variations. If we leave these influences aside, the dominant cause of long-term trends in the thermosphere is thought to be the increase in atmospheric CO 2 concentration (Akmaev et al., 2006;Cnossen, 2014;Qian et al., 2013), while changes in the Earth's magnetic field play a relatively more important role in long-term trends in the ionosphere (Cai et al., 2019;Cnossen, 2014).
The attribution of observed upper atmosphere trends is strongly based on controlled model experiments, where only the driver of interest (e.g., the CO 2 concentration) is varied, while all other factors are kept constant, including solar and geomagnetic activity levels. While this approach provides a good indication of the effect of the driver of interest, direct comparisons with observational data sets are complicated, due to differences in the analysis techniques used, differences in the time windows being studied, and the observational trends likely being a mixture of trends caused by different mechanisms. To resolve these issues, we present here for the first time a long transient simulation with the Whole Atmosphere Community Climate Model eXtension (WACCM-X) for the period 1950-2015, which includes all known drivers of long-term changes in the upper atmosphere. The data set is published separately (Cnossen, 2020), providing full access to the wider community to perform model-observation comparisons or to use for other purposes.
The data set is used here to study trends in selected global mean upper atmosphere quantities, as well as spatial variations in trends. Spatial patterns are difficult to study with observations due to uneven and limited global coverage but have also received rather little attention in modeling studies. We demonstrate that the spatial variations in the trends we obtain provide important insights in the underlying driving mechanisms, while also providing context for local trend observations. Further, two additional model experiments were conducted to test the role of the lower atmosphere (<50 km) in causing long-term trends in the upper atmosphere, as might take place through changes in wave forcing (Oliver et al., 2013). This paper is organized as follows. Section 2 describes the methodology, including a brief introduction to the model (described in detail elsewhere) and a description of the simulation setup and the data analysis procedures. Results are presented in sections 3-5, with section 3 focusing on global mean trends, section 4 on spatial variations in trends, and section 5 on the indirect effects of lower atmosphere climate change on the upper atmosphere. Each of these sections also contains some additional background material on the topic at hand. Section 6 finishes the paper with a discussion of the results together with a summary of the main conclusions.

Methodology
This study uses WACCM-X 2.0, part of the Community Earth System Model (CESM) (Hurrell et al., 2013) release 2.1.0, developed by the National Center for Atmospheric Research (NCAR). WACCM-X extends the standard WACCM, which has a model top at ∼140 km, to an altitude of ∼500 km. WACCM-X was run with the standard horizontal resolution of 1.9 • latitude × 2.5 • longitude. In free-running mode the model has 126 vertical levels, which is extended to 145 levels in "specified dynamics" (SD) mode, in which the temperature, zonal wind, and meridional wind fields up to ∼50-km altitude are relaxed toward reanalysis data.
WACCM-X is based on NCAR's Community Atmospheric Model (CAM) and includes all of the physical parameterizations of that model. In addition, WACCM-X incorporates a detailed neutral chemistry model for the middle atmosphere, a model of ion chemistry in the mesosphere-lower thermosphere (MLT) region, ion drag and auroral processes, parameterizations of extreme ultraviolet (EUV) heating and infrared transfer under nonlocal thermodynamic equilibrium, and modified parameterizations of gravity wave effects and vertical diffusion. Recently, a new electrodynamics module was implemented. This describes the interactions between electric fields, plasma motions and neutral winds self-consistently, making WACCM-X ideally suited for this study. A more detailed description of WACCM-X 2.0 is given by Liu et al. (2018), and references therein.
A transient historical simulation with WACCM-X in free-running mode for the period 1950-2015 was conducted. The simulation was initialized from a previous WACCM-X simulation and run for a year (January to December 1949) to allow the climate to reach a quasi-steady state. Lower boundary forcings and chemical emissions were specified according to the Climate Model Intercomparison Project 6 (CMIP6) historical simulation (Eyring et al., 2016). Sea surface temperatures were prescribed from a coupled atmosphere-ocean CNOSSEN 2 of 14 10.1029/2020JA028623 simulation with CESM2-WACCM6 (Gettelman et al., 2019). Solar radiative and particle forcings were prescribed following the reference scenario of the CMIP6 recommendation by Matthes et al. (2017). The main magnetic field was specified according to the International Geomagnetic Reference Field (IGRF; Thébault et al., 2015). The simulation thereby includes all known drivers of long-term change in the upper atmosphere: the increase in CO 2 concentration and other greenhouse gases, changes in ozone-depleting substances, changes in the Earth's magnetic field, and variations in solar and geomagnetic activity levels. Further, any influences of climate change in the lower atmosphere, which some have suggested could be important (e.g., Oliver et al., 2013), are also accounted for, as WACCM-X includes the lower atmosphere.
A single transient simulation was done rather than an ensemble, since the upper atmosphere, unlike the lower atmosphere, is a strongly forced system, where the climate is largely determined by external drivers (e.g., Codrescu et al., 2018;Siscoe & Solomon, 2006). Small changes in the initial state therefore do not significantly affect the calculated climate of the thermosphere and ionosphere. On the other hand, one could argue that the simulated climate change in the lower atmosphere may not be representative, as in the lower atmosphere, small differences in the initial state could result in larger variations in the climate. This could then potentially impact on our ability to simulate any effects of climatic changes in the lower atmosphere on the upper atmosphere climate correctly. In addition, it is difficult to isolate any indirect effects of climate change in the lower atmosphere on the upper atmosphere from just this single transient simulation. For these reasons, two additional model experiments were done specifically to examine the role of the lower atmosphere in causing upper atmosphere climate change.
For these two experiments the SD capability of WACCM-X was used. The simulations were set up similarly to the transient simulation described above, except that the temperature and winds in the troposphere and stratosphere were nudged with the Modern Era Retrospective Analysis for Research and Applications 2 (MERRA-2) reanalysis data (Gelaro et al., 2017). Both simulations were run for the period 1983-1987, after a spin-up period of 1 year. However, simulation N0 was nudged with MERRA-2 data for 1983-1987, while for simulation N1 the MERRA-2 data for 2015-2019 were used instead. The simulated climate in the atmosphere up to 50 km should then be relatively close to the actual climate for 1983-1987 for N0, while the simulated climate up to 50 km for N1 should be close to the climate of 2015-2019. The rest of the simulation setup was identical, so that the mean difference between the two simulations in the upper atmosphere provides an estimate of the indirect effects of climate change in the lower/middle atmosphere (<50-km altitude) between 1983-1987 and 2015-2019 on the upper atmosphere. The two intervals were chosen to have maximum separation between them within the period for which MERRA-2 data are available (1980 to present), so that the largest possible effect might be seen, while avoiding the early 1980s, when solar activity was very high. The 5 years in each interval serve as a small ensemble, similar to the approach of Solomon et al. (2018Solomon et al. ( , 2019 in their studies of long-term trends in the upper atmosphere. For all our analysis, monthly mean outputs were used. Global means were calculated using a cosine (latitude) weighting. While the model uses a vertical pressure coordinate, observational analyses are usually done with reference to geometric height, and therefore model outputs analyzed here were interpolated to geometric height to facilitate comparisons.
Removal of solar and geomagnetic activity variations from long-term data sets prior to or as part of trend analysis is notoriously difficult and one of the key problems in upper atmosphere climate research (Laštovička & Jelínek, 2019). Complications arise due to the relatively large amplitude of the ∼11-year solar cycle (e.g., Clilverd et al., 2003), as well as long-term trends in solar and geomagnetic activity levels themselves (e.g., Cnossen & Franzke, 2014). The most commonly used method to remove/suppress solar cycle variations is through a standard (multi)linear regression (e.g., Laštovička et al., 2006;She et al., 2009;Ulich & Turunen, 1997). We therefore took this same approach here, exploring several different multilinear regression models to extract long-term trends from our transient model simulation . The simplest regression model is as follows: where Y is the variable of interest, for example, the neutral temperature at 400-km altitude, and F10.7a is the 81-day average of the F10.7 index. F10.7a performed generally slightly better than the F10.7 index itself (resulting in a greater adjusted R 2 value) and therefore this was used here. More complex models, with additional predictor variables, were also tested, such as adding a geomagnetic activity term: and a further term for (F 10.7a) 2 : Where a geomagnetic activity term was included, the K P index was used, which performed better than the A P index. Inclusion of additional terms to describe seasonal variations did not alter trend estimates significantly and were therefore omitted. The F10.7, F10.7a, and K P index time series are shown in Figure 1. Note that the phase of the solar cycle is similar for the begin and end points of our time series, which should help suppress the influence of solar activity variations on the calculated long-term trends (e.g., Laštovička & Jelínek, 2019). The data record is also long (66 years, covering six solar cycles) compared to most observational data sets, which should put us in a very good position to obtain reliable trends, free of solar cycle influences. Table 1 gives an overview of the trends (± standard error) in various thermosphere-ionosphere parameters calculated with the three different multilinear regression models described above. The two most commonly used regression models for long-term trend analysis in the upper atmosphere are Model 1, which includes terms for only solar activity and a long-term trend (Equation 1), and Model 2, which includes an additional term for geomagnetic activity (Equation 2), while Model 3 includes an additional term for (F10.7a) 2 (Equation 3). We will compare the trends obtained here with those obtained by Solomon et al. (2018Solomon et al. ( , 2019, which are also summarized in Table 1. Solomon et al. (2018Solomon et al. ( , 2019 performed controlled experiments with WACCM-X to determine the magnitude of anthropogenic climate change throughout the atmosphere, comparing the intervals 1972-1976 and 2001-2005 with each other. Solar and geomagnetic activity levels were kept fixed in their simulations, which examined solar minimum and solar maximum background conditions separately. In our simulation the solar and geomagnetic activity levels varied mostly in between the extremes used by Solomon et al. (2018Solomon et al. ( , 2019. Therefore, we may expect our trends to be somewhere in between the trends they obtained for solar minimum and maximum conditions, if we assume that (1) influences of solar and geomagnetic activity variations have been properly controlled for by our linear regression analysis, and (2) the trends in the global mean upper atmosphere are primarily anthropogenic in origin, with little influence from changes in the main magnetic field (which was kept fixed by Solomon et al., 2018Solomon et al., , 2019.

Global Mean Trends
Regression Models 1 and 2 both give rather large trends for all variables and especially for the global mean neutral temperature and density at 400-km altitude. For instance, Solomon et al. (2018Solomon et al. ( , 2019 indicated trends in neutral temperature and density at 400-km altitude of, respectively, −2.8 K/decade and Note. Adjusted R 2 values are also given to indicate the quality of the fit to each regression model. The last two columns give trends obtained by Solomon et al. (2018Solomon et al. ( , 2019) (S2018, S2019) for controlled model experiments to assess the effect of the increase in CO 2 concentration under solar maximum (smax) and solar minimum (smin) background conditions, where available.
−3.9%/decade for solar minimum and −1.8 K/decade and −1.7%/decade for solar maximum conditions, while our trends obtained with regression Models 1 and 2 are about a factor 2 larger. The trends we find in the height of the peak electron density of the F 2 layer, h m F 2 , and the peak electron density of the F 2 layer, N m F 2 , appear somewhat more reasonable in comparison with Solomon et al. (2018Solomon et al. ( , 2019, but also these trends seem somewhat large. In general, there is little difference between regression Models 1 and 2.
Adding an additional term for (F10.7a) 2 (model/Equation 3) does not have a major effect on the trends in the neutral temperature at 400-km altitude or h m F 2 or the fit to these data. However, we find a marked decrease in trend magnitudes for the neutral density at 400-km altitude, N m F 2 and the total electron content (TEC): trends in these variables are about a factor 2 smaller than for regression Models 1 and 2. Further, the adjusted R 2 values for these variables increase, indicating an improved fit to the data. It therefore seems worthwhile to include an (F10.7a) 2 term. The neutral density trend obtained with regression Model 3 is also in excellent agreement with the Solomon et al. (2018Solomon et al. ( , 2019 results. The trend in N m F 2 obtained with regression Model 3, instead of somewhat large, now appears somewhat small, but considering the error margins, we can consider it consistent with Solomon et al. (2018Solomon et al. ( , 2019 as well. In contrast, the trend in neutral temperature at 400-km altitude and, to a lesser degree, the trend in h m F 2 are still larger than what we would expect. To check for any issues with appropriate removal of solar cycle influences, Figure 2 shows trends in the global mean neutral temperature and density at 400-km altitude, h m F 2 , N m F 2 , and TEC as a function of the end year of the time interval studied, with a minimum interval length of 40 years. If there are any remaining solar cycle influences that are not properly filtered out by the linear regression analysis, these would be expected to show up as oscillations in the trend magnitudes, also referred to as "ringing" (Clilverd et al., 2003). We do indeed see such oscillatory behavior that is indicative of remaining solar cycle influences. For regression Models 1 and 2 this is present for all variables shown, while the oscillations in the trends in the neutral density, N m F 2 , and TEC are much reduced for regression Model 3, which is consistent with the improved model performance for these variables. However, the trends in the neutral temperature at 400-km altitude and h m F 2 show quite large oscillations and notably weaker trends for end years beyond about 2008-2009. This could indicate trends are getting stronger over time (for instance, due to the CO 2 concentration increasing more rapidly over time), but might also be related to the low solar minimum of 2008/2009 and the subsequent weak solar cycle. If the solar cycle influence on the data is not fully accounted for by the linear regression analysis, this could have an unduly large influence on the trend magnitudes obtained for time windows ending in 2008 and after.
It is also worth noting that several studies have reported that the F10.7 proxy overestimates the true solar EUV forcing during the 2008-2009 solar minimum period, resulting in lower thermosphere and ionosphere densities than would be expected based on the observed F10.7 values (e.g., Elias et al., 2014;Solomon et al., 2013), while others have more generally reported long-term changes in the relationship between F10.7 and ionospheric parameters (e.g., Lastovicka, 2019). While this can be an additional cause for unexpectedly strong trends when the period from 2008-2009 is included in trend analyses done with real world data, CNOSSEN 5 of 14 the model's ionosphere-thermosphere system is driven by the observed F10.7 values rather than actually observed EUV emissions. Therefore, the model should not reflect any such effects.
To remove the remaining solar cycle influences as much as possible and avoid the potentially large influence of the weak solar activity since the 2008/2009 solar minimum, we may consider using the average trend values obtained for the period starting in 1950 and ending in 1997-2007, that is, the average of the trends obtained for these 11 periods, as a representative value. The standard deviation over these 11 values then gives us a measure of the uncertainty. An overview of the trends in the same parameters as listed in Table 1, but calculated with this new method, which we will refer to as Model 4, is given in Table 2. The new trends in thermosphere temperature at 400-km altitude and h m F 2 are in better agreement with Solomon et al. (2018Solomon et al. ( , 2019, although the trend in neutral temperature is still higher than expected. Further, the trend in N m F 2 is now rather weak compared to Solomon et al. (2018Solomon et al. ( , 2019. It is not immediately obvious what the cause of these discrepancies is.
However, now that we have spent considerable efforts minimizing the influence of solar cycle variations, we should examine whether other factors could play a role, such as changes in the Earth's main magnetic field. This is the other key difference between our study and Solomon et al. (2018Solomon et al. ( , 2019, who kept the Earth's magnetic field fixed. Since magnetic field changes and their effects on the upper atmosphere are location-dependent, it is appropriate to go beyond global mean quantities for this and examine the spatial variations in trends.

Spatial Variations in Trends
Observational evidence for spatial variations in long-term trends in the upper atmosphere is very limited. Observed trends in thermosphere density are based on orbital characteristics of many objects in near-Earth space (Emmert, 2015;Emmert et al., 2008;Saunders et al., 2011), and therefore represent some sort of global average trend. Thermosphere temperature trends are not directly observed, but inferred from either the decreasing trends in neutral density or from long-term ion temperature records from just a few locations (Donaldson et al., 2010;Zhang & Holt, 2013)-insufficient to establish any spatial patterns. Long-term ionosonde measurements are available from many more locations. While the distribution is far from globally uniform, several studies have examined spatial variations in trends in h m F 2 and the critical frequency of the F 2 layer, f o F 2 , based on these measurements ( note that o F 2 ∝ √ N m F 2 ). Bremer (1998) found a longitudinal variation in trends in h m F 2 and f o F 2 , with trends in both parameters being mostly negative (positive) in the European region west (east) of 30 • E. However, Upadhyay and Mahajan (1998) found no obvious spatial pattern after calculating trends from 31 stations, and a more recent analysis by Bremer et al. (2012) based on 37 stations agreed with this. Bremer et al. (2012) also specifically stated there is no significant correlation between trends in h m F 2 or f o F 2 and geomagnetic latitude, contrasting with Danilov and Mikhailov (1999) and Mikhailov (2002), who argued that trends in f o F 2 do depend on geomagnetic latitude, with stronger trends occurring at high latitudes. Mikhailov (2002) found that trends in h m F 2 do not show such latitudinal dependence.
Previous modeling studies suggest that at least the thermosphere response to increased CO 2 concentration is relatively uniform, though some latitudinal structure is present. Solomon et al. (2018Solomon et al. ( , 2019 showed slightly stronger annual mean trends at low latitudes to midlatitudes in the Northern Hemisphere. (Akmaev & Fomichev, 1998, 2000 obtained slightly stronger trends at low latitudes at equinox, but found a noticeably weaker response in the summer hemisphere at solstice. Qian et al. (2009) found that changes in h m F 2 and N m F 2 associated with an increased CO 2 concentration show considerably more spatial structure, organized to a degree by the Earth's magnetic field. For instance, belts of relatively strong noontime trends occurred in each hemisphere about 25 • away from the magnetic equator. Qian et al. (2009) argued that both changes in the O/N 2 ratio and changes in plasma transport processes contributed to these spatial variations, with the latter being most important at solar minimum. On the other hand, Cnossen (2014) found very little spatial variation in both thermosphere and ionosphere responses to increased CO 2 concentration around equinox.
Changes in the Earth's magnetic field produce changes in the ionosphere that are very location dependent, with the strongest effects occurring in the region of ∼50-50 • N and ∼100 • W to 50 • E (Cnossen & Richmond, 2008Cnossen, 2014). Effects of changes in the Earth's magnetic field on thermosphere temperature are strongest at high magnetic latitudes (Cnossen, 2014;Cnossen et al., 2016;Cnossen & Maute, 2020), although exact patterns differed between different studies. The main way in which the magnetic field can affect the thermosphere temperature is via changes in Joule heating. The magnetic field intensity affects the Joule heating directly through its influence on ionospheric conductivity, with larger intensity giving lower conductivity, as well as indirectly by affecting ion ⃗ E × ⃗ B drift and neutral wind patterns (Cnossen et al., 2011. In addition, the geographic positioning of the magnetic poles and the associated auroral zones affects the geographic distribution of the Joule heating, and to some degree also its magnitude . Figure 3 shows a map of the long-term trend in neutral temperature at 400 km, based on linear regression Model 3. Results for the other regression models, as well as Model 4 defined in the previous section, show very similar spatial variations, though the trend magnitudes may vary (especially for Model 4). Trends are clearly larger at high magnetic latitudes, especially in the Northern Hemisphere around ∼50 • N, ∼75 • W, while a belt of relatively weak trends roughly tracks the magnetic equator. This latitudinal variation does not match the expected pattern to result from the increase in CO 2 concentration: this would more likely give larger trends at low latitudes (Qian et al., 2009). The spatial variation of trends in neutral density at 400-km altitude (not shown) is less pronounced, but also features an enhanced negative trend in approximately the same region as the strong cooling patch in the Northern Hemisphere.
Since the spatial variations in the neutral temperature and density trends at 400 km cannot be explained by the effects of an increase in CO 2 concentration, let us examine whether changes in the Earth's magnetic field might play a role. Figure 4 shows the trend in the vertically integrated Joule heating power, also calculated with linear regression Model 3. This shows a strong negative trend (up to nearly -30%/decade) in roughly the same location as the strong negative trend patch in neutral temperature in the Northern Hemisphere. The negative Joule heating power trend peaks somewhat more northeastward, around 50-70 • N, 80-20 • W, and is more localized, but Joule heating effects are likely to be smeared out by neutral winds. It therefore seems likely that reduced Joule heating power is responsible for the patch of strongly negative temperature trends centered at ∼50 • N, ∼70 • W. Other high-latitude regions also show decreases in Joule heating power, but these are considerably smaller. Still, the decreases in Joule heating power at high magnetic latitudes may explain why thermosphere temperature trends are generally larger in these regions, and it might also be responsible for global mean temperature trends being larger than expected from the increase in CO 2 concentration alone. We note that there are some relatively large low-latitude changes in Joule heating power at ∼60-0 • W as well, which are clearly associated with the relatively large movement of the magnetic equator in this longitude sector. However, these are unlikely to have a significant effect on the thermosphere temperature, as the absolute magnitude of the Joule heating power at low latitudes is much smaller than at high latitudes.  Influences of solar and geomagnetic activity were removed as much as possible from the Joule heating time series by the regression analysis, so that the main source of the trend in Joule heating shown in Figure 4 is likely to be the changes in the Earth's main magnetic field. While the overall magnetic field strength has been decreasing, this is dominated by the expansion and deepening of the South Atlantic Anomaly region (Thébault et al., 2015), far away from the auroral zones. In the high-latitude regions near the magnetic poles, the main magnetic field has for the most part been strengthening slightly, especially in the Northern Hemisphere. A strengthening of the main field reduces the ionospheric conductivity, and thereby acts to reduce the Joule heating, adding to the CO 2 cooling effect. In addition to any effects of changes in field strength, changes in the geographic positioning of the magnetic poles, and therefore the auroral zones, are probably important as well, again especially in the Northern Hemisphere, where the magnetic pole has moved much more than in the Southern Hemisphere. Further evidence to support the idea that the decrease in Joule heating and the large cooling in the Northern Hemisphere region around ∼50 • N, ∼75 • W are associated with main magnetic field changes comes from controlled model experiments by Cnossen (2014) which isolated the role of main magnetic field changes from 1900 to 2000 on the neutral temperature in the thermosphere. These experiments also showed a strong cooling patch in the Northern Hemisphere high magnetic latitude region, similar to what we find here. Figure 5 shows maps of the long-term trends in h m F 2 and N m F 2 based on regression Model 3. Again, the spatial variation obtained with the other regression models is very similar. The spatial pattern of trends in TEC looks very similar in structure to that of the trends in N m F 2 and is therefore not shown. The spatial structure of trends in h m F 2 and N m F 2 indicates a superposition of CO 2 and magnetic field effects, as expected. The low latitude to midlatitude trends around ∼60 • W to 20 • E are clearly associated with main magnetic field changes (especially the movement of the magnetic equator), matching quite well with expectations from previous modeling studies (Cnossen, 2014;Cnossen & Richmond, 2008. There also appear to be small effects on h m F 2 and N m F 2 of the movement of the magnetic pole and associated auroral region, in particular in the Northern Hemisphere, where this movement has been more pronounced. The overall decrease in h m F 2 elsewhere is at least in qualitative agreement with the expected effect of the increase in CO 2 concentration, although trend magnitudes appear somewhat high compared to Solomon et al. (2018Solomon et al. ( , 2019. Trends in N m F 2 in contrast may appear very small outside the regions that are strongly affected by main magnetic field changes, but we note that the white regions on the map still correspond to trends of up to -1%/decade, which is approximately the order of magnitude expected from the increase in CO 2 concentration based on Solomon et al. (2018Solomon et al. ( , 2019. Oliver et al. (2013) suggested that an increase in gravity wave (GW) forcing could potentially contribute to cooling trends in the thermosphere. GW effects are parameterized in WACCM-X, but the implemented scheme, which includes both orographic and nonorographic waves, links the sources of the GWs, as well as their propagation and dissipation, to atmospheric quantities calculated by the model (Garcia et al., 2017;Richter et al., 2010). Both GW generation and GW propagation and dissipation are therefore able to respond to long-term changes in the simulated atmospheric conditions. GWs can affect the energy and momentum budget of the thermosphere through the vertical transport of horizontal momentum, heating associated with wave dissipation, and diffusive mixing.

Indirect Effects From the Lower Atmosphere
While GW breaking initially results in heating due to the conversion of mechanical energy into heat, the process also results in a downward heat flux so that the net effect is a cooling of the thermosphere (Walterscheid, 1981). Yigit and Medvedev (2009) simulated average GW cooling rates of up to 170 K/day, so long-term changes in GW forcing could potentially have an important effect on the thermosphere temperature. However, the temperature tendencies associated with GW breaking in our simulation tend to be much smaller and are also at least an order of magnitude smaller than the CO 2 cooling rates or Joule heating rates. While we found that long-term trends in GW-induced temperature tendencies are not necessarily small in percentage terms (up to 10-30%), their absolute magnitude should be too small to matter for the long-term trends in thermosphere temperature we find here.
The diffusive mixing induced by GW breaking can affect the thermosphere composition and density. The global mean GW-induced eddy diffusion coefficient peaks at ∼85-95-km altitude and shows a small (up to 2%/decade) long-term enhancement in that altitude range. However, it is not clear based on just our transient simulation, whether this could have contributed significantly to the calculated thermosphere density trend. Similarly, any indirect effects of changes in large-scale atmospheric circulation that might result from changes in momentum deposition by altered GW forcing cannot be readily separated from the effects of more direct, in-situ processes.
To isolate the potential effects of lower atmosphere climate change on the upper atmosphere, we analyzed the mean differences between the N1 and N0 experiments that were done specifically to quantify this aspect. No significant differences in h m F 2 , N m F 2 , and TEC were found. Differences in the neutral temperature and mass density at 400-km altitude were small (around ±1-2 K and <1%, respectively) and not significant either. Below ∼100 km, in the upper mesosphere and lower thermosphere region, differences in mass density of up to -4% (corresponding to ∼−1.3%/decade) were found that were largely significant (not shown), while neutral temperature differences in the same region were again not significant.
These results suggest that indirect effects of climate change in the lower atmosphere (<50 km) on the upper atmosphere (>100-km altitude) are negligible. However, we must bear in mind that WACCM-X is not a perfect reflection of reality and certainly in terms of representing the effects of GWs there are likely to be inaccuracies due to the simplifying assumptions made in the parameterization scheme. For instance, the scheme does not allow for any horizontal wave propagation and effects of secondary waves generated in the middle atmosphere are not included. Recent studies have shown that secondary waves could actually be quite important for the thermosphere Bossert et al., 2017;. Also, the largest GW cooling rates in WACCM-X were about a factor 3 smaller than those calculated by Yigit and Medvedev (2009). While we do not know whether their results are necessarily correct, it is certainly possible that WACCM-X could be underestimating thermal GW effects in the thermosphere. This means that the indirect, dynamical effects of climate change in the lower and middle atmosphere on the upper 10.1029/2020JA028623 atmosphere are still uncertain. Improvements in GW parameterization schemes for whole-atmosphere models such as WACCM-X are essential to make progress on this issue.

Discussion and Conclusions
We analyzed long-term trends in several upper atmosphere parameters based on a transient simulation with WACCM-X from 1950 to 2015, including all known drivers of upper atmosphere climate change. While we concentrated on trends for the full time interval of 1950-2015, trends can in principle be calculated for any subinterval at any location in the world, facilitating direct comparisons between model-based trend estimates and observational trend analyses.
Several different multilinear regression models were explored to extract trends from time series that are dominated by solar and geomagnetic activity variations as part of the ∼11-year solar cycle. We found that the most commonly used linear regression models, here referred to as models 1 and 2, struggle to adequately remove solar cycle variations from a long-term data series, even with 40+ years' worth of continuous monthly mean data. This is worrying, as most observational data records will be considerably shorter and may contain data gaps, making it even harder to suppress solar cycle influences successfully. It therefore seems especially important when analyzing long-term observational data sets to check carefully how well solar cycle influences have been removed. This can be done, for instance, by plotting the trends obtained as a function of the end year of the time series, as we did here. This offers more insight in the magnitude of any remaining solar cycle influence and provides important information on the reliability of the trends obtained. Simply reporting trends without such checks could actually be very unhelpful, as the trends might be unduly influenced by solar cycle influences, while we are trying to build a global picture of the "true" trends, that is, those trends unrelated to solar cycle variations.
For some variables, namely the neutral density at 400-km altitude, N m F 2 and TEC, an improved fit to the data was obtained by adding a term for (F10.7a) 2 to the linear regression model. This also reduced remaining solar cycle influences in the trends obtained for these variables. For the neutral temperature at 400-km altitude and h m F 2 the additional (F10.7a) 2 term did not make much difference. A few observational studies used a similar (solar activity) 2 term in their trend analyses (e.g., Emmert, 2015;Zhang & Holt, 2013), but it does not appear to be standard practice to include such a term (e.g., Laštovička & Jelínek, 2019). It would be interesting to test more widely how adding a (solar activity) 2 term in observational linear regression-based trend analyses would affect trend estimates and any spatial patterns in trends.
Spatial variations in trends can provide important clues on trend drivers. Our results confirm that CO 2 is probably the main driver of trends in the thermosphere, but at high (magnetic) latitudes, effects of changes in the Earth's magnetic field appear to be nonnegligible. Main magnetic field changes are likely responsible for a long-term decrease in Joule heating, which is especially important in the Northern Hemisphere high-latitude region. This is presumably due to the relatively large movement of the Northern Hemisphere magnetic pole and the associated auroral region. The largest thermosphere temperature trends were located in approximately the same region as the strongest trends in Joule heating and were likely the result of a superposition of CO 2 cooling effects and reduced Joule heating associated with the change in the Earth's magnetic field.
We note that Millstone Hill observatory (42.6 • N, 71.5 • W) is located within the patch of relatively strong cooling shown in Figure 3. This should make it an ideal location to examine the combined long-term effects of the increase in CO 2 concentration and main magnetic field changes. Previous work based on simulations by Cnossen and Richmond (2013) indicated that main magnetic field changes might explain ∼8% of the observed trend of ∼−4 K/decade in ion temperature at Millstone Hill (Zhang & Holt, 2013), that is, a relatively small percentage. It would be interesting to revisit this point with the new simulation results presented here. Nevertheless, it is unlikely that main magnetic field changes can explain a much larger part of the observed ion temperature trend at Millstone Hill, given that the ion temperature trend at Saint Santin in France (Donaldson et al., 2010) is of similar order of magnitude, while being located in a region of little impact from geomagnetic field changes.
The spatial patterns of trends in h m F 2 , N m F 2 , and TEC (the latter similar to the pattern in N m F 2 ) clearly point to a superposition of CO 2 and main magnetic field effects, with the latter dominating the trends in the region of ∼50 • S to 20 • N, ∼60 • W to 20 • E, in good agreement with previous studies (Cnossen, 2014; 10.1029/2020JA028623 Cnossen & Richmond, 2013). While trends associated with main magnetic field changes can be either positive or negative, depending on the location, patches of negative trends are considerably stronger and larger in size than patches of positive trends. This is not simply because of the increase in CO 2 concentration producing negative trends; it can also be seen in controlled experiments isolating the effect of main magnetic field changes on h m F 2 and f o F 2 (Cnossen, 2014;Cnossen & Richmond, 2013). This can therefore also push global mean trends to be more negative than they would be due to the increase in CO 2 concentration alone. This is true for the neutral temperature at 400 km as well and may help to explain why the trends in especially the neutral temperature and h m F 2 we found remained somewhat large compared to Solomon et al. (2018Solomon et al. ( , 2019 even after solar cycle influences were removed as much as possible. The possible role of climate change in the lower/middle atmosphere on the upper atmosphere, for instance, via altered GW forcing, is still an open question. Our simulations suggest that indirect, dynamical effects of climate change below 50-km altitude on the upper atmosphere (>100-km altitude) are small and not significant. However, due to simplifying assumptions made in the GW parameterization and the omission of secondary waves, we cannot rule out that GW effects could be more important for the upper atmosphere in the real world, and could have a larger effect on long-term trends in the thermosphere and ionosphere than we simulated here. Improvements in GW paramaterization schemes for whole-atmosphere models such as WACCM-X are essential to achieve closure on this point.