Response of atmospheric carbon dioxide to the secular variation of weakening geomagnetic field in whole atmosphere simulations

Responses of atmospheric carbon dioxide (CO2) density to geomagnetic secular variation are investigated using the Whole Atmosphere Community Climate Model‐eXtended (WACCM‐X). Our ensemble simulations show that CO2 volume mixing ratios (VMRs) increase at high latitudes and decrease at mid and low latitudes by several ppmv in response to a 50% weakening of the geomagnetic field. Statistically significant changes in CO2 are mainly found above ~90 km altitude and primarily redetermine the energy budget at ~100‐110 km. Our analysis of transformed Eulerian mean (TEM) circulation found that CO2 change is caused by enhanced upwelling at high latitudes and downwelling at mid and low latitudes as a result of increased Joule heating. We further analyzed the atmospheric CO2 response to realistic geomagnetic weakening between 1978 and 2013, and found increasing (decreasing) CO2 VMRs at high latitudes (mid and low latitudes) accordingly. For the first time, our simulation results demonstrate that the impact of geomagnetic variation on atmospheric CO2 distribution is noticeable on a time scale of decades.


Introduction
Previous studies indicate that variation in the geomagnetic field has the potential to influence climate changes in the upper atmosphere by affecting the motion of ions and electrons, ion-neutral coupling, and thus the energy budget, despite that changes in the geomagnetic field do not act directly on the neutral atmosphere. The underlying physical mechanism has been investigated through idealized modeling (Cnossen et al., 2011Zossi et al., 2018). One widely accepted view is that geomagnetic field changes lead to different ionospheric conductivities and electric fields, which results in different Joule heating rates that influence the temperature, dynamics, and composition of the neutral upper atmosphere.
However, the effect of geomagnetic field variation on the middle and lower atmosphere down to the surface remains controversial. Wollin et al. (1971) first proposed a linkage between the geomagnetic field and surface climate based on the study of deep-sea sediment cores. Researchers further found a negative correlation between the strength of geomagnetic fields and global mean surface temperature, inferred from the measured rise and fall of global sea level as well as the advance and retreat of Alpine glaciers (e.g., Gallet et al., 2005;De Santis et al., 2012). Kitaba et al. (2013) also suggested that tropospheric cooling could decrease the geomagnetic field. Although several lines of observational evidence have suggested a correlation between the geomagnetic field and climate, the mechanism is unclear. One hypothesis, the so-called "umbrella effect", postulates that a strong geomagnetic field could reduce cosmic rays flux, affecting cloud formation (e.g., Kitaba et al., 2017;Svensmark et al., 2017). Others have proposed that geomagnetic field modulation in the upper atmosphere would have downward influences on chemical and dynamic processes in the lower atmosphere, including wave propagation (e.g., Arnold and Robinson, 2001;Seppälä et al., 2013;Cnossen et al., 2016), and the descent of odd nitrogen (e.g., Randall et al., 2007). To better understand the underlying mechanism, simulations by a climate model covering the entire atmosphere region would be helpful, since they may shed light on the impacts of geomagnetic field variations on the dynamics and chemistry of the upper atmosphere to the middle and lower atmosphere. However, the interpretation of whole-atmosphere climate simulations is complicated by the relatively small climate signal in comparison to large climate variability. Hence, an ensemble with adequate samples is needed (the detailed settings in this work are discussed in Section 2). It is also noted that testing the "umbrella effect" is still challenging; although cloud microphysics have been parameterized in the climate model, the relationship to cosmic rays is still under debate.
As a radiatively active greenhouse gas, carbon dioxide (CO 2 ) plays a key role in the energy budget of the entire atmosphere. Previous studies have established that CO 2 increase produces global warming in the troposphere and cooling in the stratosphere and above, because it absorbs and radiates in the infrared band which becomes optically thin with increasing altitude (Fels et al., 1980;Dickinson, 1984). In particular, the cooling effects of CO 2 in the upper atmosphere were first studied by Roble and Dickinson (1989) and confirmed by later studies (e.g., Akmaev and Fomichev, 2000;Akmaev et al., 2006;Qian KY et al., 2006Fomichev et al., 2007;Garcia et al., 2007;Marsh et al., 2013;Solomon et al., 2015Solomon et al., , 2018. The effects of CO 2 change versus geomagnetic field change on the climate of the upper atmosphere and ionosphere have been compared previously (e.g. Cnossen, 2014;Yue XN et al., 2018). Based on the Thermosphere Ionosphere Electrodynamics General Circulation Model (TIEGCM), the two terms are treated as independent factors and then used to estimate their contributions to upper atmospheric and ionospheric climate changes. CO 2 in the model is specified as a boundary condition in the lower thermosphere and is uniformly distributed horizontally, although without consideration of the spatial and temporal variability of CO 2 in the upper atmosphere. Considering the CO 2 variability, Cnossen (2020) analyzed the long-term database of WACCM-X using a modified regression method to examine the impact of CO 2 and geomagnetic field changes on upper atmosphere climate change. Based on a whole atmosphere model, GAIA, the CO 2 doubling effect in the upper atmosphere and the dependence of geomagnetic activity upon the CO 2 -driven trend have been investigated recently (Liu HX et al., 2020, 2021. Nonetheless, to our knowledge, the interaction between CO 2 change and geomagnetic field variation has not been investigated before. CO 2 variability has been confirmed to be coupled with atmospheric dynamics (e.g. Rezac et al., 2015 and references therein), and this could also be affected by geomagnetic field variation. Thus, it is conceivable that changes in the geomagnetic field would influence CO 2 distribution and in turn, the energy budget.
In this study, we examine the potential influence of geomagnetic field variations on atmospheric CO 2 and its energy budget. We perform several sets of numerical experiments using the Whole Atmosphere Community Climate Model-eXtended (WACCM-X) v2.1. To test the CO 2 response, we artificially reduce the strength of the geomagnetic field to 50% of its original level. This helps us establish whether CO 2 would redistribute under such extreme geomagnetic field changes. By controlling for other factors, this whole-atmosphere sensitivity study can reveal the role of the geo-magnetic field in CO 2 changes and suggest the corresponding mechanism. To determine whether a measurable effect exists for realistic geomagnetic field changes, we compare the CO 2 response to different geomagnetic fields in the past three decades (from 1978 to 2013).

Overviews of WACCM-X and Simulation Setups
WACCM-X v.2.1 is one of the atmospheric components of the National Center for Atmospheric Research (NCAR) Community Earth System Model (CESM), covering Earth's entire atmosphere extending from the surface to the exobase. Detailed descriptions of WACCM-X can be found in Liu HL et al. (2018). Briefly, the upper boundary of the model is 4.1 × 10 −10 hPa (~500-700 km, depending on solar activity), with one quarter of the scale height resolution in the mesosphere and thermosphere, and a horizontal resolution of 2.5° × 1.9° (longitude × latitude). The physics package is based on CAM4 and WACCM4 (Marsh et al., 2013;Neale et al., 2013). High-latitude electric potential and auroral particle precipitation are specified by empirical models (Heelis et al., 1982;Roble and Ridley, 1987). The configurations of the geomagnetic field given by the International Geomagnetic Reference Field (IGRF) (Thébault et al., 2015) and electrodynamic processes are calculated based on the Magnetic Apex coordinate (Richmond, 1995). In WACCM-X, local thermodynamic equilibrium (LTE) and non-LTE processes are considered for different altitudes; these processes are defined in detail by Andrews et al. (1987). Below the mesosphere and lower thermosphere (MLT) region, the standard longwave formulation in CAM is used; in the region between 85− 110 km, as the atmospheric density decreases, non-LTE effects strengthen and are calculated using the recurrence formula in Fomichev et al. (1998). Above 110 km, as in TIEGCM, the coolingto-space approximation is assumed. This implies that the radiative cooling of a given atmospheric layer is controlled by the emission directly cooling to space, while the exchange from above and below is offset and therefore can be ignored (Jeevanjee and Fueglistaler, 2020). In addition, horizontal variability of CO 2 concentration is considered to be coupled with neutral dynamics. The more comprehensive consideration of the CO 2 cooling term in WACCM-X enables us to study its response to geomagnetic field variations from the thermosphere down to the mesosphere.
To investigate whether CO 2 would be affected by geomagnetic field variations, other external forces should be held constant: specifically, solar activity (F = 70 sfu.) and geomagnetically quiet conditions (Kp = 2). To exclude inter-annual variability from below, anthropogenic emissions and sea surface temperature at the lower boundary are set as the average condition around the year of 2000. To more clearly identify the signal over the large variability of the system, we artificially reduce the strength of the geomagnetic field to 50% of the original level (base case) and keep the inclination and declination angle constant by modifying the IGRF coefficients in the model. The two cases are called half-B control and original-B base case, hereafter. To enlarge the signal in the simulations, we chose to set Kp = 2 rather than use an extreme quiet condition; e.g. 0.33. 100-year continuous simulations are performed for each case to form ensembles with samples large enough for statistical analysis. Notably, although external forces are kept constant in the 100-year perpetual runs, the atmospheric state is different for each year due to internal atmospheric variability. Figure S1 depicts the variations of zonal mean temperature from January to March in the original-B case. The standard deviation of zonal mean temperature among the ensembles is 1-3 K, which is comparable to the magnitude of day-to-day variability for an individual sample. The result indicates that the runs do in fact generate statistics in the ensembles. Although this hypothetical change of geomagnetic fields is large in comparison to the geomagnetic dipole moment, which decreases at a rate of 5% on century scales (Gubbins et al., 2006), it could still be valuable for studies on the paleoclimate time scale. For instance, a Holocene geomagnetic field model, CASL10k, demonstrates that the dipole moment varies from ~6 × 10 22 to 1 × 10 23 Am 2 spanning the past 10 kys (Korte et al., 2011;Constable et al., 2016). For the past several million years, analysis of paleomagnetic records indicates that the strength of the geomagnetic dipole moment could diminish to below 2 × 10 22 Am 2 (Valet and Meynadier, 1993;Guyodo and Valet, 1999). This simulation is also beneficial for investigating how deep the impact of such extreme changes may extend from the upper atmosphere to the middle atmosphere.
To examine the effect of realistic geomagnetic condition changes during the modern era, we carry out two additional sets of simulations. These are time-slice simulations (Solomon et al., 2018) with the anthropogenic emission rates set at current-day levels (2011)(2012)(2013)(2014)(2015), compared to historic geomagnetic fields, specifically for 1978 vs. 2013 (labeled B1978 and B2013 hereafter). Figure S2 shows the geomagnetic intensity, inclination angle, and declination angle in the two configurations and their respective differences. Each case continuously runs 5 years after a 1-year spin-up. The 5-year interval forms a small ensemble to reduce the impact of inter-annual variability. Concentrations of anthropogen-ic gases are specified by the historical time-dependent lower boundary condition.

Results and Discussion
As shown in Figures 1a−1c, zonal mean CO 2 volume mixing ratios (VMRs) of the half-B case increase by ~3-10 ppmv at high latitudes for different seasons and decrease by ~2 ppmv at low latitudes and the equator in the lower thermosphere (~100 to 150 km), compared with the original-B case. In the summer hemisphere, the increase of CO 2 is larger than that in the winter hemisphere, and the former has double peaks at ~8.9 × 10 −6 hPa (~140 km) and 1.4 × 10 −4 hPa (~103 km), whereas the latter has only one peak at ~8.9 × 10 −6 hPa (~140 km). The student's t-test is used to calculate the statistical significance of these differences. The increase of CO 2 in the summer hemisphere due to the change of geomagnetic fields can extend to the mesosphere (~1.7 × 10 −3 hPa, ~88 km) relative to the winter hemisphere (~2.4 × 10 −5 hPa, ~115 km), with 99% statistical significance level (p value is 0.01). As a comparison, the anthropogenic change in CO 2 over past decades in the mesosphere and the lower thermosphere is 5.5% per decade (Qian LY et al., 2017), which is approximately 5-10 ppmv per decade at the pressure levels of 10 −4 to 10 −3 hPa.
Differences in CO 2 cooling rates between the half-B and original-B cases are depicted in Figures 1d-1f, which denotes different features from the change of CO 2 VMRs. As the geomagnetic field weakens, CO 2 cooling increases by ~2-4 K/day above ~100 km, with a larger value for high latitudes and a smaller value for middle and low latitudes. Below ~100 km, CO 2 cooling weakened slightly at high latitudes in summer, and in both hemispheres during equinoxes. The response of CO 2 cooling to the weakening geomagnetic field is also deeper in the summer hemisphere,  whereas the cooling rate decreases more in the winter hemisphere; this was statistically significant using the student's t-test.
The enhanced cooling effect in the thermosphere is stronger in the winter hemisphere than in the summer hemisphere. Notably, the CO 2 cooling is not only dependent on its abundance, but is also determined by temperature and atomic oxygen (O) abundance. In the thermosphere, the temperature significantly increases due to enhanced Joule heating. The corresponding change of O is shown in Figure S3, decreasing over 0.15 mol/mol at high latitudes and increasing ~0.02 mol/mol at low latitudes. The relative decrease in O VMRs reaches a maximum of ~20% for high latitudes at the height of ~150 km. The zonal mean CO 2 VMR and CO 2 cooling rate during the equinox for the original-B case are shown in Figures 1g and 1h. The changed CO 2 VMR at high latitudes is about 1-10% for the height of 90-200 km, and the changed CO 2 cooling rate is ~5-10% throughout the thermosphere. Figure 2a compares the magnitude change of global mean CO 2 cooling rates with other important energy sources/sinks, and the corresponding heating/cooling rates in the original case are shown in Figure 2b as a reference. Joule heating is the dominant energy source that differs between half-B and original-B cases above ~115 km. Ionospheric conductivity is inversely proportional to the strength of the geomagnetic field , and the Joule heating rate is proportional to Pedersen conductivity. Thus, a weak geomagnetic field leads to stronger Joule heating for the thermosphere. In response to increased Joule heating, nitric oxide (NO) plays a more important role in cooling the thermosphere at 5.3 μm (e.g., Kockarts, 1980), because the NO cooling rate is sensitive to temperature and quickly reacts to energy inputs, serving as a "natural thermostat" for the thermosphere (e.g., Mlynczak et al., 2003Mlynczak et al., , 2018Lu G et al., 2010). This is seen in Figure 2a: the change of NO cooling is the dominant term from 115 km to 220 km and is over -20 K/day at the peak (~100%). A detailed comparison of NO VMRs and corresponding cooling effects is shown in Figure 4. Generally, CO 2 ceases to be important in the heat budget of the lower thermosphere above ~140 km (Mlynczak et al., 2018). In this case, CO 2 cooling increases bỹ 2 K/day throughout the thermosphere above ~105 km (~5%) when the geomagnetic field weakens. Between 100−110 km, CO 2 is the primary changed cooling agent. The change in another cooling term, O(3P) cooling at 63 μm, is relatively insignificant (less than ±2 K/day for all altitudes).
As a trace species with a long lifetime, the transport of CO 2 and its v * w * spatial distribution in the mesosphere is determined by atmospheric general circulation, which is quantified by the residual mean circulation Mcintyre, 1976, Andrews et al., 1987). The meridional and vertical components of the residual mean circulation, and , are expressed as follows: where z is the log-pressure altitude, defined as .
is the constant approximate scale height of 7 km, and is the mean surface pressure. and are the zonal mean meridional and vertical velocity, respectively. is the Earth's mean radius of 6371 km, is latitude, and is the potential temperature. is the basic density, defined by , where is the gas constant of 287 J/kg/K, and is mean surface temperature. The notation of prime indicates a departure from zonal mean fields. 110 km. The primary driver of the summer to winter circulation in the mesosphere/mesopause region is gravity wave forcing. In the region between the altitudes, there are several clockwise rotation cells with weak upward motion in winter and downward motion in summer at high latitudes, which are also driven by gravity wave forcing; this flow pattern is consistent with Fig. 1 in Smith et al. (2011). The opposing circulations below and above 90 km lead to the vertical gradient of CO 2 being larger in the summer hemisphere than the winter hemisphere around the mesopause, as shown by the contour lines in Figure 3a. To illustrate the effect on the residual circulations caused by the weakening geomagnetic field, Figure 3b displays the differences in residual circulation between half-B and original-B cases, normalized to unit vector length, in the same month of June in the region where the change of CO 2 VMRs is large. It also presents the corresponding differences of the CO 2 VMRs in contour lines. Normalizing the wind vector better illustrates the flow direction at all altitudes. The strong enhanced upward/equatorward circulation above ~105 km is found at high latitudes in both hemispheres, with downward flow at mid and low latitudes. This change is caused by enhanced Joule heating rates above 115 km. Consequently, upwelling at high latitudes could bring CO 2 -rich air from below, and downwelling at mid and low latitudes transports CO 2 from the thermosphere where it is rare. The increased upward flow extends to ~90 km and 100 km in the summer and winter hemispheres, respectively. In addition, the vertical gradient of CO 2 is generally larger in the summer hemisphere (contour lines in Figure 3a) such that upwelling in this hemisphere brings more CO 2 from below, increasing CO 2 VMRs (Figures 1b and 1c). Figure 4 illustrates the differences among CO 2 VMRs, CO 2 cooling rates, and neutral temperature between the B2013 and B1978 cases with realistic changes in the geomagnetic field for the northern winter at a pressure level of 6.6 × 10 −5 hPa (~110 km). Significant enhancement of CO 2 VMRs is found at high latitudes (~3-4 ppmv), especially in the southern hemisphere where the intensity of the geomagnetic field is consistently smaller in 2013 than in 1978 ( Figure S1d). At mid and low latitudes, CO 2 VMRs decrease by ~1-2 ppmv. Correspondingly, the CO 2 cooling rate increases at high latitudes by ~2 K/day and decreases at mid and low latitudes by ~1 K/day. As for neutral temperature changes, this significantly increases above ~115 km due to enhanced Joule heating. Below this region dominated by Joule heating, the cool-ing effect becomes important. Figure 4c indicates that the neutral temperature changes at ~110 km, which decreases by ~2-3 K at high latitudes in the southern hemisphere with statistical significance. This temperature decrease should also be influenced by another important cooling effect: enhanced upwelling induced adiabatic cooling (Li JY et al., 2018, 2019. The relative contribution of the dynamic effect and CO 2 radiative effect need to be further studied in detail, although the latter seems to be small according to our analysis. Retrospective comparison between half-B and Original-B cases in Figure 2 confirms that NO cools the atmosphere above ~110 km, but is not the dominant factor of the energy budget at 100-110 km. In addition, compared with CO 2 changes between half-B and original-B cases, the magnitude of the CO 2 difference for B2013 is ~30% less than B1978, whereas the corresponding geomagnetic dipole moment decreases only by ~2%, which is ~1/25 of the 50% change found in the former experiments. This inconsistency in magnitude may be partially caused by three reasons: (a) B2013 and B1978 cases include larger atmospheric interannual variability, and half-B and original-B cases are perpetually run for the year 2000; (b) a large number of members in half-B and original-B would reduce the effect of atmospheric internal variability; and (c) changes in circulation between B2013 and B1978 cases should be more complex because the inclination and declination angles of the geomagnetic field are also different. The nonlinear response of CO 2 VMR to the change in geomagnetic field intensity is another interesting factor that prompts future investigation.
Although we did not determine whether surface climate would change due to geomagnetic field variation in the present work, according to the known physics included in WACCM-X which lacks a coupling ocean, a large geomagnetic field change can have effects down to the mesosphere, but no significant change below that (Figure 1). CO 2 changes in the middle and upper mesosphere presented here might provide a progressive insight that geomagnetic field variation could influence neutral dynamics, distribution of neutral compositions, and thermal structure of the middle atmosphere. When considering other possible factors related to the geomagnetic field (cosmic rays and particle precipitation, etc.), and more extreme geomagnetic conditions (e.g. polarity reversals), in the future the CO 2 and temperature responses might be expected to be larger and may even reach lower altitudes. . Differences in (a) VRM of CO 2 , (b) CO 2 heating rates, and (c) temperature between B2013 and B1978 cases averaged for winter (November to February) at a pressure level of 6.6×10 −5 hPa (~110 km). The gray dots indicate the area where the differences are statistically significant at a 90% level. The black lines are continental coastlines.

Conclusion
In this study, CO 2 distribution responses to changes in the geomagnetic field are investigated based on WACCM-X simulations. An artificially extreme experiment reducing geomagnetic field strength to 50% of the original level illustrates that the response exists in the atmosphere above ~90 km with 99% statistical significance. CO 2 VMRs increase at high latitudes and decrease elsewhere. CO 2 changes significantly determine the energy budget at 100-110 km. We propose that the changed transformed Eulerian mean (TEM) circulation and the associated increased Joule heating are responsible for these CO 2 changes. Another experiment based on realistic geomagnetic field changes between 1978 and 2013 demonstrated that such an effect is still tenable, indicating that the linkage persists on a time scale of decades. This paper reveals for the first time that CO 2 distribution may be influenced by geomagnetic variation, and that climate change may be induced by secular variations in geomagnetic fields. ΔQNO Figure S4. Differences of the zonal mean (top row) mixing ratio and (bottom row) heating rate of NO in (left) equinox, (middle) northern winter, and (right) northern summer. The grey dots indicate the area where the differences are statistically significant at the 99% significance level. The approximate height of corresponding pressure level is labelled on the right.