Controlling of merging electric field and IMF magnitude on storm-time changes in thermospheric mass density

The controls of merging electrical field, Em, and IMF (interplanetary magnetic field) magnitude, B on the storm-time changes in upper thermospheric mass density are statistically investigated using GRACE accelerometer observations and the OMNI data of solar wind and IMF for 35 great storms during 2002–2006. It reveals the following: (1) The correlation coefficients between the air mass density changes and the parameters of Em andB are generally larger at lower latitudes than at higher latitudes, and larger in noon and midnight sectors than in dawn and dusk. (2) The most likely delay time (MLDT) of mass density changes in respect to Em is about 1.5 h (4.5 h) at high (low) latitudes, having no distinct local time dependence, while it is 6 h at middle latitudes in all the local time sectors except for noon, which is longer than at low latitudes. A similar fact of longer delay time at mid-latitude is also seen for B. The MLDTs for B at various latitudes are all local time dependent distinctly with shorter delay time in noon/midnight sector and larger in dawn/dusk. Despite of widely spread of the delay time, IMF B exhibits still larger correlation coefficients with mass density changes among the interplanetary parameters. (3) The linear control factor of B on the density changes increases for largeB, in contrast to somewhat saturation trend for larger Em. (4) The influence ofB andEm on the mass densities shows different behavior for different types of storms. The influence intensity ofEm is much stronger for CIR-driven than for CME-driven storm, while it is not so distinct for B. On the local time asymmetry of the influence, both Em andB have largest influence at noon sector for CME-driven storms, while an obviously larger intensification of the influence is found in dawn/dusk sector during CIR storms, especially for parameterEm.


Introduction
Thermospheric total mass density is important not only for understanding the coupling process in the solar windmagnetosphere-ionosphere-thermosphere system, but also for predicting the atmospheric drag that is needed in precise orbit determination and tracking of low orbit satellites.During geomagnetic disturbances, huge amounts of energy of solar wind origin enter into geospace by magnetic field reconnection, resulting in a series of dramatic disturbances in the thermosphere (Fuller-Rowell et al., 1997;Prölss, 1997).It is recently reported that the enhancements of storm-time thermospheric total mass density reach 800 %, 300 % and 400 % relative to the quiet time values during the storms of 20-21 November 2003, 29-30 October 2003, and 7-9 November 2004(Liu and Lühr, 2005;Sutton et al., 2005;Lei et al., 2010).
Since the launch of low orbit satellites, for example CHAMP and GRACE, global coverage and high-quality observations of the thermospheric total mass density have been obtained, thanks to the high-precision STAR and SuperSTAR accelerometer (Reigber et al., 2002;Tapley et al., 2004).By using these observations, the investigation of thermospheric total mass density, especially the response of thermospheric Published by Copernicus Publications on behalf of the European Geosciences Union.Y. L. Zhou et al.: Merging electric field and IMF magnitude control on storm-time thermospheric density total mass density to solar wind/IMF (interplanetary magnetic field) and geomagnetic activity, has been widely studied.Burke et al. (2007) found that thermospheric total mass density observed by GRACE is roughly proportional to polar cap potentials and magnetospheric electric fields derived from interplanetary parameters with lag time of about 4 h.Guo et al. (2010) examined quantitatively the relationship between thermosphere density variation and solar wind energy input during intense geomagnetic storms by using mass density measurements at 72 • S, 0 • , and 72 • N latitude observed by CHAMP.Their statistical analysis shows that, out of the chosen solar wind parameters including coupling functions, the Borovsky parameter correlates best with the global-scale density variations.Kwak et al. (2009) studied the influences of the IMF B y and B z on observed thermospheric mass density using the high-latitude southern summer thermospheric mass density near 400 km altitude derived from accelerometer on board CHAMP.They found that the difference density distributions, which are obtained by subtracting values for zero IMF from those for nonzero IMF, vary strongly with respect to the direction of IMF.Subsequently, they systematically analyzed the thermospheric density response to changes in the IMF sector polarity (Kwak et al., 2011).Their results showed that the IMF sector polarity changes influence strongly the high-latitude thermospheric density variations, especially in equinox seasons.Using CHAMP and ACE data during 2002-2005, Liu et al. (2010) analyzed the dependence of thermospheric mass density at low and midlatitudes on the merging electric field during major magnetic storms.They obtained a linear empirical relation between mass density averaged over two latitudinal segments (low latitude segment and mid-latitude one) and lag-time-integrated merging electric fields.They also suggested that the dynamic pressure may also have an influence on the storm-time density enhancement by analyzing the unusual magnetic storm of 21 January 2005.
Moreover, the recent results also indicated that the existing models, for example MSIS90 and NRLMSISE-00, not only underestimate the magnitude of storm-time mass density, but also cannot adequately reproduce the temporal and spatial distribution of it (Bruinsma et al., 2004;Sutton et al., 2005;Liu et al., 2005;Liu and Luehr, 2005;Forbes et al., 2005).In recent years, much attempt has been made to improve the capability of the existing model to predict the stormtime mass density.With CHAMP and GRACE air drag data, JB2008 model was developed, using the hourly Dst (disturbance storm time) index as the driver of storm-time density changes (Bowman et al., 2008).Liu et al. (2011) developed an orbit-averaged mass density prediction model for both altitude of CHAMP and GRACE with preconditioned merging electric field as input, which could reproduce fairly well the storm-time orbit-averaged mass density, but not efficient for different latitude in detail.With the help of CHAMP air drag data, we have previously established an empirical relation of storm-time mass density changes at about 400 km with both SYM-H index and total Joule heating power (Ma et al., 2006;Zhou et al., 2009), which could predict mass density changes at different latitude every 5 • for dayside and nightside separately.Using the empirical relation, the prediction quality of the NRLMSISE-00 model during storm-time is improved greatly.However, we also find that the calculation and prediction of the total Joule heating power from solar wind and IMF parameters are not so convenient in practice.Thus, we try to find the direct connection of storm-time mass density changes with solar wind and IMF parameters.
In order to find effective control parameters in practice for predicting storm-time changes in thermospheric mass densities from near-Earth solar wind/IMF data, we have firstly made analyses on the relationships between the mass density and various interplanetary parameters such as Akasofu coupling function, solar wind dynamic pressure and electric field, various IMF components and so on.It is recognized that, besides the SYM-H index, the merging electric field E m and IMF magnitude B are much better correlated with storm-time mass density changes among the interplanetary parameters cited above.In the present study, the controls of the selected parameters of E m and B on the storm-time mass density changes in the upper thermosphere are investigated by using GRACE accelerometer observations, emphasizing on the varying of the control factors and delay times versus latitude and local time mainly for mid-and low latitude thermosphere.

