New method in computer simulations of electron and ion densities and temperatures in the plasmasphere and low-latitude ionosphere

A new theoretical model of the Earth's low- and mid-latitude ionosphere and plasmasphere has been devel- oped. The new model uses a new method in ionospheric and plasmaspheric simulations which is a combination of the Eu- lerian and Lagrangian approaches in model simulations. The electron and ion continuity and energy equations are solved in a Lagrangian frame of reference which moves with an in- dividual parcel of plasma with the local plasma drift velocity perpendicular to the magnetic and electric fields. As a re- sult, only the time-dependent, one-dimension electron and ion continuity and energy equations are solved in this La- grangian frame of reference. The new method makes use of an Eulerian computational grid which is fixed in space co- ordinates and chooses the set of the plasma parcels at every time step, so that all the plasma parcels arrive at points which are located between grid lines of the regularly spaced Eule- rian computational grid at the next time step. The solution values of electron and ion densities Ne and Ni and temper- atures Te and Ti at the Eulerian computational grid are ob- tained by interpolation. Equations which determine the tra- jectory of the ionospheric plasma perpendicular to magnetic field lines and take into account that magnetic field lines are "frozen" in the ionospheric plasma are derived and included in the new model. We have presented a comparison between the modeled NmF2 and hmF2 and NmF2 and hmF2 which were observed at the anomaly crest and close to the geomagnetic equator simultaneously by the Huancayo, Chiclayo, Talara, Bogota, Panama, and Puerto Rico ionospheric sounders during the 7 October 1957 geomagnetically quiet time period at solar maximum. The model calculations show that there is a need to revise the model local time dependence of the equatorial upward E ◊ B drift velocity given by Scherliess and Fe- jer (1999) at solar maximum during quiet daytime equinox conditions. Uncertainties in the calculated Ni , Ne, Te, and Ti resulting from the difference between the NRLMSISE- 00 and MSIS-86 neutral temperatures and densities and from the difference between the EUV97 and EUVAC solar fluxes are evaluated. The decrease in the NRLMSISE-00 model (O)/(N2) ratio by a factor of 1.7-2.1 from 16:12 UT to 23:12 UT on 7 October brings the modeled and measured NmF2 and hmF2 into satisfactory agreement. It is shown that the daytime peak values in Te, and Ti above the ionosonde stations result from the daytime peak in the neutral tem- perature. Our calculations show that the value of Te at F2- region altitudes becomes almost independent of the electron heat flow along the magnetic field line above the Huancayo, Chiclayo, and Talara ionosonde stations, because the near- horizontal magnetic field inhibits the heat flow of electrons. The increase in geomagnetic latitude leads to the increase in the effects of the electron heat flow along the magnetic field line on Te. It is found that at sunrise, there is a rapid heat- ing of the ambient electrons by photoelectrons and the differ- ence between the electron and neutral temperatures could be increased because nighttime electron densities are less than those by day, and the electron cooling during morning con- ditions is less than that by day. This expands the altitude region at which the ion temperature is less than the electron temperature near the equator and leads to the sunrise electron temperature peaks at hmF2 altitudes above the ionosonde sta- tions. After the abrupt increase at sunrise, the value of Te de- creases, owing to the increasing electron density due to the increase in the cooling rate of thermal electrons and due to the decrease in the relative role of the electron heat flow along the magnetic field line in comparison with cooling of ther- mal electrons. These physical processes lead to the creation of sunrise electron temperature peaks which are calculated above the ionosonde stations at hmF2 altitudes. We found that the main cooling rates of thermal electrons are electron- ion Coulomb collisions, vibrational excitation of N2 and O2, and rotational excitation of N2. It is shown that the increase in the loss rate of O + ( 4 S) ions due to the vibrational excited N2 and O2 leads to the decrease in the calculated NmF2 by a factor of 1.06-1.44 and to the increase in the calculated hmF2, up to the maximum value of 32 km in the low-latitude ionosphere between -30 and +30 of the geomagnetic lati-

the difference between the EUV97 and EUVAC solar fluxes are evaluated.The decrease in the NRLMSISE-00 model [O]/[N 2 ] ratio by a factor of 1.7-2.1 from 16:12 UT to 23:12 UT on 7 October brings the modeled and measured NmF2 and hmF2 into satisfactory agreement.It is shown that the daytime peak values in T e , and T i above the ionosonde stations result from the daytime peak in the neutral temperature.Our calculations show that the value of T e at F2region altitudes becomes almost independent of the electron heat flow along the magnetic field line above the Huancayo, Chiclayo, and Talara ionosonde stations, because the nearhorizontal magnetic field inhibits the heat flow of electrons.The increase in geomagnetic latitude leads to the increase in the effects of the electron heat flow along the magnetic field line on T e .It is found that at sunrise, there is a rapid heating of the ambient electrons by photoelectrons and the difference between the electron and neutral temperatures could be increased because nighttime electron densities are less than those by day, and the electron cooling during morning conditions is less than that by day.This expands the altitude region at which the ion temperature is less than the electron temperature near the equator and leads to the sunrise electron temperature peaks at hmF2 altitudes above the ionosonde stations.After the abrupt increase at sunrise, the value of T e decreases, owing to the increasing electron density due to the increase in the cooling rate of thermal electrons and due to the decrease in the relative role of the electron heat flow along the magnetic field line in comparison with cooling of thermal electrons.These physical processes lead to the creation of sunrise electron temperature peaks which are calculated above the ionosonde stations at hmF2 altitudes.We found that the main cooling rates of thermal electrons are electronion Coulomb collisions, vibrational excitation of N 2 and O 2 , and rotational excitation of N 2 .It is shown that the increase in the loss rate of O + ( 4 S) ions due to the vibrational excited N 2 and O 2 leads to the decrease in the calculated NmF2 by a factor of 1.06-1.44 and to the increase in the calculated hmF2, up to the maximum value of 32 km in the low-latitude ionosphere between -30 and +30 • of the geomagnetic lati-

Introduction
At low-and mid-latitudes, the Earth's magnetic field can be represented, to a good approximation, by a dipole.The horizontal orientation of the geomagnetic field at the geomagnetic equator is known to be the basic reason for the active nature of the low-latitude ionosphere, which is characterised by the equatorial electrojet, equatorial plasma fountain, equatorial (Appleton) anomaly, additional layers, plasma bubbles, and spread-F.These equatorial characteristic properties of the ionosphere have been studied observationally and theoretically for many years (see, for example, review papers presented by Moffett, 1979;Anderson, 1981;Walker, 1981;Rishbeth, 2000;Abdu, 1997, 2001, andreferences therein).Many theoretical models of the plasmasphere and low-latitude ionosphere were constructed and have been applied to study a wide variety of equatorial ionosphere characteristic properties.Among these models, it is necessary to point out the following major sophisticated plasmaspheric and low-latitude ionospheric models: the Sheffeld University plasmasphere-ionosphere model (Bailey and Sellec, 1990;Bailey and Balan, 1996), the coupled thermosphere-ionosphere-plasmasphere model (CTIP) (Fuller-Rowell et al., 1988;Millward et al., 1996), a coupled thermosphere-ionosphere model (CTIM) (Fuller-Rowell et al., 1996), the global theoretical ionospheric model (GTIM) (Anderson, 1973;Anderson et al., 1996), and the global numerical self-consistent and time-dependent model of the thermosphere, ionosphere, and protonosphere (Namgaladze et al., 1988).These models include transport of plasma by geomagnetic field-aligned diffusion and neutral windinduced plasma drift of ions and electrons, and plasma motion perpendicular to the geomagnetic field, B, due to an electric field, E, which is generated in the E-region.This electric field affects F-region plasma causing both ions and electrons to drift in the same direction with an drift velocity, V E = E × B/B 2 .In a Lagrangian method, the finite-difference grid moves with the local plasma drift velocity V E perpendicular to the magnetic and electric fields.The rate of change of electron and ion number densities and temperatures in a moving frame of reference is much easier to compute because the convective terms in the continuity and energy equations are absent in the moving frame.As a result, it is needed to solve only one-dimensional, time dependent ion and electron continuity and energy equations along magnetic field lines in this moving frame of reference.
Contrary to a Lagrangian computational grid, an Eulerian computational grid is fixed in space coordinates.The main aim of our work is to elaborate a new approach which includes the advantages of both approaches in solving electron and ion continuity and energy equations in the ionosphere and plasmasphere.Our new approach is a combination of the Eulerian and Lagrangian approaches in model simulations.This new method is used to construct a new model of the plasmasphere and ionosphere which will be used to calculate electron and ion densities and temperature in the plasmasphere and ionosphere at low and middle latitudes.
In the present work we investigate the equatorial anomaly using the constructed new model of the plasmasphere and ionosphere and the progress in understanding the F2-layer physics that has come from the development of models of the thermosphere and ionosphere.Our purpose is to discuss the models' success in reproducing the equatorial anomaly phenomenon.In contrast to previous studies of the equatorial anomaly, the model of the ionosphere and plasmasphere used in this work includes the fundamental laboratory rate coefficient measurements of O + ( 4 S) ions with vibrationally excited N 2 and O 2 given by Hierl et al. (1997), the quenching rate coefficients for O + ( 2 D) and O + ( 2 D) by N 2 measured by Li et al. (1997), the updated Einstein coefficients for the O + ( 2 P) → O + ( 4 S) + hν and O + ( 2 P) → O + ( 2 D) + hν transitions given by Kaufman and Sugar (1986), and the updated photoionization and photoabsorption cross sections for the N 2 , O 2 , and O photoionization reactions which form N + 2 , O + 2 , O + ( 4 S), O + ( 2 D), O + ( 2 P), O + ( 4 P), and O + ( 2 P * ) ions (Richards et al., 1994;Schaphorst et al., 1995;Berkowitz, 1997).
There is a strong dependence of the equatorial anomaly characteristics (i.e.crest latitudes and magnitudes) on the vertical drift velocity of the equatorial F-layer, and the theoretically modeled low-latitude distributions of the electron density are very sensitive to input drift velocities (Klobuchar et al., 1991).The present work reports the attempt to study some features of this relationship in the case study in which NmF2 electron densities are observed at the anomaly crest and close to the geomagnetic equator simultaneously, near approximately the same geomagnetic meridian by the Panama, Bogota, Talara, Chiclayo, and Huancayo ionospheric sounders during the 7 October 1957 time period.The model wishes to look at the effects of changing V E .
The model of the ionosphere and plasmasphere uses the solar EUV flux and the neutral temperature and densities as the model inputs.As a result, the model/data discrepancies arise due to uncertainties in EUV fluxes and a possible inability of the neutral atmosphere model to accurately predict the thermospheric response to the studied time period in the upper atmosphere.Over the years, testing and modification of the MSIS neutral atmosphere model has continued, and it has led to improvements through several main versions of this neutral atmosphere model: MSIS-77 (Hedin et al., 1977 a, b), MSIS-86 (Hedin, 1987), and NRLMSISE-00 (Picone et al., 2000(Picone et al., , 2002)).In the present work we investigate how well the Panama, Bogota, Talara, Chiclayo, and Huancayo ionospheric sounder measurements of electron densities taken during the geomagnetically quiet period of 7 October 1957 agree with those calculated by the model of the ionosphere and plasmasphere using the MSIS-86 or NRLMSISE-00 neutral temperature and densities.The model of the ionosphere and plasmasphere has an option to use the solar EUV fluxes from the EUVAC model (Richards et al., 1994) or the EUV97 model (Tobiska and Eparvier, 1998).As a result, the model/data agreement can be better or worse when we use NRLMSISE-00, as opposed to MSIS-86 and EUVAC, as opposed to EUV97.One objective of the model/data comparison which is carried out in this work is to present an evaluation of uncertainties in model calculations of electron and ion densities and temperatures from the comparison between neutral atmosphere models and between the solar flux models as input model parameters.
The O + ( 4 S) ions that predominate at ionospheric F2region altitudes are lost in the reactions of O + ( 4 S) with unexcited N 2 (v = 0) and O 2 (v = 0) and vibrationally excited N 2 (v) and O 2 (v) molecules at vibrational levels, v > 0. Vibrationally excited N 2 and O 2 react more strongly with O + ( 4 S) ions in comparison with unexcited N 2 and O 2 (Schmeltekopf et al., 1968, Hierl et al., 1997).As a result, an additional reduction in the electron density is caused by the reactions of O + ( 4 S) ions with vibrationally excited N 2 and O 2 .Numerical simulations of the ionosphere show that the daytime mid-latitude electron density of the F2-region should be reduced by a factor of 1.5-2.5, due to enhanced vibrational excitation of N 2 at high solar activity during geomagnetically quiet and storm periods (see Pavlov and Foster, 2001, and references therein).The reduction to two-thirds of its value due to vibrationally excited N 2 is found in the low-latitude F-region electron density at the location of the equatorial trough at solar maximum (Jenkins et al., 1997).The increase in the O + + O 2 loss rate due to vibrationally excited O 2 decreases the simulated daytime F2 peak density by up to a factor of 1.7 at high solar activity (Pavlov, 1998b;Pavlov et al., 1999Pavlov et al., , 2000;;2001;Pavlov and Foster, 2001).In this paper we examine the latitude dependence of the effects of vibrationally excited N 2 and O 2 on the electron density and temperature at solar maximum during geomagnetically quiet conditions of 7 October 1957, to investigate the role of vibrationally excited N 2 and O 2 in the formation of the observed electron density equatorial anomaly variations.

