Effect of the solar activity variation on the Global Ionosphere Thermosphere Model (GITM)

Abstract

The accuracy of global atmospheric models used to predict the middle/lower thermosphere characteristics is still an open topic. Uncertainties in the prediction of the gas properties in the thermosphere lead to inaccurate computations of the drag force on space objects (i.e. satellites or debris). Currently the lifetime of space objects and therefore the population of debris in low Earth orbit (LEO) cannot be quantified with a satisfactory degree of accuracy. In this paper, the Global Ionosphere Thermosphere Model (GITM) developed at the University of Michigan has been validated in order to provide detailed simulations of the thermosphere. First, a sensitivity analysis has been performed to investigate the effect of the boundary conditions on the final simulations results. Then, results of simulations have been compared with flight measurements from the CHallenging Minisatellite Payload (CHAMP) and Gravity Recovery and Climate Experiment (GRACE) satellites and with existing semi-empirical atmospheric models (IRI and MSIS). The comparison shows a linear dependency of the neutral density values with respect to the solar activity. In particular, GITM shows an over-predicting or under-predicting behaviour under high or low solar activity respectively. The reasons for such behaviour can be attributed to a wrong implementation of the chemical processes or the gas transport properties in the model.

Ann. Geophys., 34, 725736, 2016 www.ann-geophys.net/34/725/2016/ doi:10.5194/angeo-34-725-2016 Author(s) 2016. CC Attribution 3.0 License.

Davide Masutti1, Gnther March2, Aaron J. Ridley3, and Jan Thoemel1

1Von Karman Institute for Fluid Dynamics, Sint-Genesius-Rode 1640, Belgium

2Delft University of Technology, Delft 2629 HS, the Netherlands

3University of Michigan, Ann Arbor, MI 48109, USA

Correspondence to: Davide Masutti ([email protected])

Received: 14 January 2016 Revised: 4 July 2016 Accepted: 5 July 2016 Published: 1 September 2016

Abstract. The accuracy of global atmospheric models used to predict the middle/lower thermosphere characteristics is still an open topic. Uncertainties in the prediction of the gas properties in the thermosphere lead to inaccurate computations of the drag force on space objects (i.e. satellites or debris). Currently the lifetime of space objects and therefore the population of debris in low Earth orbit (LEO) cannot be quantied with a satisfactory degree of accuracy. In this paper, the Global Ionosphere Thermosphere Model (GITM) developed at the University of Michigan has been validated in order to provide detailed simulations of the thermosphere. First, a sensitivity analysis has been performed to investigate the effect of the boundary conditions on the nal simulations results. Then, results of simulations have been compared with ight measurements from the CHallenging Minisatellite Payload (CHAMP) and Gravity Recovery and Climate Experiment (GRACE) satellites and with existing semi-empirical atmospheric models (IRI and MSIS). The comparison shows a linear dependency of the neutral density values with respect to the solar activity. In particular, GITM shows an over-predicting or under-predicting behaviour under high or low solar activity respectively. The reasons for such behaviour can be attributed to a wrong implementation of the chemical processes or the gas transport properties in the model.

Keywords. Ionosphere (modelling and forecasting) magnetospheric physics (solar windmagnetosphere interactions) meteorology and atmospheric dynamics (thermospheric dynamics)

Effect of the solar activity variation on the Global Ionosphere Thermosphere Model (GITM)

1 Introduction

The thermosphere is the layer of the Earths atmosphere directly above the mesosphere and directly below the exo-sphere (90600 km altitude). The thermosphere can be described in terms of density, chemical composition, temperature and wind velocity. The thermosphere, largely driven by the solar and geomagnetic activity, shows a strong coupling with the ionosphere and the magnetosphere (Lu et al., 1995). This atmospheric layer is characterised by a high degree of variability. Thermosphere density varies over a wide range of spatial and temporal scales under the inuence of the complex interactions between the Earth system and the solar processes. Unfortunately, the measurements of the thermosphere have always been quite sparse and their uncertainties have often introduced a relatively high level of ambiguity. The accurate description of the thermosphere can be a key aspect for the design of a space mission in low Earth orbit (LEO). In particular, the aerodynamic properties impact the lifetime and the attitude control of LEO satellites. In order to accurately predict the satellite orbital decay, for example, a more detailed modelling of the atmosphere density is necessary. Indeed, the neutral air density is affected signicantly by the incoming solar radiation and its value can even increase by 2 orders of magnitude during geomagnetic storms as shown by Sutton et al. (2005) and Krauss et al. (2015). During a solar storm event in March 1991 as reported by Sagalyn and Bowhill (1993), a high number of satellites were lost from the NORAD catalogue because their trajectory changed substantially as a result of the increased neutral density and consequent drag increase. This prevented the necessary conjunction analysis between LEO objects, as well as possible eva-

Published by Copernicus Publications on behalf of the European Geosciences Union.

726 D. Masutti et al.: Effect of the solar activity variation on GITM

sive manoeuvres leading to a risk to space assets and astronauts.

Currently, the atmospheric forecasts are based on semi-empirical models for the gas neutral composition; among those models there are MSIS-90 (Hedin, 1991), NRLMSISE-00 (Picone et al., 2002), DTM (Bruinsma, 2015) and JB2008 (Bowman et al., 2008). These models are nely tuned to match a database of ight measurements, but when applied outside the interpolated ranges they are subjected to large uncertainties in atmosphere density and composition. These uncertainties affect the orbital propagator calculation by leading to an error of several kilometres in the estimated position of LEO satellites (Vallado and Finkleman, 2014). This suggests that physics-based modelling should be used to forecast the state of the upper atmosphere. Such modelling inherently needs a thorough validation.

Among the different physical models available, the general circulation models (GCMs) cover an important eld of atmospheric modelling. GCMs like CMAT2 (Harris, 2001), TIEGCM (Richmond et al., 1992), CTIPe (Millward et al., 1996) or GITM (Ridley et al., 2006) require an accurate implementation of the thermosphere phenomena in order to provide validated results. The last of these models, the Global Ionosphere Thermosphere Model (GITM), has been used for this investigation. Neutral density and electron density data available from CHAMP (CHallenging Minisatellite Payload) and GRACE (Gravity Recovery and Climate Experiment) have been compared with the simulations in order to investigate the effect of the solar activity on GITM.

