Thermospheric mass density variations during geomagnetic storms and a prediction model based on the merging electric field

With the help of four years (2002–2005) of CHAMP accelerometer data we have investigated the dependence of low and mid latitude thermospheric density on the merging electric field,Em, during major magnetic storms. Altogether 30 intensive storm events (Dst min < −100 nT) are chosen for a statistical study. In order to achieve a good correlationEm is preconditioned. Contrary to general opinion, Em has to be applied without saturation effect in order to obtain good results for magnetic storms of all activity levels. The memory effect of the thermosphere is accounted for by a weighted integration of Em over the past 3 h. In addition, a lag time of the mass density response to solar wind input of 0 to 4.5 h depending on latitude and local time is considered. A linear model using the preconditioned Em as main controlling parameter for predicting mass density changes during magnetic storms is developed: ρ = 0.5Em+ρamb, where ρamb is based on the mean density during the quiet day before the storm. We show that this simple relation predicts all storm-induced mass density variations at CHAMP altitude fairly well especially if orbital averages are considered.


Introduction
The variation of the thermospheric mass density during geomagnetic storms is a rather complex phenomenon.Both the density and the composition experience a series of dramatic changes.The thermospheric response comprises the effects of Joule/particle heating, Lorentz force, thermal expansion, upwelling, and horizontal wind circulation.Studying the Correspondence to: R. Liu (ruosi@gfz-potsdam.de)thermospheric density distribution during storm-time is of great importance not only for the precise orbit determination of satellites flying at several hundred kilometers altitude, but also for the interpretation of magnetosphere-ionosphericthermosphere (MIT) coupling.The relationship between geomagnetic input and ionospheric and thermospheric density has been studied for decades (e.g.Matuura, 1972;Prölss, 1980;Crowley, 1991;Forbes et al., 1996).Meanwhile many atmospheric models are developed with the purpose of predicting the mass density as close as possible.Pätzold's model (Pätzold, 1963) is one of the first that contains the geomagnetic heating as a cause for the density enhancement.Signatures of particle energy flow into the upper atmosphere were first documented by Jacchia in 1959.In 1964 a K p or a p dependent exospheric temperature was included in the Jacchia model (Jacchia and Slowley, 1964).Moe and Moe developed a global density model in 1975, based on measurements of Spades and Logacs satellites (Moe et al., 1977).More recent models like the MSIS series (e.g.Hedin, 1991;Picone et al., 2002), based on incoherent scatter radar measurements of temperature and in situ satellite measurements of density and composition, try to consider the magnetospheric and solar influence as good as possible.
With the launch of CHAMP and GRACE satellites measurements from the accelerometers on board enabled systematic studies of the thermospheric mass density on a global scale (Bruinsma et al., 2004;Liu et al., 2005).The global measurements of these satellites have been compared with predictions of models in many studies (e.g.Bruinsma et al., 2004Bruinsma et al., , 2006;;Forbes et al., 2005;Sutton et al., 2005;Liu and Lühr, 2005).All the studies have shown that the model predictions are in general agreement with CHAMP/GRACE measurements during quiet times, but they significantly underestimate the density during enhanced geomagnetic activity.The open issue to be addressed is, how to predict the thermospheric density as close as possible, especially during geomagnetic storms.This question is further extended to another: is there a geophysical quantity to which the density responds most efficiently or to which it is related most closely?
In recent years many efforts have been made to study the density dependence on various controlling parameters both during geomagnetically quiet and active times.Lathuillère et al. (2008) have studied the global thermosphere response to auroral magnetic activity forcing at mid and low latitudes using a method based on a singular value decomposition of the satellite data.They separated the spatial variation from the time variation, while the former is captured by the singular variation and the latter is captured by the projection coefficient.The projection coefficient is then used to define a disturbance coefficient, which is found to correlate better with the density than the a p index during storms.Zhou et al. (2008) have studied the thermospheric density response to the total global Joule heating power, Q j and the high resolution ring current index, Sym-H.A statistical study based on CHAMP mass density data from 19 great magnetic storms during 2001 to 2004 has been carried out.With Q j and Sym-H as the main controlling parameters they improved the prediction quality of the NRLMSISE-00 model during storm-time.Müller et al. (2009) performed a statistical analysis using CHAMP density data from low and equatorial regions of the years 2002 to 2005.They provide a comprehensive study of density dependences on the main drivers of the thermosphere, such as solar flux (P10.7),local time, season and magnetic activity.They have provided functional relations for all of these drivers.A systematic analysis of interplanetary magnetic field (IMF) B y and B z influences on high-latitude thermospheric density is provided by Kwak et al. (2009).In their study density observations from CHAMP during 17 October 2001 through 24 February 2002 are statistically investigated as a function of the direction and strength of the interplanetary magnetic field (IMF) for the southern summer hemisphere.Their study shows that the thermospheric density poleward of 60 • magnetic latitude tend to be strongest when IMF Bz is negative, and weakest when Bz is positive; for positive IMF B y the density changes appear in opposite time sectors to those for negative B y .Burke et al. (2007) have investigated thermospheric density correlation with Dst, polar cap potential, PC , and magnetospheric electric field, E VS , by making use of GRACE accelerometer measurements during two geomagnetic storm periods in 2004.They pointed out that A p (or a p as well), which reflects the responses of ground magnetometers to ionospheric currents at mid latitudes, significantly underestimates the thermospheric energy budgets during severe geomagnetic storms.Moreover, it has been shown in their Fig. 2 that during storm-time the thermosphere responds faster than the approximation provided by the a p index.These two points can explain to some extent why the standard models fail to reproduce storm-time thermospheric density distribution with reasonable fidelity.
In this study we aim to investigate, how the thermospheric density responds to solar wind inputs as quantified by the merging electric field, E m , during magnetic storms.Based on a period of four years starting from 2002 density data derived from CHAMP accelerometer measurements of 30 magnetic storms are investigated along with their correlation to the preprocessed E m .At the end we will bring forward a linear model using the merging electric field as the dominant controlling parameter to predict storm-time mass density distribution at 400 km altitude.
In the next section the datasets are introduced, which include CHAMP accelerometer measurements needed for the density and the ACE observations needed for deriving the merging electric field.Section 3 discusses first the correlation between density and E m , after that an empirical formula is introduced for reproducing the storm-time density variations along the CHAMP orbit.Section 4 shows an example of prediction using the formula obtained in Sect.3. Then we discuss the results and draw conclusions in Sects.5 and 6.