Theoretical model
We present a new model of the middle-and low-latitude ionosphere and plasmasphere.This model uses a dipole approximation to the Earth's magnetic field and takes into account the offset between the geographic and geomagnetic axes.The horizontal components of the neutral wind, which are used in calculations of the wind-induced plasma drift velocity along the magnetic field, are specified using the HWW90 wind model of Hedin et al. (1991).In the model, time-dependent ion continuity equations for the three major ions, O + ( 4 S), H + , and He + and for the minor ions, NO + , O + 2 , and N + 2 are solved by taking into account the production and loss rates of ions, transport of plasma by geomagnetic fieldaligned diffusion and neutral wind-induced plasma drift of ions and electrons and plasma motion perpendicular to the geomagnetic field due to an electric field which is generated in the E-region.The approach of the local chemical equilibrium is used to calculate steady-state number densities of O + ( 2 D), O + ( 2 P), O + ( 4 P), and O + ( 2 P * ) ions.Timedependent electron and ion energy balance equations are solved in the model.These equations include heating and cooling rates of electron and ions and a term due to the E × B drift of electrons and ions.Modelled electron heating caused by collisions between thermal electrons and photoelectrons is provided by a solution of the Boltzmann equation for photoelectron flux along a centered -dipole magnetic field line, the same field line used for solving for number densities of electron and ions, and electron and ion temperatures at the same grid point.The model uses Boltzmann distributions of N 2 (v) and O 2 (v) to calculate [N 2 (v)] and [O 2 (v)] which are included in the model loss rate of O + ( 4 S) ions and cooling rates of thermal electrons, due to vibrational excitation of N 2 and O 2 .The chemistry, physics, and solution procedure have been described in detail in Appendix A.
The coordinate system considered in this work is presented in Appendix A. Orthogonal curvilinear coordinates are: q, U , and a geomagnetic longitude, .The important properties of these coordinates are that q is aligned with, and U and are perpendicular to, the magnetic field, the U and coordinates are constant along a dipole magnetic field line, and the McIlwain parameter L can be presented as L = U −1 .The model takes into account that the plasma E × B drift velocity can be presented as the dipole coordinate system, e and e U are unit vectors in and U directions, respectively.
The zonal component, V E , of the E × B drift is not included in the our model calculations (see Appendix A) as it is believed that this E × B drift component has a negligible effect on the electron density profiles (Anderson, 1981).It should be noted that, as far as the author knows, possible effects of the zonal component of the E × B drift on electron and ion densities and temperatures are not included in the published model calculations of the ionospheric equatorial anomaly variations (see, for example, Bailey and Sellec, 1990;Bailey and Balan, 1996;Su et al., 1997).
The equatorial magnitude of the meridional component of the E × B drift velocity has been found to vary greatly from day to day, and these drift velocities have large seasonal and solar cycle variations (Woodman, 1970;Fejer et al., 1989Fejer et al., , 1995;;Scherliess and Fejer, 1999).It is also known to be longitude dependent (Schieldge et al., 1973;Fejer et al., 1995).8 of Scherliess and Fejer (1999) for high solar activity and equinox conditions was used to find the equatorial value of E (dash-dotted line).Solid line shows the empirical equatorial electric field which was modified in the time range between 07:27 LT and 11:00 LT and between 15:00 LT and 18:30 LT by the use of the comparison between the measured and modelled values of NmF2 and hmF2 at 00:00 UT and 16:00 UT.The average quiet time value of E at the F-region altitudes over Arecibo (dashed line) is found from the average quiet time perpendicular/northward F-region plasma drifts for equinox conditions presented in Fig. 2 of Fejer (1993).
There is evidence that the vertical E × B drift velocity varies with altitude at the geomagnetic equator (Pingree and Fejer, 1987).In the present study, the simplistic approach is used to calculate the dependence of E on t ge , where t ge is the local time at the geomagnetic equator for the magnetic longitude of each ionosonde station.
In the model, the value of E (t ge ) over the geomagnetic equator given by the dash-dotted line in Fig. 1 is obtained from the empirical F-region quiet time equatorial vertical drift velocity presented in Fig. 8 of Scherliess and Fejer (1999) for high solar activity and equinox conditions.As it will be discussed later in Sect.4, this empirical equatorial electric field is modified in the time range between 07:27 LT and 11:00 LT and between 15:0 LT and 18:30 LT by the use of the comparison between the measured and modelled nighttime values of hmF2.The resulting equatorial magnitude of E (t ge ), which is used in the model calculations is shown by the solid line in Fig. 1.The average quiet time value of E at the F-region altitudes over Arecibo (dashed line in Fig. 1) is found from Fig. 2 of Fejer (1993), where the average quiet time perpendicular/northward F-region plasma drifts for high solar activity and equinox conditions is presented.
Equations (A11), (A12), and (A15) determine the trajectory of the ionospheric plasma perpendicular to magnetic field lines and the moving coordinate system.It follows from Eq. (A11) that time variations of U caused by the existence of the E component of the electric field are determined by time variations of the component, E eff , of the effective electric field given by Eq. (A12).We have to take into account that the magnetic field lines are "frozen" in the ionospheric plasma (see Sect.A2.5.1 of Appendix A).As a result, E eff (t) is not changed along magnetic field lines (see Eq. A15).The equatorial and Arecibo values of E (t ge ) are used to find the equatorial and Arecibo values of E eff (t ge ) from Eqs. (A12) and (A15).The equatorial value of E eff (t ge ) (the equatorial E (t ge ) is given by the solid line in Fig. 1) is used for magnetic field lines with an apex altitude, h ap = R eq − R E , less than 600 km, where R eq is the equatorial radial distance of the magnetic field line from the Earth's center and R E is the Earth's radius.The Arecibo value of E eff (t ge ) (the Arecibo E (t ge ) is given by the dashed line in Fig. 1) is used if the apex altitude is greater than 2 126 km.Linear interpolation of the equatorial and Arecibo values of E eff (t ge ) is employed at the intermediate apex altitudes.
The model starts at 15:12 UT on 5 October 1957.This UT corresponds to 10:00 LT at the geomagnetic equator and 351.9 • of the geomagnetic longitude (see explanations of the value of the geomagnetic longitude in Sect.3).First of all, the steady-state N i , N e , T i , and T e are found by the use of the model of the ionosphere and plasmasphere with E = 0 (i.e.without the E × B drift velocity).It means that the onedimensional time dependent Eqs.(A1), (A6), and (A7) of Appendix A are solved along each computational grid dipole magnetic field line at 10:00 UT on 5 October 1957, to obtain the N i , N e , T i , and T e initial conditions.These steadystate daytime values of N i , N e , T i , and T e are used as initial conditions to solve the two-dimensional, time dependent  Fejer (1999) given by the dash-dotted line in Fig. 1 is used.Right panels (e), (f), (g), (h) show the model results when the empirical equatorial electric field found the equatorial perpendicular plasma drift velocity of Scherliess and Fejer (1999) was modified in the time range between 07:27 LT and 11:00 LT and between 15:00 LT and 18:30 LT (this modified equatorial electric field is shown by a solid line in Fig. 1

Solar geophysical conditions and data
The value of the geomagnetic K p index was between 0 and 2 for the studied time period of 7 October 1957.It should be noted that when the thermosphere is disturbed, it takes time for it to relax back to its initial state, and this thermosphere relaxation time determines the time for the disturbed ionosphere to relax back to the quiet state.It means that not every time period with K p ≤ 3 can be considered as a magnetically quiet time period.The characteristic time of the neutral composition recovery after a storm impulse event ranges from 7 to 12 hours, on average (Hedin, 1987), while it may need up to days for all altitudes down to 120 km in the atmosphere to recover completely back to the undisturbed state of the atmosphere (Richmond and Lu, 2000).The value of K p was between 0 and 3 for the previous 4-6 October 1957 time period, i.e. the studied time period of 7 October 1957 can be considered as a magnetically quiet time period.The F10.7 solar activity index was 254 on 7 October 1957, while the 81-day averaged F10.7 solar activity index was 234.
Our study is based on hourly critical frequencies fof2, and foE, of the F2 and E-layers, and the maximum usable frequency parameter, M(3000)F2, data from the Huancayo, Chiclayo, Talara, Bogota, Panama, and, Puerto Rico ionospheric sounder stations available on the Ionospheric Digital Database of the National Geophysical Data Center, Boulder, Colorado.Locations of these ionospheric sounder stations are shown in Table 1.The first five sounders at low latitude are within ±3.5 • geomagnetic longitude of one another.As a result, all model simulations are carried out for the geomagnetic longitude of 351.9 • .To complete the picture of the latitude dependence of NmF2 and hmF2 variations, we compare the modeled NmF2 and hmF2 at the geomagnetic longitude of 351.9 • with Nmf2 and hmF2 measured by the Puerto Rico ionosonde station with geomagnetic longitude 2.8 • .The Puerto Rico sounder departs slightly from the near conjugacy of the Huancayo, Chiclayo, Talara, Bogota, and Panama ionospheric sounder stations, but this geomagnetic longitude deviation is nonsignificant for our study because the equatorial anomaly effects are less pronounced at the Puerto Rico sounder in comparison with those at the first five sounders.The values of the peak density, NmF2, of the F2-layer is related to the critical frequencies fof2 as NmF2 = 1.24 • 10 10 fof2 2 , where the unit of NmF2 is m −3 , the unit of fof2 is MHz.To determine the ionosonde values of hmF2, we use the relation between hmF2 and the values of M(3000)F2, fof2, and foE recommended by Dudeney (1983) from the comparison of different approaches as hmF2 = 1490/[M(3000)F2 + M]-176, where M = 0.253/(fof2/foE−1.215)−0.012.If there are no foE data, then it is suggested that M = 0, i.e. the hmF2 formula of Shimazaki (1955) is used.

Equatorial perpendicular electric field modification
In Fig. 2, geomagnetic latitude plots are shown of hmF2 and NmF2 at 00:00 UT (panels (a), (b), (e), and (f)) and 16:00 UT (panels (c), (d), (g), and (h)) on 7 October 1957 from the ionospheric sounder station measurements (squares) and model calculations (solid, dotted, and dashed lines).Four left panels (a), (b), (c), and (d) show the model results when the original equatorial perpendicular plasma drift of Scherliess and Fejer (1999), given by the dash-dotted line in Fig. 1, is used.Four right panels (e), (f), (g), and (h) show the model results when the empirical equatorial electric field, found from the equatorial perpendicular plasma drift velocity of Scherliess and Fejer (1999), was modified in the time range between 07:27 LT and 11:00 LT and between 15:00 LT and 18:30 LT (this modified equatorial electric field is shown by a solid line in Fig. 1).Dashed lines show the model results when the original NRLMSISE-00 neutral temperature and densities are used.Solid lines show the model results when the NRLMSISE-00 model [O] was decreased by a factor of 1.7 from 16:12 UT to 23:12 UT (from 11:00 LT to 18:00 LT, where LT is the local time at the geomagnetic equator and 351.9 • of the geomagnetic longitude) during all model simulation periods.Dotted lines show the model results when the NRLMSISE-00 model [N 2 ] and [O 2 ] were increased by a factor of 2.1 from 16:12 UT to 23:12 UT during all model simulation periods.The vibrationally excited N 2 (v > 0) and O 2 (v > 0) are included in the model, and the EUVAC solar flux model is used as the input model parameter in all model calculations presented in Fig. 2.
The comparison between the results shown in the two upper panels (a) and (e) of Fig. 2 clearly indicates that there is a large disagreement between the measured and modelled hmF2 at 00:00 UT on 7 October 1957, if the equatorial upward E × B drift given by Scherliess and Fejer (1999) is used.The results presented in panels (a) and (b) of Fig. 2 provide evidence that we can match the measured and modeled NmF2 using the corrected neutral densities.However, the corrections of the NRLMSISE-00 model do not bring the measured and modeled hmF2 into agreement.We conclude that this disagreement in hmF2 is caused by the long time duration of the pre-reversal strengthening of the equatorial upward E × B drift given by Scherliess and Fejer (1999).The high estimate of this pulse duration in E leads to unreal, high-modeled F2 peak altitudes at 00:00 UT.Our calculations presented in the panels (e) and (f) of Fig. 2 provide evidence that, to bring the measured and modeled F2-region main peak altitudes into agreement, the magnitude of E has to be approximately constant in the time range between 15:00 LT and 18:00 LT with the following peak in E , which has a shorter time width in comparison with the time duration of the pre-reversal strengthening of the original equatorial perpendicular plasma drift given by Scherliess and Fejer (1999).Fejer et al. (1989) show solar maximum ion vertical drifts over Jicamarca, which is very close to Huancayo, near the magnetic equator.The pre-reversal enhancements during quiet equinox periods can begin as late as 18:00 LT and peak after 19:00 LT, although the average enhancements are earlier and broader.In addition, Batista et al. (1986) estimate the pre-reversal enhancement at Huancayo to peak between 18:00 LT and 19:00 LT based on hmF2 changes during equinox solar maximum conditions.Thus, our delay of the pre-reversal enhancement until 18:00 LT is in agreement with the observed day-to-day variability at Jicamarca, and previous estimates for Huancayo.
The principal feature of the equatorial anomaly is that the crest-to-trough ratio is increased with increasing upward E × B drift (Dunford, 1967;Su et al., 1997;Rishbeth, 2000).The measurements show that, by mid-afternoon (15:00 UT), the equatorial anomaly crests are forming away from the geomagnetic equator, while the model calculations with the equatorial E (t ge ), given by the dash-dotted line in Fig. 1, produce the onset of the equatorial anomaly crest formation close to 16:00 UT (see lines in panels (c) and (d) of Fig. 2).The disagreement between the sizes of the equato- Our calculations show that a strengthening of the equatorial upward E × B drift before 17:00 UT on 7 October 1957 leads to an increase in the northern and southern depths of the equatorial NmF2 trough (these depths can be expressed as ratios of NmF2 at F2-region northern and southern crests to an equatorial NmF2) at 17:00 UT on 7 October 1957.The modification of E (t ge ) is shown by a solid line in Fig. 1.This modification, which was carried out in the time range between 07:27 LT and 11:00 LT, includes the strengthening of E and the time shift of the peak in E (t ge ) relative to the peak in E (t ge ), shown by the dash-dotted line in Fig. 1.The first maximum (E = 1.1 mVm −1 ) of the modified E (t ge ) occurs at 08:30 LT while the first maximum (E = 0.66 mVm −1 ) of E (t ge ), given by the dash-dotted line in Fig. 1, is located between 10:00 LT and 11:00 LT.It should be noted that the revised magnitude of the first peak in E is close to the magnitude of the second peak in E , although the second peak is expected to be about double the size of the first peak in solar maximum (Fesen et al., 2000).However, the Fejer et al. (1995) quiet-time equinox model from the AE-E satellite observations during moderate to high solar flux conditions at 260 • E suggest a morning peak at 09:30 LT, which is nearly equal to the sharper evening reversal peak, similar to our proposed changes.The comparison between the squares and solid line in panel (h) of Fig. 2 shows that the northern depth of the equatorial NmF2 trough in the calculated NmF2 is approximately consistent with the measured depth, if the modified E (t ge ) is used.The model NmF2 is higher than the observations with the anomaly crest shifted poleward, if the original NRLMSISE-00 model is used (see dashed line in panel (h) of Fig. 2).Panel (g) of Fig. 2 shows that the agreement between measured and modeled hmF2 is somewhat worse.It should be noted that the model with the modified value of E (t ge ) produces the onset of the equatorial anomaly crest formation close to 15:00 UT, in agreement with the measured onset of the equatorial anomaly crest formation given by ionosonde stations.
As a result, the equatorial electric field, shown by the solid line in Fig. 1, and the Arecibo value of E (t ge ), shown by the dashed line in Fig. 1, are used in all subsequent model calculations presented in this paper, as described in Sect. 2.

Evaluation of uncertainties in model calculations of
NmF2 and hmF2 from the comparison between neutral atmosphere models and between the solar flux models as input model parameters The measured (squares) and calculated (lines) NmF2 and hmF2 are displayed in the two lower panels of Figs.3-8 for the 7 October 1957 time period above the Huancayo (Fig. 3), Chiclayo (Fig. 4), Talara (Fig. 5), Bogota (Fig. 6), Panama (Fig. 7), and Puerto Rico (Fig. 8) ionosonde stations.The results obtained from the model of the ionosphere and plasmasphere, using the combinations of the original NRLMSISE-00 or MSIS-86 neutral temperature and density models and the EUVAC or EUV97 solar flux models as the input model parameters, are shown by solid lines (the NRLMSISE-00 model in combination with the EUVAC model), dotted lines (the NRLMSISE-00 model in combination with the EUV97 model), dash-dotted lines (the MSIS-86 model in combination with the EUVAC model), and dashed lines (the MSIS-86 model in combination with the EUV97 model).
The differences in original neutral densities and temperatures from the NRLMSISE-00 and MSIS-86 models result in the differences between solid and dash-dotted lines (the EUVAC solar flux model is used) and between dotted and dashed lines (the EUV97 solar flux model is used).We found that the use of the NRLMSISE-00 model, as opposed to the MSIS-86 model, leads to the highest possible increase in the calculated NmF2 by a maximum factor of 1. 29, 1.23, 1.22, 1.28, 1.30, and 1.39 and to the highest possible decrease in the calculated hmF2 by 38, 42, 40, 22, 22 and 23 km above the Huancayo, Chiclayo, Talara, Bogota, Panama, and Puerto Rico ionosonde stations, respectively.The use of the EUV97 solar flux model as opposed to the EUVAC solar flux model leads to the increase in the calculated NmF2 by a factor of 1.13-1.34and to the highest possible variations in calculated hmF2 of 12 km.Our calculations clearly show that the best agreement between the measured and modeled electron densities is obtained if the MSIS-86 neutral densities and temperature in combination with the EUVAC solar flux (dash-dotted lines in Figs.3-8) are used as the input model parameters.At the same time, the NRLMSISE-00 model is the outgrowth of the MSIS-86 model, and we have a right to expect that the NRLMSISE-00 model describes real neutral temperature and densities variations more accurately in comparison to the MSIS-86 model.Therefore, the NRLMSISE-00 neutral temperature and density model of Picone et al. (2000Picone et al. ( , 2002)), and the EUVAC solar flux model of Richards et al. (1994) are used in the further model calculations presented in this work.
The results obtained from the model of the ionosphere and plasmasphere, using the combinations of the NRLMSISE-00 or MSIS-86 neutral temperature and density models and the EUVAC or EUV97 solar flux models as the input model parameters, are shown by solid lines (the NRLMSISE-00 model in combination with the EUVAC model), dotted lines (the NRLMSISE-00 model in combination with the EUV97 model), dash-dotted lines (the MSIS-86 model in combination with the EUVAC model), and dashed lines (the MSIS-86 model in combination with the EUV97 model).It is evident from Figs. 3-8 that the electron and ion temperature changes created by the difference between the NRLMSISE-00 and MSIS-86 neutral temperatures and number densities or by the difference between the EUV97 and EUVAC solar fluxes are negligible.
The electron and ion temperatures start to increase from their night-time values close to 10:12 UT.The electron temperature reaches a morning peak at about 10:52-11:12 UT, above the ionosonde stations of Table 1, while the ion temperatures above all the ionosonde stations presented in Table 1 have no morning peaks at hmF2.Following the morn- ing peak, there is a rapid decrease in the electron temperature, which reaches a minimum at around 11:32-13:12 UT.After the electron temperature minimum, the electron temperature increases again above the ionosonde stations.The peak values in the electron and ion temperatures above all the ionosonde stations presented in Table 1 occur at about 20:32-21:22 UT.Our calculations show that the magnitudes of the electron and ion temperatures at hmF2 are close to the neutral temperature at hmF2 during most of the daytime conditions.As a result, the peak values in the electron and ion temperatures result from the peak in the neutral temperature at hmF2, which occurs very close to the time of the peaks in the electron and ion temperatures above the ionosonde stations presented in Table 1.
It is well known that in the ionospheric F-region, there is an inverse relationship between electron temperature and electron density, i.e. greater electron densities produce lower electron temperatures.As a result, the electron temperature is close to the neutral temperature during most of the daytime period at hmF2 altitudes at solar maximum, due to high magnitudes of electron cooling rates in comparison with the input of the electron heat flow along the magnetic field line in the daytime electron temperature variations.
Electron and ion temperatures profiles measured at Jicamarca close to the geomagnetic equator between June 1965 and November 1966 are such that T e > T i between about 200 km and 300 km during daytime conditions, T e ≈ T i ≈ const between about 300 km and 500 km, and the value of T e is close to the neutral temperature in this daytime isothermal altitude region (McClure, 1969;Schunk and Nagy, 1978).This daytime electron and ion isothermal region can come up to 600 km (Bailey et al., 1975).Throughout this region the electron thermal conduction term in the thermal balance equation of electrons is negligible in comparison with cooling of electrons due to collisions of thermal electrons with ions and neutral gases (Bailey et al., 1975;Schunk and Nagy, 1978).Our calculations show that the values of the electron temperatures at F2-region altitudes become almost independent of the electron heat flow along the magnetic field line above the Huancayo (Fig. 3), Chiclayo (Fig. 4), and Talara (Fig. 5) ionosonde stations because the near-horizontal magnetic field inhibits this heat flow of electrons.The increase in geomagnetic latitude leads to the increase in the effects of the electron heat flow along the magnetic field line on T e .
It follows from the electron and ion temperatures profiles measured at Jicamarca that the enlargement of the altitude region with T e > T i occurs at sunrise at all heights to at least 600 km (McClure, 1969).Our calculations show that at sunrise, there is a rapid heating of the ambient electrons by photoelectrons, and the difference between the electron and neutral temperatures could be increased because nighttime electron densities are less than those by day, and the electron cooling during morning conditions is less than that by day.This expands the altitude region at which T e > T i near the equator, and leads to the sunrise electron temperature peaks at hmF2 altitudes above the ionosonde stations.After the abrupt increase at sunrise, the electron temperature decreases, owing to the increasing electron density due to the increase in the cooling rate of thermal electrons and due to the decrease in the relative role of the electron heat flow along the magnetic field line in comparison with cooling of thermal electrons.As a result, the morning electron temperature peaks which are found above the ionosonde stations at hmF2 altitudes are explained by these physical processes.
Early studies have pointed out that the radar T e measured at Jicamarca are lower than T e measured by using probes on satellites, and there was a problem with unreal night-time radar T e < T i (McClure et al., 1973;Aponte et al., 2001).This problem was solved by Sulzer and Gonzalez (1999) and Aponte et al. (2001), who found that electron-electron and electron-ion Coulomb collisions are responsible for the additional incoherent backscatter spectral narrowing above Jicamarca, leading to the change in the measured T e /T i ratio.Specifically, T e = T i at night, at F-region altitudes above Jicamarca (Aponte et al., 2001).In agreement with this conclusion, our calculations show that T e = T i at night, at F-region altitudes close to the geomagnetic equator.
The relative magnitudes of the cooling rates are of particular interest for understanding the main processes which determine the electron temperature.The model of the iono- sphere and plasmasphere uses the electron cooling rates due to electron-ion Coulomb collisions and elastic collisions of electrons with N 2 , O 2 , O, He and H, the thermal electron impact excitation of O 2 (a 1 g ), O 2 (b 1 g + ), and the fine structure levels of the ground state of atomic oxygen, the rates of electron cooling through vibrational and rotational excitation of N 2 and O 2 , and the electron energy loss arising from electron-impact-induced transitions 3 P → 1 D for atomic oxygen (see Sect.A2.2 of Appendix A).The relative role of the electron cooling rates was evaluated.We found that the main cooling rates of thermal electrons on 7 October 1957, are electron-ion Coulomb collisions, vibrational excitation of N 2 and O 2 , and rotational excitation of N 2 .The relative role of the cooling rates of thermal electrons by low-lying electronic excitation of O 2 (a 1 g ) and O 2 (b 1 g + ), from rotational excitation of O 2 , in collision of O( 3 P) with thermal electrons with the O( 1 D) formation, and by the atomic oxygen fine structure excitation, is negligible in comparison with the effects of the main cooling rates on the electron temperature for the geomagnetically quiet period on 7 October 1957.

Effects of corrections in the NRLMSISE-00 model [O]
or [N 2 ] and [O 2 ] on the ionosphere.
Figures 3-8 show that the calculated NmF2 is systematically higher than the measured one during most of the studied time period.We can expect that the neutral models have some inadequacies in predicting the number densities with accuracy, and we have to change the number densities by correction factors at all altitudes, to bring the modeled electron densities into agreement with the measurements.
The comparison between the measured (squares) and modeled (lines) NmF2 and hmF2 latitude variations are shown in Fig. 9, at 17:00 UT (two upper panels) and 19:00 UT (two lower panels) on 7 October 1957 and in Fig. 10   It is found from numerical simulations of the mid-latitude ionosphere that the daytime magnitude of NmF2 should be reduced by about a factor of 2-3, due to enhanced vibrational excitation of N 2 and O 2 at high solar activity during geomagnetically quiet and storm periods (see Pavlov and Foster, 2001 and references therein).It is apparent from the results of our calculations that the NmF2 decrease caused by the reactions of O + ( 4 S) ions with vibrationally excited N 2 and O 2 is less at low geomagnetic latitudes, in comparison with that at middle geomagnetic latitudes at solar maximum.
The excitation of N 2 and O 2 by thermal electrons pro-vides the main contribution to the values of N 2 and O 2 vibrational excitations, if the electron temperature is higher than about 1600-1800 K at F-region altitudes (Pavlov, 1988(Pavlov, , 1997(Pavlov, , 1998b;;Pavlov and Namgaladze, 1988;Pavlov and Buonsanto, 1997;Pavlov and Foster, 2001).The difference between the vibrational temperature, T N 2V , of N 2 and the neutral temperature, T n , and the difference between the vibrational temperature, T O2V , of O 2 and the neutral temperature increase with increasing electron temperature, and the noticeable differences T N2V T n > 50 − 200 K and T O2V − T n > 50 − 200 K, are realized at hmF2, if T e >1700- .Observed (squares) and calculated (lines) hmF2 and NmF2 at 21:00 UT (two upper panels) and 23:00 UT (two lower panels) on 7 October 1957.The difference between the universal time and the local time at the geomagnetic equator is 05:12.The curves are the same as in Fig. 9.
1800 K at F-region altitudes.Thus, as a result of low electron temperatures at F-region altitudes of the low latitude ionosphere on 7 October 1957, the values of T N2V and T O2V are close to T n , while the differences between T n and the middle latitude T N2V and T O2V are noticeable during daytime conditions at solar maximum (see Pavlov et al., 1999;Pavlov and Foster, 2001;Pavlov et al., 2001, and references therein).This is the first reason which explains the weaker decrease in the low-latitude NmF2 due to N 2 (v) and O 2 (v), in comparison with that at middle geomagnetic latitudes.The effects of N 2 (v) and O 2 (v) on middle-latitude NmF2 were usually evaluated by comparing the measured and modeled NmF2 using the original (Richards, 1991) or modified (Pavlov and Buonsanto, 1997) Richards method, when model plasma drift velocities, caused by neutral winds, are found from the agreement between measured and modeled hmF2.As a result, the middle-latitude model, including vi-brationally excited N 2 and O 2 in the loss rate of O + ( 4 S) ions produces hmF2 values very close to hmF2 produced by the middle-latitude model without including vibrationally excited N 2 and O 2 .The low-latitude model described in Appendix A uses the HWW90 thermospheric wind model of Hedin et al. (1991) to calculate the thermospheric wind components and the corresponding plasma drift velocities along magnetic field lines.As Figs. 9 and 10 show, the lowlatitude model, including N 2 (v) and O 2 (v) in the loss rate of O + ( 4 S) ions produces hmF2 with values higher than hmF2 produced by the low-latitude model without including N 2 (v) and O 2 (v).As a result of including N 2 (v) and O 2 (v) in the loss rate of O + ( 4 S) ions, the equatorial F2-layer is lifted to great heights, where the loss rate of O + ( 4 S) ions is decreased, and this leads to an increase in NmF2, which is masked by the general decrease in NmF2 due to vibrationally excited N 2 and O 2 .In other words, this additional increase in hmF2 de-creases the effects of vibrationally excited N 2 and O 2 on the electron density in the low-latitude ionosphere.

Conclusions
A new theoretical model of the Earth's low-and middlelatitude ionosphere and plasmasphere has been developed.The new model uses a new method in ionospheric and plasmaspheric simulations which is a combination of the Eulerian and Lagrangian approaches in model simulations.The electron and ion continuity, and energy equations, are solved in a Lagrangian frame of reference which moves with an individual parcel of plasma, with the local plasma drift velocity perpendicular to the magnetic and electric fields.As a result, only the time-dependent, one-dimensional electron and ion continuity and energy equations are solved in this Lagrangian frame of reference.The new method makes use of an Eulerian computational grid, which is fixed in space coordinates and chooses the set of the plasma parcels at every time step, so that all the plasma parcels arrive at points which are located between grid lines of the regularly spaced Eulerian computational grid at the next time step.The solution values of electron and ion densities, and temperatures at the Eulerian computational grid, are obtained by interpolation.
Dipole orthogonal curvilinear coordinates q, U , and are used, where q is aligned with, and U and are perpendicular to, the magnetic field, and the U and coordinates are constant along a dipole magnetic field line.Equations A11-A14, and Eqs.A15-A17, which determine the trajectory of the ionospheric plasma perpendicular to magnetic field lines and the moving coordinate system, are derived.It follows from these equations that time variations of U , caused by the existence of the E component of the electric field, are determined by time variations of the component, E eff , of the effective electric field and time variations of , caused by the existence of the E U component of the electric field, are determined by time variations of the U component, E U eff , of the effective electric field.It is shown that the magnetic field lines are "frozen" in the ionospheric plasma, if the values of E eff and E U eff are not changed along magnetic field lines, and there is the interdependency given by Eq.A17 between changes in E U eff in the direction and changes in E eff in the U direction.
The Eulerian computational grid used consists of a distribution of the dipole magnetic field lines in the ionosphere and plasmasphere.One hundred dipole magnetic field lines are used in the model for each fixed value of the geomagnetic longitude.The number of the fixed nodes taken along each magnetic field line is 191.For each fixed value of the geomagnetic longitude, the region of study is a (q, U ) plane which is bounded by two dipole magnetic field lines.The low boundary dipole magnetic field line has the apex altitude of 150 km.The upper boundary dipole magnetic field line has the apex altitude of 4 491 km and intersects the Earth's surface at two geomagnetic latitudes: ±40 • .The Eulerian computational grid dipole magnetic field lines are distributed between these two boundary lines.
The model takes into account the offset between the geographic and geomagnetic poles.The horizontal components of the neutral wind, which are used in calculations of the wind-induced plasma drift velocity along the magnetic field, are specified using the HWW90 wind model of Hedin et al. (1991).In the model, time-dependent ion continuity equations for the three major ions, O + ( 4 S), H + and He + and for the minor ions, NO We have presented a comparison between the modeled NmF2 and hmF2, and NmF2 and hmF2 which were observed at the anomaly crest and close to the geomagnetic equator simultaneously by the Panama, Bogota, Talara, Chiclayo, and Huancayo ionospheric sounders during the 7 October 1957 geomagnetically quiet time period at solar maximum, near approximately the same geomagnetic meridian of 351.9 • .To complete the picture of the latitude dependence of NmF2 and NmF2 variations, we compare the modeled NmF2 and hmF2 at the geomagnetic longitude of 351.9 • with NmF2 and hmF2 measured by the Puerto Rico ionosonde station with geomagnetic longitude of 2.8 • .
A two-peaked structure in the time dependence of the equatorial vertical E × B drift velocity is given by the model of Scherliess and Fejer (1999) at solar maximum during quiet daytime equinox conditions.It leads to a two-peaked structure in the time dependence of the equatorial value of E .The model results highlight the relationship between local time variations of the low-latitude electron densities and the equatorial value of E .The model calculations show that there is a need to revise the model dependence of the equatorial E in local time by elevating and displacing the morning peak to earlier times, and by compressing the time of the pre-reversal peak.It is found that the large disagreement between the measured and modelled hmF2 at 00:00 UT on 7 October 1957 (at 18:48 LT on 6 October 1957) is caused by the long time duration of the prereversal strengthening of the equatorial upward E × B drift given by Scherliess and Fejer (1999).The long period of the pre-reversal enhancement in E leads to unreal high modeled F2 peak altitudes at 00:00 UT.Our calculations provide evidence that, to bring the measured and modeled F2-region main peak altitudes into agreement, the magnitude of E has to be approximately constant in the time range between 15:00 LT and 18:00 LT with the following peak in E , which has a shorter time width in comparison to the time duration of the pre-reversal strengthening of the original equatorial perpendicular plasma drift given by Scherliess and Fejer (1999).The modification of E (t ge ) was carried out in the time range between 07:27 LT and 11:00 LT, to mimic the depth of the anomaly crest and trough in the modeled and measured NmF2.This modification includes the strengthening of E and the time shift of the first peak in E (t ge ) relative to the first peak in the original equatorial E (t ge ) found from the equatorial perpendicular plasma drift of Scherliess and Fejer (1999).The first maximum (E = 1.1 mVm −1 ) of the modified E (t ge ) occurs at 08:30 LT, while the first maximum (E = 0.66 mVm −1 ) of the original equatorial E (t ge ) is located between 10:00 LT and 11:00 LT.The northern depth of the equatorial NmF2 trough in the calculated NmF2 is approximately consistent with the measured depth, if the modified E (t ge ) is used.The model with the modified value of E (t ge ) produces the onset of the equatorial anomaly crest formation close to 15:00 UT, in agreement with the measured onset of the equatorial anomaly crest formation given by the ionosonde stations.
Electron and ion densities and temperature uncertainties resulting from the difference between the NRLMSISE-00 and MSIS-86 neutral temperatures and densities, and from the difference between the EUV97 and EUVAC solar fluxes, are evaluated.Our calculations show that the best agreement between the measured and modeled electron densities is obtained if the MSIS-86 neutral densities and temperature model is used, in combination with the EUVAC solar flux model.We found that the electron and ion temperature uncertainties caused by these differences are negligible.
The thermal electron energy budget in the low-latitude ionosphere at solar maximum was examined.It is shown that the daytime peak values in the electron and ion temperatures, which occur at about 20:32 UT-21:22 UT above the ionosonde stations, result from the peak in the neutral temperature at F2-region altitudes, which occurs very close to the time of the peaks in the electron and ion temperatures.Our calculations show that the values of the electron temperatures at F2-region altitudes become almost independent of the electron heat flow along the magnetic field lines above the Huancayo, Chiclayo, and Talara ionosonde stations because the near-horizontal magnetic field inhibits this heat flow of electrons.The increase in geomagnetic latitude leads to an increase in the effects of the electron heat flow along the magnetic field line on T e .
At sunrise, there is a rapid heating of the ambient electrons by photoelectrons, and the difference between the electron and neutral temperatures could be increased because nighttime electron densities are less than those by day, and the electron cooling during morning conditions is less than that by day.This expands the altitude region at which T e > T i near the equator and leads to the sunrise electron temperature peaks at hmF2 altitudes above the ionosonde stations.After the abrupt increase at sunrise, the electron temperature decreases, owing to the increasing electron density, due to the increase in the cooling rate of thermal electrons and due to the decrease in the relative role of the electron heat flow along the magnetic field line, in comparison with cooling of thermal electrons.As a result, the morning electron temperature peaks, which are found above the ionosonde stations at hmF2 altitudes, are explained by these physical processes.
The relative role of the electron cooling rates was evaluated.We found that the main cooling rates of thermal electrons on 7 October 1957 are electron-ion Coulomb collisions, vibrational excitation of N 2 and O 2 , and rotational excitation of N 2 .The relative role of the cooling rates of thermal electrons by low-lying electronic excitation of O 2 (a 1 g ) and O 2 (b 1 g + ), from rotational excitation of O 2 , in collision of O( 3 P) with thermal electrons with the O( 1 D) formation, and by the atomic oxygen fine structure excitation, is negligible, in comparison with the effects of the main cooling rates on the electron temperature for the geomagnetically quiet period on 7 October 1957.
The model of the ionosphere and plasmasphere was able to reproduce F-region main peak electron densities and altitudes observed on 7 October, if modified NRLMSISE-00 [O] or [N 2 ] and [O 2 ] are used.We found that it is necessary to decrease the NRLMSISE-00 model [O]/[N 2 ] ratio by a factor of 1.7-2.1 from 16:12 UT to 23:12 UT, to bring the modeled and measured NmF2 and hmF2 into satisfactory agreement.This result indicates that the NRLMSISE-00 model may need improvements during geomagnetically quiet periods in equinox at solar maximum at low latitudes.
The increase in the loss rate of O + ( 4 S) ions due to the vibrationally excited N 2 and O 2 leads to the decrease in the calculated NmF2 by a factor of 1.06-1.44 and to the increase in the calculated hmF2, up to the maximum value of 32 km in the low-latitude ionosphere between -30 • and +30 • of geomagnetic latitude.Inclusion of vibrationally excited N 2 and O 2 brings the model and data into better agreement.

A1.1 Ion continuity equations
The model used is a three-dimensional, time-dependent model of the ionosphere and plasmasphere that uses a dipole approximation to the Earth's magnetic field and takes into account the offset between the geographic and geomagnetic axes.The model includes the coupled ion continuity equations for the three major ions, O + ( 4 S), H + , and He + , and the ion continuity equations for the minor ions, NO + , O + 2 , and N + 2 , which can be written as: where N i is the ion concentration, B=|B| is the absolute value of the geomagnetic field whose magnitude can be calculated as B=B O (R E /R) 3 (1 + 3cos 2 ) 1/2 , =90 • -ϕ is the geomagnetic co-latitude, ϕ is the geomagnetic latitude, R E is the Earth's radius, R is the radial distance from the Earth's center, B 0 is the equatorial value of B for R=R E and =0, t is a local time, S is the distance along the magnetic field line, positive in the direction north to south, P i and L i are the production rates and the ion loss rates by the chemical reactions, Q i and Q i are the production rates of ions by photoionization and due to photoelectrons, C i = V i + W , V i is the field-aligned diffusion velocity, W is the field-aligned windinduced plasma drift velocity, with the total time-derivative is defined by The model also includes O + ( 2 D), O + ( 2 P), O + ( 4 P), and O + ( 2 P * ) ions whose number densities are obtained from local chemical equilibrium (see Sect. 2.2).The value of the electron number density is calculated as the sum of the ion number densities: N e = i N i .
The ion diffusion velocities are calculated by solving the system of equations given by Pavlov (1997).The windinduced plasma drift velocity along the magnetic field is determined as W = U cos I , where U is the magnetic meridianal component of the thermospheric wind in spherical polar geomagnetic coordinates (R, , ), is the geomagnetic longitude, I is the magnetic field dip angle, cos I = sin (1 + 3 cos 2 ) −1/2 .The model takes into account the offset between the geographic and geomagnetic poles and calculates U = U (g) cos D − U (g) sin D (see Bailey and Balan, 1996), where D is the magnetic declination angle, (g) = 90 • − ϕ(g) is the geographic co-latitude, ϕ(g) is the geographic latitude, (g) is the geographic longitude, U (g) and U (g) are the horizontal components of thermospheric wind in spherical polar geographic coordinates, which are positive in the southward and eastward directions, respectively.The magnitudes of U (g) and U (g) are obtained from the thermospheric wind components given by the HWW90 thermospheric wind model (Hedin et al., 1991).To calculate the magnetic declination angle, we use the approach described in detail by Bailey and Balan (1996).
The model includes the production and loss rates of O + ( 4 S), NO + , O 2 + , N 2 + , H + , and He + ions by the chemical reactions described in detail by Pavlov and Foster (2001), except for the dissociative recombination rate coefficient for O + 2 ions, whose value is taken from Peverall et al. (2001) Pavlov, 1998b) using the photoionization and photoabsorption cross sections of N 2 , O 2 , and O compiled by Richards et al. (1994) and the photoionization cross section for He presented by Samson et al. (1994).The total and partial atomic oxygen photoionization cross sections of Richards et al. (1994) were updated using the measurements of these cross sections given by Schaphorst et al. (1995), andBerkowitz (1997) and compiled by Avakyan et al. (2000).Additional production rates of O + ( 4 S), O + ( 2 D) and O + ( 2 P) ions are obtained in the model by inclusion of O + ( 4 P) and O + ( 2 P * ) ions using photoionization cross sections given by Richards et al. (1994).The O + ( 4 P) state decays to O + ( 4 S) ions and O + ( 2 P * ) ions decay to either O + ( 2 D) or O + ( 2 P) with a branching ratio of 2:1 for O + ( 2 D): O + ( 2 D) (Pavlov and Foster, 2001).Ionization of O by photoelectrons produces little extra O + ions.The electronically excited oxygen ions are converted to unexcited O + ( 4 S) ions, and N + 2 and O + 2 ions, by chemical reactions that are included in the model of the ionosphere and plasmasphere (see Table A1).As a result, the total production rate of unexcited O + ( 4 S) ions (Q tot ( 4 S) is calculated in the model of the ionosphere and plasmasphere as where Q( 4 S) and Q( 4 S) are the production rates of Q( 4 S) ions by photoionization and due to photoelectrons, Q( 4 P) and Q( 4 P) are the production rates of Q( 4 P) ions by photoionization and due to photoelectrons, K 1 , K 3 , K 5 , K 7 , K 10 are the rate coefficients of the chemical reactions displayed in Table A1 McFarland et al. (1974) and Park and Banks (1974), respectively.

A1.2 Number densities of O
The difficulty of preparing O + ( 2 D) and O + ( 2 P) metastable states leads to the difficulty in the determination of the rate coefficients of the reactions involving these ions.In the experiment by Johnsen and Biondi (1980a, b), metastable O + ( 2 D) and O + ( 2 D) ions are prepared by dissociative charge transfer: He + + O 2 → He + O+O + ( 4 S, 2 D, 2 P).Johnsen and Biondi (1980 a, b) reported rate constants of 8 • 10 −10 cm 3 s −1 and 7 • 10 −10 cm 3 s −1 for chemical reactions 8 and 9 of Table A1, respectively, assuming that the reaction He + + O 2 produces only O + ( 4 S) and O + ( 2 D) ions with the predominance of O + ( 2 D) ions.The branching ratios for O + ( 4 S), O + ( 2 D), and O + ( 2 P) ions formed in the reaction He + + O 2 have been measured as a function of the center-of-mass kinetic energy, E, by Bischof and Linder (1986) and by Gerlich (1991).Contrary to the assumption of Johnsen and Biondi (1980a, b), more than 60% of the O + ions formed at thermal energies are found to be in the O + ( 2 D) state (Bischof and Linder, 1986;Gerlich, 1991).As a result, the rate coefficients of the O + ( 2 D)+N 2 and O + ( 2 D) + O 2 reactions given by Johnsen and Biondi (1980 a, b) include a large error.Li et al. (1997) have measured the dependence of the cross sections for reactions 2 and 8 of Table A1 on E, in the energy range 0.006-40 eV, employing the differential retarding potential method.They found that these cross sections are nearly independent of E, and estimated the values of the charge transfer rate constants of reactions 2 and 8 of Table A1 at thermal energies as 2.0 • 10 −10 cm 3 s −1 and 1.5 • 10 −10 cm 3 s −1 , respectively.
The rate coefficient for an ion-neutral chemical reaction can be calculated from the measured cross section, σ (E), using the relation (St.-Maurice and Torr, 1978) where u is the relative speed between the reactants, the Maxwellian distribution of relative speeds is calculated as and m n denote the masses of the ion and neutral reactants, respectively, a relative drift velocity between ion and neutral reactants.
It follows from Eq. A2, that if the cross section for a reaction is approximately constant, then this reaction rate coefficient is approximated by the relation of K(T ) = σ {8kT (π µ) −1 } 0.5 .The measured cross sections of reactions 2 and 8 of Table A1 are nearly independent of E in the energy range 0.006-40 eV (Li et al., 1997).As a result, we can find the rate coefficients K 2 and K 8 of reactions 2 and 8 of Table A1 as K 2 = 2.0 • 10 −10 (T /300) 0.5 cm 3 s −1 , K 8 = 1.5 • 10 −10 (T /300) 0.5 cm 3 s −1 .These values of K 2 and K 8 were used by Pavlov and Foster (2001) without the explanations presented above.
The rate coefficients K 8 and K 9 of the O + ( 2 D) + N 2 and O + ( 2 D) + O 2 reactions (the reactions 8 and 9 of Table A1) given by Johnsen and Biondi (1980a, b) include the same error due to the same source of O + ( 2 D) ions in the reaction He + + O 2 , and, therefore, it is possible to assume that the K 9 /K 8 ratio given by Johnsen and Biondi (1980a, b) as 7/8 is close to the correct value.As a result, we believe that K 9 = 1.3 • 10 −10 (T /300) 0.5 cm 3 s −1 .A1.2.2 Quenching of O + ( 2 D) and O + ( 2 P) by electrons The metastable O + ( 2 D) and O + ( 2 P) electron quenching rate coefficients have not been measured in the laboratory.The rate coefficients for these chemical reactions can be calculated as K ij (T e ) = ∞ 0 σ ij (u)uf (u)du , where σ ij is the cross section for the transition i → j , u is the relative speed between O + ( 2 D) or O + ( 2 D) ions and electrons, the Maxwellian distribution of relative speeds is calculated as , where E = m e u 2 /2, g i = 2i + 1 is the statistical weight of the i-th level, a 0 is the Bohr radius, and R y is the Rydberg constant.The effective collision strength, Q ij (T e ), is determined as Q ij (T e ) = ∞ 0 ij (x)exp(−x)dx, where x = E(kT e ) −1 .As a result, the rate coefficients for quenching of O + ( 2 D) and O + ( 2 P) by electrons may be obtained from the effective collision strengths as The values of Q ij (T e ) are nearly independent of T e in the thermal en-ergy range of electrons (Henry et al., 1969, McLaughlin andBell, 1998), and we conclude that K ij (T e ) ≈ const T −0.5 e .The updated rate coefficients can be calculated from the updated effective collision strengths of McLaughlin and Bell (1998) as K 5 = 2.5 • 10 −8 (300/T e ) 0.5 cm 3 s −1 , K 6 = 7.0 • 10 −8 (300/T e ) 0.5 cm 3 s −1 , K 10 = 4.0 • 10 −8 (300/T e ) 0.5 cm 3 s −1 .
It should be noted that the rate coefficients K 5 , K 6 , and K 10 of reactions 5, 6 and 10 of Table A1, which are used in the current models of the ionosphere and plasmasphere (e.g.Pavlov and Foster, 2001), were presented by Torr and Torr (1982).These rate coefficients were found using the effective collision strengths of metastable O + ions calculated by Henry et al. (1969).The updated rate coefficients are less than those given by Torr and Torr (1982) by a factor of 1.9-2.1.A1.2.3 Steady-state number density of O + ( 4 P), O + ( 2 P * ) and O + ( 2 P) ions The O + ( 4 P) state decays very promptly to O + ( 4 S) ions with the measured radiative lifetime of O + ( 4 P) ions as τ ( 4 P) = (1.26± 0.10) • 10 −9 s (Smith et al., 1971).The radiative lifetime of O + ( 2 P * ) ions is found by Tayal and Richardson (2000) as τ ( 2 P * ) = 1.4 • 10 −10 s.Using the measured oscillator strengths of transitions among states of O + ions and measured excitation energies of O + states, compiled by Tayal and Richardson (2000), Pavlov and Foster (2001) found that O + ( 2 P * ) ions decay to either O + ( 2 D) or O + ( 2 D) with a branching ratio of 2:1 for O + ( 2 D): O + ( 2 P).The value of τ ( 4 P) or τ ( 2 P * ) is much less than the studied characteristic times of changes in densities of ions and electrons in the ionosphere and plasmasphere.As a result, the model of the ionosphere and plasmasphere uses the steady-state number densities of O + ( 4 P) and O + ( 2 P * ) as The loss rate of O + ( 2 D) ions is determined by the reactions (1)-( 6) of Table A1.
The characteristic time, τ ( 2 P), of O + ( 2 D) number density changes in these chemical reactions is given as , where K 1 − K 6 are the rate coefficients of the chemical reactions (1)-( 6), respectively, displayed in Table A1.The value of τ ( 2 P) is less than 2.78 s, i.e. this time constant is much less than the studied characteristic times of changes in densities of ions and electrons in the ionosphere and plasmasphere.As a result, the steady-state number density of O + ( 2 D) ions is used in the model calculations as where Q( 2 P) and Q( 2 P) are the production rates of O + ( 2 P) ions by photoionization and due to photoelectrons, respec- -K 10 are the rate coefficients of the chemical reactions ( 4) and (6)-10), respectively, displayed in Table A1.

A2.2. Energy balance equations
To determine the temperatures T i and T e of ions and electrons, we use the energy balance equations for ions and electrons given by Pavlov (1997) and add in these equations the terms which take into account the drift of plasma perpendicular to the magnetic field line (see, for example, Bailey and Balan, 1996): where k is the Boltzmann's coefficient, m i denotes the mass of the i-th ion, mn denotes the mass of the nth neutral component of the atmosphere, ν ie and ν in are the collision frequencies for momentum transfer between ions and electrons and between ions and neutrals, , ω i is the ion cyclotron frequency, E ⊥ is the perpendicular component of the electric field with respect to the geomagnetic field, λ i and λ e are the thermal conductivities of ions and electrons, L el is the electron cooling rate in the process "l", P e is the heating rate of the electron gas by photoelectrons, P rc is an additional heating rate of the electron gas due to Coulomb collisions between ring current ions and plasmaspheric electrons and wave-electron interactions (the value of P rc = 0 is used in the calculations presented in this work), the total timederivatives are defined by d dt T i = ∂ ∂t T i + V E • grad T i and by • grad T e , the value of C i is the same as in Eq. (A1), and to calculate the field-aligned electron velocity, C e , we assume that there are no field-aligned currents, i.e.C e = i C i N i /N e .
We use the same equations for ν ie and ν in as given by Bailey and Balan (1996), except for the O + − O frequency, whose value is taken from Pesnell et al. (1993).The model uses the generally accepted electron cooling rates due to electron-ion Coulomb collisions and elastic collisions of electrons with N 2 , O 2 , O, He and H presented by Schunk and Nagy (1978), and the thermal electron impact excitation of O 2 (a 1 g ) and O 2 (b 1 + g ) given by Prasad and Furman (1973).The revised electron cooling rates by vibrational and rotational excitation of O 2 and N 2 of Pavlov (1998 a, c), the expression for atomic oxygen fine structure cooling rate of thermal electrons derived by Pavlov and Berrington (1999), and the thermal electron cooling rate due to electron-impactinduced transitions 3 P → 1 D for atomic oxygen (Lobzin et al., 1999) are included in the model.The model uses the updated expression for λ e given by Pavlov et al. (2000Pavlov et al. ( , 2001) ) and the expression for λ i given by Hochstim (1969).
The heating rate of electrons by photoelectrons is calculated by the use of the approach of Hoegy (1984) (the same approach was derived from the kinetic equation by Krinberg and Tachilin (1984)).The value of P e is a function of the photoelectron flux.
We split the studied ionosphere and plasmasphere region into two regions.The first region includes the plasmasphere and ionosphere with the magnetic field lines which intersect the point S=0 above 1500 km, i.e. the magnetic field lines have the apex altitudes h ap ≥ 1500 km in the first region.Modelled electron heating due to photoelectrons is provided by a solution of the Boltzmann equation for photoelectron flux along a centered-dipole magnetic field line, using the method of Krinberg and Tachilin (1984) on the same field line grid which is used in solving for the electron and ion temperatures.In the altitude range 130-700 km in the Northern and Southern Hemispheres, the model solves the Boltzmann equation for photoelectron flux using the method of Krinberg and Tachilin (1984) and the updated elastic and inelastic cross sections of the neutral components of the atmosphere described by Lobzin et al. (1999) and calculates the value of P e .In the approach of Krinberg and Tachilin (1984), photoelectron transport and loss processes due to the elastic and inelastic collisions of electrons with neutral components of the atmosphere and Coulomb electron-electron collisions are taken into account.Their technique is based on the solution of simpler transport flux equations derived from the Boltzmann equation and determines the ionospheric electron heating rate with an accuracy of about 10%, in comparison with that obtained by the numerical solution of the kinetic equation (Krinberg and Tachilin, 1984).Above 700 km, the energy lost by photoelectrons in heating the plasma is calculated using the analytical equation for the plasmaspheric transparency, P (E), (Krinberg and Matafonov, 1978;Krinberg and Tachilin, 1984) that determines the probability of the magnetically trapped photoelectrons with an energy, E, of entering the magnetically conjugated ionosphere.The transparency depends mainly on a single parameter proportional to the Coulomb cross section and the total content of electrons in the plasmasphere magnetic flux tube (the transparency approaches unity as photoelectrons pass through the plasmasphere without significant absorption, and P(E) = 0, if photoelectrons are absorbed by the plasmasphere).This analytical transparency approach determines the plasmaspheric heating rate with an accuracy of about 10-20% in comparison with that obtained by the numerical solution of the kinetic equation (Krinberg and Matafonov, 1978;Krinberg and Tachilin, 1984).Below 700 km, the photoelectron transfer equations, which were derived by Krinberg and Taschilin (1984) from the kinetic equation for the photoelectron flux, are solved in both hemispheres.This approach is used in the models of the plasmasphere and ionosphere (Krinberg and Taschilin, 1984;Pavlov and Foster, 2001;Korenkov et al., 1996).
The second region is the ionosphere region with the magnetic field lines that intersect the apex altitudes below 1500 km.To calculate the value of P e in the second region, we solve the equations which were derived by Krinberg and Taschilin (1984) from the kinetic equation for the photoelectron flux below the points S = 0 in both hemispheres, i.e. the transparency approach is not used.
Time dependent continuity and energy equations, which determine [N 2 (v)] and [O 2 (v)], are presented by Pavlov (1997) and Pavlov (1998b).The model uses Boltzmann distributions of N 2 (v) and O 2 (v) as where E 1 = 3353 K is the energy of the first level of N 2 given by Radzig and Smirnov (1980) v)] are the total number densities of N 2 and O 2 .It follows from Eq. (A25) that ) where the values of [N 2 (0)] and [O 2 (0)] are calculated from [N 2 ] and [O 2 ] (whose values are given by the MSIS-86 model of Hedin (1987) or by the NRLMSISE-00 model of Picone et al. (2000Picone et al. ( , 2002))) by The values of α and β are determined by solving the time-dependent, one dimensional quanta energy equations given by Pavlov (1997Pavlov ( , 1998b)).The quanta energy equation solution procedure is described in detail by Pavlov (1997Pavlov ( , 1998b)).