Because of the increasing number of objects in low Earth orbit, the need for an accurate real-time forecasting tool has become critical. In 2016, the QB50 mission will launch one of the biggest CubeSat constellations in order to perform multi-point measurements of the predominant species down to 200 km altitude (Thoemel et al., 2014). These scientic measurements will enhance the understanding of atmospheric characteristics in the middle and lower thermosphere and they will provide data for a more accurate validation of GCMs.

In Sect. 2, it is presented a brief description of GITM and the ight data used for the validation. The results obtained by comparing the computations and the ight data over a complete range of solar ux values are discussed in Sect. 3. Section 4 is devoted to the conclusions.

2 Methodology

2.1 GITM

GITM is a three-dimensional NavierStokes equation solver that models the Earths thermosphere and ionosphere system using a spherical coordinate grid. The GITM code solves the coupled continuity, momentum and energy equations of neutrals and ions in a time-marching fashion. For the ions, the

Longitude [deg]

10-10

3

2.9

2.8

-3 2.7

Density [kg m]

2.6

2.5

2.4

2.3

2.2

2.1

-180 -150 -120 -90 -60 -30 0 30 60 90 120 150 180

Figure 1. Grid sensitivity analysis at 0 latitude and 200 km altitude (long lat cells).

velocity variation as a function of time is ignored, so the steady-state ion ow velocity is solved. Moreover, for this model the ion temperature consists of a mixture of the electron and neutral temperatures, and the neutrals are solved for using the NavierStokes equations. GITM explicitly solves for the neutral densities of O, O2, N(2D), N(2P), N(4S), N2 and NO, as well as ion species O+(4S), O+(2D), O+(2P), O+2, N+, N+2 and NO+. The simulations are initiated using the neutral/ion densities and temperatures from the MSIS-90 and IRI (Bilitza, 2001) semi-empirical models. In general, the results of the simulations become independent of the initial conditions after 72 h of simulated time. In contrast to other general circulation models of the thermosphere, GITM solves the vertical and horizontal advection in separate equations without assuming any hydrostatic equilibrium. This allows the code to take into account the non-constant gravity in the vertical direction and the Coriolis forces (Deng et al., 2008).

Because GITM 3-D computations are based on a structured spherical grid with adjustable spatial discretisation in longitude and latitude, the accuracy of the results needs to be veried as a function of the grid renement. A reference case simulation on 19 April 2004 has been run for 4 days (3 days of spin-up + 1 day of simulation). From the grid analysis

results shown in Figs. 13, it is observed how the grid renement increases the accuracy of the simulations. This is illustrated, for example, in Fig. 1, where the neutral density values of the atmosphere at an altitude of 200 km at 0 latitude are displayed as a function of the longitudinal direction.Even if the density values are conned in a very narrow range between 21010 and 31010 kg m3, the different longi

tudinal proles are scattered as a function of the grid renement. By increasing the number of cells in the mesh, the density proles tends to converge to similar results. In fact, the relative error between the density proles with a renement of 144 72 cells (2.5 long 2.5 lat) and 104 7252 cells

Ann. Geophys., 34, 725736, 2016 www.ann-geophys.net/34/725/2016/

D. Masutti et al.: Effect of the solar activity variation on GITM 727

Longitude [deg]

10-12

10-12

3.5

1.6

36x18 cells 52x26 cells 72x36 cells 104x52 cells

3

1.4

-3 2.5

Density [kg m]

Density [kg m]

-3

1.2

2

1

1.5

1

0.8

144x72 cells

0.5

0.6

-180 -150 -120 -90 -60 -30 0 30 60 90 120 150 180

-90 -60 -30 0 30 60 90

Figure 2. Grid sensitivity analysis at 0 latitude and 400 km altitude (long lat cells).

(3.46 long 3.46 lat) is less than 0.5 %, whereas the rela

tive error between the coarser (3618 cells or 10 long 10

lat) and the ner (144 72 cells or 2.5 long 2.5 lat) grid

is 7 %. In both cases, the density proles seem to have peaks and valleys at similar longitudinal locations. With the same grid analysis repeated at an altitude of 400 km (Fig. 2), the relative error between the coarser and ner grid is increased to almost 19 %. However, the relative error between the two most rened grids is as low as 1.2 %. The analysis suggests that the quality of the results obtained from GITM simulations are in general grid-dependent below a certain renement threshold. Because GITM uses a spherical grid for its computations, the cell dimensions change as a function of the altitude for a constant grid renement. This implies that, for a grid with 14472 cells (5 long 5 lat), the side length of

a squared cell at 200 km is about 286 km, whereas this value linearly increases to about 295 km length for an altitude of 400 km.

Figure 3 again shows a grid sensitivity analysis of the density proles over a meridian at 0 longitude and 400 km altitude. It can be seen from also looking at the meridian proles that the relative difference between the two ner grids (144 72 and 104 7252 cells) is less than 1.5 %. Because

of the low relative difference observed between the two ner grids in Figs. 1, 2 and 3, it can be assumed that the accuracy of the results will not benet from a further renement of the mesh in both longitude and latitude. In order to nd a tradeoff between the accuracy of the simulations and the computation time, it was decided to select a grid with 7236 cells in

longitude and latitude respectively. This grid, used hereafter, allows for results to be obtained with only 3 % difference from the most rened grid in 400 CPU hours instead of 1500 CPU hours (for 4 simulated days). The GITM results, discussed in Sect. 3, have been obtained from 3-D runs with a total of 15 simulated days each (3 days of spin-up + 12 days

Latitude [deg]

Figure 3. Grid sensitivity analysis at 0 longitude and 400 km altitude (long lat cells).

of simulation) in order to be fully representative of the model

performances.

2.2 Validation data