Dataset
The Challenging Minisatellite Payloads (CHAMP) spacecraft was launched on 15 July 2000.It cycles the Earth in a near-polar orbit with an inclination of 87.25 • .Within 131 days the satellite orbit covers all the local times.In this study a time period from 2002 to 2005 is considered.During these four years the orbital altitude decayed from 425 km to 360 km.CHAMP is very suitable for long-term monitoring of thermospheric characteristics due to its orbit height and the coverage of all latitudes and local times.
The mass density data are derived from measurements of the STAR (Space Triaxial Accelerometer for Research missions) sensors on CHAMP.Along-track axis data are used to calculate the density.The data are preprocessed by omitting spurious spikes and accelerations caused by solar radiation pressure or attitude maneuvers.The time resolution of the processed data is 10 s.As the spacecraft is flying, the aerodynamic acceleration vector, a, can be expressed as: where m is the spacecraft mass (≈ 500 kg), v r the spacecraft velocity relative to the atmosphere.C a is the aerodynamic force coefficient vector, the area A ref is a fixed reference value used to make C a dimensionless.A more detailed description of the calculating algorithms can be found in Doornbos et al. (2010).The total mass density of the air, ρ, is determined from the projection of the aerodynamic acceleration on the x-axis of the spacecraft body-fixed (SBF) coordinates.The density is deduced from the along track (x) component of the vectors in Eq. ( 1): Ann. Geophys., 28,[1633][1634][1635][1636][1637][1638][1639][1640][1641][1642][1643][1644][1645]2010 www.ann-geophys.net/28/1633/2010/Neglect of in-track wind in the density calculation can lead to errors.The meridional wind speeds at mid and low latitudes is in the range of 0-100 m/s under most conditions, which is small compared to the orbital velocity of 7.6 km/s.Another possible error source is the corotation of the atmosphere contributing to the relative velocity of the spacecraft.The corotation velocity over the equator is in the range of 490 m/s at CHAMP altitude.This has been taken into account in the analysis.We will not discuss the errors in detail because in this study only major storms are investigated.Given the density enhancement of several hundred percent during storms, the measurements allow for a significant level of physical interpretation even with a few percent uncertainty.
Based on the conclusion of Burke et al. (2007) that the polar cap potential PC correlates with density fairly well during all storm phases, PC could be used to predict the density changes during storms.Nevertheless, the availability of PC measurements from DMSP is poor.As a matter of fact, the merging electric field, E m , is a physical quantity which closely correlates with PC and can be easily obtained from solar wind and IMF data.From this point of view we assume that E m might correlate with density well during all storm phases, and due to its easy accessibility E m could be a more suitable parameter than PC to predict the density changes during storms.According to Kan and Lee (1979) PC is due to the component of E m perpendicular to the geomagnetic field lines in the reconnection region.E m can be written as: where v SW is the solar wind speed, is the transverse component of IMF in Geocentric-Solar-Magnetospheric (GSM) coordinates, θ is the angle between the interplanetary and the geomagnetic field in the reconnection, or so-called "merging" region.It approximately equals the polar angle between B T and the GSM Z-axis, i.e. tan(θ The solar wind and IMF data are taken from the Advanced Composition Explorer (ACE) spacecraft.ACE is an explorer mission that was launched on 25 August 1997 and positioned near the L1 Lagrangian point which is a point of Earth-Sun gravitational equilibrium about 1.5 million km from Earth and 148.5 million km from the Sun.It takes typically 20-60 min for the solar wind observed by ACE to arrive at the Earth's magnetopause.The data used are from the magnetic field instrument, MAG (Smith et al., 1998), at 16 s resolution and from the plasma instrument SWEPAM (McComas et al., 1998) at 64 s resolution.Both are resampled to 1min resolution and then time shifted using the so-called phase front propagation technique (Weimer et al., 2003) to represent the solar wind and IMF conditions at the front side magnetopause, which is assumed to be located at a distance of 10 R E at the sub-solar point.Finally, an additional average delay of 15 min has been added as travel time to the ther-mosphere for solar wind input at the magnetopause (Vennerström et al., 2002).

Methodology
In this section the storm-time thermospheric density dependence on the E m is studied, then a linear empirical model is proposed using the E m as the primary controlling parameter to predict density distribution during storms

Density processing
We concentrate our investigations at mid and low latitude, i.e. within ±60 • geomagnetic latitude (MLAT).This area plays a dominant role because it covers 90% of the Earth's surface.The high latitude area is omitted because during storm-time the heating in this area causes thermal expansion and thermospheric composition change.Therefore, at high latitude regions we may find density enhancements or depletions in response to Joule heating (e.g., Lei et al., 2010).
Before we analyse the density measurements, each CHAMP orbit is first divided into an ascending and a descending half, which are subdivided into three latitudinal segments, namely −60 • to −30 • , −30 • to 30 • , and 30 • to 60 • .Thus every orbit is divided into 6 MLT-latitude bins.By this means the effects of latitude, season and local time are separated.For each bin one average density value per overflight is calculated.We use averaged density data along the orbit that eliminates local features to a great extent, since we are interested in the large-scale and global features.

E m processing
The processing of the merging electric field, E m , in this study includes two steps -truncation and integration.In the following we will describe the two procedures in detail.
The Hill model (Hill et al., 1976) predicts that PC saturates for strong solar wind driving.Nagatsuma (2002) concluded that the polar cap potential tends to be saturated when the value of E m exceeds 5 mV/m, and this saturation only depends on the intensity of E m .In this study we show that the saturation effect can be presented by the E m , such that E m has an approximately linear relationship with PC .The saturation equation we use for E m is expressed as The relation between the transpolar potential of PC and the merging electric field E m .PC is calculated using the method described in Ober and Maynard (2003).The blue curve shows that PC saturates when E m increases, while the red curve shows a almost linear relation between PC and E m(sat) .
where the value of 8 mV/m has been found as a suitable E m truncation value.To test the efficiency of this equation, we calculated PC according to the Hill model and make a comparison of its respective variations with E m and E m(sat) .The transpolar potential, PC , given by the Hill model can be written as it assumes that the transpolar potential consists of an unsaturated magnetospheric convection potential M and a saturated transpolar potential S .
where E SW , the solar wind electric field, is defined as v SW B T sin(θ/2).P SW = nm p v 2 SW is the solar wind ram pressure, n is the solar wind number density, and m p is the proton mass, P is the ionospheric Pedersen conductance.A value of P between 6 to 10S fits best to DMSP measurements of PC .The constant 30 kV in Eq. ( 6) is the baseline potential not attributable to E SW (Ober and Maynard, 2003).
In our calculation of PC using the approach above, we select P SW = 5 nPa, P = 6S, and E m is used in Eq. ( 6) instead of E SW .With those inputs we obtain S ≈ 450 kV.We treat the calculation result as the geo-effective component of PC .What we are really concerned of is the correlation between PC and E m .The variation of PC as function of E m is shown by the blue solid curve in Fig. 1.The black dashed line marks the saturation value of the merging electric field E = 8 mV/m.PC increases rapidly as the E m increases when E m is below 8 mV/m, but when E m reaches higher values, the increasing rate of PC begins to drop off, and PC will finally saturate somewhere above 400 kV.The red curve reveals the variation of PC versus E m(sat) .Apparently, PC is related much more linearly to the E m(sat) .If PC reaches 350 kV or higher, the linear relation will break down, but as long as the PC is below 250 kV, we can treat the relation between PC and E m(sat) as linear.
Figure 2 shows one example of the saturated electric field E m(sat) in contrast to the original E m computed from the ACE measurements for the storm in 22-28 July 2004.The black curve reflects the original E m , exhibiting three peaks around the storm times 22 UT, 80 UT and 130 UT, respectively.The first and second peaks reach to about 12 mV/m, the last one reaches even 21 mV/m.Between the peaks the value drops below 3 mV/m.The E m(sat) , plotted by the blue curve, duplicates the E m perfectly at values below 4 mV/m, but is dramatically truncated nearby the peaks, so that all the peaks are confined under the saturation level.The higher the peak, the more remarkable the truncation.For example, the first and second peaks are truncated to 60%, while the third peak is truncated to 35%.
It is pointed out by Richmond et al. (2003) that, due to the inertia of the air, the changes of thermospheric wind lag behind changes in the IMF.In their study densities are correlated with lagged and time-averaged IMF values.This point of view holds also true for the merging electric field.Therefore we use the same procedure as they do for deducing effective E m values.Based on their Eq.(1), the lagged, timeaveraged E m can be defined as where E m is treated as continuous function of time t, t = 0 is chosen 24 h before the first E m data is used, τ is the effective averaging and lag time for the exponential weighting function in the integrands.The correlation between the mass density and E m differs only weakly when τ ranges from 1 h to 10 h.So the exact choice of τ is not critical.Richmond et al. (2003) used τ = 3 h for their calculation of IMF By and Bz.Since this value gives a relative good correlation in our study, it is also used for our case.The green curve in Fig. 2 shows the E m for 22-28 July 2004 calculated with τ = 3 h.E m follows the variation of E m quite closely with a 3 h delay.All the positive and negative spikes seen in E m are now smoothed out by the integration, by which the electric field is better correlated with the thermospheric density.
The red curve in Fig. 2 shows E m(sat) as derived from Eq. ( 8) when applied to the truncated E m(sat) as done in the first step.The E m(sat) shows smoothing and a 3h-delay effect compared to the E m(sat) , and a truncation effect compared to the E m .In the following E m(sat) is used to develop a predictive model for storm-time density evolution.As well, E m(sat) is calculated for each passage over the 6 latitude bins mentioned in Sect.3.1.

A linear model
In this section we want to present an empirical model of mass density around 400 km altitude versus E m(sat) .This model will later serve to predict the density variations during geomagnetic activities.A linear relation is selected under the assumption that the density variations are linearly correlated with the merging electric field.
The empirical linear equation used in this study is: where ρ is the storm-time mass density at 400 km, a and b are the fitting coefficients acquired from the line regression.The fitting is done for each storm individually.Here we would like to mention that the fitting also takes into account a delay time between the density and the merging electric field, which will be fully discussed in the next subsection.In this subsection we only want to mention the conception and focus on the discussion of forming the linear equation.The values of a obtained from the super storms become outliers of the scatter plots, which is caused by the truncation done in Sect.3.2.During the super storms the value of E m increases greatly above the truncation threshold, sometimes the enhancement reaches 300% or even more.Therefore the great reduction of E m(sat) incurred by truncation needs to be balanced by a much higher factor a. Considering the significance of the super storms, it is unreasonable to exclude them as exceptions, when we are trying to build up a density model for storm-time.Therefore the truncation of E m might not reflect the reality too well.We have redone the analysis with the E m by skipping the truncation step, as is shown by the green curve in Fig. 2. The new results are shown in Fig. 4. The right two panels for b remains quite the same as in Fig. 3.In the left two panels for a, now the red dots from the super storms fit quite well into the main scatter area.Moreover, the scatter in both panels is reduced by 20%, showing a more concentrated distribution between 0.2 to 0.6.A median value a ≈ 0.5 can be obtained from the two panels in Fig. 5.The standard deviation of a drops from 0.25 in Fig. 3 to 0.15 in Fig. 4.
To testify our assumption that the relation between thermospheric density and E m is linear, we have done the scatter plots of density versus E m for the 6 latitude bins for the storm during 22-28 July 2004, as shown in Fig. 5.The top panels are for the day-side (12:00 MLT) and the bottom panels for the night-side (00:00 MLT).From left to right, the three columns are for the northern mid latitude, low latitude and southern mid latitude, respectively.A linear fitting is shown by the red straight line in each panel, with the line equation at the top left corner.For each of the latitude bins the density displays a clear linear relation with E m .The linear relation can be seen more clearly in the nightside, where it is less affected by sunlight.The correlation coefficient for the nightside is above 0.9, compared to coefficients around 0.8 for the dayside.The slopes are approximately the same on the dayand nightside.
Since we would like to deduce a more or less general predicting formula for all storm events, a constant value of 0.5 is used for a in Eq. ( 9).For b it is hard to deduce an universal number from Fig. 4. From the physical point of view b in Eq. ( 9) represents the ambient density before the storm which is rather dependent on the prevailing solar flux level.So we use the ambient density ρ amb to replace b in Eq. ( 9), where ρ amb can be individually obtained from a geomagnetically quiet day (reference day) prior to each storm.Müller et al. (2009) have deduced the mass density dependence on solar flux during quiet days.We use this to account for changing solar flux level (P10.7)during the storm.From the fitting equations provided in their Fig.3a, b the changes with respect to a reference day can be expressed as: where ρ amb is the ambient density value for a storm day, ρ r is the mean density derived from the reference day, P10.7 and P10.7 r represent the corresponding daily solar flux index for the storm day and the reference day, respectively.The final linear model for thermospheric mass density at 400 km is rewritten as: ρ(10 −12 kg/m 3 )=0.5Em (mV/m)+ρ amb (10 −12 kg/m 3 ) (11)

Delay time of the density
It is known that the mass density at 400 km altitude reacts after a delay time, which is expected to depend on latitude and MLT.The best delay time of density response is determined by cross-correlating with E m .We tested delays from 0 to 10.5 h, which correspond to 0 to 7 orbital periods.In order to investigate the dependence of MLT, the delays are sorted into 4 different MLT sectors: 05:00 MLT to 09:00 MLT as the morning sector, 10:00 MLT to 16:00 MLT as the noon sector, 17:00 MLT to 20:00 MLT as the evening sector, and 21:00 MLT to 04:00 MLT as the night sector.
Figure 6 shows bar diagrams of determined optimal delay times of all 30 storm events in the four MLT sectors for low (top) and mid (bottom) latitudes separately.Basically, the time lag of mass density changes behind E m varies from storm to storm.More than 95% of the delays are within 6 h.For the low latitude, no clear MLT dependence can be found.Except for the morning sector, which gives a relative smaller delay time of 1.5 h, the other three sectors give a delay time of 3 h.The low latitude band shows a good correlation with delayed E m .In 90% of the events a correlation coefficient over 0.8 is obtained (not shown here).At mid latitudes the delay time tends to be smaller at morning (1.5 h) and noon (0 h) sectors and larger at evening and night sectors (both 4.5 h).A correlation coefficient over 0.8 can be found in 73% of the events.A few low correlation coefficients are found in the morning and evening sectors, indicating relative bad correlation in these sectors.
Based on the results presented in Fig. 6, we use hereafter a constant delay of 3 h for low latitude.For mid latitude, variable delays of 0 h, 1.5 h, 4.5 h are used for noon, morning and evening/night, respectively.

Prediction results
The empirical model derived in the above section is used to reproduce the density variation at 400 km altitude during geomagnetic storms.We take the storm series in 22-28 July 2004 to test the prediction capability of this model.The six panels of Fig. 7 presents the prediction (red curves) for the six latitude classes along with the corresponding CHAMP measurements (blue curve).Also shown at the top left corner in each panel are the correlation coefficients, R, between measurements and predictions, the mean value E( ) and standard deviation σ ( ) of the relative error , where = prediction-measurement measurement × 100%.The CHAMP measurements show 3 prominent peaks coinciding with the 3 E m peaks in Fig. 2. The measurements show day-night and summer-winter asymmetry in magnitude.The nightside density reacts 2-3 h later than the dayside, which is consistent with our analysis of density delay time.The main features of the density measurements are correctly reproduced by the model.Correlation coefficients between 0.83 to 0.95 indicate good agreement between measurements and predictions.The correlation is (1) better on the nightside (bottom panels) than on the dayside, (2) better at low latitudes than at mid latitudes.And amplitudes are predicted better in the summer than in the winter hemisphere.One may notice that, the predictions in general overestimate the measurements since we choose a constant value of a for all E m , the overestimation level is different from peak to peak.For example, the first peak is generally overestimated by 20%, while the third peak The four columns show respectively the delay times of density at mid latitudes for morning, noon, evening and night sectors, where the morning sector represents 05:00 to 09:00 MLT, the noon sector 10:00 to 16:00 MLT, the evening sector 17:00 to 20:00 MLT, and the night sector 21:00 to 04:00 MLT.11).The six panels represent the 6 latitude/MLT bins.The blue curve in each panel represents the density measurements from CHAMP and the red curve denotes the prediction.The correlation coefficient, R, between the measurement and prediction, as well as the mean and standard deviation of relative error E( ) and σ ( ), are given in each panel.by 15%, therefore the standard deviation of has a relatively large value between 20% to 30%.Moreover, the bottom right panel shows an overall overestimation throughout the whole event.This is probably due to the combined effect of the season (winter) and MLT (midnight).
It is known that the mass density has local time and seasonal asymmetry as seen in Fig. 7.If we average the data over one orbit, the density from the day-and night-side, summer and winter hemisphere effects tend to compensate each other, hence the orbit-averaged density is expected to correlate better with the geophysical forces than the segmented density average.The orbit-averaged density is derived by averaging the mean density of the 6 segments in one orbit, where the two low latitude bins are weighed by a factor of 2, because they have a latitude range of 60 • rather than 30 • .The merging electric field is processed in the same way as described in Sect.3, here a constant delay time of 3h derived from the mean delay time of orbit-averaged density is taken into account.The linear model remains the same.Figure 8 presents the prediction of orbit-averaged density.The orbit-averaged density measurement is much smoother than the segment-averaged ones.The prediction result is quite similar to that in Fig. 7.As expected, the correlation coefficient, R, has a high value of 0.93.The first and third density peak are overestimated by only 7%, and the valley between the first and second peak is underestimated by 20%, but both the mean value and standard deviation are smaller than the results in Fig. 7.The mean relative error is quite small (2.6%), this means the over-predicted and under-predicted parts largely counterbalance.