A2.4 Neutral temperature and densities and solar EUV fluxes
The model of the ionosphere and plasmasphere includes an option to use two sets of the models of the neutral temperature and densities.First, one uses the NRLMSISE-00 neutral temperature and densities given by Picone et al. (2000Picone et al. ( , 2002)), and second, one uses the MSIS-86 neutral temperature and densities model of Hedin (1987).Both neutral temperature and density models are run using 3-h geomagnetic Ap indices.To calculate the density of NO, the model given by Titheridge (1997) is used.The model of the ionosphere and plasmasphere can use the solar EUV fluxes from the EUVAC model (Richards et al., 1994) or the EUV97 model (Tobiska and Eparvier, 1998).At night our model includes the neutral ionization by scattered solar 121.6, 102.6 and 58.4 nm fluxes, as described by Pavlov (1997).
A2.5 Solution of continuity and energy equations A2.5.1 Dipole coordinates The coordinate system considered here is similar to that of Anderson (1973) and Rasmussen et al. (1993).Orthogonal curvilinear coordinates are: q = (R E /R) 2 cos , U = (R E /R) sin 2 , and a geomagnetic longitude, .The important properties of these coordinates are that q is aligned with, and U and are perpendicular to, the magnetic field, the U and coordinates are constant along a dipole magnetic field line, and the McIlwain parameter L = R/(R E sin 2 ) can be presented as L = U −1 .Also, for a dipole magnetic field line Bailey and Balan, 1996).We take into account that the plasma E × B drift velocity can be presented as is the component of E in the dipole coordinate system, E U is the U component of E in the dipole coordinate system, e and e U are unit vectors in and U directions, respectively.For the plasma trajectory tracking, the definition of the velocity vector can be used as follows: where the coordinate scale factors, h U and h , are given by Rasmussen et al. (1993) as h U = RL cos I and h = R sin , I is the magnetic field dip angle, cos I = sin (1 + 3cos 2 ) −1/2 .Taking into account that B = B 0 (R E /R) 3 (1 + 3 cos 2 ) 1/2 , Eq. (A10) which determines the trajectory of the ionospheric plasma and the moving coordinate system can be rewritten as If the E eff and E eff U components of the effective electric field are changed along magnetic field lines (i.e. the values of E sin 3 and E U sin 3 (1 + 3 cos 2 ) −1/2 are changed along magnetic field lines), then the values of ∂ ∂t U (i.e. ∂ ∂t L) and ∂ ∂t are not constants, and magnetic field lines are not "frozen" in the ionospheric plasma.However, the ionospheric plasma is assumed to be with a magnetic field "frozen" in above about 150 km, where the drift velocities of ions and electrons perpendicular to the geomagnetic and electric fields are approximately the same (Ratcliffe, 1956).To overcome this difficulty, it is necessary to prove that the values of E h and E U h U are not changed along magnetic field lines.To understand the physics, there must first be a clear understanding of the concept of "magnetic line" preservation.The concept of magnetic line preservation or moving magnetic field lines gives rise to the mnemonic of magnetic field lines as entities with an integrity that is locally convected by plasma flow fields, which are locally orthogonal to magnetic field lines (Alfven and Fälthammer, 1963).Magnetic field lines are carried about with any plasma flow fields whose motions are perpendicular to magnetic field lines.For the magnetic field line to be "frozen" to the perpendicular motion of plasma, the magnetic field diffusion is assumed to be negligible, and the evolution of the magnetic field is governed by the induction equation as (Alfven and Fälthammer, 1963) ∂ ∂t B = rotB × V E .We have to take into account that As a result, the electric field restrictions which follow from the fact that the magnetic field lines are frozen in the ionospheric plasma can be formulated from rotB × V E = 0 as These equations yield the frozen-in-field conditions in the ionosphere and plasmasphere above about 150 km.It follows late the values of N i (q, U, , t + t), T i (q, U, , t + t) and T e (q, U, , t + t) simultaneously for all computational grid dipole magnetic field lines.The subsequent strategy for the case when E >0 is different from that for the case when E <0.
It is well known that during most of the daytime conditions, the value of E >0, and this electric field leads to an upward drift of plasma above the geomagnetic equator, producing an equatorial plasma fountain (see, for example, Moffett, 1979;Anderson, 1981;Rishbeth, 2000;Bailey and Balan, 1996).In this case, the plasma moves along geomagnetic field lines and perpendicular to magnetic field lines along the U coordinate from the low boundary line (the number of this grid line is k = 1) to the upper boundary line (let us assume that the number of this grid line is k = kk).We define the (q, U) coordinate at the k-th magnetic field grid line as the ((q k , U k ) coordinates.Using the values of N i (q k , U k , , t), T i (q k , U k , , t) and T e (q k , U k , , t) at the k-th grid line, we calculate the values of N i (q k , U k + U k , , t + t), T i (q k , k + U k , , t + t) and T e (q k , U k + U k , , t + t) in the moving Lagrangian frame by solving the one-dimensional, time-dependent Eqs.(A1), (A6) and (A7) in this frame.Now we want to recalculate these results to our Eulerian computational grid.The value of U k at each time step, t, is determined from Eqs. (A11), (A12), and (A15).It follows from these equations that the value of U k is the same for different departure points of the computational grid dipole magnetic field line, i.e. the arriving points follow within the single grid magnetic field line.In addition to the first set of the calculated values of N i (q k , U k + U k , , t + t), T i (q k , k + U k , , t + t) and T e (q k , U k + U k , , t + t), we have the second set of the calculated values of N i (q k−1 , U k−1 , , t + t), T i (q k−1 , U k−1 , , t + t) and T e (q k−1 , U k−1 , , t + t) at the neighboring underlying computational grid dipole magnetic field line.The value of q k−1 does not coincide with the value of q k and, thus, a search-interpolation procedure is needed to calculate the third set of N i (q k , U k−1 , , t + t), T i (q k , U k−1 , , t + t) and T e (q k , U k−1 , , t + t) (these values of N i , T i and T e correspond to the neighboring, underlying computational grid dipole magnetic field line as well) from the second set of the calculated N i , T i and T e .Using the first and the third sets of the calculated N i , T i and T e and the interpolation procedure, we calculate the desired quantities of N i (q k , U k , , t + t), T i (q k , U k , , t + t) and T e (q k , U k , , t + t).
To put into practice this strategy, it is necessary to have the desired quantities of N i (q k , U k , , t + t), T i (q k , U k , , t + t) and T e (q k , U k , , t + t) for k = 1 (the low boundary dipole magnetic field line).For this computational grid dipole magnetic field line with the apex altitude of 150 km, the calculations are carried out without E × B drift velocity by solving the one-dimensional, timedependent Eqs.(A1), (A6), and (A7).It should be noted that there are no differences in the strategy of calculations of N i , T i and T e for the grid line with k = kk and for the grid lines with 1<k<kk.
For each Eulerian computational grid dipole magnetic field line, we can find that q k (min) ≤ q k ≤ 0, where q k (min) is a minimum value of q k .The magnitude of | q k (min) | is increased if the value of k is increased.It means that we cannot use our method to find the value of N i , T i and T e close to both ends of each Eulerian computational grid dipole magnetic field line with k > 1.These grid points are located below 150 km altitude for the Eulerian computational grid point distribution which is used in this study.The calculations of N i , T i and T e for these Eulerian computational grid points are carried out without E × B drift velocity by solving the one-dimensional, time-dependent Eqs.(A1), (A6) and (A7).
If the value of E < 0, then the plasma moves along geomagnetic field lines and perpendicular to magnetic field lines along the U coordinate from the upper boundary grid line, with k = kk to the low boundary grid line with k = 1.Like the previous case, we use the values of N i (q k , U k , , t), T i (q k , U k , , t) and T e (q k , U k , , t) at the k-th grid line, to calculate the values of N i (q k , U k + U k , , t + t), T i (q k , U k + U k , , t + t) and T e (q k , U k + U k , , t + t) in the moving Lagrangian frame by solving the one-dimensional, time-dependent Eqs.(A1), (A6), and (A7) in this frame.After that, the model recalculates the results to our Eulerian computational grid.
Equations (A11), (A12), and (A15) are used to find the value of U k .Like the previous case, the value of U k is the same for different departure points of the computational grid dipole magnetic field line, i.e. the arriving points follow within the single grid magnetic field line.
In addition to the first set of the calculated values of N i (q k , U k + U k , , t + t), T i (q k , U k + U k , , t + t), and T e (q k , U k + U k , , t + t), we have the second set of the calculated values of N i (q k+1 , U k+1 , , t + t), T i (q k+1 , U k+1 , , t + t) and T e (q k+1 , U k+1 , , t + t) at the neighboring overlying computational grid dipole magnetic field line.The value of q k+1 does not coincide with the value of q k and, thus, a search-interpolation procedure is carried out to calculate the third set of N i (q k , U k+1 , , t + t), T i (q k , U k+1 , , t + t) and T e (q k , U k+1 , , t + t) (these values of N i , T i and T e correspond to the neighboring overlying computational grid dipole magnetic field line as well) from the second set of the calculated N i , T i and T e .The first and the third sets of the calculated N i , T i and T e are used to obtain the desired quantities of N i (q k , U k , , t + t), T i (q k , U k , , t + t) and T e (q k , U k , , t + t) by interpolation.
To put into practice this part of our method, it is necessary to have the sought out quantities of N i (q k , U k , , t + t), T i (q k , U k , , t + t) and T e (q k , U k , , t + t) for k = kk (the upper boundary dipole magnetic field line).This computational grid dipole magnetic field line intersects the Earth's surface at middle geomagnetic latitudes (two geomagnetic latitudes of ± 40 • are used in this study).It is well known that, unlike the auroral and equatorial ionosphere, electric fields have little effect on the mid-latitude ionosphere, and even relatively strong electric fields measured by the Millstone Hill radar (43 • N, 288 • E) during a January 1997 magnetic storm had little effect on electron and ion densities (Richards et al., 2000).Therefore, we can suggest that the effects of the equatorial electric field on N i , T i and T e are negligible at middle geomagnetic latitudes.As a result, the model calculations are carried out without E × B drift velocity for k = kk by solving the one-dimensional, time-dependent Eqs.(A1), (A6) and (A7).It is necessary to point out that, in the case of E <0, there are no differences in the strategy of the calculations of N i , T i and T e for the grid line with k = 1 and for the grid lines with 1<k<kk.
It is necessary to determine the boundary conditions to solve the one-dimensional, time-dependent Eqs.(A1), (A6) and (A7) in the first and second parts of our method.At the lower ends of each Eulerian computational grid dipole magnetic field line (which are located at the 130 km altitudes in the Northern and Southern Hemispheres), the diffusion and drift processes of ions and the processes of transfer of thermal energy of electrons and ions by the thermal conductivity and drift are neglected in the model calculations (C i = 0, V E = 0, λ i = 0 and λ e = 0).
The numerical simulations of N i , N e , T i and T e presented in Sect. 4 of this work give practically the same results, if the lower boundary grid line has more a lower apex altitude of 140 km and the upper boundary grid line intersects the Earth's surface at two more high middle-latitude geomagnetic latitudes: ±45 • .This may confirm the validity of the present method for simulations for geomagnetically quiet time periods.If there are regions with E <0 and E >0 and grid lines with E = 0 at the same time then the direction splitting technique is used.The calculations of N i , T i and T e are carried out without E × B drift velocity for both boundary grid lines with k = 1 and k = kk or (and) for grid lines with E = 0, and two different strategies for solving the continuity and energy equations described above are employed at the same time.
A2.6 Comparison between the presented new solution procedure and the CTIP and CTIM model solution procedures in the ionospheric code The semi-Lagrangian approach has long been used in meteorology for numerical weather prediction, where the use of a large time step is essential for efficiency (Smolarkiewicz and Pudykiewicz, 1992).This approach has been introduced by Robert (1981), and the basic idea is to discretize the Lagrangian derivative of the solution in time, instead of the Eulerian derivative.The extension of the semi-Lagrangian method to the solution of Navier-Stokes equations was presented in the pioneering work of Pironneau (1982).As an efficient and accurate approach to computing the advection process, semi-Lagrangian schemes have been extensively studied and widely incorporated into many numerical models for atmospheric flows (Smolarkiewicz and Pudykiewicz, 1992).A sort of semi-Lagrangian technique was developed and employed in Sect.A2.5.2 conformably to the ionosphere and plasmasphere.
Lagrangian and Eulerian frames were used by the early version of the CTIP model (Fuller-Rowell et all., 1988).The CTIP model integrated the Global Thermospheric Model and the Sheffield University High-Latitude Ionospheric Convection Model.Unfortunately, the solution procedure used by the CTIP model for the ionospheric code is described very briefly by Fuller-Rowell et al. (1988) and, as a result, it is not possible to go into details.For geomagnetic latitudes equatorward of 65 degrees, the empirical ionospheric model of Chiu (1975) was used in the CTIP model, while our model is the theoretical model of the ionosphere and plasmasphere for geomagnetic latitudes equatorward of ±40 degrees.For geomagnetic latitudes northward of 65 degrees, the CTIP model uses the ionospheric code with the backward convection path integration algorithm, while the forward integration is performed in the model presented in this work.There are two fundamentally different strategies for solving the continuity and energy equations in our model for E >0 and E <0 in the ionospheric and plasmaspheric code, with the use of the low or upper boundary grid line conditions (see Sect.A2.5.2), while the CTIP model ionospheric code does not depend on the sign of E .It should be noted that an attempt to use low boundary grid line conditions for E <0 leads to the divergence of the numerical scheme presented in this work, and this generates the need of using upper boundary grid line conditions.It is also necessary to emphasise the newness and importance of the derived Eqs.(A15)-(A17) in the ionospheric code.
The CTIP model was enhanced by including the theoretical model of the low-and middle-latitude ionosphere and plasmasphere (Millward et al., 1996).Each plasma flux tube of the updated CTIP model circulates under the influence of the E × B drift such that, over a 24-hour simulation, each plasma flux tube returns exactly to its starting position (Millward et al., 1996, pages 239 and 254).It means that a Lagrangian approach is used in the CTIP model ionospheric code of Millward et al. (1996), and this determines the differences between the solution procedure presented in Sect.A2.5.2 and the solution procedure of the CTIP ionospheric code described by Millward et al. (1996).The CTIM model uses the technique of Fuller-Rowell et al. (1988) in calculating the ionospheric parameters (Fuller-Rowell et al., 1996, page 224).As a result, the differences between the solution procedure used by the CTIM model in the ionospheric code and the solution procedure presented in Sect.A2.5.2 are the same as those described above.
Topical Editor M. Lester thanks B. Emery and two other referees for their help in evaluating this paper.

Fig. 7 .Fig. 8 .
Figures 3-8show that the calculated NmF2 is systematically higher than the measured one during most of the studied time period.We can expect that the neutral models have some inadequacies in predicting the number densities with accuracy, and we have to change the number densities by correction factors at all altitudes, to bring the modeled electron densities into agreement with the measurements.The comparison between the measured (squares) and modeled (lines) NmF2 and hmF2 latitude variations are shown in Fig.9, at 17:00 UT (two upper panels) and 19:00 UT (two lower panels) on 7 October 1957 and in Fig.10at 21:00 UT (two upper panels) and 23:00 UT (two lower panels) on 7 October 1957.Dashed lines show the model results when the original NRLMSISE-00 neutral temperature and densi-

Fig. 9 .
Fig.9 Fig.10 A5) where Q( 2 D) and Q( 2 D) are the production rates of O + ( 2 D) ions by photoionization and by photoelectrons, and τ ( 2 D) is the characteristic time of O + ( 2 D) number density changes in chemical reactions that is determined as divV E ) + P e + P re − l L el , (A7)