CHAMP is a German small-satellite mission for geoscientic and atmospheric research and applications, managed by the GeoForschungsZentrum (GFZ). Between 2000 and 2010, CHAMP simultaneously generated, for the rst time, high-precision gravity and magnetic eld measurements. These measurements allowed for the detection of not only the spatial variations in both elds but also their variability with time (Reigber and Schwintzer, 1995). The accelerations on CHAMP are measured by the STAR (Space Three-axis Accelerometer for Research) instrument developed by ONERA (Touboul et al., 1998). Its objective is to measure all non-gravitational accelerations of the satellite (drag, solar and Earth radiation pressure) in order to determine the Earths gravity eld from purely gravitational orbit perturbations.The accelerometer measurement principle is based on electrostatic suspension of a proof mass in a cage. The instantaneous position of the proof mass is measured by three capacitive sensors, which permit a determination of the acceleration vector and resulting drag and density. The instrument has a dynamic range of 1104 m s2, a resolution of better than 3109 m s2, and a frequency range of 101 to 104 Hz.

The STAR instrument is positioned at the centre of gravity of CHAMP to minimise the effect of rotational accelerations and gravity gradients on the measurements.

Data derived from STAR can be used to calculate the neutral density with an uncertainty level of 1015 % (1) for moderate geomagnetic activity (Bruinsma et al., 2004). Sutton (2009) shows an overestimation of the CHAMP density values of about 20 % when compared to the HASDM semi-empirical model. Doornbos (2012) observed a similar overestimation when comparing the same semi-empirical model

www.ann-geophys.net/34/725/2016/ Ann. Geophys., 34, 725736, 2016

728 D. Masutti et al.: Effect of the solar activity variation on GITM

Figure 4. Daily (grey line) and 81-day-averaged (black line) solar ux received on Earth at 10.7 cm wavelength between 1980 and 2015. The grey region represents the selected investigation period.

data with CHAMP densities processed with a different retrieval algorithm. Both Sutton (2009) and Doornbos (2012) conclude that the bias in the density values is caused by the approximations in the modelling of the external satellite geometry and the aerodynamic interaction of the gas with the spacecraft walls. The design of the CHAMP satellite, with its complex elongated shape and protruding instruments, antennae and bafes, is not ideal for the purpose of density measurements.

CHAMP has also a planar Langmuir probe (PLP) in order to measure the atmosphere electron density with an uncertainty within 10 % (1) (Liu et al., 2007). The PLP provides measurements of the electron density and temperature every 15 s for 1 s in the voltage sweep mode. Thus, CHAMP data from the Sutton (2009) database v2.2 are available every 15 s, which corresponds to a spatial resolution of about 115 km considering a satellite velocity of 7.6 km s1. Further details of the PLP design and performance are summarised by Cooke et al. (2003) and McNamara et al. (2007). In the present paper, novel comparisons have been performed between numerical models and electron density data from CHAMP.

On top of CHAMP data, GITM simulations have also been compared with GRACE ight data. GRACE was launched in 2002 and consists of two almost identical satellites which follow each other on the same orbit at a distance of about 220 km (Kang et al., 2006). The orbits are circular polar, with an initial altitude of 485 km at launch (near solar maximum23) and an inclination of 89 . The design of the GRACE satellites is similar to the CHAMP satellite. Indeed, these twin satellites contain a GPS receiver, an extremely precise accelerometer and star tracker cameras. In addition, each of the GRACE satellites carries a K-band horn directed at the other satellite for inter-communication. Like CHAMP, these satellites have the main objective of investigating Earths gravity eld and its small variations.

Figure 5. Effect of the geomagnetic storminess on GITM and CHAMP neutral densities from 10 to 30 July 2004.

Table 1. Boundary conditions of the simulations.

Date (DD/MM/YYYY) F10.7 Ap index

(W m2 Hz1)

Case a 18/12/2001 233.1 24.5 9 7

Case b 15/06/2002 145.6 5.2 7 4

Case c 28/08/2003 112.9 8.2 12 8

Case d 10/04/2004 100.3 9.9 8 5

Case e 10/10/2005 77.4 1.8 5 4

Case f 24/10/2008 66.7 0.8 5 4

3 Results

3.1 Sensitivity to solar and geomagnetic variability

The solar ux F10.7 is a good proxy to describe the daily solar activity; Fig. 4 shows the cyclic behaviour of this index between the years 1980 and 2015. The grey line in Fig. 4 shows the daily variation in the F10.7 proxy, whereas the black line represents the 81-day-averaged variation. The grey area indicates the corresponding solar activity period during the descending phase of the 23rd solar cycle, where the F10.7 values are ranging from 66.7 to 233.1 W m2 Hz1. Table 1 summarises the boundary conditions of the cases under investigation in this paper. These specic simulation windows have been selected in periods of low geomagnetic activity (Ap < 30) in order to avoid any geomagnetic storm that could interfere with the study. Figure 5 shows the impact of a solar storm on CHAMP ight data and GITM simulations. Three distinct storms (Ap >= 150) impacted the thermosphere be

tween 22 and 27 July 2004. Both the data and simulations show an increase of the neutral density by almost 1 order of magnitude.

The neutral density of the GITM simulations from Table 1 has been averaged over a period of 12 days and presented as a function of altitude, latitude and local time in order to extract

Ann. Geophys., 34, 725736, 2016 www.ann-geophys.net/34/725/2016/

D. Masutti et al.: Effect of the solar activity variation on GITM 729

Figure 6. Logarithmic variation in the atmosphere density at 0 and 80 N latitude and 02:00 LT as a function of altitude and adjusted solar ux (F10.7).

Figure 7. Logarithmic variation in the atmosphere density at 0 and 80 N latitude and 14:00 LT as a function of altitude and adjusted solar ux (F10.7).

the long term effect of the solar ux input. Figures 6 and 7 show the effect of the F10.7 values on the simulated density in the thermosphere at 0 and 80 N latitude as a function of altitude at 02:00 and 14:00 local time (LT) respectively. As can be observed from Fig. 6, the sensitivity of the atmosphere density to the solar ux is increasing with the increasing altitude. At about 150 km altitude and below there is a small effect of the Suns activity, and the air density can be considered almost independent of the energy radiated by the Sun. However, at an altitude of 400 km, the air density value can increase by 2 orders of magnitudes between low and high solar activity (i.e. 66.7 to 233.1 W m2 Hz1 respectively).