The air mass density data
The total mass density data used in this study are derived from the accelerometer measurements onboard the Gravity Recovery and Climate Experiment (GRACE) satellites by an ESA-authorized group in Delft Institute for Earth-Oriented Space Research (Doornbos et al., 2009).The GRACE mission (Tapley et al., 2004) is jointly implemented by NASA and the German Aerospace Center (DLR), and its objective is to map the global gravity field.The GRACE mission consists of two identical satellites separated from each other by approximately 220 km along track.The twin GRACE satellites were launched on 17 March 2002 into a near-circular polar orbit with an initial altitude of about 525 km.It has decayed down to 449 km in January 2012.The orbit inclination is 89 • , covering all the local times every 160 days.The satellite mass at launch is 487.2 kg.Both satellites are equipped with a three-axis SuperSTAR accelerometer with a resolution of 10 −10 m s −2 to observe the non-gravitational forces, which was developed by Centre National d'Etudes Spatiales, CNES, France.Although the GRACE mission is designed to improve our understanding of Earth's gravity field, the highly sensitive accelerometers on each satellite provide valuable data that can be analyzed to obtain thermospheric

Selection of storm events
For the present study, the air mass density data from

Solar wind and IMF data
In this paper, the solar wind and IMF data are obtained from OMNI database (ftp://omniweb.gsfc.nasa.gov/omni/high res omni/), which have been lagged to the nose of the Earth's bow shock from original solar wind data onboard multiple spacecrafts of ACE, WIND, Geotail and so on.The merging electric field can be written as (Kan and Lee, 1979) In Eq. (1) v sw is solar wind speed, B t = B 2 y + B 2 z is the magnitude of IMF component in the yz plane in GSM coordinates, and θ is the angle between the z-direction and the projection of the IMF vector on the yz plane in GSM coordinates.The magnitude of the IMF vector in GSM coordinate system is denoted by B, and B = B 2 x + B 2 y + B 2 z .All the OMNI data used in this study have the same temporal resolution of 1 min.
Figure 2 gives two samples to show the temporal variation of GRACE-A observed total mass density, E m and B during 19-22 November 2003 and6-11 November 2004.It can be seen from Fig. 2 that the GRACE-A observed mass density is closely correlated to E m and B with delay time of several hours.

Calculation of storm-time mass density changes
The storm-time changes in mass density are calculated as the deviations of GRACE-ACC derived mass density from a reference that characterizes the quiet-time temporal and spatial distribution of the thermospheric mass density.This reference baseline is obtained with NRLMSISE-00 model by using a fixed daily Ap index at a low level as the model input while the 3-hourly geomagnetic activity index ap for current time and 57 h prior to current time are switched off.The fixed daily Ap is selected artificially as 18 nT, which corresponds to a Kp level of 3 1991), the high limit of an un-disturbed geomagnetic condition.When calculating the mass density reference, daily solar flux F10.7 for previous day and its moving average over a window of 81 days centered on the current day are normally used as the model input parameters to account for the influence of the solar radiance.By using such a reference baseline, the quiet time variation of mass density with solar activity, latitude and local time etc. could be roughly removed and the density changes induced mainly by great magnetic storm would then be obtained.About the fixed daily Ap, other lowlevel values (say 9 nT, etc.) rather than 18 nT can also be used for the present study.The specific values of the low-level daily Ap have negligible influences on the results of correlation analyses of mass density changes with interplanetary parameters, excluding the intercept of linear fitting that is not of concern in this paper.Before making statistical analysis, the storm-time density changes are sorted into grids of latitude by local time (LT).The geographic latitude ranging from 87.5 • S to 87.5 • N is divided into 35 groups with an interval of 5 • .The local time is divided into four groups of midnight, dawn, noon and dusk sectors, centered separately at 00:00, 06:00, 12:00 and 18:00 LT.Thus there are 140 grids globally.For each orbit the storm-time changes calculated by the method cited in the previous paragraph are averaged over each grid.On account of the GRACE satellite circling the Earth about 16 times every day, a time series of storm time mass density changes with time resolution of about 95 min (i.e., the Kepler period of the GRACE satellite orbit) can be obtained in each grid for each of the selected 35 storm events.
Figure 3 shows some examples of the time series of storm time mass density changes for different local time grids at equatorial latitude during two specified superstorm events, along with the corresponding time series of interplanetary parameters of E m and B. In Fig. 3 both E m and B have been moving averaged within a time window of one orbit period of the GRACE satellite.

