An empirical model of the thermospheric mass density derived from CHAMP satellite

In this study, we present an empirical model, named CH-Therm-2018, of the thermospheric mass density derived from 9-year (from August 2000 to July 2009) accelerometer measurements from the CHAllenging Mini-satellite Payload (CHAMP) satellite at altitudes from 460 to 310 km. The CHAMP dataset is divided into two 5-year periods with 1-year overlap (from August 2000 to July 2005 and from August 2004 to July 2009) to represent the high-to-moderate and moderate-to-low solar activity conditions, respectively. The CH-Therm-2018 model describes the thermospheric density as a function of seven key parameters, namely the height, solar flux index, season (day of year), magnetic local time, geographic latitude and longitude, as well as magnetic activity represented by the solar wind merging electric field. Predictions of the CH-Therm-2018 model agree well with CHAMP observations (within 20 %) and show different features of thermospheric mass density during the two solar activity levels, e.g., the March–September equinox asymmetry and the longitudinal wave pattern. From the analysis of satellite laser ranging (SLR) observations of the ANDE-Pollux satellite during August–September 2009, we estimate 6 h scaling factors of the thermospheric mass density provided by our model and obtain the median value equal to 1.267±0.60. Subsequently, we scale up our CH-Therm-2018 mass density predictions by a scale factor of 1.267. We further compare the CH-Therm-2018 predictions with the Naval Research Laboratory Mass Spectrometer Incoherent Scatter Radar Extended (NRLMSISE-00) model. The result shows that our model better predicts the density evolution during the last solar minimum (2008–2009) than the NRLMSISE-00 model.

1 Introduction

The thermosphere is the top layer of the gravitationally bound part of the atmosphere, which is partly ionized and extends from about 90 to over 600 km (Lühr et al., 2004). Its density variations are mainly driven by the extreme solar ultraviolet (EUV) irradiance, the energetic particles and electrical energy from the magnetosphere and solar wind, as well as by waves originating in the lower atmosphere that propagate upward into the thermosphere. The thermospheric mass density in general falls off exponentially with increasing altitude, with scale heights of about 25 to 75 km in the upper atmosphere, depending on altitude and solar flux levels. In addition to the vertical variation, the mass density also varies horizontally (latitude and longitude) as well as with solar flux, geomagnetic activity, season and local time (Emmert, 2015).

The thermosphere plays a crucial role for near-Earth space operations, as the total mass density is the key parameter for orbit perturbation of low Earth orbit (LEO) satellites. Therefore, knowledge of the thermospheric density is critical in the planning of LEO missions, such as their orbital altitudes, lifetime and re-entry prediction. As the ionosphere is embedded in the thermosphere, the knowledge of thermospheric density will also help to improve our understanding of the coupling between the thermosphere, ionosphere and lower atmosphere (Liu et al., 2013; Emmert, 2015).

There are several tools for measuring the thermospheric mass density. The atmospheric drag provides the most direct means, which can be measured by onboard accelerometers (e.g., Champion and Marcos, 1973; Lühr et al., 2004; Doornbos et al., 2010) or estimated from the changes of LEO object trajectories (e.g., King-Hele, 1987; Emmert et al., 2004). Other instruments, such as neutral mass spectrometers (e.g., von Zahn, 1970; Hedin, 1983), ultraviolet remote sensing (e.g., Meier and Picone, 1994; Christensen et al., 2003) and the pressure gauge mounted on rockets (e.g., The Rocket Panel, 1952; Clemmons et al., 2008), can also be used for inferring the mass density. The details of these techniques have been reviewed by several earlier studies (e.g., Osborne et al., 2011; Clemmons et al., 2008; Emmert, 2015). Various empirical models have also been developed to describe the thermospheric mass density variability. The most widely used are the Mass Spectrometer Incoherent Scatter Radar Extended (MSISE) model family (Hedin, 1991; Picone et al., 2002), the Drag Temperature Model (Bruinsma et al., 2003, 2012) and the Jacchia-Bowman 2008 (JB2008) model series (Bowman et al., 2008a, b). Liu et al. (2013) and Yamazaki et al. (2015) reported two empirical models derived from recent LEO missions, such as the CHAllenging Minisatellite Payload (CHAMP, Reigber et al., 2002) and the Gravity Recovery and Climate Experiment (GRACE, Tapley et al., 2004). These two models represent well the prominent thermospheric structures at low latitudes like the equatorial mass density anomaly (EMA) and the wave-4 longitudinal pattern, as well as the solar wind influence on the high-latitude thermosphere, respectively.