The comparison between Figs. 6 and 7 shows that the increment of density between low and high solar activity is dependent on both local time and latitude. In particular, the density increase is slightly bigger at 80 N with respect to 0 (at constant local time) and on the nightside with respect to the dayside (at constant latitude). Because of this large variation

Figure 8. Comparison of neutral density as a function of time from GITM simulations (solid black line), CHAMP ight data (solid grey line) and semi-empirical models (dashed line).

Figure 9. Comparison of electron number density as a function of time from GITM simulations (solid black line), CHAMP ight data (solid grey line) and a semi-empirical model (dashed line).

in density in the mid-lower thermosphere, satellites in low Earth orbit must increase their orbital altitude to take into account the increased aerodynamic drag during periods of high solar activity. An example is the International Space Station (Zhukov et al., 2012).

3.2 Comparison with ight data

The results of the GITM simulations are compared with neutral densities and electron densities measured by CHAMP and GRACE. As an example, the trend of the simulations and the ight data is presented in Figs. 8 and 9 as a function of time (only 6 h) for the neutral and the electron density respectively. As can be observed in Fig. 8, the density shows a sinusoidal-like behaviour typical of the daynight sequence.The higher peaks correspond to a high density over the Equator on the dayside, whereas the valleys correspond to the pas-

www.ann-geophys.net/34/725/2016/ Ann. Geophys., 34, 725736, 2016

730 D. Masutti et al.: Effect of the solar activity variation on GITM

Figure 10. Density ratio (GITM / CHAMP) as a function of time and latitude for case a (local time of ascending node, LTAN = 15:10; local time of the descending node, LTDN = 03:10;

hapogee = 469 km; hperigee = 405 km).

Figure 11. Density ratio (GITM / CHAMP) as a function of time and latitude for case b (LTAN = 22:50; LTDN = 10:50;

hapogee = 449 km; hperigee = 400 km).

sage over the South Pole on the nightside. The agreement between the ight data and the semi-empirical models is very good, especially with the NRLMSISE-00 model. However, a marked overestimation of the order of 50 % can be observed from the comparison between the GITM simulations and the ight data.

In Fig. 9 the electron densities obtained from CHAMP, the GITM simulations and the semi-empirical model are compared. The agreement with the ight data and the IRI-2012 model is good, and again less good agreement with the GITM simulations is observed. In particular, on the nightside, the simulations around the Equator show a highly over-predicted electron density with respect to ight data. The reason for such an overestimation can be attributed to the missing divergence term in the electron continuity equation in GITM. Such a missing term together with the missing electrodynamic model does not allow GITM to account for the proper electron downwelling which causes accumulations of electrons on the nightside in the equatorial region.

Figure 12. Density ratio (GITM / CHAMP) as a function of time and latitude for case c (LTAN = 06:44; LTDN = 18:44;

hapogee = 423 km; hperigee = 392 km).

Figure 13. Density ratio (GITM / CHAMP) as a function of time and latitude for case d (LTAN = 10:01; LTDN = 22:01;

hapogee = 410 km; hperigee = 378 km).

Figures 1015 show a comprehensive comparison between GITM simulations and CHAMP neutral densities as a function of latitude and time for all cases presented in Table 1. All gures represent the ratio between GITM-simulated density and CHAMP ight data. A density ratio value of 1 corresponds to a good validation between the simulations and the experiments, whereas values larger or smaller than 1 indicate that the GITM density is over-predicted or under-predicted with respect to CHAMP data. The ratio values are shown as a function of time and latitude for the ascending/descending branch of the orbit. Local times of ascending/descending nodes and orbit altitude are given for each case under investigation. In general, Figs. 1015 show that GITM simulations are consistent in time and they are representative of the behaviour of the model under the different solar activity periods.

For all cases (af), GITM simulations present a different prediction behaviour over the poles than over the Equator.Over the poles the computed neutral density tends always to under-predict the ight data with respect to the predictions at the Equator. This is true both for the ascending/descending

Ann. Geophys., 34, 725736, 2016 www.ann-geophys.net/34/725/2016/

D. Masutti et al.: Effect of the solar activity variation on GITM 731

Figure 14. Density ratio (GITM / CHAMP) as a function of time and latitude for case e (LTAN = 07:38; LTDN = 19:38;

hapogee = 379 km; hperigee = 350 km).

Figure 15. Density ratio (GITM / CHAMP) as a function of time and latitude for case f (LTAN = 01:06; LTDN = 13:06;

hapogee = 353 km; hperigee = 323 km).

branch of the orbit. Figures 10, 11, 13 and 15 show that, around 30 latitude, there are two distinct zones of high

relative ratio values with respect to the rest of the nightside orbit branch. These high relative ratio values correspond to the under-prediction or the over-prediction of the equatorial mass anomaly (EMA) by GITM. This phenomenon, generated by the equatorial electrojet (EEJ) and recently described by Liu et al. (2005) using CHAMP data, creates an accumulation of density between 15 and 30 latitude poleward with respect to the magnetic equator. The interaction between the Earth magnetic eld and the ions transported by the EEJ causes a Lorentz lift-off mechanism that pushes ions and electrons to higher altitudes over the magnetic equator. When the particles reach the maximum altitude, they fall down and generate an accumulation of density with a fountain-like shape. Because of the experimental nature of the electrodynamics model in GITM (i.e. dynamo modelling), it was decided to disable such a feature during the investigation so that the EMA phenomenon cannot be correctly modelled in the current simulations. As can be seen from Figs. 12 and 14, when the orbit crosses the ascending/descending nodes at

Figure 16. Electron density ratio (GITM / CHAMP) as a function of time and latitude for case a (LTAN = 15:10; LTDN = 03:10).

Figure 17. Electron density ratio (GITM / CHAMP) as a function of time and latitude for case b (LTAN = 22:50; LTDN = 10:50).

dawn and dusk respectively, the ratio values of the ascending/descending branch of the orbit are quite similar to each other. Moreover, for both cases the under/over-prediction of the EMA is visible during the ascending branch as well as the descending branch of the orbit.