Comparison of correlations between density changes and various interplanetary parameters
In order to find the effective and practical control parameters composed by solar wind and IMF observation data for predicting storm-time changes in thermospheric mass densities, we firstly examined the relationships between the mass density changes and various interplanetary parameters including merging electric field, various IMF components, Akasofu coupling function, solar wind speed and plasma density, solar wind dynamic pressure, and so on.
When making cross-correlation analyses, the length of the time series of density changes and solar wind parameters should properly be determined.Firstly, the start and end time of the time series of mass density changes are fixed by looking around their excursions, referring to the development progress of the magnetic storm represented by index of Dst or SYM-H.The moment when the density is just about to increase significantly is taken as the start time.As for the ending time of the time series, it has somewhat multiform relating to the magnetic storm profile and type.In general, the time when the mass density changes recovery to the quiettime level before the storm (usually for single and multi-main phases of larger CME storms) or approaching but more and   calculated using Ap = 18 nT, the time series calculated from Ap = 9 nT are also shown in Fig. 4 for a comparison that convinces that the selection of 18 nT induces few effects on the correlation analysis results in the present study.
According to the fixed time series of the mass density changes, the time series of various interplanetary parameters are then determined by taking their lengths equal to that of the mass density changes and the start times leading the start times of mass density changes from 12 to −4 orbits step by step.
The cross-correlation between the storm-time mass density change and the specific interplanetary parameter (E m or IMF B) is calculated in each grid of latitude by local time for each storm.In total 17 linear correlation coefficients corresponding to different leading times are calculated one by one.Out of the calculated coefficients, the maximum one and its corresponding leading time of interplanetary parameters relative to density changes are taken in most cases as the resultant correlation coefficient and the delay time.Positive delay time means that the mass density changes lag behind interplanetary parameters.Meanwhile, negative delay time implies an opposite mean.In case there appear more than one maximum, the physically reasonable one of delay time is chosen.
Figure 5 presents the statistical distribution of the correlation coefficients for the parameters cited above at mid-and low latitudes, taking all the selected 35 storms into account in statistics.It can be seen from Fig. 5 that density changes generally have good correlation with a group of parameters such as merging electric field E m , Akasofu coupling function ε, IMF magnitude B, and magnitude of IMF component in the GSM yz plane B t .Among them, the merging electric field shows nearly 90 % probability for correlation coefficients larger than 0.6 and nearly 50 % probability for those exceeding 0.8, ranking the best one with respect to the correlation coefficient with mass density changes.The correlations between the density changes and the interplanetary parameters of B, B t and ε are somewhat comparable, showing coefficients exceeding 0.6 with nearly 86 %, 84 % and 85 % probability, respectively.However, the parameter of Akasofu coupling function, ε, shows much less probability than B for coefficients larger than 0.8.Parameters of B, B t and ε exhibit respectively about 48 %, 44 % and 45 % probability for coefficients exceeding 0.8.The main purpose of this study is to find such interplanetary parameters that would serve as effective and practical controller to predict storm-time changes in mass density.Considering the advantage of being more easily available in practice for IMF magnitude B, we chose it also as a candidate control parameter in addition to the merging electrical field, which has the best correlation with mass density changes.So we focus on the relationships between the mass density changes with E m and B in the present study.

Control of E m and B on the density changes
In this section, we investigate in detail the controlling of E m and B on the storm-time mass density changes.Figure 7 gives the distribution of correlation coefficients of the storm-time density changes with E m and B in different LT sectors for all the 35 storms.It is indicated that the most likely correlation coefficient (abbreviated to MLCC) is 0.8 in all four LT sectors.The probability that the correlation coefficient exceeds 0.8 is greater in noon and midnight than in dawn and dusk sectors.