As reported by previous studies, the height and solar activity are the two most important factors that affect the thermosphere mass density (Liu et al., 2005; Guo et al., 2008; Lei et al., 2012). The CHAMP altitude decreased coincidentally within the declining phase of solar cycle 23. Therefore, it is difficult to fully separate the height and solar activity effects on mass density from CHAMP observations. By assuming a linear dependence on height variation, Liu et al. (2013) used the dataset from 2002 to 2005 when CHAMP was at the altitude of 420 to 350 km to construct a model, focusing on low altitudes and mid-latitudes. They argued that a linear approximation is applicable within an error of about 3.5 % over one scale height. To reduce the height variation effects on the model, Yamazaki et al. (2015) used the MSISE-00 model to normalize the CHAMP and GRACE densities to a common height of 450 km, focusing on high latitudes. However, as the MSISE-00 model was not accurate during the extreme solar minimum of 2008 to 2009 (Thayer et al., 2012; Liu et al., 2014), it would possibly affect their height correction during the solar minimum period; therefore, they also used the dataset from 2002 to 2006. Both models mentioned above considered only the dataset from high to moderate solar activity, while the dataset from the solar minimum (2008 to 2009) has not been included.

Different to Liu et al. (2013) and Yamazaki et al. (2015), we take into account in this study the dataset from August 2000 to July 2009 for constructing our empirical models of the thermospheric mass density, to make more efficient use of the CHAMP observation. This period includes high and low solar activities and the CHAMP satellite altitude varies from 450 to 310 km. Both these dependences had not been considered in the aforementioned models. Furthermore, we compare the density results from CHAMP with estimates from a spherical calibration satellite, ANDE-Pollux, which allows us to scale the obtained values to quasi-absolute levels. The rest of the paper is organized as follows. In Sect. 2, we first briefly introduce the CHAMP satellite and its accelerometer measurements, then describe our model construction approach and present the CH-Therm-2018 itself. Our model predictions are presented and discussed in Sect. 3. Section 4 presents a validation of our model using satellite laser ranging (SLR) measurements to the spherical satellite ANDE-Pollux. In Sect. 5 we provide the comparison between our model and the NRLMSISE-00 model. The relevant discussion and summary are given in Sect. 6.

Figure 1 The satellite altitude (a) and thermospheric mass density (b) measured by the CHAMP satellite for the whole mission period.

2 Data and model construction

2.1 CHAMP satellite and its accelerometer measurements

The CHAMP spacecraft was launched on 15 July 2000 into a near-circular polar orbit (inclination: 87.3^∘) with an initial altitude of 456 km. By the end of the mission, 19 September 2010, the orbit had decayed to about 250 km. For covering all local times, CHAMP needs 131 days. The thermospheric mass density measurements were deduced from the accelerometer onboard CHAMP, which aimed to measure the nonconservative forces exerted on the satellite with a resolution of < 10⁻⁹ m s⁻² in along-track and cross-track directions (Reigber et al., 2002). The basic equations for deriving the thermospheric mass density from accelerometer measurements have been described by Lühr et al. (2004) and Liu et al. (2005). And by means of an improved approach the mass density is provided with a resolution of less than 10⁻¹⁴ kg m⁻³ (Doornbos et al., 2010). For this study we used the dataset analyzed with the new approach by the Delft group and made available at http://thermosphere.tudelft.nl/acceldrag/data.php (last access: December 2014).

2.2 The approach for constructing an empirical model

To give an overview of the CHAMP mission, Fig. 1a shows the satellite altitude variations for the whole mission period. Its mean value decayed from about 460 km in July 2000 to 260 km in September 2010. We see that the satellite was lifted four times (twice in 2002, once in 2006 and 2009) to higher altitude where the air drag is smaller, for extending the lifetime. The thermospheric mass density derived from the onboard accelerometer is presented in the bottom panel, which shows decreasing density from 2002 to 2009, coinciding with the reducing solar flux. But from August 2009 to the end of mission, the derived mass density has increased dramatically from about $\mathrm{5}\times {\mathrm{10}}^{-\mathrm{12}}$ to $\mathrm{40}\times {\mathrm{10}}^{-\mathrm{12}}$ kg m⁻³, which is mainly caused by the rapid decrease in satellite altitude during the last mission year but also influenced by the rising activity of the solar cycle 24.

Figure 2 The variations of solar flux index (P10.7, a) and solar wind merging electric field (E_m, b) from 2000 to 2010. The mean values of two parameters, P10.7_ref and $E_{{\mathrm{m}}_{\mathrm{ref}}}$, during two 5-year periods (from August 2000 to July 2005 and from August 2004 to July 2009, respectively) are given in the upper part of each panel.