Figures 1621 represent the comparison between GITM simulations and CHAMP electron number densities as a function of latitude and time for all cases presented in Table 1. The values shown in the gures are again the ratios between the simulations and the ight data. As discussed above for Fig. 9, in the comprehensive analysis shown in Figs. 1621 it can be observed that GITM simulations are unable to correctly predict the electron density on the night-side at 30 over the Equator. All cases indicate that the ratio

between the simulations and the ight data is equal to 5 (or higher) on the nightside. Again, in Figs. 1621 it can be ob-served that the results are consistent over time and therefore representative of the behaviour of the model. The comparison between the simulated and measured electron densities shows similar features to those in the analysis carried out for the neutral densities. In particular, on the dayside the EMA phenomenon is still visible at 30 poleward over the mag

netic equator, where the typical double peak is either over-

www.ann-geophys.net/34/725/2016/ Ann. Geophys., 34, 725736, 2016

732 D. Masutti et al.: Effect of the solar activity variation on GITM

Figure 18. Electron density ratio (GITM / CHAMP) as a function of time and latitude for case c (LTAN = 06:44; LTDN = 18:44).

Figure 19. Electron density ratio (GITM / CHAMP) as a function of time and latitude for case d (LTAN = 10:01; LTDN = 22:01).

or underestimated by the simulations because of the missing electrodynamic model. The analysis over the poles and over the Equator shows again a different behaviour of the model in the two regions. In general, it can be observed that the ratios over the poles are closer to a value of 1 for all cases under investigation, whereas the ratio values over the Equator are more scattered.

3.3 Effect of the solar activity

Figure 22 summarises the behaviour on GITM for a wide spectrum of solar ux values. The ratios between the measured and simulated neutral densities are given as a function of the F10.7 index and for CHAMP and GRACE satellites. Every point shown in Fig. 22 is the result of an orbital average of the ratios between the simulations and the ight data. Figure 22 indicates that there is a linear dependency between the accuracy of the GITM simulations and the solar activity. In particular, the simulated neutral densities are underestimating the ight data by a factor of 0.4 during low solar activity (85 W m2 Hz1), whereas they are overestimating the ight data by a factor of 1.7 during high solar activity (260 W m2 Hz1). The density simulated by GITM

Figure 20. Electron density ratio (GITM / CHAMP) as a function of time and latitude for case e (LTAN = 7:38, LTDN = 19:38).

Figure 21. Electron density ratio (GITM / CHAMP) as a function of time and latitude for case f (LTAN = 1:06, LTDN = 13:06).

agrees well (1 : 1 ratio) with the ight measurement only for F10.7 values between 150 and 170 W m2 Hz1 (i.e. during medium solar activity).

Figure 23 shows the results obtained from the comparison between the simulated and measured electron densities as a function of the solar ux. The results shown are affected by a high degree of scattering, in particular on the nightside. As explained in Sect. 3.2, a wrong implementation of the electron continuity equation in GITM is the reason for a high density of electrons on the nightside. Because of such an implementation, electrons are not correctly downwelling to lower altitudes on the nightside, which is the reason for their accumulation and the extremely high ratio on the nightside.However, the results from the comparison on the dayside, even if they are still affected by a marked scattering, again show a linear dependency with the solar ux. In conclusion, the sensitivity of GITM to the solar ux is higher for the neutral density (solid line) than for the electron density (dashed line).

An investigation of the reasons for such sensitivity is currently ongoing and some preliminary causes are presented here.

Ann. Geophys., 34, 725736, 2016 www.ann-geophys.net/34/725/2016/

D. Masutti et al.: Effect of the solar activity variation on GITM 733

Figure 22. Effect of the solar activity on a large dataset of GITM simulations in comparison with CHAMP and GRACE neutral densities. The solid line represents the linear interpolation of all neutral density ratios.

Figure 23. Effect of the solar activity on a large dataset of GITM simulations in comparison with CHAMP electron number density. The solid line represents the linear interpolation of all neutral density ratios (from Fig. 22). The dashed line represents the linear interpolation of dayside-only electron number density ratios.

The transport properties, namely the dynamic viscosity and the thermal conductivity, are not correctly implemented in GITM. For example, the viscosity model presents a wrong trend as a function of temperature if compared to a simple Sutherland formulation. Nevertheless, a preliminary comparison by (Giulioni, 2016) between neutral densities computed with GITM using the original viscosity model and a Sutherland-based model shows a change of less than 5 % in density. The latter comparison suggests that the effect of the viscosity model in GITM is negligible; nevertheless, the effect of a different thermal conductivity model is yet to be investigated.

Figure 24. Effect of the new and old allocation of the exothermic excess energy on the simulated density. The density is simulated along the CHAMP trajectory as an example.

In the version of GITM used in this paper the photoelectron heating is modelled through a semi-empirical law controlled by a xed photoelectron heating efciency.This heating coefcient is known to have an important bias effect on the simulations as explained by Burrell et al. (2015). Its value is in fact obtained from data comparison because of the lack of ight measurements. In the last version of GITM, described in the paper by Zhu et al. (2016), the photoelectron efciency is now treated with an improved formulation that takes into account a heating efciency as a function of the altitude.

The coefcients for the chemical reaction rates implemented in GITM are mainly taken from Rees (1989) and Schunk and Nagy (2009). A number of these rates are given as constant, while others are specied as a function of the electron temperature or the mean between the ion and the neutral temperature. Recently, Richards and Voglozin (2011) re-examined the photochemical reactions in the ionosphere and presented a literature survey of the latest chemical reactions with updated coefcient rates in the upper atmosphere. An appreciable difference can be found by comparing the reactions in the paper of Richards and Voglozin (2011) and chemical processes implemented in GITM.

In GITM, all of the excess energy from the exothermic chemical processes is allocated into the neutrals and none into the ions. This causes an overheating of the model especially in times of large ionisation rates (i.e. with large F10.7 values). As explained by Rees (1989), this approach is not correct because the excess energy should be partitioned between the products of the chemical reaction with an inverse proportion to their masses.Figure 24 presents the results of a preliminary investigation on the effect of the different energy allocations on