LT and latitude dependence of correlation degree
Figure 8 shows the distribution of correlation coefficients of the storm-time density changes with E m and B in dif- ferent latitude regions for all the 35 storms.The high latitude, mid-latitude and low latitude are for the geographic latitude regions of ±87.5 • ∼ ±52.5 • , ±52.5 • ∼ ±27.5 • and ±27.5 • ∼ 0 • , separately.We can clearly see from Fig. 8 that the MLCC (most likely correlation coefficient) for E m and B is 0.8 at middle and low latitudes.The MLCC is 0.7 for E m and 0.8 for B at high latitudes.At low and middle latitudes, the probability that the correlation coefficients equal and exceed 0.8 is much greater than at high latitudes.
In short, the correlation degree characterized by MLCC of density changes with E m and B and the probability for correlation coefficients equal and exceeding 0.8 are larger at low and middle latitudes than at high latitudes, and larger in noon and midnight sectors than in dawn and dusk ones.Figure 10 gives the statistic distribution of delay times of the storm-time density changes in respect to E m at different (low, middle, and high) latitudes in different 4 LT sectors for all the 35 storms studied.It shows that the storm-time mass density changes almost absolutely lag behind merging electric field, except for at high latitudes where a few zero delays occurred that may not imply real zero delay but be due to low temporal resolution of one orbit number of 95 min.As a whole, the delay times that have occurrence rates larger than 5 % range from 0-7 orbits (about 0-10.5 h).The maximum occurrence rate as shown in each histogram composing Fig. 10 is thereafter called the most likely delay time (abbreviated to MLDT).When the second largest occurrence rate adjoins to and is comparable with MLDT, say the ratio to MLDT being larger than 80 %, we will consider the weighted mean of the first two largest occurrence rates and name it as weighted MLDT.

Delay time for E m
At low latitudes, the delay times seem not so remarkably dependent on local time, having the same MLDT about 4.5 h (3 orbits) in 4 various sectors.In view of weighted MLDT, the delay time measured by orbit number is 3.5 in midnight, 2.5 in noon/dawn sector, and a midway of 3 in dusk.Besides, the distributions of delay time in the noon/midnight sectors are more concentrated, while they are more dispersive in dawn/dusk.
At mid-latitudes, the delay time distributions have the same MLDT of 4 orbits in all the local time sectors except for noon.For the dawn sector, we can see a slightly shorter weighted MLDT of 3.5 orbits.On the other hand, it has much shorter MLDT of 1 orbit in noon sector with weighted MLDT about 2 orbits.
At high-latitudes, the delay time distributions have the same MLDT of 1 orbit in all four sectors without exception.However, except for noon, the distributions show two discrete peaks in all the other 3 sectors.One peak is at the most likely delay time (MLDT) of 1 orbit, another one at delay time of 4 orbits.Such a structured pattern of the delay time distribution may be attributed to the storm-time thermospheric heating source region located within high latitudes, which we allowed a rather wide range of 52.5 • to 87.5 • .
The distribution features of delay time to E m are shown in Table 2, giving the MLDT at different local time sectors for high, middle and low latitudes.Inside the parentheses in Table 2 is the second maximum of delay time in each specified histogram that is comparable with or equal to MLDT.
If we consider the latitude dependence of the delay time containing all 4 local time sectors and the local time dependence for all the latitudes, it resulted in the following distributions as shown in Fig. 11a and b.It can be seen from Fig. 11a that the most likely delay time (MLDT) is about 4.5 h (3 orbits) at low and middle latitudes when making statistics for all local time sectors, while it takes 1 orbit of MLDT for high latitudes.On the other hand, the delay times Low latitudes Middle latitudes High latitudes noon 4.5 h (3.0 h) 1.5 h (4.5 h) 1.5 h midnight 4.5 h (6.0 h) 6.0 h 1.5 h (6.0 h) dawn 4.5 h (3.0 h) 6.0 h (4.5 h) 1.5 h (6.0 h) dusk 4.5 h 6.0 h 1.5 h (6.0 h) Table 3.The MLDT of the storm-time density changes with B at different local time sectors for high, middle and low latitudes.
Low latitudes Middle latitudes High latitudes noon 4.5 h 4.5 h (1.5 h) 1.5 h (7.5 h) midnight 3.0 h 6.0 h (4.5 h) 6.0 h (0-1.5 h) dawn 7.5 h (6.0 h) 6.0 h 7.5 h (0 h) dusk 6.0 h 7.5 h (4.5 h) 7.5 h (0 h) have a shortest one of 1 orbit (95 min) at noon sector, while longest delay of 4 orbits (about 6 h) at midnight and dusk sectors.In general, the delay time of density changes behind E m characterized by MLDT (the most likely delay time) is longer at low and middle latitudes than at high latitudes, and shorter in noon than in other sectors.Figure 13 gives the statistic distribution of delay times of the storm-time density changes in respect to IMF B at different latitudes in different LT sectors for all the 35 storms.As a whole, the delay times of mass density changes behind B spread over a quite large range from −4 to 12 or more orbits, showing much more dispersion than E m .