Discussion
In this paper we have investigated how the thermospheric density at low and mid latitude responds to the solar wind input during geomagnetic storms.As a controlling parameter we chose the merging electric field, E m .To our knowledge the correlation between E m and density has so far not yet been studied thoroughly.interpreted for 3 latitude ranges on both dayside and nightside, in order to separate the effects of local time and geomagnetic latitude.The applied E m(sat) is preconditioned by considering a saturation effect and a lagged response of density features to solar wind input.However, the saturation effect is later removed to obtain more reasonable predictive results.For each storm event the delay time is investigated, and the results are interpreted by latitude range as function of MLT.Linear regressions have been carried out for each storm to obtain an empirical model using E m as the main input to predict thermospheric density distribution at 400 km during storm-time.The final empirical relation is written as ρ = 0.5E m + ρ amb .An example for predicting the mass density response of the storm event on 22-28 July 2004 is shown at the end.

The delayed response of density
During magnetic storms, the density responds to E m with a delay time, which varies somewhat from one storm to another.Based on the statistical study we conclude that the delay time at low latitude is about 3 h, at mid latitude it varies with local time, from 0 h at noon to 4.5 h at night-time.With respect to the orbit-averaged density, the mean delay time is 3 h, which is quite similar to that of low latitude.The delay time is difficult to determine precisely, for it shows dependence on various factors.Three kinds of dependences are herewith summarized: 1.The stronger the storm, the smaller the delay time.This trend was also found by Zhou et al. (2008) in their study on density response to the Sym-H index and total Joule heating.We find that during the 3 super storms the density reacts to E m almost immediately or within 1 to 2 h.
2. The delay time is latitude dependent.At the morning and noon sectors, the low latitude shows on average a 1-2 h larger delay than the mid latitude.This can be explained by the equatorward propagation of the solar energy impinging into the high latitude during storms.Interestingly, at the evening and night sectors, the mid latitude reacts approximately 1.5 h later than the low latitude.We can not offer an explanation for this so far.
3. The delay time is local time dependent.This point of view is only valid at mid latitude, where the morning and noon sectors have a 3 h smaller delay than the evening and night sectors.At low latitude no clear local time dependence is observed.Wang et al. (2006) provided evidence that during storms the latitude of fieldaligned currents (FAC) followed the IMF B z variation quite closely on the dayside, while the FAC latitude traces the change of the Dst index better on the nightside.The Dst is lagging behind IMF B z changes by about 2 h.Müller et al. (2009) also pointed out that the low latitude thermospheric density responds to changes in a p by 1 to 2 h earlier on the dayside than on the nightside.Although in our case the dependence is more clear at mid latitudes.
4. The delay time of orbit-averaged density is almost the same as that of low latitude.
5. The delay time shows no seasonal dependence.