www.ann-geophys.net/34/725/2016/ Ann. Geophys., 34, 725736, 2016

734 D. Masutti et al.: Effect of the solar activity variation on GITM

Figure 25. Effect of the new (solid lines) and old (dashed lines) allocation of the exothermic excess energy on the neutral, ion and electron temperature respectively (date: 10 March 2002, 00:00; 0

latitude and 0 longitude).

the simulated density. Figure 24 shows the comparison of the simulated densities along the CHAMP trajectory using the new (following Rees, 1989) and the old (original GITM implementation) allocation of the exothermic excess energy. Results show that by allocating the excess energy only in the neutral temperature energy equation, an articial heating of the model is generated, increasing the density by almost 70 % with respect to the implementation suggested by Rees (1989). This articial heating is also shown in Fig. 25, where the vertical proles of all temperatures (neutral, ion and electron) obtained with the original energy allocation are in general 1025 % higher with respect to the same temperatures obtained using the new energy allocation.

It is found by Giulioni (2016) that the NO cooling rates computed in GITM are at least 2 times higher with respect to the measured cooling rates by the TIMED/SABER instrument. The reasons for such discrepancy could be attributed to not only a wrong concentration of NO or O in the simulations but also incorrect modelling of the cooling mechanism and the related reaction rate coefcients. To further conclude on the NO cooling, the implementation in GITM must be compared to the most recent work in this eld presented by Hwang et al. (2003).

In general, Pawlowski and Ridley (2009) showed that, using the broad published values of thermal conductivity, NO cooling rates and other parameters, there is a signicant variation in the temperatures and densities derived from the model. This study shows that the default parameters in GITM cause the model to be too sensitive to F10.7 when compared to the mass density. The results presented by Giulioni (2016) using an improved chemistry module (which includes both

the updated reaction rates and the corrected excess energy allocation) show that the neutral densities computed with the new module are lower with respect to the same densities simulated with the original implementation of the chemistry module. Giulioni (2016) also shows the performances of the new chemistry module over a limited range of solar activities and observes that the sensitivity of GITM to the F10.7 solar ux is not yet xed.

4 Conclusions

From a grid sensitivity study on GITM a trade-off between the accuracy of the simulations and the computation time has been found. A grid with 72 36 cells (5 5 ) in lon

gitude and latitude respectively allows for results to be obtained with only 3 % difference from the most rened grid (144 72 cells or 2.5 2.5 ) in a quarter of the computa

tion time (1500 vs. 400 CPU hours).

The comprehensive validation of the GITM simulations and the ight data for both neutral densities and electron number densities conrms the consistency of the results during the different solar activity periods investigated here. The comparison shows the predicting capabilities of GITM, indicating that in general there is a relative under-prediction of the neutral density over the poles. The same comparison also shows that, without the electrodynamic model enabled(i.e. dynamo model), GITM is unable to correctly model the EMA, and this is visible especially on the nightside, where the comparison shows two distinct zones of under/over-prediction at 30 poleward with respect to the magnetic

equator.

It is also found that there is a clear dependency between the GITM bias and the solar activity. In particular, this dependency shows that the ratio between the simulations and the ight measurements is a linear function of the solar ux F10.7. From the data presented it is possible to observe that GITM is underestimating the ight data by a factor of 0.4 during periods of low solar activity(85 W m2 Hz1), whereas it is overestimating the ight data by almost a factor of 2 during periods of high solar activity (260 W m2 Hz1). The simulated neutral densities agree well with the ight measurements only for F10.7 values between 150 and 170 W m2 Hz1. The results also show that

GITM has a different sensitivity to the solar ux variation for the neutral density and the electron density respectively.In particular, the electron number density bias is less sensitive to a change in the solar ux, and the results have a lower bias with respect to the ight data on the dayside.

There are many reasons that could explain the high sensitivity of GITM to the solar activity. An ongoing study at the University of Michigan and the von Karman Institute is addressing this topic, and the answer is not currently known.Schemes underlying GITM such as the transport properties

Ann. Geophys., 34, 725736, 2016 www.ann-geophys.net/34/725/2016/

D. Masutti et al.: Effect of the solar activity variation on GITM 735

modelling and the chemical reaction rates are the suspected causes.

5 Code availability

The GITM code is a property of the University of Michigan and cannot be distributed.

6 Data availability

The results from the GITM simulations are publicly available. Because of the considerable amount of data (several gigabytes), the results can only be obtained on request from the main author. All the ight data used to compare with the simulations are third-party data (CHAMP Density v2.2: http://sisko.colorado.edu/sutton/data/ver2.2/champ/density/

Web End =http://sisko.colorado.edu/sutton/ http://sisko.colorado.edu/sutton/data/ver2.2/champ/density/

Web End =data/ver2.2/champ/density/ ; CHAMP Electron Density: http://isdc-old.gfz-potsdam.de/

Web End =http: http://isdc-old.gfz-potsdam.de/

Web End =//isdc-old.gfz-potsdam.de/ http://isdc-old.gfz-potsdam.de/

Web End = ; GRACE Density v2.2: http://sisko.colorado.edu/sutton/data/ver2.2/grace/density/

Web End =http:// http://sisko.colorado.edu/sutton/data/ver2.2/grace/density/

Web End =sisko.colorado.edu/sutton/data/ver2.2/grace/density/ ).

Author contributions. Davide Masutti prepared the manuscript with the contribution of all co-authors. In particular, Gunther March collected the data and prepared the preliminary simulations. Aaron Ridley provided continuous support with the GITM code and the simulations. Jan Thoemel reviewed the work.

Acknowledgements. The research leading to these results has received funding from the European Unions Seventh Framework Programme (FP7) under REA grant agreement no. 284427. Support for the University of Michigan is from the National Science Foundation (grant ATM1242839).

The topical editor, C. Jacobi, thanks two anonymous referees for help in evaluating this paper.

References

Bilitza, D.: International Reference Ionosphere 2000, Radio Sci.,36, 261275, doi:http://dx.doi.org/10.1029/2000RS002432