Most important for the variation of thermospheric density is the altitude. In the CH-Therm-2018 model, we consider an exponential dependence on height with a constant scale height for the variation of the mass density. However, as seen in Fig. 1, the CHAMP-measured density has dramatically increased by almost a factor of 8 when its altitude goes below 310 km, which also indicates that a constant scale height is not appropriate for the whole altitude range down to 250 km. Therefore, in this study we consider the 9-year dataset from August 2000 to July 2009 when the satellite was above 310 km, and divide the dataset into two 5-year periods with a 1-year overlap. The two sets of results represent the high-to-moderate and moderate-to-low solar activity conditions, and the altitude of CHAMP decayed from about 460 to 370 km and from 390 to 310 km during the two periods, respectively.

The second most important parameter for the mass density variation is the solar flux level. According to Guo et al. (2008), the solar flux index P10.7 is more suitable than F10.7 for characterizing thermospheric density variations. P10.7 is defined as P10.7 = (F10.7 + F10.7A) ∕ 2, where F10.7A is the 81-day averaged value of the daily F10.7. Figure 2a shows the P10.7 variations from 2000 to the end of 2010, which decreases from over 250 sfu (solar flux unit) in 2002 to below 70 sfu in 2008–2009, and then slightly increases back to 75 sfu at the end 2010. The mean values of P10.7 during the two 5-year periods considered, hereafter referred to as P10.7_ref, are 144.7 and 79.7 sfu, respectively. The bottom panel in Fig. 2 shows the variations of the solar wind merging electric field, E_m. Liu et al. (2010, 2011) and Zhou et al. (2013) found that E_m is an appropriate parameter to describe the disturbance of the thermospheric mass density by magnetic activities. Considering the memory effect of the magnetosphere–ionosphere–thermosphere system to the solar wind input (Werner and Prölss, 1997; Liu et al., 2010), E_m can be defined as follows:

$$\begin{array}{}\text{(1)}& E_{\mathrm{m}}(t,\mathit{\tau })={\displaystyle \frac{{\int }_{t_{\mathrm{1}}}^t{E}_{\mathrm{m}}{}^{\prime }\left(t^{\prime }\right)e^{(t^{\prime }-t)/\mathit{\tau }}\mathrm{d}{t}^{\prime }}{{\int }_{t_{\mathrm{1}}}^t{e}^{(t^{\prime }-t)/\mathit{\tau }}\mathrm{d}{t}^{\prime }}},\end{array}$$

where E_m^′ represents a continuous function of time t^′ of the actual merging electric field at the magnetopause, t₁ is chosen 3 h before the actual epoch (t), and τ (here 0.5 h) is the e-folding time of the weighting function in the integrands. For calculating E_m^′, we use the solar wind to magnetosphere coupling functions, as defined by Newell et al. (2007), and to make E_m^′ values comparable with the solar wind electric field, the function has been rescaled as follows:

$$\begin{array}{}\text{(2)}& E_{\mathrm{m}}{}^{\prime }={\displaystyle \frac{\mathrm{1}}{\mathrm{3000}}}{V}_{\mathrm{SW}}^{\mathrm{4}/\mathrm{3}}(\sqrt{B_y^{\mathrm{2}}+B_z^{\mathrm{2}}}{)}^{\mathrm{2}/\mathrm{3}}{\sin }^{\mathrm{8}/\mathrm{3}}\left({\displaystyle \frac{\mathit{\theta }}{\mathrm{2}}}\right),\end{array}$$

where V_SW (km s⁻¹) is the solar wind velocity and B_y and B_z (nT) are the interplanetary magnetic field (IMF) components in geocentric solar magnetospheric (GSM) coordinates, θ is the clock angle of the IMF ($\tan \left(\mathit{\theta }\right)=\mathrm{|}{B}_y\mathrm{|}/B_z$). With these units the value of the merging electric field will result in millivolts per meter (mV m⁻¹). This approach for calculating the merging electric field has also been used by Xiong and Lühr (2014) and Xiong et al. (2016). From Fig. 2 we see that the values of the merging electric field are below 5 mV m⁻¹ most of the time (slightly higher during higher solar activity years), with mean values hereafter referred to as $E_{{\mathrm{m}}_{\mathrm{ref}}}$ of 1.6 and 1.1 mV m⁻¹ for the two 5-year periods, respectively.

Lei et al. (2012) investigated the annual and semi-annual variations of thermospheric density observed by the CHAMP and GRACE satellites, based on the empirical orthogonal function (EOF) analysis. However, the EOF method does not consider the physical characteristics, and the basic functions of an EOF-derived model can change significantly when a different dataset is used. Therefore, in this study we use the multivariable least-square fitting method for constructing our empirical model. A similar approach has been applied by Marinov et al. (2004) and Liu et al. (2013). In our model, we consider the dependences on height (h), solar flux (P10.7), season (DoY, day of year), magnetic local time (MLT), geographic latitude (θ) and longitude (ϕ), as well as magnetic activity (E_m). We use a set of parameters for fitting the coefficient matrix to the CHAMP measurements, which is expressed as follows:

$$\begin{array}{ll}{\displaystyle }& {\displaystyle }\mathit{\rho }=f_{\mathrm{1}}({\mathit{\rho }}_{\mathrm{0}},h,H_{\mathrm{d}})\cdot f_{\mathrm{2}}\left(P\mathrm{10.7}\right)\cdot f_{\mathrm{3}}\left(\mathrm{DoY}\right)\cdot f_{\mathrm{4}}\left(\mathrm{MLT}\right)\\ \text{(3)}& {\displaystyle }& {\displaystyle }\cdot f_{\mathrm{5}}\left(\mathit{\theta }\right)\cdot f_{\mathrm{6}}\left(\mathit{\phi }\right)\cdot f_{\mathrm{7}}\left(E_{\mathrm{m}}\right),\end{array}$$

where ρ₀ (10⁻¹² kg m⁻³) is the mass density at the reference height (310 km, the lowest height of CHAMP during the 9-year period considered), and H_d denotes the mass density scale height (km). Both parameters are valid for the reference environmental conditions during the two periods (see below). More discussion of these parameters will follow in Sect. 4. The seven subfunctions are defined as follows.

$$\begin{array}{ll}\text{(4)}& {\displaystyle }& {\displaystyle }{f}_{\mathrm{1}}({\mathit{\rho }}_{\mathrm{0}},h,H_{\mathrm{d}})={\mathit{\rho }}_{\mathrm{0}}\cdot e^{\left(-(h-\mathrm{310})/H_{\mathrm{d}}\right)}{\displaystyle }& {\displaystyle }{f}_{\mathrm{2}}\left(\mathrm{P}\mathrm{10.7}\right)=a_{\mathrm{0}}+a_{\mathrm{1}}\cdot \left(\mathrm{P}\mathrm{10.7}-\mathrm{P}{\mathrm{10.7}}_{\mathrm{ref}}\right)\\ \text{(5)}& {\displaystyle }& {\displaystyle }+a_{\mathrm{2}}\cdot {\left(\mathrm{P}\mathrm{10.7}-\mathrm{P}{\mathrm{10.7}}_{\mathrm{ref}}\right)}^{\mathrm{2}}{\displaystyle }& {\displaystyle }{f}_{\mathrm{3}}\left(\mathrm{DoY}\right)=b_{\mathrm{0}}+\sum _{i=\mathrm{1}}^{\mathrm{3}}\mathit{\left\{}{b}_{\mathrm{1}}\right(i)\cdot \cos \left({\displaystyle \frac{i\cdot \mathrm{2}\mathit{\pi }\cdot \mathrm{DoY}}{\mathrm{365.25}}}\right)\\ \text{(6)}& {\displaystyle }& {\displaystyle }+b_{\mathrm{2}}\left(i\right)\cdot \sin \left({\displaystyle \frac{i\cdot \mathrm{2}\mathit{\pi }\cdot \mathrm{DoY}}{\mathrm{365.25}}}\right)\mathit{\}}{\displaystyle }& {\displaystyle }{f}_{\mathrm{4}}\left(\mathrm{MLT}\right)=c_{\mathrm{0}}+\sum _{j=\mathrm{1}}^{\mathrm{4}}\mathit{\left\{}{c}_{\mathrm{1}}\right(j)\cdot \cos \left({\displaystyle \frac{j\cdot \mathrm{2}\mathit{\pi }\cdot \mathrm{MLT}}{\mathrm{24}}}\right)\\ \text{(7)}& {\displaystyle }& {\displaystyle }+c_{\mathrm{2}}\left(j\right)\cdot \sin \left({\displaystyle \frac{j\cdot \mathrm{2}\mathit{\pi }\cdot \mathrm{MLT}}{\mathrm{24}}}\right)\mathit{\}}{\displaystyle }& {\displaystyle }{f}_{\mathrm{5}}\left(\mathit{\theta }\right)=d_{\mathrm{0}}+\sum _{k=\mathrm{1}}^{\mathrm{6}}\mathit{\left\{}{d}_{\mathrm{1}}\right(k)\cdot \cos \left({\displaystyle \frac{k\cdot \mathrm{2}\mathit{\pi }\cdot \mathit{\theta }}{\mathrm{180}}}\right)\\ \text{(8)}& {\displaystyle }& {\displaystyle }+d_{\mathrm{2}}\left(k\right)\cdot \sin \left({\displaystyle \frac{k\cdot \mathrm{2}\mathit{\pi }\cdot \mathit{\theta }}{\mathrm{180}}}\right)\mathit{\}}{\displaystyle }& {\displaystyle }{f}_{\mathrm{6}}\left(\mathit{\phi }\right)=g_{\mathrm{0}}+\sum _{l=\mathrm{1}}^{\mathrm{4}}\mathit{\left\{}{g}_{\mathrm{1}}\right(l)\cdot \cos \left({\displaystyle \frac{l\cdot \mathrm{2}\mathit{\pi }\cdot \mathit{\phi }}{\mathrm{360}}}\right)\\ \text{(9)}& {\displaystyle }& {\displaystyle }+g_{\mathrm{2}}\left(l\right)\cdot \sin \left({\displaystyle \frac{l\cdot \mathrm{2}\mathit{\pi }\cdot \mathit{\phi }}{\mathrm{360}}}\right)\mathit{\}}{\displaystyle }& {\displaystyle }{f}_{\mathrm{7}}\left(E_{\mathrm{m}}\right)=m_{\mathrm{0}}+m_{\mathrm{1}}\cdot \left(E_{\mathrm{m}}-E_{{\mathrm{m}}_{\mathrm{ref}}}\right)\\ \text{(10)}& {\displaystyle }& {\displaystyle }+m_{\mathrm{2}}\cdot {\left(E_{\mathrm{m}}-E_{{\mathrm{m}}_{\mathrm{ref}}}\right)}^{\mathrm{2}}\end{array}$$