The linear relation between density and geomagnetic input
In this study a linear relation between E m and storm-time thermospheric mass density is obtained as ρ = 0.5E m +ρ amb , where ρ amb implies the ambient density at quiet time.The density correlation with geomagnetic inputs are also studied in other papers.Burke et al. (2007) have investigated the thermospheric density correlation with the Dst index and the polar cap potential index, PC .Coincidently, they have chosen the same storm event as we did in Sect. 4. The orbitaveraged density is used in their study.They concluded that Dst follows thermospheric density changes through the main and early recovery phases of geomagnetic storms rather well, but in the later recovery phase the density returns much faster to quiet time values than Dst; whereas the PC estimated by the Siscoe-Hill model correlates with thermospheric mass density fairly well during all storm phases.In our study, the merging electric field, E m , also shows good correlation with density throughout the storm phases, including the recovery phase.Our result is quite consistent with the PC based result of them.Interestingly, we have to omit the saturation effect of the solar wind input.Only if we employ the merging electric field without truncation, storms of all levels can be predicted reasonably well.We have no simple explanation for this observation.Zhou et al. (2008) have utilized CHAMP accelerometer measurements during major storms to study the density dependence on total Joule heating power, Q j , and the highresolution ring current index Sym-H.They assume that the thermospheric density at 400 km is proportional to a linear combination of Q j and Sym-H.A multiple linear regression with a proper time shift of Q j and Sym-H is used.A linear empirical relation of storm-time density changes at 400 km altitude versus Sym-H index and Q j is worked out for different latitudinal segments and sunlight conditions.Prime purpose of their work is to find an empirical relation that can account for the deficiencies of the NRLMSISE-00 model, leading to a better prediction of thermospheric mass density variations during storm-time.Müller et al. (2009) also studied the dependence of CHAMP density measurements on the a m index for geomagnetically active days (A p > 15).When the effect of solar flux and seasonal variation is removed, an additive linear increase of the air density with a m is observed.The functional dependence for both day and night sides is δρ(10 −12 kg/m 3 ) = 0.012a m .The form of their function is quiet similar to our result if we rewrite Eq. ( 11) as ρ − ρ amb = 0.5E m .Both equations indicate that the effect of solar wind input is linear and additive to the background density.Opposed to that solar EUV flux causes a relative change of the mass density, i.e. regions of high density are affected stronger than low density regions (Müller et al., 2009).