Web End =10.1029/2000RS002432 http://dx.doi.org/10.1029/2000RS002432

Web End = , 2001.

Bowman, B. R., Tobiska, W. K., Marcos, F., Huang, C. Y., Lin, C.S., and Burke, W. J.: A New Empirical Thermospheric Density Model JB2008 Using New Solar and Geomagnetic Indices, in: Astrodynamics Specialist Conference, Honolulu, Hawaii, 1821 August 2008, AIAA 2008-6438, 2008.

Bruinsma, S.: The DTM-2013 thermosphere model, J. Space Weather Space Clim., 5, A1, doi:http://dx.doi.org/10.1051/swsc/2015001

Web End =10.1051/swsc/2015001 http://dx.doi.org/10.1051/swsc/2015001

Web End = , 2015.

Bruinsma, S., Tamagnan, D., and Biancale, R.: Atmospheric densities derived from CHAMP/STAR accelerometer observations, Planet. Space Sci., 52, 297312, doi:http://dx.doi.org/10.1016/j.pss.2003.11.004

Web End =10.1016/j.pss.2003.11.004 http://dx.doi.org/10.1016/j.pss.2003.11.004

Web End = , 2004.

Burrell, A., Goel, A., Ridley, A., and Bernstein, D.: Correction of the photoelectron heating efciency within the global ionosphere-thermosphere model using Retrospective Cost Model Renement, J. Atmos. Sol.-Terr. Phy., 124, 3038, doi:http://dx.doi.org/10.1016/j.jastp.2015.01.004

Web End =10.1016/j.jastp.2015.01.004 http://dx.doi.org/10.1016/j.jastp.2015.01.004

Web End = , 2015.

Cooke, D. L., Turnbull, W., Roth, C., Morgan, A., and Redus, R.: Ion Drift-Meter Status and Calibration, in: First Champ Mission Results for Gravity, Magnetic, and Atmospheric Studies, Springer, Berlin, 212219, doi:http://dx.doi.org/10.1007/978-3-540-38366-6_31

Web End =10.1007/978-3-540-38366-6_31 http://dx.doi.org/10.1007/978-3-540-38366-6_31

Web End = , 2003.

Deng, Y., Richmond, A. D., Ridley, A. J., and Liu, H.-L.: Assessment of the non-hydrostatic effect on the upper atmosphere using a general circulation model (GCM), Geophys. Res. Lett., 35, L01104, doi:http://dx.doi.org/10.1029/2007GL032182

Web End =10.1029/2007GL032182 http://dx.doi.org/10.1029/2007GL032182

Web End = , 2008.

Doornbos, E.: Thermospheric Density and Wind Determination from Satellite Dynamics, Springer-Verlag, Berlin Heidelberg, doi:http://dx.doi.org/10.1007/978-3-642-25129-0

Web End =10.1007/978-3-642-25129-0 http://dx.doi.org/10.1007/978-3-642-25129-0

Web End = , 2012.

Giulioni, L.: Accurate modeling of the gas properties in the thermosphere under high and low solar activity conditions, Project report rm 2016-12, von Karman Institute, Sint-Genesius-Rode, Belgium, 2016.

Harris, M. J.: A new coupled middle atmosphere and thermosphere circulation model: Studies of dynamic, energetic and photochemical coupling in the middle and upper atmosphere, PhD thesis, University College London, UK, 2001.

Hedin, A. E.: Extension of the MSIS Thermospheric Model into the Middle and Lower Atmosphere, J. Geophys. Res., 96, 1159, doi:http://dx.doi.org/10.1029/90JA02125

Web End =10.1029/90JA02125 http://dx.doi.org/10.1029/90JA02125

Web End = , 1991.

Hwang, E. S., Castle, K. J., and Dodd, J. A.: Vibrational relaxation of NO(v = 1) by oxygen atoms between 295 and 825 K, J. Geo

phys. Res.-Space, 108, 1109, doi:http://dx.doi.org/10.1029/2002JA009688

Web End =10.1029/2002JA009688 http://dx.doi.org/10.1029/2002JA009688

Web End = , 2003. Kang, Z., Tapley, B., Bettadpur, S., Ries, J., Nagel, P., and Pastor,R.: Precise orbit determination for the GRACE mission using only GPS data, J. Geodesy, 80, 322331, doi:http://dx.doi.org/10.1007/s00190-006-0073-5

Web End =10.1007/s00190- http://dx.doi.org/10.1007/s00190-006-0073-5

Web End =006-0073-5 , 2006.

Krauss, S., Temmer, M., Veronig, A., Baur, O., and Lammer,H.: Thermospheric and geomagnetic responses to interplanetary coronal mass ejections observed by ACE and GRACE: Statistical results, J. Geophys. Res.-Space, 120, 88488860, doi:http://dx.doi.org/10.1002/2015JA021702

Web End =10.1002/2015JA021702 http://dx.doi.org/10.1002/2015JA021702

Web End = , 2015.

Liu, H., Luehr, H., Henize, V., and Koehler, W.: Global distribution of the thermospheric total mass density derived from CHAMP, J. Geophys. Res., 110, A04301, doi:http://dx.doi.org/10.1029/2004JA010741

Web End =10.1029/2004JA010741 http://dx.doi.org/10.1029/2004JA010741

Web End = , 2005. Liu, H., Stolle, C., Forster, M., and Watanabe, S.: Solar activity dependence of the electron density at 400 km at equatorial and low latitudes observed by CHAMP, J. Geophys. Res., 112, A11311, doi:http://dx.doi.org/10.1029/2007JA012616

Web End =10.1029/2007JA012616 http://dx.doi.org/10.1029/2007JA012616

Web End = , 2007.

Lu, G., Richmond, A. D., Emery, B. A., and Roble, R. G.: Magnetosphere-ionosphere-thermosphere coupling: Effect of neutral winds on energy transfer and eld-aligned current,J. Geophys. Res., 100, 1964319659, doi:http://dx.doi.org/10.1029/95JA00766

Web End =10.1029/95JA00766 http://dx.doi.org/10.1029/95JA00766

Web End = , 1995.