The height variation of mass density is described by an exponential function, i.e., Eq. (4), and normalized to the altitude at 310 km. To better use the linear and quadratic fitting, P10.7 and E_m have been centered to their mean values (144.7 and 79.7 sfu and 1.6 and 1.1 mV m⁻¹, respectively) of the two 5-year periods as seen in Eqs. (5) and (10), respectively. The dependences of the other parameters, such as season, magnetic local time, geographic latitude and longitude, have been approximated by trigonometric functions including harmonics from 3 to 6 orders, as shown in Eqs. (6)–(9). In this way 46 parameters are needed to construct the model, and all the bias values in the Eqs. (5) to (10), namely a₀, b₀, c₀, d₀, g₀ and m₀, have been set to 1.

Table 1 The derived values of parameters as defined in Eqs. (4) to (10) for constructing the CH-Therm-2018 empirical model.

3 CH-Therm-2018 model results

As described above, by using each 5-year period of CHAMP measurements we have derived empirical models based on 46 free parameters. The values of these parameters are listed in Table 1. Taking all interrelations into account it results in a number of $\mathrm{3}\times \mathrm{3}\times \mathrm{7}\times \mathrm{8}\times \mathrm{12}\times \mathrm{8}\times \mathrm{3}$ = 145 152 coefficients in our empirical models, both for the high and low solar activity periods. On top we find the reference density at 310 km altitudes. For the first more active period (mean P10.7 = 144.7 sfu and E_m=1.6 mV m⁻¹) we get a value for ρ₀ of $\mathrm{7.65}\times {\mathrm{10}}^{-\mathrm{12}}$ kg m⁻³ and for the second low activity period (mean P10.7 = 79.7 sfu and E_m= 1.1 mV m⁻¹) we get $\mathrm{3.37}\times {\mathrm{10}}^{-\mathrm{12}}$ kg m⁻³. This decrease by a factor of 2.2 reflects primarily the effect of the change in solar flux level. Next in line in Table 1 is the scale height. The derived values of 94 and 80 km for the two activity periods are quite large. For comparison, Liu et al. (2011) estimated from comparisons of CHAMP and GRACE density measurements scale heights of 83 and 60 km for solar flux levels of 200 and 80 sfu, respectively. A more detailed discussion of our constant scale height is given in Sect. 6.

Figure 3 The altitude versus solar activity variations of model-predicted thermospheric mass density around noon at (a) high and (b) low solar activity conditions. The longitude has been chosen at Greenwich meridian. Panels (c) and (d) are the altitude versus geographic latitude variations of model-predicted mass density for high and low solar activity conditions, respectively. Panels (e) and (f) show the dependence of model-predicted mass density on merging electric field for both periods.

The obtained dependence of mass density on solar flux level is twice as high during the low solar flux period as during the solar maximum years. This result has to be seen in connection with the obtained scale height. The harmonically varying dependences on season, local time, latitude and longitude show no pronounced dependence on the activity level when combining the two amplitudes (cosine and sine) of the fundamental oscillations. Different from that, the relative dependence of the mass density on magnetic activity (parameter at bottom) is significantly higher for low solar activity. In the following we are going to present the main features captured by the two different model solutions.

Figure 4 Similar as Fig. 3, but for the distribution of (a) and (b): geographic latitude versus magnetic local time; (c) and (d): geographic latitude versus day of year; (e) and (f): geographic latitude versus longitude.

