A model of high‐latitude thermospheric density

We present an empirical model of the high‐latitude air density at 450 km, derived from accelerometer measurements by the CHAllenging Minisatellite Payload and Gravity Recovery and Climate Experiment satellites during 2002–2006, which we call HANDY (High‐Latitude Atmospheric Neutral DensitY). HANDY consists of a quiet model and disturbance model. The quiet model represents the background thermospheric density for “zero geomagnetic activity” conditions. The disturbance model represents the response of the thermospheric density to solar wind forcing at high latitudes. The solar wind inputs used are the following: (1) solar wind electric field ESW, (2) interplanetary magnetic field (IMF) clock angle CSW, and (3) solar wind dynamic pressure PSW. Both quiet and disturbance models are constructed on the basis of spherical harmonic function fitting to the data. Magnetic coordinates are used for the disturbance model, while geographical coordinates are used for the quiet model. HANDY reproduces main features of the solar wind influence on the high‐latitude thermospheric density, such as the IMF By effect that produces a hemispheric asymmetry in the density distribution.


Introduction
The magnetosphere-ionosphere-thermosphere system is under significant influence of the solar wind at high latitudes. Deposition of the energy and momentum from the solar wind results in heating of the upper atmosphere, causing the thermosphere to expand. The density of the thermosphere increases at a fixed altitude as the thermal expansion brings molecular-rich air to higher levels. During geomagnetic storms, the thermospheric density can increase by several hundred percent in comparison with quiet periods [e.g., Sutton et al., 2005;Bruinsma et al., 2006], which is an obvious concern for operators of low Earth orbit satellites.
Various empirical models have been proposed to describe temperature and density variability in the thermosphere. The following are examples: Mass Spectrometer Incoherent Scatter Radar Extended (MSISE) models [Hedin, 1991;Picone et al., 2002]; Drag Temperature Model [Bruinsma et al., , 2012; Jacchia-Bowman (JB) models [Bowman et al., 2008a[Bowman et al., , 2008b; and CHAllenging Minisatellite Payload (CHAMP) empirical model [Liu et al., 2013]. Those models evaluate the air density as a function of altitude, latitude, longitude, solar time, solar and geomagnetic activities, and day of year. The models are useful not only for satellite operations (i.e., orbital tracking and prediction) but also for characterizing the nature of thermospheric variability.
The previous models have been focused at low latitudes. Low latitudes account for substantial area of the globe, and thus, it has been the primary interest for scientists and satellite operators. Besides, in situ satellite measurements have been sparse at high latitudes, as it requires the satellite to be in a high-inclination orbit. Consequently, the three-dimensional models mentioned above do not include high-latitude density features. This issue was brought to light when  compared thermospheric densities measured by the CHAllenging Minisatellite Payload (CHAMP) satellite [Reigber et al., 2002] with the MSISE-90 model [Hedin, 1991].  showed that the MSISE-90 model fails to reproduce high-latitude density structures observed by the CHAMP satellite. The discrepancy was especially evident at magnetic high latitudes around the noon sector and premidnight sector where the MSISE-90 model significantly underestimates the air density. The MSISE-90 model, like other global models of the thermosphere, is based on fitting of low-order global spherical functions to observations, which tend to smooth out relatively small structures at high latitudes.
Air drag measurements by CHAMP have revealed high-latitude thermospheric density distributions in great detail. The high-inclination near-circular orbit of CHAMP (I = 87 ∘ ) enabled pole-to-pole measurements during its operation from July 2000 through September 2010. Lühr et al. [2004] found a region of enhanced 10.1002/2015JA021371 thermospheric density at high latitudes around the dayside cusp. It was later realized that the amplitude of this density bulge is largely controlled by the magnitude of the solar wind electric field [e.g., Rentz and Lühr, 2008;Kervalishvili and Lühr, 2013] and the orientation of the interplanetary magnetic field (IMF) [Kwak et al., 2009;Yamazaki et al., 2015]. Numerical studies have shown that the cusp region density enhancement could arise from local heating due to soft electron precipitation and Poynting flux [Crowley et al., 2010;Li et al., 2011;Zhang et al., 2012;Deng et al., 2013].
Extensive high-latitude measurements by CHAMP also revealed a region of relatively enhanced thermospheric density in the nighttime sector around 22-01 magnetic local time, which is associated with substorms [Ritter et al., 2010]. The simulation study by Zhang et al. [2012] showed that soft electron precipitation increases the thermospheric density not only in the dayside cusp region but also in the premidnight region.
The Gravity Recovery and Climate Experiment (GRACE) satellite [Tapley et al., 2004] has also collected a large quantity of thermospheric density data in a high-inclination (I = 89 ∘ ) near-circular orbit since March 2002. The GRACE altitude was higher than CHAMP's by some 100 km during the period when both GRACE and CHAMP were operative. Lei et al. [2012] explored density data from the two satellites, normalizing those data into a single fixed height at 400 km. This way, they were able to improve the horizontal spatial coverage of the data, which facilitates the separation between spatial and temporal variability. In the present study, we also use thermospheric density data from CHAMP and GRACE. The main focus of this paper is to present the first model of the high-latitude thermospheric density at 450 km that uses solar wind parameters as inputs. The model results are discussed in comparison with previous results in the literature.

Data and Model Construction
We analyze thermospheric total mass densities derived from the CHAMP and GRACE accelerometer data. The density retrieval procedures and error evaluations were detailed in Sutton [2008]. Briefly, the overall accuracy of the data is ∼11%. Errors mainly come from the drag coefficient, neglect of winds, and instrument calibration. The error from neglecting winds is typically 2-10%, but at high latitudes during storms, the error can be up to ∼25%. The error due to the precision of the accelerometer is less than 1%.
We use the data for 4 years from August 2002 to July 2006, when both satellites were operating. During this period the CHAMP satellite gradually descended from approximately 420 km to 350 km, and the GRACE altitude was around 500-480 km, as shown in Figure 1 (top). All the data from the two satellites at different altitudes were normalized at a single height of 450 km, which is approximately in the middle of the two satellites. The MSISE-00 model was used for the altitude corrections. As discussed in previous studies [Forbes et al., , 2011Lei et al., 2012], it is necessary to intercalibrate CHAMP and GRACE densities before the two data sets are combined. Following these studies, we first computed the mean ratio between CHAMP measurements and MSISE-00 model, as well as the mean ratio between GRACE measurements and MSISE-00 model. A correction factor was, then, determined so as to adjust one of the obtained ratios to the other. An analysis indicated the following relationship: GRACE/MSISE = 0.897 × CHAMP/MSISE. The period of our data analysis was limited to 4 years (August 2002to July 2006. This choice was made as a compromise between including as much data as possible and eliminating measurements in which the CHAMP altitude was too low compared to the target height of 450 km. Besides, studies showed that the MSISE-00 model was not accurate during the extreme solar minimum of 2008-2009 owing to an unexpectedly large amount of helium [Thayer et al., 2012;Liu et al., 2014a], which would affect our height corrections based on MSISE-00. The period we investigate (i.e., 2002-2006) is in the declining phase of solar cycle 23. Figure 1 (middle) shows monthly mean values of the M 10.7 index, which is the Mg II core-to-wing ratio [Viereck et al., 2001] scaled to the F 10.7 index. We derived M 10.7 based on a linear least squares fit of daily F 10.7 to the corresponding Mg II values, following Solomon et al. [2011]. Daily values of the M 10.7 index will be used later as an input parameter of our models. Monthly mean ap index is plotted in Figure 1 (bottom), representing geomagnetic activity during 2002-2006. As is known, the declining phase of a solar cycle provides relatively high geomagnetic activity [e.g., Lockwood et al., 1999], which is favorable for the purpose of our investigation.
Our model, HANDY (High-latitude Atmospheric Neutral Density), was constructed in terms of the logarithm of the air density (not the absolute density), as is the case for most empirical models of the thermosphere. The analysis in terms of log density (or relative density changes) is more appropriate for empirical modeling than the analysis in terms of absolute density. This is primarily owing to the fact that variance of absolute density varies significantly through a solar cycle. The variance of absolute density increases with increasing solar activity as the background density undergoes an order-of-magnitude increase from solar minimum to solar maximum. Since ordinary least squares fitting assumes uniform variance, if absolute density were chosen to use for the fitting, the solar maximum data would receive much greater weight in the fit than the solar minimum data. This problem can be avoided by using log density, as its variance is much more uniform through a solar cycle. Further discussion on the advantage of log density over absolute density in climatological studies can be found in the paper by Emmert and Picone [2010, section 2.2].
The construction of the model involves two steps. First, we created a global quiet model, which represents the background density for "zero geomagnetic activity" conditions. The quiet model was then subtracted from the original data. The residuals were analyzed focusing on high latitudes in order to construct a disturbance model, which represents the response of the high-latitude thermosphere density to solar wind forcing. Letting M represent the model density, where log q is the quiet model and log d is the disturbance model. It is obvious, from equation (1), that the construction of the disturbance model is, in effect, to model the relative density perturbation from the quiet-time background density, i.e., d ∼ ∕ q .
We use geographical coordinates for the quiet model and magnetic coordinates for the disturbance model. The quiet-time background density is primarily controlled by solar heating, and thus, geographical coordinates are suitable to describe the density distribution and variability. On the other hand, density perturbations during geomagnetically disturbed periods arise mainly from high-latitude Joule heating and other processes

10.1002/2015JA021371
of magnetosphere-ionosphere-thermosphere coupling, which are better organized in magnetic coordinates than in geographical coordinates. Specifically, we use magnetic apex coordinates described by Richmond [1995] and .

Quiet Model
The quiet model was constructed in the following way. First, we collected all the CHAMP and GRACE data when the average of the ap index during the last 24 h is less than 9 nT (equivalently, Kp ≤ 2). Next, these quiet-day measurements were evaluated as a function of geographic latitude, longitude, universal time, local time, day of year, solar activity, and geomagnetic activity. The equations used for the quiet model are basically the same as that of the MSISE-00 model, which were detailed in Hedin [1983,1987]. The differences are that we use the M 10.7 index instead of the F 10.7 index and that the dependence of log densities on the daily ap index is simplified considering only the linear dependence. (The exact formula for the quiet model is shown in the supporting information.) Previous studies found that the M 10.7 index is able to represent the solar activity influence on the thermosphere, often better than F 10.7 [e.g., Guo et al., 2007]. The dependence of the density on geomagnetic activity is fairly small for ap < 9 nT, which can be assumed to be linear [e.g., Vickers et al., 2013].
Fitting coefficients were determined using the least squares technique. The quiet model of HANDY is then obtained by normalizing the reconstructed densities to ap = 0 nT, which gives density estimates for the zero geomagnetic activity condition. It is noted that the zero geomagnetic activity condition defined above does not mean that the energy input from the solar wind to the upper atmosphere is actually zero. Generally, the high-latitude atmosphere is subject to solar wind disturbances even when ap = 0 nT. The zero geomagnetic activity condition merely gives an objective reference level of geomagnetic activity for the density perturbation.
As is clear from the analysis procedure, our quiet model is designed to be consistent with the MSISE-00 model at ap = 0 nT. We demonstrate in Figure 2 that seasonal and solar activity variations in the density are in good agreement between our quiet model and the MSISE-00 model (at ap = 0 nT). The plots are limited to poleward of ±60 ∘ latitude, as we are interested in high latitudes only. Both models show a strong dependence on solar activity. The enhanced solar EUV heating during solar maximum increases the thermospheric temperature, which leads to an increase of the density at a fixed altitude. The density tends to be greater in the summer hemisphere than in the winter hemisphere, which is due to higher solar insolation in summer. High-latitude densities are comparable between the Northern and Southern Hemispheres under the equinoctial condition. The effect of solar insolation also explains the daily variation of the density, causing greater densities during daytime than nighttime.
Figure 3 (top row) compares the quiet model with the original measurements at high magnetic latitudes (poleward of ±60 ∘ magnetic latitude). As expected, the quiet model accounts a substantial part of the density but slightly underestimates the measurements as it does not include contributions of geomagnetic activity. A study by Emmert and Picone [2010] showed that more than 90% of global density perturbations can be attributed to the effect of solar activity and solar insolation (i.e., solar time and season).
It is known that the composition of the high-latitude thermosphere varies significantly with solar activity and season. According to MSISE-00, the neutral population at 450 km is dominated by atomic oxygen (O), which accounts for ∼90% of the total mass on average. The other 10% is mainly by helium (He) and molecular nitrogen (N 2 ). The contribution of He can be over 50% in the winter hemisphere under low solar and geomagnetic activity conditions owing to the existence of the winter helium bulge [see Liu et al., 2014b, and references therein]. Different constituents have different response to geomagnetic activity. Thus, although we do not have composition data for the present study, the change in the composition is expected to add complexity to the thermospheric response to solar wind forcing described in the following section. (See also Thayer et al. [2012] for the composition effect on the total mass density at CHAMP and GRACE altitudes.)

Disturbance Model
We analyzed the residuals between the measurements and quiet model estimates for the disturbance model. The model was constructed separately for the Northern Hemisphere and Southern Hemisphere in the apex magnetic coordinates. Our approach is somewhat similar to that of Weimer [1995,1996], who constructed empirical models of the high-latitude potential electric field on the basis of a spherical cap harmonic analysis. Following Weimer [1995Weimer [ , 1996, the lower boundaries were set to ±45 ∘ magnetic apex latitudes, and spherical harmonic functions were fitted to the data only poleward of the boundaries. The model formulation is as follows: where MLAT is apex magnetic latitude in degrees, and MLT is magnetic local time in hours. A lm and B lm are functions of solar wind parameters, which will be discussed later. P lm is the Schmidt normalized associated Legendre function with order l and degree m, expressed as follows: where s = 1 for m = 0, s = Thus, varies from 1 at the North Pole to −1 at the lower boundary of 45 ∘ in the Northern Hemisphere, and from 1 at the South Pole to −1 at the lower boundary of at −45 ∘ in the Southern Hemisphere. The spherical harmonic expansion was truncated at l = 6 and m = 3. The inclusion of higher-order/higher-degree terms does not add any new steady structure but tends to exaggerate noise. It is noted that the quiet model also uses spherical harmonics truncated at l = 6 and m = 3. An important difference between the quiet model and disturbance model is that the quiet model is based on a global fitting while the disturbance model is based on a regional fitting (i.e., poleward of ±45 ∘ magnetic apex latitudes). By limiting the area of fitting, we ensure that the disturbance model can properly capture density structures in the high-latitude regions.
The A lm and B lm terms in equation (2) are functions of the following: solar wind electric field magnitude E SW , solar wind dynamic pressure P SW , IMF clock angle C SW , solar activity M, and day of year DoY. For the solar activity parameter M, we use the daily value of the M 10.7 index 1 day prior to the observed density. For the solar wind parameters, we used 1 min OMNI solar wind data adjusted to the Earth's bow shock nose. An additional time delay of 15 min was added to account for the magnetosphere-ionosphere distance [Vennerstrøm et al., 2002;Rentz and Lühr, 2008]. We tested various combinations of averaging window widths and time lags using the CHAMP and GRACE data for the polar regions above ±60 ∘ magnetic latitude. The best fit was obtained when the solar wind data are averaged for the past 10 h from the present time. Any time lag from the 15 min adjustment did not improve the fitting. The suitable time lag may be different at other latitudes, as the time lag for the thermospheric density response to solar wind disturbances is generally larger at lower latitudes [Sutton et al., 2009].
Mathematical expressions for A lm and B lm were determined on a trial and error basis. So far, the best model performance was obtained using the following: where C lm is either A lm or B lm . We found that the use of √ E SW and √ P SW improves the fitting, compared to the results with E SW and P SW . As will be shown later, log d does not linearly increase with E SW and P SW .
The coefficients lm , lm , lm , lm , lm , lm , lm , lm , and lm were determined for A lm and B lm , separately in the Northern Hemisphere and the Southern Hemisphere, using the least squares method. The model coefficients are included in the supporting information, with corresponding 1 errors evaluated using the bootstrap method. Figure 3 (middle row) shows good agreement between the final model estimates (i.e., the quiet model plus disturbance model) and CHAMP/GRACE measurements at high magnetic latitudes. The goodness of fit, evaluated as the square of the correlation coefficient, is 0.90 in the Northern Hemisphere and 0.91 in the Southern Hemisphere. The root-mean-square error, defined here as √ ( log − log M ) 2 , is plotted in Figure 3 (bottom row). The results indicate that the average root-mean-square error for the HANDY model is 15.0% in the Northern Hemisphere and 15.6% in the Southern Hemisphere during the period investigated. As can be seen in Figure 3 (bottom), the performance of HANDY does not significantly depend on solar activity or season. On the other hand, the root-mean-square error for the MSISE-00 model shows seasonal variations, indicating poor performance during local summer. The average error for MSISE-00 is 21.5% and 24.1% in the Northern Hemisphere and Southern Hemisphere, respectively. The larger error for MSISE-00 than HANDY is partly due to the fact that MSISE-00 does not include CHAMP/GRACE data in fitting while HANDY does. Thus, we also calculated the root-mean-square errors for the JB2008 model [Bowman et al., 2008b], which is constrained by CHAMP andGRACE data (2001-2005). The results are shown in Figure 3 (bottom). The average error for JB2008 is 19.0% in the Northern Hemisphere and 19.4% in the Southern Hemisphere. Previous studies have reported that the JB2008 model generally performs slightly better than the MSISE-00 model [Bowman et al., 2008b;Shim et al., 2012].  Figure 4 shows the high-latitude density response to geomagnetic/solar wind disturbances in the Northern Hemisphere (Figure 4, top) and in the Southern Hemisphere (Figure 4, bottom) for the September equinox condition. Figure 4 (left column) is for the MSISE-00 model, and Figure 4 (right column) is for the HANDY model. For the MSISE-00 results, the color indicates the density ratio between ap = 27 nT (equivalently Kp = 4) and ap = 0 nT for a moderate solar activity condition (F 10.7 = 135 sfu). The HANDY results are produced for the following solar wind inputs: E SW = 3.3 mV/m, P SW = 3.1 hPa, and C SW = 180 ∘ . These E SW and P SW values correspond to geomagnetic activity of ap = 27 nT on average. It is immediately obvious in Figure 4 that there are significant discrepancies between the MSISE-00 and HANDY results. The HANDY results indicate regions of locally enhanced and reduced response around the noon sector and predawn sector, respectively, which are completely missing from the MISISE-00 results. The weak response around the predawn sector is interesting, as it is indeed where Joule heating from the closure of magnetic field-aligned currents peaks [Weimer, 2005]. Yamazaki et al. [2015] discussed that the production rate of nitric oxide, which is a strong radiative coolant, may be locally enhanced due to precipitation of hard electrons. This mechanism needs to be validated by independent measurements.

Comparison With Other Models
Another clear discrepancy is that HANDY predicts the largest increase in the relative density in the premidnight sector, while the MSISIE-00 shows it after the midnight sector. It should be noted that the MSIS-00 model uses geographical coordinates while the disturbance model of HANDY uses magnetic coordinates. However, the difference in the coordinates does not seem to explain all the differences between the MSISE-00 and HANDY results. The premidnight sector is where substorm onset is often observed [e.g., Frey et al., 2004]. Heating due to energetic particle precipitation during substorms may be a reason for the density enhancement in this region. The solar wind response in HANDY is in general agreement with other climatological studies [e.g., Rentz and Lühr, 2008;Kervalishvili and Lühr, 2013].

Dependence on Solar Wind Drivers
In order to provide insight into how the relative density d varies with the solar wind electric field magnitude and solar wind dynamic pressure, we have run the HANDY model for various E SW and P SW conditions. The calculations started with E SW = 1.7 mV/m, P SW = 1.9 hPa, C SW = 180 ∘ , M 10.7 = 120 sfu, and DoY = 264. These E SW and P SW values roughly correspond to ap = 4 nT (or Kp = 1). Then, in one of the runs, the E SW value was increased up to 15.1 mV/s without changing the other parameters. In the other run, the P SW value was increased up to 10.5 hPa while keeping the other parameters the same. These upper values for E SW and P SW correspond to approximately ap = 208 nT (or Kp = 8), i.e., a severe storm condition. The results are shown in Figure 5, where the relative density is averaged for poleward of ±60 ∘ magnetic latitude. As expected, the relative density increases with increasing E SW in both the Northern and Southern Hemispheres. Since the equinox condition was assumed, the E SW response is comparable in the two hemispheres. Later, we will show how the relative density changes with the IMF clock angle, solar activity, and season. It can be seen in Figure 5 that the relative density approximately linearly changes with E SW . This is why a nonlinear function √ E SW fitted better than E SW to the log densities: See our model parameterization in equation (5). The solar wind electric field is directly related to high-latitude ionospheric electric fields [e.g., Shepherd et al., 2002] and thus is a good indicator of the amount of the energy deposited from the solar wind to the magnetosphere-ionosphere-thermosphere system. The mechanisms by which the thermosphere density responses to the solar wind energy input were discussed in detail by Lei et al. [2010].
It is obvious from the results in Figure 5 that the HANDY attribution of density perturbations to P SW is much smaller than that to E SW . This indicates that the solar wind density makes a relatively small contribution to thermospheric density variability. It is noted that high-speed solar wind increases both E SW and P SW .
The dependence of the relative density on the IMF clock angle is depicted in Figure 6 for the Northern Hemisphere. The density response to the solar wind becomes greater as the IMF B z changes from positive to negative. The overall pattern of the relative density remains the same for different C SW conditions. Figure 7 shows the results for the Southern Hemisphere. Again, the IMF B z control of the relative density is evident. The IMF B y also modulates the high-latitude density response to the solar wind. As we reported in the previous paper [Yamazaki et al., 2015], the influence of the positive IMF B y in the Northern Hemisphere resembles the influence of the negative IMF B y in the Southern Hemisphere. In Figures 6 and 7, this effect is most evident in the density response in the dawn sector. In the Northern Hemisphere, the negative IMF B y reduces the density response at dawn, while in the Southern Hemisphere the positive IMF B y reduces the density response at dawn. Figure 8 illustrates, more clearly, the hemispheric asymmetry of the high-latitude thermospheric density due to IMF B y . The density ratios for the negative to positive IMF B y results show approximately the opposite pattern between the Northern Hemisphere and the Southern Hemisphere. It is noted that the results in Figure 8 vary little with solar activity and season. It is interesting to note that the pattern of the IMF B y effect on the high-latitude thermosphere bears some resemblance with the average high-latitude electric potential pattern [e.g., Weimer, 2005;Cousins and Shepherd, 2010]. That is, the regions of the large IMF B y effect in Figure 8 roughly correspond to the regions of the maximum and minimum electric potential. The mechanism by which the IMF B y affects the high-latitude density is not known, and hence, further investigation will be necessary. One possible way is through the action of vertical winds. It is known that at high latitudes, the pattern of horizontal winds depends on the IMF B y [Förster et al., 2008]. Changes in the horizontal wind system possibly alter vertical winds, which would affect the density by causing adiabatic heating or cooling and by altering the composition.  Figure 9 shows the high-latitude thermospheric density response to the solar wind under different solar EUV activity conditions. It can be seen that relative density perturbations are slightly smaller when solar activity is higher. As we showed in Figure 2, the background thermospheric density changes by an order of magnitude due to changes in solar activity. The fact that the dependence of the relative density on solar activity is small indicates that solar activity affects not only the background density but also the magnitude of density perturbations by a similar (but slightly smaller) rate. Emmert and Picone [2010] obtained similar results based on the analysis of the global average thermospheric density. They discussed that the reduced density response during high solar activity period may result from mitigated temperature perturbations due to increased scale height.

Dependence on the Background Density
The seasonal variation in the relative density is illustrated in Figure 10. The relative density tends to be small during local summer when the background density is relatively high. Previous studies have shown that the response of the absolute density to geomagnetic activity tends to be comparable in the Northern Hemisphere and the Southern Hemisphere [e.g., Yamazaki et al., 2015].

Summary and Conclusions
Using extensive accelerometer data by the CHAMP and GRACE satellites, we have constructed the first model to reproduce the response of the high-latitude thermospheric density to solar wind forcing. We call the model HANDY. HANDY is composed of a quiet model and disturbance model. Both are created at an altitude of 450 km. The quiet model represents the background density, and it was designed to be consistent with the MSISE-00 for the zero geomagnetic activity (ap = 0 nT) condition. The disturbance model represents relative density variations, from the quiet-time background, due to forcing by the solar wind. HANDY has new features that other global thermosphere models do not have, such as the use of solar wind parameter inputs (in contrast to conventional geomagnetic activity inputs) and the use of magnetic apex coordinates along with geographical coordinates.
The root-mean-square error for the HANDY model is approximately 15% during the period of investigation. A comparison with other commonly used models indicates that the HANDY performance is reasonably good.
(The corresponding error is ∼23% for the MSISE-00 model and ∼19% for the JB2008 model.) It should be noted, however, that the data we analyzed for the HANDY construction were limited for August 2002 to

10.1002/2015JA021371
July 2006, and thus, the model is not constrained by the measurements around solar maximum or solar minimum. The representation of the solar cycle variation in the background density may be improved by including more data especially during high or low solar activity.
The HANDY results for the high-latitude density response to solar wind forcing are consistent with previous studies. The magnitude of relative density perturbations depends strongly on the solar wind electric field intensity and the orientation of the IMF. The response of the relative density to the solar wind electric field is most significant in the premidnight sector, where substorm onset often takes place. The solar wind response is relatively strong in the noon sector and relatively weak in the predawn sector. These features are completely missing from the MSISE-00 model. Also, the HANDY results clearly demonstrated the IMF B y effect, which we found in the previous study [Yamazaki et al., 2015]. The effect of the IMF B y in the Northern Hemisphere is similar to that in the Southern Hemisphere for the opposite sense of the IMF B y .
Although HANDY has been shown to be useful in revealing features of the high-latitude thermospheric density, the model restriction to a single altitude is an obvious disadvantage for practical purposes. One possible approach for the vertical extension of the model is to construct similar models at different altitudes and interpolate them, which would require a more extensive data set from various satellites at different heights. Another approach is to introduce an assumption on the height profile of the air density. Most global empirical models assume diffusive equilibrium, which is a good approximation above ∼200 km. However, this approach requires knowledge of temperature and composition for at least two heights. We will leave the extension of HANDY for a topic of a future study.
Finally, the supporting information contains a document describing the model formats, along with lists of quiet and disturbance model coefficients with corresponding 1 errors.