McNamara, L. F., Cooke, D. L., Valladares, C. E., and Reinisch, B. W.: Comparison of CHAMP and Digisonde plasma frequencies at Jicamarca, Peru, Radio Sci., 42, RS2005, doi:http://dx.doi.org/10.1029/2006RS003491

Web End =10.1029/2006RS003491 http://dx.doi.org/10.1029/2006RS003491

Web End = , 2007.

Millward, G. H., Moffett, R. J., Quegan, S., and Fuller-Rowell,T. J.: A coupled thermosphere ionosphere plasmasphere model (CTIP), in: STEP Handbook on Ionospheric Models, edited by: Schunk, R. W., Utah State University, 239279, 1996. Pawlowski, D. J. and Ridley, A. J.: Quantifying the effect of thermospheric parameterization in a global model, J. Atmos. Sol.-Terr. Phy., 71, 20172026, doi:http://dx.doi.org/10.1016/j.jastp.2009.09.007

Web End =10.1016/j.jastp.2009.09.007 http://dx.doi.org/10.1016/j.jastp.2009.09.007

Web End = , 2009.

www.ann-geophys.net/34/725/2016/ Ann. Geophys., 34, 725736, 2016

736 D. Masutti et al.: Effect of the solar activity variation on GITM

Picone, J., Hedin, A., Drob, D., and Aikin, A.: NRLMSISE-00 empirical model of the atmosphere: Statistical comparisons and scientic issues, J. Geophys. Res., 107, 1468, doi:http://dx.doi.org/10.1029/2002JA009430

Web End =10.1029/2002JA009430 http://dx.doi.org/10.1029/2002JA009430

Web End = , 2002.

Rees, M. H.: Physics and chemistry of the upper atmosphere, Cambridge University Press, Cambridge, 1989.

Reigber, C. and Schwintzer, P.: A Challenging Microsatellite Pay-load for Geophysical Research and Application, in: Small Satellites for Remote Sensing, Space Congress 95, Bremen, 1995.Richards, P. G. and Voglozin, D.: Reexamination of ionospheric photochemistry, J. Geophys. Res.-Space, 116, A08307, doi:http://dx.doi.org/10.1029/2011JA016613

Web End =10.1029/2011JA016613 http://dx.doi.org/10.1029/2011JA016613

Web End = , 2011.

Richmond, A. D., Ridley, E. C., and Roble, R. G.: A thermosphere/ionosphere general circulation model with coupled electrodynamics, Geophys. Res. Lett., 19, 601604, doi:http://dx.doi.org/10.1029/92GL00401

Web End =10.1029/92GL00401 http://dx.doi.org/10.1029/92GL00401

Web End = , 1992.

Ridley, A. J., Deng, Y., and Toth, G.: The Global Ionosphere-Thermosphere Model, J. Atmos. Sol.-Terr. Phy., 68, 839864, doi:http://dx.doi.org/10.1016/j.jastp.2006.01.008

Web End =10.1016/j.jastp.2006.01.008 http://dx.doi.org/10.1016/j.jastp.2006.01.008

Web End = , 2006.

Sagalyn, R. C. and Bowhill, S. A.: Progress in Geomagnetic Storm

Prediction, in: Geophysical Monograph Series, 01/1993, American Geophysical Union, doi:http://dx.doi.org/10.1029/GM073p0157

Web End =10.1029/GM073p0157 http://dx.doi.org/10.1029/GM073p0157

Web End = , 1993.Schunk, R. and Nagy, A.: Ionospheres: Physics, Plasma Physics, and Chemistry, Cambridge Atmospheric and Space Science Series, 2nd Edn., Cambridge University Press, 2009.

Sutton, E. K.: Normalized Force Coefcients for Satellites with Elongated Shapes, J. Spacecraft Rockets, 46, 112116, doi:http://dx.doi.org/10.2514/1.40940

Web End =10.2514/1.40940 http://dx.doi.org/10.2514/1.40940

Web End = , 2009.

Sutton, E. K., Forbes, J. M., and Nerem, R. S.: Global thermospheric neutral density and wind response to the severe 2003 geomagnetic storms from CHAMP accelerometer data, J. Geophys. Res., 110, A09S40, doi:http://dx.doi.org/10.1029/2004JA010985

Web End =10.1029/2004JA010985 http://dx.doi.org/10.1029/2004JA010985

Web End = , 2005.

Thoemel, J., Singarayar, F., Scholz, T., Masutti, D., Testani, P., Asma, C., einhard, R., and Muylaert, J.: Status of the QB50 cubesat constellation mission, in: 65th International Astronautical Congress, Toronto, Canada, 29 September3 October 2014, IAC-14.B4.2.11, 2014.

Touboul, P., Foulon, B., and LeClerc, G. M.: STAR, The Accelerometer of the Geodesic Mission CHAMP, in: Proceedings of the 49th IAF Congress, Melbourne, Australia, 28 September2 October 1998, IAF-98-B.3.07, 1998.

Vallado, D. A. and Finkleman, D.: A critical assessment of satellite drag and atmospheric density modeling, Acta Astronaut., 95, 141165, doi:http://dx.doi.org/10.1016/j.actaastro.2013.10.005

Web End =10.1016/j.actaastro.2013.10.005 http://dx.doi.org/10.1016/j.actaastro.2013.10.005

Web End = , 2014.

Zhu, J., Ridley, A. J., and Deng, Y.: Simulating electron and ion temperature in a global ionosphere thermosphere model: Validation and modeling an idealized substorm, J. Atmos. Sol.-Terr. Phy., 138139, 243260, doi:http://dx.doi.org/10.1016/j.jastp.2016.01.005

Web End =10.1016/j.jastp.2016.01.005 http://dx.doi.org/10.1016/j.jastp.2016.01.005

Web End = , 2016. Zhukov, V. N., Melnikov, E. K., and Smirnov, A. I.: Solving the

International Space Station (ISS) motion control long-term planning problem, in: 23rd International Symposium on Space Flight Dynamics, Jet Propulsion Laboratory, 2012.

Ann. Geophys., 34, 725736, 2016 www.ann-geophys.net/34/725/2016/