The prediction capability
Our linear model based on E m can reproduce the density distribution during major magnetic storms as demonstrated by the example given in Sect. 4. The main features are successfully predicted by the model, but some under/overpredictions are found in the 6 latitude bins.This is largely due to the constant factor 0.5 in Eq. ( 11).For all latitudes and local times a = 0.5 is found as an average for all storms, but as a matter of fact, a standard deviation of 0.15 exists for both low and mid latitudes.The prediction improves if orbital averages are considered.For the application in orbit predictions these average densities are most relevant.For that reason we regard the suggested constant 0.5 as a proper general value for all the storms.Burke et al. (2007) concluded that, during geomagnetic storms, local density varies widely but orbit-averaged density evolve systematically.An example of the prediction for orbit-averaged density is given in our Fig. 8. Consistent with Burke et al. (2007), the result shows better correlation coefficient than the prediction for the segmented density.
In this study we have described and tested a procedure for predicting the thermospheric mass density changes during geomagnetic storm.In a follow-on paper we will perform detailed comparisons of our prediction with observed density changes for various types of storms at two observation heights.

Conclusion
Based on four years (2002)(2003)(2004)(2005) of CHAMP accelerometer data we have studied the dependence of mass density at low and mid latitudes on the merging electric field, E m , during major magnetic storms.A linear empirical relation between the density and the E m is obtained.We have proven that this relation can be used for prediction of density changes during storm-time.The main results of this study can be summarized as follows: 1.The properly preprocessed merging electric field tracks the density changes in all phases of magnetic storms sufficiently well.It is possible to utilize it to predict density variations during a storm.In order to account for the memory effect of the thermosphere a 3 h integrated E m value is applied.
2. A truncation of E m for strong solar wind driving, according to the saturation of the cross-polar cap potential, PC , does not reflect the variation of thermospheric density well.Only if the full E m swing is considered the simple linear relation holds also for super storms.
3. The delay time of density changes behind E m depends on various factors.Basically, it depends on geomagnetic latitude, local time and the storm intensity.No seasonal dependence is found.In this study, a constant delay of 3 h is used for low latitude, and a local time dependent delay between 0 to 4.5 h is used at mid latitudes for day and nightside, respectively.The mean delay of orbitaveraged density is almost the same as the low latitude delay.
4. The storm-induced mass density perturbation at constant altitude has been found to be an additive enhancement on top of the quiet-time density.It can be expressed as linear relation between the density and E m : ρ = 0.5E m +ρ amb , independent of local time and for all the storms during 2002 to 2005.The ambient density, ρ amb , is determined from the quiet day before the storm, and the solar flux influence on ρ amb during the storm is taken into account.
5. The linear relation can reproduce the storm-time density changes, especially the orbit-averaged density sufficiently well.For that reason the proposed model can, for example, be used for calculating the storm effect on satellite drag.
. The time interval of interest is 2002 to 2005, four years on the declining phase of Solar Cycle 23.During the investigated period the yearly-averaged solar flux, F10.7, decays from 179 in 2002 to 92 in 2005.Altogether 30 major geomagnetic storms are chosen with their minimum Dst < −100 nT.There are 3 super storms among the 30 storms: the storm in 29-31 October 2003 with Dst min = −383 nT, the storm in 20-22 November 2003 with Dst min = −422 nT and the storm on 7-11 November 2004 with Dst min = −373 nT.More features of the 30 storms can be found in Zhang et al. (2007).
Fig. 1.The relation between the transpolar potential of PC and the merging electric field E m .PC is calculated using the method described inOber and Maynard (2003).The blue curve shows that