Table 1 .
Ionosonde station names and locations ge ) on 7 October 1957.The results presented in panel (d) of Fig. 2 provide evidence that we cannot match the measured and modeled NmF2 and the sizes of the measured and modeled equatorial troughs using the corrections of the NRLMSISE-00 model [O], [N 2 ], or [O 2 ].
+ , O 2 + , and N 2 + , are solved, and the approach of local chemical equilibrium is used to calculate steady-state number densities of O + ( 2 D), O + ( 2 P), O + ( 4 P), and O + ( 2 P * ) ions.The model includes the solution of time-dependent electron and ion energy balance equations.The model uses Boltzmann distributions of N 2 (v) and O 2 (v) to calculate [N 2 (v)] and [O 2 (v)], which are included in the model loss rate of O + ( 4 S) ions and cooling rates of thermal electrons due to vibrational excitation of N 2 and O 2 .
. The reactions of photoionization and ionization by photoelectrons of N 2 , O 2 , O, and He, which form N + 2 , O + 2 , O + ( 4 S), O + ( 2 D), O + ( 2 P), and He + ions, are included in the model (see details in Appendix A of , K 11 and K 12 are the rate coefficients of the chemical reactions N + 2 + O → O + ( 4 S) + N 2 and O + H + → O + ( 4 S) + H given by

Table A1 .
Chemistry of O + ( 2 D) and O + ( 2 D) ions ( a The effective temperature T = (m 1, where m i and m n denote the masses of ion and neutral reactants, respectively, V d is a plasma drift velocity) ) and Q( 2 P * ) are the production rates of O + ( 2 P * ) ions by photoionization and due to photoelectrons, respectively. * Radzig and Smirnov (1980)rgy of the first level of O 2 given byRadzig and Smirnov (1980), T N2v and T O2v are the vibrational temperatures of N 2 and O 2 .The vibrational quanta, α and β, of N 2 and O 2 are determined as α