Delay time for IMF B
The statistic distributions of delay times of storm-time density changes with B for high latitudes (the right column in Fig. 13) usually have two peaks locating near zero or 5 orbits, which is somewhat similar to E m delay feature at high latitudes.At low latitudes the most likely delay time (MLDT) is obviously shorter in noon (3 orbits)/midnight (2 orbits) sectors than in dawn (5 orbits)/dusk (4 orbits) sectors, as shown in Table 3.This is also true at mid-and high latitudes though not so distinctly as at low latitudes.In the parentheses in Table 3 is the second maximum of delay time that is comparable with or equal to MLDT.It can be seen from Table 3 that except for dawn sector the MLDTs at middle latitudes are larger than or comparable to those at low latitudes in various sectors.We should recall that, as stated in Sect.3.2.2, the most likely delay times of mass density to E m at mid-latitude are also longer than at low latitudes, but except for noon sector rather than dawn sector here.
When considering the latitude dependence of the delay time mass density changes to IMF, containing all 4 local time sectors and the local time dependence for all the latitudes, it resulted in the following distributions as shown in Fig. 14a  and b.It can be seen from Fig. 14a that the MLDT is about 4.5 h (3 orbits) at low and middle latitudes, while it takes zero orbit of MLDT for high latitudes.It should be noted that there are two discrete comparable peaks of delay time at high latitude: one is the mentioned MLDT of zero; other one is at 5 orbits (∼7.5 h).As for the local time dependence, the most likely delay time at noon sector is at 3 orbit (∼4.5 h) with a comparable second largest occurrence rate at 2 orbit (∼3 h), while it is longer in other 3 local time sectors, being respectively 5 orbits (∼7.5 h) for dusk sector and 4 orbits (∼6 h) for midnight and dawn.So, we can say in general the delay time of density changes behind IMF B characterized by MLDT (the most likely delay time) is longer at low and middle latitudes than at high latitudes, and shorter in noon than in other sectors.

Linear control factors
By the aforementioned cross correlation analysis, it is found that the storm-time mass density changes have good correlation with E m and B. Now we examine the linear control of E m and B on the mass density changes.A linear relation, ρ = ax +b, is used , where ρ is the storm-time mass density changes, x is the parameters E m or B, which are timeshifted properly according to the delay time described in the previous section, and a and b are respectively the linear control factor and intercept.By means of linear regression fitting, the control factor and intercept can be obtained.When making the linear fitting, the E m data samples are divided into two groups: one is larger than 8 mV m −1 and the other one less than it.For parameter of B, they are similarly classified into two groups according to being or not exceeding 20 nT.
Figure 15 shows the scatter point plots of the storm-time mass density changes versus the time-shifted E m and B in different LT sectors at low and middle latitudes for the 35 storms studied.The fitted regression coefficients, a and b, are given in the figure.It is interesting that the linear control factors for B on the density changes are larger for larger B (>20 nT) in all four LT sectors.In contrast, the control factors for E m on the mass density changes are smaller for larger E m (>8 mV m −1 ), which has been found to be a suitable truncation value of E m relating to saturation of the polar cap potential (Liu et al., 2010;Ober and Maynard, 2003).

Nonlinear control factors
In the following context we focus on the influence of E m and B on the mass density changes for different storm types.On account of the fact that the density variations are not exactly linearly correlated with both parameters according to the analysis results cited above, a nonlinear relation, ln ρ = a • x + b, is selected.Here ln ρ is the natural logarithm of storm-time mass density changes, a and b are the fitting coefficients which can be obtained by means of linear regression, and x is the parameter of E m or B.
Figures 16 and 17 present the scatter point plots of the natural logarithm of storm-time mass density changes versus the time-shifted E m and B in different local time sectors at low and middle latitudes for 26 CME and 9 CIR events under consideration, respectively.Table 4 gives the fitted factor of a at four different local time sectors for different storm types.
Here, the unit of the factor a is (g cm −3 )/(mV m −1 ) for parameter of E m and it is (g cm −3 )/(nT) for parameter of B. It reveals that the influence of B and E m on the storm-time mass densities shows different behavior for different types of storms.Firstly, the magnitudes of the influence factors of E m for CIR storm are larger than that for CME storm, showing an increase by 2.7 times from 0.0214 for CME storm to 0.0576 for CIR storm at dawn sector as shown in Table 4.For IMF B, however, the magnitude difference in the influence factors between CME and CIR storm is not so distinct as E m .Secondly, the the local time dependencies of the influence factors show different profiles during the two different types of storms.For CME-driven storms both E m and B have largest influence (characterized by linear fitting factor of natural logarithm of storm-time mass density changes, ln ρ, versus E m or B) at noon sector and smallest one at dawn (for E m ) or midnight (for IMF B) sector, while for CIR-driven storms, E m and B have larger influence at dawn sector and smaller at midnight sector.In terms of the ratio of influence factor for CIR over CME listed in Table 4, it shows clearly an intensification larger in dawn/dusk sector than in noon/midnight, especially for parameter of E m .

A brief discussion
In this section a few points of new findings and unexpected results obtained in this study are briefly discussed.