Fig. 2 .
Fig. 2. Comparison of the various types of E m during the storm time 22-28 July 2004.The black curve shows E m directly derived from ACE measurements, The blue curve shows E m(sat) truncated to 8 mV/m, the red and green curves show the 3 h-integrated E m(sat) and E m with and without truncation, respectively.
Figure  3shows the fitting results separated by latitude classes as a function of MLT.The top two panels are for low latitudes and bottom for mid latitudes.For low latitudes the distribution of a shows no MLT dependence, 85% of the values are concentrated between 0.2 to 0.6, a mean value of 0.5 can be obtained.A few high values larger than 1.0 are shown in red color.They are from the 3 super storms mentioned at the beginning of Sect. 3. The distribution of b shows a pronounced peak in the afternoon sector and is more peaceful in the other MLT sectors.The mid latitudes show very similar outcomes as the low latitude, again the dots from the super storms are marked by red color.

Fig. 3 .
Fig. 3. Coefficients derived from the linear regression of Eq. (9) as the function of MLT.Results are derived separately for low and mid latitudes.The distribution of parameters a and b are displayed in the left and right columns, respectively.Red dots mark the a-values of 3 super storms.

Fig. 5 .Fig. 6 .
Fig.5.The scatter plots of density versus E m for the 6 latitude bins for the storm 22-28 July 2004.The top panels are for the day-side (12:00 MLT) and the bottom panels for the night-side (00:00 MLT).From left to right, the three columns are for the northern mid latitude, low latitude and southern mid latitude, respectively.A linear fitting is shown by the red straight line in each panel, with the fitting equation at the top left corner.

Fig. 7 .
Fig. 7.The density prediction results of the storm series in 22-28 July 2004 using Eq.(11).The six panels represent the 6 latitude/MLT bins.The blue curve in each panel represents the density measurements from CHAMP and the red curve denotes the prediction.The correlation coefficient, R, between the measurement and prediction, as well as the mean and standard deviation of relative error E( ) and σ ( ), are given in each panel.