The panels (a) and (b) of Fig. 3 show the altitude versus solar activity variations from the two periods, over an altitude range from 310 to 470 km. As the level of solar activity is quite different for the two periods, the range of P10.7 has been limited to 100–280 and 65–125 sfu, respectively. The model-predicted mass density shows generally similar variations for both periods, which increases with larger solar activity but decreases with altitude. The borders between different colors can be interpreted as constant pressure levels. Panels (c) and (d) of Fig. 3 show the altitude versus geographic latitude variation of the mass density around noon hours. The P10.7 values for the two periods have been set to 150 and 80 sfu, respectively. The mass density generally decreases from low to high latitudes during both periods. For the higher solar activity conditions, the equatorial mass density anomaly, which was earlier described by Liu et al. (2005, 2007), can be seen, with the peak mass density appearing around ±20^∘ latitude. The panels (e) and (f) of Fig. 3 show the dependence of the model-predicted mass density on merging electric field during both periods. We see that the mass density increases roughly linearly with the merging electric field with hardly any indication of a saturation effect, which is consistent with results published by Müller et al. (2009) and Liu et al. (2011).

In Fig. 4 the dependence on periodically varying parameters is shown. The panels (a) and (b) present the MLT versus latitude distribution of the mass density. The solar activity has been set again to 150 and 80 sfu for the two periods and the altitude has been set to 400 and 340 km, respectively. During both solar activity periods, the mass density reaches its maximum and minimum around 14:00 and 03:00 MLT, respectively. The EMA feature is more evident at higher solar activity conditions, as shown in panel (a) of Fig. 4, with larger crest density in the Northern Hemisphere, as we have chosen predictions for the September equinox. Additionally, a clear density trough is seen around −75^∘ in the Southern Hemisphere during the lower solar activity conditions. The panels (c) and (d) of Fig. 4 present the seasonal versus latitude variations, showing the mass density peaks at the two equinox seasons and a pronounced minimum around the June solstice, which is a well-known feature (e.g., Emmert, 2015). An interesting detail is that the mass density exhibits larger amplitudes during the March equinox than during the September equinox for high solar activity conditions, while it exhibits an opposite ratio for lower solar activity conditions. This equinox asymmetry of thermospheric mass density is consistent with the findings of Liu et al. (2013), who reported that the equinox asymmetry weakens or disappears when the solar flux level falls to below P10.7 = 110 sfu. Guo et al. (2008) argued that the March–September equinox asymmetry can partly be attributed to the interannual variability of the thermosphere mass density. Another interesting feature seen from the model-predicted result is that at all latitudes the thermospheric mass densities are lower during the June solstice than those during the December solstice, while the expected hemispheric asymmetry between high-latitude densities during solstice seasons is not evident in our model outputs. We checked the mean annual variations of CHAMP density measurements at various latitude bands and confirm the dominance of the July minimum at all latitudes with a deeper trough in the Southern Hemisphere.

The coupling between the lower atmosphere and upper atmosphere–ionosphere has been widely investigated in relation to longitudinal wave patterns of different thermospheric and ionospheric parameters (e.g., Immel et al., 2006; Häusler et al., 2007; Liu et al., 2009). The tides excited by latent heat release in tropospheric deep convection tropical clouds can propagate vertically upward (Hagan and Forbes, 2003). These tides vary with season, causing longitudinal patterns with varying wave numbers over the course of a year. Best known are the wave number 4 (WN4) pattern during the months around August and wave number 3 (WN3) pattern around solstice seasons, corresponding to the diurnal eastward-propagating DE3 and DE2 tidal components, respectively (e.g., Forbes et al., 2006; Lühr et al., 2008; Wan et al., 2010; Xiong and Lühr., 2013). The panels (e) and (f) of Fig. 4 show the global distribution of the mass density around noontime for the two considered conditions. Here we find again the EMA signature. Some tidal features, a mixture of longitudinal wave-3 and wave-4 patterns, are found at EMA crest regions in particular during the higher solar activity period, while for the lower solar activity, wave-2 and wave-3 patterns are more prominent. The difference in longitudinal wave patterns may be due to their different wavelengths and their relative susceptibility to molecular dissipation at different solar flux conditions (Bruinsma and Forbes, 2010).

Figure 5 Scaling factors of thermospheric density derived from the analysis of SLR data from the ANDE-Pollux during 16 August to 30 September 2009 for two models: JB2008 and CH-Therm-2018.

4 Density validations by SLR measurements to calibration satellites