Longer delay time at middle latitudes
It is very interesting to note that, except for noon sector, the delay times of mass density changes in respect to E m at low latitudes are shorter than at mid-latitudes, and similarly occur also for IMF B except for dawn sector.This phenomenon seems in contravention of high latitude origin of mass density changes and consequent propagation equatorward.In fact it was noticed by Liu et al. (2010) at CHAMP altitudes but puzzled the author to interpret.In this study, additional independent evidence drawn from GRACE data at about 490 km altitude has been found to confirm this natural phenomenon.We surmise that this may imply some additional heating or disturbance sources rather than high latitude origin at working for low and/or mid-latitude mass density changes.Considering the recognized fact is not true at noon (for E m ) or dawn (for IMF B) sector, we suppose that it may be associated with some coupling processes in the magnetosphere-ionosphere/thermosphere system occurring mainly in the night-and duskside geospace.One possible source may be energetic neutral atom (ENA) precipitation of storm-time ring current origin, caused by charge exchange between energetic RC (ring current) ions and cold atoms of geo-corona.During the storm main phase, the RC is a partial ring with strong local time asymmetry, located mainly in nighttime favoring duskside near equator.Observational evidence has been issued by DeMajistre et al. (2005) to present direct connection of ENA precipitation and enhanced thermospheric airglow emissions at mid-latitude during storm times.In some cases, intensive large-scale travelling atmospheric disturbances (TAD) propagating from high to lower latitudes (e.g., Forbes et al., 2005) may modulate the latitude profile of storm-time mass density changes, affecting the latitude distribution of the density delay time.Be-sides, the coupling between low-latitude thermosphere and ionosphere associated with prompt penetration of interplanetary electric field that has larger penetration efficiency during night may be another possible candidate.This topic in detail is beyond the scope of this paper and will be pursued in our future study.

Discrete peaks of delay time at high latitude
As cited in Sect.3.2.2, at high latitudes, the distribution histograms of mass density delay times in respect to E m have the same MLDT of 1 orbit in all four sectors without exception.However, except for noon, the distributions show two  discrete peaks in all the other three sectors.One peak is at the most likely delay time (MLDT) of 1 orbit, another one at delay time of 4 orbits.The delay times of mass density changes to IMF B at high latitude have also two discrete peaks: one is at zero orbit and other one at 4-6 orbits.Such a structured pattern of the delay time distribution may be attributed to the storm-time thermospheric heating source region located within high latitudes, which we allowed a rather wide range of 52.5 • to 87.5 • .The shorter delay time peak may probably relate to the heating source region, while the longer one may be a little far from the source region at highest latitudes as seen on the right in Fig. 9.

Uniform LT distribution of MLDT for E m
In view of the most likely delay times as shown in Table 2 and Fig. 10, the mass density changes lag behind E m by the same time length when making the statistics for all storms over a specified latitude extent except for mid-latitude one.That is, the delay time that occurs the most likely at a given latitude segment is almost independent of local time, although the distribution patterns of occurrence rate vary with LT and the delay times at different LT sectors for individual storm seen by GRACE are always local time dependent.This result is different from that obtained by Liu et al. (2010, Fig. 6) at  In addition, it is also noticed that the most likely delay times to E m at GRACE altitude are systemically one orbit longer than that obtained by Liu et al. (2010) at CHAMP altitude.Due to low time resolution of an orbit period, this difference can be biased to 47.5 min.This longer delay time could completely be attributed neither to different interplanetary data from ACE used by Liu and OMNI in this paper nor to different analysis methods, which could cause Liu's result to be only 15 min in advance of that by the present work.Perhaps, it implies somewhat an upward propagation (or heat conducting) delay effect.

Delay time for IMF B
As seen from Fig. 13, the distribution of delay time for B is much more dispersive than for E m .Despite of the wide spread of the delay time, IMF B exhibits larger correlation coefficients with mass density changes among the interplanetary parameters.Why it is and what it means physically is a puzzle.We cannot interpret it properly yet.

Different influences on CME and CIR storms
The CME (coronal mass ejection) and CIR (co-rotating interaction region)-type storms have different solar and interplanetary origin.The CME storm driver includes multi-forms such as CME sheath, magnetic clouds and ejecta, while CIR storm is driven by recurring high-speed streams associated with large coronal holes (Tsurutani et al., 2006a, b, and therein).The former occur dominantly during solar maximum, and almost all the superstorms are generated by this kind of driver, while the latter at declining solar cycle phase and their intensity indicated by Dst index can rarely reach below −100 nT.Moreover, in some cases CIR storm has recovery phase lasting as long as a few weeks.Borovsky and Denton (2006) summarized systematically the differences between the two types of storms.They conclude that CMEdriven storms have denser magnetospheric plasma sheets, stronger ring currents as well as great auroras and can produce harmful solar energetic particle events and new radiation belts and so on.CIR-driven storms are of longer duration, have hotter plasma sheets and hence stronger spacecraft charging, and produce higher fluxes of relativistic electrons.They also suggest that CME-and CIR-driven storms should be studied separately when geomagnetic storms are studied.In Sect.3.2.4 of this paper we examined preliminarily the influences of E m and IMF B on the mass density changes for different types of CME-driven and CIR-driven storms, characterized by linear fitting factor of natural logarithm of storm-time mass density changes, ln ρ, versus E m or B. It reveals that the magnitudes of the influence factors for CIR storm are in general larger than that for CME storm.This seems much more remarkable for E m parameter, showing an increase by 2.7 times at dawn sector and ∼ 1.4 − −2.1 times at other 3 sectors.In addition, the local time dependencies of the influence factors show different profiles.For parameter of E m , CME-driven storms have larger influence at noon/midnight sector and smaller at dawn/dusk sector, while CIR-driven storms have largest influence at dawn sector and smallest at midnight.For parameter of IMF B, it shows somewhat similar LT-dependent trend, showing the largest influence factor at noon for CME-driven storm, while it maximized at dawn and minimized at noon/ midnight for CIR-driven storm.Recently, Chen et al. (2012) has compared the effects of the CIR-and CME-driven geomagnetic activity on thermospheric density and spacecraft orbits by using CHAMP satellite observations.They concluded that the larger changes in thermospheric density during CIR storms are caused by the longer duration of CIR-storms.In the present study, the larger influence for CIR concerns the instantaneous effect, not the accumulative effect of E m and B on the thermospheric mass density.Thus the longer duration time could not explain our results.About the stronger influence of E m on mass density changes for CIR storm, Liu et al. (2011) have found this phenomenon and attribute it to storm-type dependence of cross polar cap potential (CPCP) saturation.According to Borovsky and Denton (2006), saturation of CPCP occurs rarely for high-speed stream-driven storms but commonly for CME-driven storms.Besides this reason, we guess there may be other additional mechanism responsible for larger effects in mass density changes for CIR storm, e.g., the nonlinear Alfvén wave within the high streams proper (Tsurutani et al., 2006a).The different pat-terns in local time dependence of influence factors during CME-and CIR-driven storms enfolded in this study may also be associated with the Alfvén wave mechanism, but this is not certain yet.The exact mechanism concerning the geoeffectiveness of CME and CIR storm is in debate at this time.It is greatly worthy of our further study but beyond the scope of the present investigation.1.The merging electric field E m and IMF magnitude B are much better correlated with storm-time mass density changes among considered various interplanetary parameters.The correlation coefficients between density changes with E m and B are in general larger at lower latitudes than at higher latitudes, and larger in noon and midnight sectors than in dawn and dusk ones.
2. Usually, the storm-time mass density changes lag behind E m and B several hours, varying with latitude and local time.The following was found: (a) The most likely delay time (MLDT) for E m is about 1 (3) orbits at high (low) latitudes having no distinct local time dependence.However, it is 4 orbits at middle latitudes in all the local time sectors except for noon, which is longer than at low latitudes.
A similar fact of longer delay time of mass density at mid-latitudes than at low latitudes is also seen for parameter of IMF B. (c) At high latitudes, the distribution of delay time for both E m and B has two peaks in all local time sectors except noon for E m .One is 0-1 orbits associated with polar heating source region.The other peak is about 4 to 5 orbits relating to mass density changes at latitudes a little far from source region.
3. The linear control factors of E m and B parameters on the storm-time mass density changes are examined respectively for lower and higher levels of the two parameters.It is found that the control factors of B on the density changes are larger for larger B (>20 nT).In contrast, the control factors of E m on the mass density changes become smaller for larger E m (>8 mV m −1 ).
4. The influences of B and E m on the storm-time mass densities characterized by nonlinear control factors show different behavior for different types of storms.
The influence intensity of E m on mass density is stronger for CIR-driven than for CME-driven storms, manifested as 2.7 times for CME-driven over CIRdriven storm at dawn sector.In terms of the ratio of influence factor for CIR over CME, it is discovered that there is a larger intensification in dawn/dusk sector than in noon/midnight sector clearly for both E m and IMF B.
The results suggest that merging electric field of E m and IMF magnitude of B could be selected as candidate parameters to predict storm-time mass density changes from interplanetary parameters if the delay time of storm-time mass density changes behind E m and B can be specified reasonably.
Parameter E m is expected to be of better prediction performance at the expense of all the solar wind speed and IMF vector data available, while IMF B would be a good selection in practice when only IMF magnitude data are at hand.What practical factors and how they influence the delay time of mass density changes behind the solar wind/IMF parameters is worthy of further study.
Fig. 1.(a) The local time and latitude coverage of the GRACE-A satellite orbits for the considered storms; (b) histogram of the local time sector (noon/midnight or dawn/dusk) for CME and CIR.

Fig. 2 .
Fig. 2. Temporal variations of GRACE-A observed total mass density, E m and B during 19-22 November 2003 and 6-11 November 2004.The black line is for mass density, red for E m and light blue for B.

Fig. 3 .
Fig. 3. Orbit (time) series of storm-time mass density changes, E m and B around equator at midnight and noon sector.(a) During 19-22 November 2003; (b) during 6-11 November 2004.Blue line is for mass density changes, red for E m and green for B.

Fig. 4 .
Fig. 4. Examples of specified time series of mass density changes.The vertical dotted lines indicate the beginning and ending time (orbit number) of the time series.In addition to the time series calculated using Ap = 18 nT (pink color for dusk or noon sector, green for dawn or midnight), the time series calculated from Ap = 9 nT (red color for dusk or noon sector, blue for dawn or midnight) are also given for a comparison.

Fig. 5 .
Fig. 5. Distribution of correlation coefficient of the storm-time density changes with various interplanetary parameters at mid-and low latitudes for all local time sectors taking all the 35 storms into account in statistics.

Fig. 6 .
Fig. 6.The latitudinal variation of correlation coefficients of mass density changes with E m (left) and B (right) for different LT sectors during magnetic storm of 19-22 November 2003 (top) and 6-12 November 2004 (bottom).

Fig. 7 .
Fig. 7. Distribution of correlation coefficients of the storm-time density changes with E m (left) and B (right) in different LT sectors for the 35 storms studied.