So far we have presented density results derived entirely from the CHAMP air drag measurements. Atmospheric drag is the major nongravitational force acting on LEO satellites, and it causes orbital decay. Since the atmospheric drag depends primarily on the mass density, SLR measurements of spherical LEO satellites can be used to estimate mass density at their altitude. Because of their simple geometry, so-called cannonball satellites can be used for quasi-absolute calibrations. This is not an easy task since, on the one hand, it requires precise modeling of all other gravitational and nongravitational perturbations acting on the satellites, and on the other hand, the amount of SLR observations contributing globally to LEO satellites observations is low. However, the derived density values can either be used to validate empirical models locally or provide scaling factors for these models (Panzetta et al., 2018).

As an example, we analyzed the SLR observations to the cannonball LEO satellite ANDE-Pollux between 16 August and 3 October 2009, and we derived from 6–12 h time series of estimated scaling factors for the thermospheric density predictions for the CH-Therm-2018 models. Figure 5 shows the comparison between SLR results and CHAMP estimates in terms of scaling factors. The mean and median values of the derived scaling factors are 1.4 and 1.267, respectively. Also included in the figure is the comparison with the JB2008 model. These values infer that the CH-Therm-2018 model underestimates the thermospheric density at least during the time interval used. In fact, the underestimation of CHAMP density estimates has earlier been suggested by Doornbos (2012), who reported that the CHAMP-derived densities were systematically lower by about 25 % than those from GRACE when normalized to a common altitude with the help of an atmospheric model like NRLMSISE-00. Some uncertainty may be introduced by the fact that the ANDE-Pollux observations we compared here are taken from August and September 2009, while the CHAMP dataset we used for the CH-Therm-2018 model ends in July 2009. By taking advantage of the obtained median factor we scaled up all the CH-Therm-2018 predicted mass density values by 1.267.

In addition we also compared the SLR-derived densities with four different empirical models: CIRA86 (Hedin et al., 1988), NRLMSISE-00 (Picone et al., 2002), DTM2013 (Bruinsma, 2015) and JB2008 (Bowman et al., 2008b). The corresponding mean values of the estimated scaling factors are 0.65±0.26 for CIRA86, 0.65±0.25 for NRLMSISE-00, 0.79±0.24 for DTM2013 and 0.89±0.27 for JB2008, respectively. It indicates that all these models clearly overestimate the thermospheric density during the period of the low recent solar minimum.

Figure 6 Panel (a) shows the CH-Therm-2018 model-predicted mass density (red) and CHAMP observations (black) from August 2000 to July 2009. Panel (b) shows the same density but for the NRLMSISE-00 model (green) and CHAMP observations (black). Panel (c) gives the annual average relative differences between the model estimates and CHAMP observations.

Figure 7 The linear regression between CHAMP observations and our model-predicted results during 131-day period (a) for high (centered on 1 March 2002) and (b) low (centered on 1 November 2008) solar activity conditions, respectively. (c) The red color shows the slope (top panel) and mean ratio (bottom panel) of the linear regression for each 2-month period from 2000 to 2010. The green color shows the results from NRLMSISE-00 model.

5 Comparison with the NRLMSISE-00 model

For reproducing the CHAMP observations with our empirical model, we have combined the results derived from both periods. For the results from August 2000 to July 2004 we use the model predictions from the first 5-year period, while for the results from August 2005 to July 2009 we use the model predictions from the second 5-year period. For the 1-year overlapping period from August 2004 to July 2005, we consider the model predictions from both periods, but use a linearly weighted combination for the time of overlap. Figure 6a presents our model predictions (red) and CHAMP observations (black) from August 2000 to July 2009. In general, our model follows the measurements quite well, and even the spikes (corresponding to high magnetic activity) are reasonably well reproduced. For comparison, the panel (b) also shows the predictions from the NRLMSISE-00 model (green), which has been divided by the scale factor of 1.267 as derived from Fig. 5. Compared to our model, the NRLMSISE-00 model is clearly overestimating the CHAMP observation during solar minimum years. The bottom panel presents quantitatively the relative differences between the model predictions and CHAMP observations:

$$\begin{array}{}\text{(11)}& \mathrm{\Delta }\mathit{\rho }={\displaystyle \frac{{\mathit{\rho }}_{\mathrm{model}}-{\mathit{\rho }}_{\mathrm{CHAMP}}}{{\mathit{\rho }}_{\mathrm{CHAMP}}}}\cdot \mathrm{100}.\end{array}$$

The annual average differences between our model and observations are within the range ±20 % for all nine years, while NRLMSISE-00 overestimates the observations by about 5 % for high and moderate solar activity years, and reaches as high as 40 % for the extremely low solar activity years. It is no surprise that our model predicts the observations better than the NRLMSISE-00 model, because our model is derived from CHAMP data, which have not been included in the NRLMSISE-00 model.