Figure 6
Figure 6 presents some examples of the latitude dependence of correlation coefficients of mass density changes with E m (left) and B (right) for different LT sectors during the magnetic storm of 19-22 November 2003 and 6-12 November 2004.It can be seen from Fig. 6 that the mass density changes are correlated closely with E m and B during the two storm events, showing correlation coefficients above 0.7 at almost all latitudes (except for 2004 storm event at high latitudes in the Southern Hemisphere).Figure7gives the distribution of correlation coefficients of the storm-time density changes with E m and B in different LT sectors for all the 35 storms.It is indicated that the most likely correlation coefficient (abbreviated to MLCC) is 0.8 in all four LT sectors.The probability that the correlation coefficient exceeds 0.8 is greater in noon and midnight than in dawn and dusk sectors.Figure8shows the distribution of correlation coefficients of the storm-time density changes with E m and B in dif-

Figure 9
Figure 9 presents some examples of the latitude dependence of delay time of mass density changes behind E m for different LT sectors during the magnetic storm of 19-22 November 2003 and 22-28 July 2004.For an individual storm event, the delay time varies with latitude and local time sectors.Figure10gives the statistic distribution of delay times of the storm-time density changes in respect to E m at different (low, middle, and high) latitudes in different 4 LT sectors for all the 35 storms studied.It shows that the storm-time mass density changes almost absolutely lag behind merging electric field, except for at high latitudes where a few zero delays occurred that may not imply real zero delay but be due to low temporal resolution of one orbit number of 95 min.As a whole, the delay times that have occurrence rates larger than 5 % range from 0-7 orbits (about 0-10.5 h).

Fig. 9 .
Fig. 9. Latitudinal variation of delay time of mass density changes with E m for different LT sectors during magnetic storms of 19-22 November 2003 (left) and 22-28 July 2004 (right).

Figure 12
Figure 12 shows a few examples of the latitude dependence of delay time of mass density changes in respect to IMF B for different LT sectors during magnetic storm of 19-22 November 2003 and 22-28 July 2004.Figure13gives the statistic distribution of delay times of the storm-time density changes in respect to IMF B at different latitudes in different LT sectors for all the 35 storms.As a whole, the delay times of mass density changes behind B spread over a quite large range from −4 to 12 or more orbits, showing much more dispersion than E m .The statistic distributions of delay times of storm-time density changes with B for high latitudes (the right column in Fig.13) usually have two peaks locating near zero or 5 orbits, which is somewhat similar to E m delay feature at high latitudes.At low latitudes the most likely delay time (MLDT) is obviously shorter in noon (3 orbits)/midnight (2 orbits) sectors than in dawn (5 orbits)/dusk (4 orbits) sectors, as shown in Table3.This is also true at mid-and high latitudes though not so distinctly as at low latitudes.In the parentheses in Table 3 is the second maximum of delay time that is comparable with or equal to MLDT.It can be seen from Table3that except for dawn sector the MLDTs at middle latitudes are larger than or comparable to those at low latitudes in various sectors.We should recall that, as stated in Sect.3.2.2, the most likely delay times of mass density to E m at mid-latitude are

Fig. 10 .
Fig. 10.Distribution of delay times of the storm-time density changes with E m at different latitudes in different LT sectors for the 35 storms studied.

Fig. 12 .
Fig. 12. Examples of latitudinal variation of delay time of mass density changes in respect to IMF B for different LT sectors during magnetic storms of 19-22 November 2003 (left) and 22-28 July 2004 (right).

Fig. 15 .
Fig. 15.Scatter point plot of the storm-time mass density changes versus the time shifted E m and B in different LT sectors at low and middle latitudes for the 35 storms studied.The linearly fitted lines are drawn separately over larger (>8 mV m −1 for E m and >20 nT for B) and smaller (<8 mV m −1 for E m and <20 nT for B) values of interplanetary parameters.

Fig. 16 .
Fig. 16.Scatter point plot of the natural logarithm of storm-time mass density changes versus the time shifted E m and B in 4 different LT sectors at low and middle latitudes for the 26 CME storms studied.
of GRACE accelerometer observations and the solar wind and IMF OMNI data, we have statistically investigated the relationships of merging electrical field, E m , and IMF magnitude, B, with the storm-time changes in the upper thermospheric mass density for 35 great storms during 2002-2006.The linear control factors of E m and B on the mass density changes and their algorithm are examined, along with the delay times of density changes behind E m and B. The dependences of the control factors on the latitude and local time are investigated for different storm types.The main results of this study can be concluded:

( b )
In comparison with E m , the delay time of mass density in respect to IMF B is much more dispersive and local time dependent.The MLDTs for B at various latitudes are generally shorter in noon/midnight sectors than in dawn/dusk, especially at low latitudes.Despite of the wide spread of the delay time, IMF B exhibits still larger correlation coefficients with mass density changes among the interplanetary parameters.

Table 2 .
The MLDT of the storm-time density changes in respect to E m at different local time sectors for high, middle and low latitudes.

Table 4 .
Linear control factor of E m and B on ln ρ, a E m and a B , at different local time sectors for CME and CIR storms and their ratio.aE m (g cm −3 )/(mV m −1 ) a B (g cm −3 )/(nT)