For a more quantitative inspection of the CHAMP model, we have divided the 9-year dataset into 2-month bins of overlapping 131-day intervals. This time period is required for covering all 24 h of local time in each bin. For the 2-month bins, we calculate the linear regression slope and the mean ratio between the CHAMP observations and model predictions. The mean ratio is defined as the ratio between the mean values of the observations divided by the model predictions during the 131 days. Examples of this analysis during high (centered on 1 March 2002) and low (centered on 1 November 2008) solar activities are presented in Fig. 7a and b, respectively. The correlation coefficients between the model predictions and observations reach 0.89 and 0.86, the slopes of the linear fitting are 1.03 and 1.07, and the mean ratios are 1.11 and 1.04. Panel (c) of Fig. 7 presents the slope (a, b) and mean ratio (bottom panel) between the observations and our empirical model (red) as well as the NRLMSISE-00 model (green), respectively. Here again the NRLMSISE-00 model (green) has been downscaled by a factor of 1.267.

The slope CH-Therm-2018 model result varies within the range of 0.6 to 1.2 and the mean ratio varies between 0.9 and 1.2 during almost all the nine years, which are better than those of the NRLMSISE-00 model during the solar minimum (2008–2009). An exception makes the excursion of the slope around 0.6 at the end of 2003. This means both our model and NRLMSISE-00 overestimate the mass density during October and November 2003 (see Fig. 6) periods of very strong magnetic storms.

It is worth noting that we have extended the model prediction to the last year of the CHAMP mission, as shown in Fig. 7c. We see that the slope and the mean ratio between observations and our empirical model have increased dramatically, reaching values of more than 4.0 and 2.0 at the end of the mission, respectively. This is a consequence of the quite low altitude of the CHAMP satellite. Therefore, we have to note that our model is suitable for the altitude range from 310 to 470 km. And the large increase of the CHAMP-measured mass density during the last mission year (see Fig. 1) might be an indication of a smaller scale height due to a composition change at altitudes below 310 km.

6 Discussions and summary

We have constructed a new model of thermospheric neutral density, called CH-Therm-2018, from the CHAMP accelerometer measurements over a 9-year period from August 2000 to July 2009, covering both high and low solar activity conditions (solar flux index P10.7 ranges from over 250 to below 70 sfu). The CHAMP altitude changed from 460 down to 310 km within this period. Good fits between model and observation are achieved when a constant scale height over this range is assumed. But in addition solar flux level and magnetic activity dependent scaling factors are introduced. This is from the physical point of view not justified because neither the solar flux nor the magnetic activity increases the amount of air particles. Both these parameters change the height distribution of neutral particles and thus modify the scale height. During the CHAMP mission the orbital altitude decreased simultaneously with the reduction of solar flux level. For that reason it is impossible to reliably determine the dependence of the scale height on solar flux from this dataset. For this modeling purpose this deficiency can be mitigated by a piecewise approximation of the real scale height relation by an exponential function with fixed scale height, and a reference density at 310 km altitude scaled by a solar flux and magnetic activity functions. The two considered periods are 5 years long.

Conventional atmospheric models often have problems with representing the magnetic activity dependence. From Table 1 (bottom rows) it is obvious that the relative dependence on magnetic activity increases significantly when the solar activity goes down. This fact has been noted frequently before. But it is also worth mentioning that the absolute change in mass density with magnetic activity is fairly independent of the solar flux background (see Fig. 3e and f). This confirms earlier claims by Müller et al. (2009) and Liu et al. (2011).

An independent validation of the model-predicted mass densities was performed by comparing with SLR observations on the spherical satellite ANDE-Pollux. Because of the simple geometry of this spacecraft, obtained density estimates can be considered as quasi-absolute. Comparisons performed during the period of low solar activity (16 August to 30 September 2009) reveal that the density values of the CH-Therm-2018 model should be upscaled by a factor of 1.267 to fit the SLR observations. This factor has been applied to all model values.

The comparison between our adjusted model predictions with the NRLMSISE-00 model shows that the thermospheric density predicted by the CH-Therm-2018 model agrees well (within ±20 %) with the CHAMP observations over the whole period, while the NRLMSISE-00 model overestimates the observations by about 40 % at the period of low solar activities.

The CH-Therm-2018 model shows quite different features of thermospheric mass density at different solar activity conditions. For example, the EMA feature is more prominent at higher solar activity. The larger density at the March equinox than that at the September equinox is only seen at higher solar activity, while this seasonal asymmetry exhibits an opposite sense during lower solar activity conditions. Concerning the tidal signatures at low and equatorial latitudes the thermospheric mass density presents mainly longitudinal wave-4 and wave-3 patterns at higher solar activity, changing to wave-3 and wave-2 patterns at lower solar activity period.

A pending issue for the future studies is a better representation of the mass density height dependence. For this, it would be helpful to take simultaneous measurements from at least two satellites into account. Also the extension of the model to lower altitudes, down to the GOCE orbit is planned for a follow-up study.