A revised method to extract thermospheric parameters from incoherent scatter observations

Height distribution of ionospheric plasma parameters in the F2-region is closely related to height distribution of the main thermospheric parameters. Therefore, they can be extracted from ionospheric observations solving an inverse problem of aeronomy. A self-consistent approach to the Ne(h) modeling at the F2-region heights has been applied to solve the problem. Using routine incoherent scatter radar observations (Ne(h), Te(h), Ti(h), Vi(h) profiles) the method yields a self-consistent set of main aeronomic parameters responsible for the F2-region formation. The list of derived parameters includes: neutral temperature profile Tn(h) depending on the exospheric temperature Tex, the temperature at 120 km T120 and the shape-parameter S, which determine the temperature profile, concentration of neutral species [O], [O2], [N2], vertical plasma drift W, which may be converted to the meridional thermospheric wind Vnx, total solar EUV flux and ion composition (O, O2, NO, N2, N) as a result of Ne(h) fitting. Therefore, the method gives a complete description of the upper atmosphere condition in the vicinity of incoherent scatter facility for the periods of observation. Analysis of all available EISCAT (CP-1, CP-2) observations has shown wide deviations from MSIS-86 model predictions for geomagnetically disturbed conditions while the retrieved parameters are close to the model ones for quiet periods. The approach turns out to be very useful for physical analyses of the F2-layer disturbance mechanisms giving a complete picture of the phenomenon in question. Limitations and problems related to method application are discussed. Under existing conditions when thermospheric observations are not conducted currently the proposed method may be considered a real tool for thermosphere investigation and monitoring at least for the periods of ISR observations.

Usually a ionospheric parameter being fitted in model calculations depends on many input aeronomic parameters -some of them are more critical, others are less.It is always possible by varying one of the main aeronomic parameters to fit the ionospheric one, meanwhile the other parameters are taken from outside (e.g., empirical models or observations).Therefore, one obtains a partial solution and this is a special problem to estimate its validity.The more ionospheric variables are simultaneously fitted, the more reliable solution (as a set of aeronomic parameters) is obtained.In this case the extracted aeronomic parameters are not arbitrary but constitute a self-consistent set.Thus we come to the idea of a self-consistent approach to the ionospheric modeling proposed by Stubbe (1973).
Fitting the whole Ne(h) profile opens many possibilities to extract aeronomic parameters as it contains a large amount of information.But this is a complicated problem which requires for its solution an adequate physical model as a tool for analysis, reliable and consistent ionospheric parameter observations and special methods to extract a consistent set of main aeronomic parameters responsible for the observed N e(h) distribution.Nevertheless, attempts are being made in this direction and some progress has been achieved.Zhang et al. (2001) using ISR observations found one to three of five thermospheric parameters by fitting Ne(h) distribution.A self-consistent approach to the F-region modeling by Mikhailov and Schlegel (1997) with later modifications (Mikhailov and Förster, 1999;Mikhailov and Schlegel, 2000) also deals with Ne(h) fitting to retrieve basic aeronomic parameters responsible for the Ne(h) formation in the mid-latitude F2-region.The choice of Ne(h) as a fitted parameter is due to the fact that this quantity is the most reliable and easily measured with the IS technique parameter while other plasma characteristics are subjected to larger uncertainties.For instance, observed Te and much stronger Ti (compared to Ne) depend on the ion composition model applied in the fit of the theoretical to the measured autocorrelation function during the incoherent data analysis.The uncertainty in the actual ion composition may lead to considerable uncertainties in the derived T i(h) and Te(h) profiles especially during disturbed conditions.Plasma velocity observation is a complicated technical problem which is not successfully solved at some ISR facilities (Kharkov, Millstone Hill, for instance).
Recent modifications of the self-consistent method improving its numerical stability, an extension of the method to non-stationary (twilight and nighttime) conditions, its limitations as well as problems encountered are discussed in the paper.Thermospheric parameters extracted from EISCAT (CP-1 and CP-2) and Millstone Hill ISR observations both for quiet and disturbed conditions are presented in comparison with the MSIS-86 thermospheric model.The role of disturbed thermospheric parameters in generating F2-layer disturbances is discussed.

PRESENT DAY STATE-OF-THE-ART
A widespread method in the F2-region physical modeling to obtain observed NmF2 and hmF2 variations is to fit hmF2 by varying vertical plasma drift or meridional thermospheric wind Vnx.While NmF2 is fitted by changing O/N2 ratio (e.g., Richards et al., 1994a;Mikhailov et al., 1995;Pavlov and Foster, 2001) or plasmaspheric flux during nighttime hours (e.g., Mikhailov and Förster, 1999).This gives some useful information about the thermosphere state especially during disturbed periods when empirical thermospheric models turn out to be inefficient.Unfortunately, the obtained thermospheric parameters present only a partial solution as NmF2, hmF2 or TEC also depend on other aeronomic parameters which are fixed to default values during this fitting.Without an independent control of extracted thermospheric parameters the validity of obtained results remains questionable.
As an example recently obtained in this direction let us consider the results by Lilensten and Blelly (2002), who modeled EISCAT observations for June 9, 1994.It was a quiet day (Ap=8) during solar minimum (F10.7=83)although some splashes of electric field up to 20 mV/m took place both at the beginning and at the end of the period.The TRANSCAR first principle model (Blelly et al., 1996) was applied in that analysis.Observed Ne(h) profiles were used to specify NmF2 and hmF2 variations and Ne(h) integration in the 80-425 km height range gave ITEC.It was shown that [O] varying within ± 20% with respect to the MSIS model can provide the observed ITEC and NmF2 variations (fig.4.1).The latter is the expected result as TEC and NmF2 are known to exhibit a pretty good correlation (Das Gupta and Basu, 1973;Kane, 1975;Jakowski et al., 1991)  results can be obtained varying [O2] and [N2], but a correction factor should be larger (1.5-2) in this case, and this looks unreal for a quiet time period.It is stressed that in reality one can hardly expect solely atomic oxygen or molecules variations and a multi-parameter fit is required.So, in quiet geomagnetic conditions when empirical thermospheric models like MSIS give reasonable neutral composition and temperature it is possible to describe fairly well F2-region plasma parameter variations and this has been demonstrated repeatedly (e.g., Pavlov and Buonsanto, 1997;Richards and Wilkinson, 1998;Pavlov et al., 1999Pavlov et al., , 2000;;Pavlov and Oyama, 2000 and references therein), but the situation is more complicated during disturbed periods when large deviations from the MSIS model take place (Buonsanto et al., 1992;Mikhailov and Foster, 1997;Mikhailov and Schlegel, 1998;Litvin et al., 2000;Pavlov and Foster, 2001).However, in the frame of space weather this is a promising and perhaps the only acceptable approach.Indeed, space weather requires real-time (or close to) and global control of the ionosphere and this problem cannot be provided with IS radars.Such service should be based on the use of cheap ionosondes for the F2-layer parameters and the positioning system receivers to estimate the TEC.As mentioned earlier, ISR observations are widely used to obtain information on the thermosphere.There are two generally used methods based on: i) the use of the energy equation for O + ions to find atomic oxygen concentration and neutral temperature in the F2-region (e.g., Evans et al., 1979;Oliver, 1979Oliver, , 1990;;Alcaydé et al., 1982;Blelly et al., 1992;Schoendorf and Oliver, 1998;Litvin et al., 2000), and ii) the use of the momentum equation for O + ions to obtain meridional thermospheric wind (e.g., Lathuillère et al., 1997;Witasse et al., 1998).
Based on the experience obtained since the first publication of the self-consistent method by Mikhailov and Schlegel (1997), some crucial points should be mentioned which are to be taken into account extracting thermospheric information from Ne(h) profiles.
1) Simultaneously measured basic plasma parameters (Ne, Te, Ti, Vi profiles) should be internally consistent and demonstrate a sufficient accuracy.For instance, EISCAT does provide such observations, while not all Millstone Hill ISR data can be used for such analysis due to rare observations or problems with Vi(h) measurements.
2) An optimal set of searched aeronomic parameters should be specified, which is sufficient for describing the observed Ne(h) distribution in the chosen geophysical conditions (day or nighttime, high, middle or equatorial latitudes, etc.).Only these key parameters can be retrieved with an acceptable accuracy from the observed N e(h) profile as their contributions are the largest.
3) Unlike direct F2-layer model calculations when all input parameters can be changed arbitrarily, in case of solving an inverse problem of aeronomy the searched parameters constitute some sort of clusters in which they are self-consistent varying within confined limits.Therefore, the main parameters should be found simultaneously within the same algorithm and this strongly complicates the problem.Usually this is not taken into account, and searching for an aeronomic parameter by fitting an observed ionospheric one, other important parameters are taken from empirical modelsthat is the set of aeronomic parameters turns out to be not self-consistent.Such approach is acceptable to some extent in quiet geomagnetic conditions when thermospheric empirical models like MSIS provide reasonable neutral temperature and concentrations, but it can not be justified for disturbed periods.

METHOD DESCRIPTION
A self-consistent approach to the F-region modeling using ISR observations is described by Mikhailov andSchlegel (1997, 2000), Mikhailov and Förster (1999).The method was successfully used to analyze F2-layer storms (Mikhailov and Foster, 1997;Mikhailov and Schlegel, 1998), ion composition (Mikhailov and Kofman, 2001), and F1-layer disturbances (Mikhailov and Schlegel, 2003).Nevertheless the method is still under development and some modifications to improve its numerical stability as well as its extension for non-stationary conditions are considered in this part of the paper.The sketch illustrates the idea of the method.Routinely observed and specially prepared (see later) Ne(h), Te(h), Ti(h), and Vi(h) profiles are used in the F2-region model to fit calculated Ne(h) to the observed one.The variable parameters are: Tex, T120, shape parameter S for the Tn(h) profile, [N2] concentrations as well as for the total EUV flux.Vertical plasma drift W(h), which can be converted to the meridional thermospheric wind Vnx, is found from the observed Vi(h) (see later).The whole ion composition is available as a result of model calculations.All parameters are found simultaneously constituting a self-consistent set.Let us consider the three items shown in the sketch.

Input ISR data reduction
ISR observations need a special reduction before being used in calculations.The EISCAT CP-1 observations provide range profiles of Ne, Te, Ti and Vi every 5 min (CP-2 every 6 min) with the antenna beam directed along the local geomagnetic field line.Each particular height profile is not smooth exhibiting fluctuations for some reasons and cannot be used for solving the inverse problem.Therefore, median (not average) profiles with Standard Deviations (SD) at each height step are calculated over 1.0-1.5 h time interval, which is close to the e-fold time of the daytime F2-layer.It gives 13-17 values to find the median.At Millstone Hill in some experiments the observations are rare (three per hour) and a 3-4 h time interval has to be used to obtain such median profiles.This is not convenient especially during disturbed periods when changes in the F2-region are pretty fast.These median vertical profiles are then smoothed by a polynomial of up to 5th degree before being used in our model calculations.
Experimental T e(h), Ti(h) and Ne(h) profiles depend on the ion composition used in the fit of the theoretical to the measured Auto-Correlation Function (ACF).An uncertainty in ion composition may lead to considerable uncertainties in the derived Te(h) and Ti(h) profiles and to somewhat smaller uncertainties in N e(h) (e.g., Waldteufel, 1971;Lathuillère et al., 1983).The effect of varying ion composition is most noticeable during disturbed periods, but an appreciable effect may also take place for quiet periods as well (e.g., Mikhailov and Schlegel, 2003).The largest uncertainties take place at the F1-region heights where the ion composition changes from molecular to atomic.Therefore, a correction of the experimental Te(h), Ti(h) and Ne(h) profiles is required.A simple correction proposed by Waldteufel (1971) may be applied when the deviations in ion composition from the model (used in the ISR data analysis) are not large, but in case of strong perturbations this simple correction results in unreal Te(h) and Ti(h) profiles.In such cases a more sophisticated iterative method considered by Mikhailov and Schlegel (1997) should be applied, which provides the proper fit to the measured ACF.Unfortunately, this iterative method requires considerable calculations and cannot be applied routinely.
An example of good EISCAT observations for a disturbed period of August 05, 1992 is given in fig.4.3.Initial median (before polynomial smoothing) profiles calculated over 11:00-12:00 UT time interval are shown.Note small SD in the whole height range telling us about the reliability of observations.

Ionospheric model
A mid-latitude F2-region model includes transport process for O + ( 4 S) and photo-chemical processes only for O + ( 2 D), O + ( 2 P), O2 + (X 2 ′), N + , N2 + and NO + ions in the 150-550 km height range.A two-component model of the solar EUV from Nusinov (1992) is used to calculate the photoionization rates in 48-wavelength intervals (10-1050 Å).To extend the method to twilight and dark hours a nighttime ionization source has been added to the model.Nighttime scattered HLyβ (1026 Å), HeI (584 Å), HeII (304 Å) and radiative recombination emission of O + ions (910 Å) are taken according to the Kashirin (1986) model.The photo-ionization and photo-absorption crosssections are obtained from Torr et al. (1979), Richards and Torr (1988), Richards et al. (1994b); and from Ivanov-Kholodny and Nikoljsky (1969) for X-ray emission.The list of chemical processes used in the model is given in table 4.I.As long as we work above 150 km the contribution of neutral NO can be neglected and corresponding reactions are not in the list.Atomic nitrogen being used in some reactions (table 4.I) is taken from the MSIS-86 model.The most important for the F2-region chemistry is the O + +N2 reaction.Recent flowing afterglow laboratory measurements by Hierl et al. (1997) of this reaction rate are used in the model.These measurements were made at Tn=Ti = Tv (where Tn -neutral, Ti -ion temperatures and Tv -vibrational temperature of the excited N2) in a wide temperature range and take into account the effects of vibrationally excited N2 * .This may be important for summer high solar activity conditions (e.g., Pavlov and Buonsanto, 1997;Pavlov et al., 1999 and references therein).A comparison of different O + +N2 reaction rate constants using EISCAT observations has shown that the Hierl et al. (1997) rate coefficient for this reaction may be recommended for aeronomic calculations (Mikhailov and Schlegel, 2000).The experimental dependence from Hierl et al. (1997) was approximated by a polynomial to be used in the model.In the auroral zone strong horizontal E × B drifts of the ions increase the effective temperature for this reaction.According to Schunk et al. (1975) we accept Teff = (miTn + mnTi)/(mn + mi) + 0.329 E 2 where E is in mV/m.Concentration of O + ( 4 S) ions is calculated from the continuity equation where Vi is measured total vertical velocity of plasma (Vi = Vl sinI at EISCAT) and W is vertical plasma drift due to thermospheric winds and electric fields.The first right hand term of eq.(4.2) is the diffusion velocity for O + ions (≠ Ne) from Banks and Kockarts (1973).Observed and smoothed (see earlier) Te and T i values are used in calculations.Expression (4.2) specifies W, which in principle can be con-verted to the meridional thermospheric wind.This is a standard approach usually applied to find meridional thermospheric winds from ISR observations (e.g., Buonsanto and Wittasse, 1999 and references therein).Neglecting vertical neutral winds (possible under small Joule heating conditions), the meridional neutral wind Vnx can be found from .sin cos cos sin cos cos = -At EISCAT due to small magnetic declination (D=1.2°), the contribution of zonal neutral wind Vny may be neglected.Diffusion collision frequencies νij for O + are related to momentum transfer collision frequencies ν* by the expression (see eq. 19.13 in Banks and Kockarts, 1973) , where i applies to O + ions and j applies to other neutral or ionized gas species.Collisions of O + ions with neutral O, O2, N2 and only with NO + ,O2 + ions are taken into account.Unlike the previous version of the method where all ions were included to νij, here we leave only major ones thereby increasing the numerical stability of calculations.All O + ion collision frequencies are taken from Banks and Kockarts (1973).A correction factor 1.2 for ν(O+-O) was applied in accordance with the results of analyses by Pesnell et al. (1993), Oliver and Glotfelty (1996), Buonsanto et al. (1997), Litvin et al. (2000).The ν ij , q i , β values in eqs.(4.1) and (4.2) depend on the aeronomic parameters searched for.
Normally the results are rather insensitive to the choice of the upper boundary height within 400-600 km range and any calculation can be used as a solution (Mikhailov and Kofman, 2001).Therefore, calculations are made with the upper boundary at 550 km where [O + ] = Ne obs − ∑ ni is specified.Ion concentrations ni are known at each iteration of fitting calculated N e (h) to the experimental one.
Usually the low boundary for [O + ] is specified at 160 km during daytime, but it may be set at 230-260 km during nighttime at EISCAT location to avoid particle precipitation effects at lower heights.Corpuscular ionization is not included in the model and periods with intensive precipitation cannot be developed with the method (see later).

Fitting procedure
An important point is the list of aeronomic parameters to be searched for.These should be the main critical parameters responsible for the F-region formation.An analytical solution of the continuity equation for electron concentration in the F2-region (e.g., Ivanov-Kholodny and Mikhailov, 1986) can help select the parameters.In the mid-latitude daytime F2-region they are: neutral temperature Tn, atomic oxygen [O], linear loss coefficient β = γ1[N2] + γ2[O2], efficiency of solar ionization, and vertical plasma drift W due to thermospheric winds and electric fields.
All aeronomic parameters may be divided into two groups: rate constants and cross-sections of processes, which are supposed to be known from laboratory experiments and those which vary with geophysical conditions.The latter present our main concern, although the former can be also analyzed using ionospheric observations (e.g., Ivanov-Kholodny et al., 1976;Oliver and Glotfelty, 1996;Buonsanto et al., 1997;Mikhailov and Schlegel, 2000).
Atomic oxygen is the main neutral constituent at the F-region heights and it specifies the total production rate as well as the diffusion rate of O + ions.The recombination rate of O + ions is specified by linear loss coefficient β, so molecular nitrogen and oxygen should be in the list, moreover they strongly contribute to the total production rate.There is still a pretty large uncertainty with the solar EUV fluxes in particular with their day-to-day variations due to varying solar activity.Total EUV flux rather than spectral distribution of the energy is important for F2-layer ionization (Ivanov- Kholodny and Nikoljski, 1969), therefore this parameter should also be in the list.Neutral temperature height profile is a crucial structural parameter for the whole thermosphere.Within the MSIS model formalism (we are working with), the Tn(h) is specified by three parameters: Tex, T120, and the shape parameter S.Although we are not dealing with heights below 160 km it was found that the method works better if the MSIS T 120 value is freed.So, T120 was formally added to the list of searched parameters.But it should be stressed that this is just a technical step and the extension of Tn, [O], [O2], [N2] down to 120 km height is just an extrapolation as we do not fit any Ne(h) profile below 160 km height.The last very important parameter is vertical plasma drift velocity W due to thermospheric winds and electric fields.Unfortunately the success of N e(h) fitting strongly depends not only on the absolute value of W at the F2-layer heights, but also on the height profile W(h).Therefore, any attempts to use height-independent W normalized by empirical model values (e.g., Hedin et al., 1991) should give unsatisfactory results.The only acceptable approach according to our experience is to use simultaneous ISR observations of total vertical plasma velocity subtracting diffusion velocity (see eq. (4.2)).This confines the application of the method to high-quality ISR data like the EISCAT facility provides.Due to technical problems with Vi measurements some of Millstone Hill observations cannot be developed with the method.
So we have seven parameters to be found: three parameters (Tex, S, T120) specifying Tn(h) profile, factors for the MSIS model [O], [O2], [N2] concentrations as well as for the total solar EUV flux with λ <1050 Å. Vertical plasma drift W(h) is calculated at each step of fitting using eq.(4.2).The whole ion composition is also available at each step of calculations.Using standard multi-regressional methods we fit the calculated Ne(h) profile to the observed one and find a self-consistent set of aeronomic parameters listed.It is important to stress that all parameters should be found simultaneously.All the attempts to divide the fitting process by steps gave unsatisfactory results.
During quiet and slightly disturbed conditions the initial values of searched parameters are set equal to the model ones with a ± 20% corridor for their variations.Usually this is enough to obtain an acceptable solution.But in disturbed conditions the thermospheric parameters are beyond the ± 20% corridor and some steps are required to localize the area of possible solution, that is to find approximate values of the searched parameters.After this the standard method is applied to find the final values of the parameters.Two criteria are used to select the solution: i) a fit of N e(h) should be good in the whole height range, that is the calculated profile should lie within the experimental ± SD corridor, the goodness of fit being estimated with ii) The extracted parameters should be within the ± 20% corridor with respect to the average values.There are stationary and non-stationary versions of the method.In case of the stationary solution the left hand term of eq.(4.1) is set to zero, while [O + ] (≠ Ne) height profile from the previous time step should be available to obtain a non-stationary solution.The stationary form of eq.(4.1) is used during sunlit conditions.A comparison of stationary and non-stationary solutions shows close results for solar zenith angles χ ≤ 95°.During twilight and in the night (χ > 95°) a non-stationary approach gives better results.Such behavior is explainable.Median profiles of observed characteristics used in calculations (see earlier) are obtained over 1.0-1.5 h time step in case of EISCAT observations, but this time interval may be larger for Millstone Hill observations.Therefore, the minimal time step may be 1 h for solving eq.(4.1).This time step ∆t is close to the e-fold time (≈ 1.5 h) of the daytime F2-region.In this case it is better to use a stationary solution (Waldman, 1973).During nighttime the e-fold time is much larger than 1.5 h and a non-stationary solution is preferable.Unfortunately, a non-stationary solution needs much handy work to prepare median smoothed Ne, Te, Ti, Vi profiles for all time steps, and can hardly be applied for routine calculations.An example for one quiet day is given later.
The quality of Ne(h) fitting is illustrated in fig.4.4 where EISCAT daytime observations for quiet (top) and disturbed (bottom) are compared with calculations.Different levels of solar activity (minimum, June 26, 1997;medium, August 02, 1992;and maximum, July 02, 1990) and different levels of geomagnetic activity are considered (Ap indices for this and previous days are shown next to curves in the bottom panel).Simple quiet time Ne(h) profiles as well as much complicated disturbed profiles with well-pronounced F1-layer (August 05, 1992) or profiles without F2-layer maximum at all (April 10, 1990) can be fitted in the whole height range with a good accuracy.

COMPARISON WITH MSIS-86 MODEL
The method was applied to EISCAT (CP-1 and CP-2) daytime observations (28 magnetically disturbed and 30 quiet days) and the results were compared to the MSIS-86 thermospheric model predictions.The results for exospheric temperature Tex, shape parameter S for Tn(h), and total neutral density ρ are given in fig. 4 Total neutral gas density -The MSIS-86 model overestimates ρ at 300 km height for disturbed days with an average ρcal/ρMSIS = 0.52 (the difference is significant at the 99% confidence level).The deviation from MSIS-86 is insignificant for quiet time conditions.The result may be important for satellite orbital characteristics.
Exospheric temperature -The MSIS-86 model underestimates Tex by 5-7% on average both for quiet and disturbed conditions.The differences are significant at the 99% level.

Shape parameter S for Tn(h) profile -
The MSIS-86 model strongly overestimates S both for disturbed (average Scal/SMSIS = 0.48) and quiet (Scal/SMSIS = 0.89) days.Both deviations are significant at the 99% level.
Therefore, there is a problem with MSIS-86 predictions for disturbed conditions and this has been stressed repeatedly (see earlier), while in quiet conditions it gives the results which in general are close to our estimates.This is not surprising as any empirical model MSIS-86 presents fairly well only the conditions corresponding to the majority of observations, that is to quiet and slightly disturbed condi-   tions which are the most probable, while strong disturbances are relatively rare.Although the dependence on 3-hour Ap (or daily Ap) global indices is formally included in the model, it gives a smoothed averaged pattern of the disturbed thermosphere and cannot describe the peculiarity of an individual disturbance.On the other hand, the closeness of our results to the quiet time MSIS predictions may tell us about the efficiency of the method developed.But for a convincing conclusion simultaneous ISR and satellite probe in-situ observations are required, which are not available at present.

DIURNAL VARIATIONS
Calculation of diurnal variation of thermospheric parameters is another testing of the method to check whether the found parameters demonstrate consistent variations during the day.All fits are independent in case of the stationary approach normally used during sunlit conditions, while they are slightly related when a non-stationary scheme is applied.EISCAT observations for quiet days of April 1-2, 1992 were used for this analysis.The calculations were made for 24 UT moments with median N e, Te, Ti, Vi profiles prepared over 1 h time intervals around the UT moments in question.Figure 4.7 gives the observed NmF2 and hmF2 variations as well as ratios of calculated parameters to the MSIS-86 values.Solar zenith angle variation for the period in question is shown in the bottom panel.Stationary solutions (as mentioned earlier) were used for sunlit ( χ ≤ 95°) period and a non-stationary scheme -for nighttime hours.The results are seen to be consistent for the sunlit hours.Average ratios along with SD deviations for the parameters considered are given in table 4.II for the sunlit hours.The difference from MSIS for average values is also estimated using the Student criterion.
The standard deviations with respect the average values are seen not to be large -about 3.5% for Tex and S and 15-21% for neutral species.The difference of thermospheric parameters from MSIS values for this particular quiet day in general agrees with the statistical estimates given earlier.
The results for nighttime hours are less satisfactory (fig.4.7).Large ratios for neutral species perhaps reflect large Tex and S values, but they look as unreal although the quality of Ne(h) fitting is excellent.The calculated Tex turn out to be larger than observed Ti values and this is impossible during nighttime hours when both temperatures should be close.Reasons for such behavior are not clear right now and additional analyses are required.One of the possible factors is the electron precipitation taking place in the auroral nighttime F2-region, but this effect is not included in the model.Changes in Ne(h) profile resulted from additional corpuscular ionization, the method ascribes to thermospheric parameters which accept unreal values.

A SEVERE STORM (CASE STUDY)
The efficiency of the method can also be illustrated analyzing the periods of severe geomagnetic disturbances when large perturbations of thermospheric parameters are expected.The April 10, 1990 severe storm with Ap = 124 presents a good example.Both Millstone Hill ISR (American sector, middle latitudes) and  EISCAT (European sector, auroral zone) provided observations for that period, which was considered by Mikhailov and Foster (1997), Mikhailov and Schlegel (1998).Observed and calculated Ne(h) profiles as well as ion composition (O + /Ne ratio) are shown in fig.4.8 for EISCAT (top panels) and Millstone Hill (bottom panels).The severe geomagnetic disturbance resulted in a complete disappearance of F2-layer at usual heights and the ionospheric maximum was presented by F1-layer peak -so called G-conditions.The disturbance picture is seen to be similar at two locations, but the storm effect was larger at Millstone Hill where electron concentration decreased by a factor of ten at the F2-layer heights.Large changes in ion composition (O + /Ne ratio) took place during that event.The ionosphere was strongly enriched with molecular ions even at the F2-layer heights, while this ratio was close to the standard model in quiet conditions.Large deviations of O + /Ne ratio from the standard model being used during the development of ISR observations required a correction of the experimental Te(h), Ti(h) and Ne(h) profiles for the disturbed day using the iterative method by Mikhailov and Schlegel (1997).Strong perturbations of the thermospheric parameters were revealed for April 10 relative to the previous quiet days.At EISCAT the calculated Tex increased by 550 K The calculated temperature profile Tn(h) also strongly differs from the MSIS one for that disturbed day. Figure 4.9 gives three daytime temperature profiles at the EISCAT location on April 10, 1990.A very low value of the shape parameter, S=0.0055 km − −1 instead of 0.016 km − −1 predicted by MSIS resulted from Ne(h) fitting.MSIS Tn(h) values are seen to be larger than observed ion temperatures below 300 km and this is impossible in reality.Therefore, the method seems to give reasonable results in strongly disturbed conditions as well.A small value of S (0.011 km − −1 with respect to the MSIS predicted S = 0.017 km -1 ) was obtained for the same day, April 10, 1990 from the Millstone Hill data analysis.Therefore, the results of calculations show a pronounced decrease in the neutral temperature in the low thermosphere for this severe geomagnetic storm.It was proposed (Mikhailov and Schlegel, 1998) that this Tn decrease resulted from the enhanced atmospheric cooling through nitric oxide.The increase in NO concentration during disturbed periods is well documented (e.g., Barth, 1989Barth, , 1995;;Ridley et al., 1999;Solomon et al., 1999).

DISCUSSION
The proposed method demonstrates high efficiency allowing us to infer basic thermospheric parameters as well as ion composition from routine ISR observations both in quiet and disturbed conditions.Along with observed plasma parameters this gives a complete picture of the ionosphere and thermosphere state in the vicinity of the ISR facility.ISR have been operating for many years and this opens an opportunity to analyze thermosphere-ionosphere interaction in various geophysical conditions.Among the problems which are the most interesting and challenging both from scientific and practical points of view is the F-region storm effects, and the developed method was widely used for such analyses (see the reference list).
Our analysis of daytime EISCAT observations has shown that MSIS-86 strongly and systematically overestimates [O] at high latitude both in quiet and especially during disturbed periods (fig. 4.6).Similar results were obtained by Litvin et al. (2000) analyzing with the energy equation Millstone Hill observations for a disturbed period of June 5-11, 1991.During the most active phase of the disturbance they found [O] to be lower by a factor of 2 than MSIS-86 predictions.On the other hand, their Tex are lower than MSIS values by 200 K during daytime.Although the authors stress that the shape parameter S strongly affects the results, they used the MSIS-86 values in their calculations, which are strongly overestimated for disturbed conditions (fig.4.5).The Tex * S product demonstrates a sort of invariant in the inverse problem considered, therefore large MSIS-86 values of S should result in lower Tex values in accordance with their conclusion.Optical observations by Grossmann et al. (2000) revealed a systematic (by 40% on average below MSIS-86) difference in atomic oxygen densities in the 130-175 km height range.
Therefore, independent results on [O] variations obtained with other methods seem to agree with our conclusions.But these are partial comparisons on one-two parameters while the proposed method provides information simultaneously on many thermospheric parameters.Only simultaneous ISR and satellite (similar to AE-C mission) observations of the whole set of thermospheric characteristics could be used for the overall testing of the method.Unfortunately, there is no such opportunity at present.
Let us consider problems and limitations related to the usage of the method.Firstly this concerns the ISR observations.The quality of the experimental material is different at different facilities.For instance, due to rare (usually 3 per hour) observations at Millstone Hill a 3-4 h time interval is needed to calculate median profiles.It is not always possible (especially for disturbed days) to find a 3-4 h period of relative stability in NmF2 and hmF2 variations.This criterion is applied to decrease the scatter in the observations and increase the reliability of median profiles.Unlike EISCAT observations there is a problem with using routine V i(h) data at Millstone Hill as they may need an additional correction due to technical reasons.Therefore, not all available routine Millstone Hill observations can be used for our analysis.The same difficulty takes place with Kharkov ISR facility which does not provide Vi(h) observations at all.
On the other hand, EISCAT observations (each particular experiment) are not normalized by foF2 values as is done at Millstone Hill.Although the declared uncertainty in the measured electron concentration at EISCAT is not large −10 ÷12% (Farmer et al., 1984;Kirkwood et al., 1986) in fact in some cases the difference in NmF2 (compared to nearby ionosonde observations) is much larger as the analysis by Mikhailov and Schlegel (2003) has revealed.The absence of such normalization also results in a shift between long pulse and multi pulse N e(h) profiles and this shift should be taken into account before the profiles are used for analysis.
The problem with dependence of experimental Te(h), Ti(h) and Ne(h) on ion composition was mentioned earlier.But despite the problems encountered the ISR data are the most informative and consistent being the only observations which can be used with such a method to retrieve the aeronomic parameters.
According to numerous testings of the method acceptable results can be obtained for sunlit (χ ≤ 95 o ) conditions only when the F2-region is mainly controlled by photo-chemical processes and the contribution of each parameter is distinct, therefore it can be retrieved.During nighttime hours dynamical processes dominate and the contribution of individual parameters to the N e(h) distribution is not that distinct.This is seen in fig.4.7 when after 22 UT all parameters are highly correlated and the reliability of each extracted parameter is not very high in this case.Additional problems with EIS-CAT nighttime observations are related to horizontal E × B drifts moving plasma via the location as well as with particle precipitation producing fresh ionization.Neither effect is included in the model and cannot be taken into account properly.These additional factors are absent at middle latitudes and further analyses of Millstone Hill nighttime observations are required to clear up the possibility of the method in nighttime conditions.
The accuracy of calculated aeronomic parameters depends on many input quantities and can hardly be estimated strictly.Only direct independent and simultaneous satellite observations of neutral composition and temperature could be used for such a strict comparison.Such observations are absent at the moment.But the lowest estimate of expected errors resulting from the uncertainty in measured ionospheric parameters can be made (Mikhailov and Schlegel, 1997) The uncertainty in measured electron concentration according to Farmer et al. (1984) and Kirkwood et al. (1986) gives ∆logNmF2 about ± 0.05.The hmF2 scatter around its median value usually is about ±10 km over the chosen one hour period of observations.The uncertainty in measured field aligned plasma velocity is ±2 m/s at the height of 300 km (Jones et al., 1986).This gives the uncertainty of ≈16% in [O] and ≈18% in the linear loss coefficient β.As β in the F2-region is mainly determined by molecular nitrogen, this estimate may be applied to [N 2].But it should be kept in mind that the above estimates are the least expected errors for the calculated [O] and [N2].A comparison with an empirical (i.e.averaged) thermospheric model like MSIS-86 in quiet geomagnetic conditions can also give an estimate of the accuracy for the retrieved parameters.According to results given in figs.4.5 and 4.6 the overall uncertainty in modelled Tex is less than ±10% and in neutral species is about ± (20-30)%.This is an acceptable result keeping in mind that the accuracy of many input parameters is of the same order.But the possibility of the method to extract thermospheric parameters for disturbed conditions is the most valuable result.The empirical models like MSIS are not efficient during storm periods and the proposed approach can provide useful information at least in the vicinity of the ISR facility.

CONCLUSIONS
The results of our analysis may be summarized as follows.1) Further development of the self-consistent approach to the F2-region modeling has been made in three directions: 1.1)An increase in numerical stability has been achieved by taking into account only major ions O2 + and NO + in the expression (4.2) for the O + diffusion velocity, while a complete list of pertinent chemical reactions is used to calculate ion concentrations.
1.2) An inclusion of nighttime ionization by scattered radiation in 1026, 584, 304, and 910 Å lines allowed us to extend the area of possible solutions up to χ ≈ 95-100°at least from the evening side.
1.3)An extension of the method for non-stationary conditions allowed us to consider twilight and nighttime hours when the F2-region is not stationary.
2) Application of the method to analysis of EISCAT daytime observations (28 magnetically disturbed and 30 quiet days) has shown that MSIS-86 strongly overestimates atomic oxygen and the shape parameter S in particular for disturbed conditions (average [O]cal/[O]MSIS = 0.43 and Scal/SMSIS = 0.48), while Tex is close (the difference is 5-7% on average) to MSIS-86 predictions both in quiet and disturbed conditions.With regard to molecular species MSIS-86 strongly overestimates O2 and N 2 concentrations for disturbed periods, while these concentrations are close to MSIS-86 values for quiet days.The results are statistically significant.A comparison on the total neutral gas density, which is important for practical applications, has shown that MSIS-86 overestimates ρ at 300 km for disturbed days with an average ρ cal /ρ MSIS = 0.52 (significant at the 99% confidence level), while the deviation from MSIS-86 is in-significant for quiet time conditions.But for a convincing testing of the method simultaneous ISR and satellite (similar to AE-C mission) observations of the whole set of thermospheric parameters is required.Unfortunately, there is no such opportunity at present.
3) The method can be efficiently used in sunlit ( χ ≤ 95-100°) conditions when the F2-region is under essential control of photo-chemical processes.The contribution of various thermospheric parameters to the Ne(h) distribution is distinct in this case and they can be retrieved with an acceptable accuracy.During nighttime conditions when dynamical processes dominate the role of individual parameters is not that distinct and found solutions may be unsatisfactory.Further analyses in this direction are required.
4) F2-region storm effects which are of great importance both for ionosphere physics and applications can be successfully analyzed with the method.Retrieved neutral and ion composition, temperature and meridional wind give a complete picture of the thermosphere state in the vicinity of the ISR facility for the conditions in question.Keeping in mind that most satellite mass-spectrometric observations were made in the past and are not conducted currently, the method may be considered a real tool for investigation of the thermosphere-ionosphere interaction under various geophysical conditions.
, in particular for this case when only the ionospheric part of TEC below 425 km was considered.With the same correction of [O] it was possible to obtain N e, Te, Ti, and Vi height profiles close to EISCAT observations (fig.4.2), [O2] and [N2] were fixedto default MSIS values.It was also shown that acceptable

Fig. 4 . 1 .
Fig. 4.1.Total electric field measured by EISCAT (top panel).The NmF2 and ITEC variations are shown in the third and fourth panels.Bold lines are the measurements and the thin lines are the fits with a correction factor for the atomic oxygen density.Second panel gives the relative errors for NmF2 (upper line) and for ITEC (lower line).The bottom panel shows the values of the time varying correction factor for atomic oxygen.

Fig. 4 . 2 .
Fig. 4.2.The four upper panels are the ionospheric profiles at 6 LT, the middle panel at 12 LT and the bottom panels at 16 LT.From left to right, one gets the electron density, electron temperature, ion temperature, ion velocity parallel to B. The broken lines are the measurements with instantaneous values in dashed lines and the surrounding extreme values in full lines.The modeling outputs are in smoothed lines.The dashed lines stand for no correction factor in the neutral atmosphere.The dotted lines are with a correction factor of 0.7.The full line (in the electron density only) for a time varying correction factor.This factor's values are 0.915 at 6 LT, 0.99 at 12 LT and 0.65 at 16 LT.
Fig. 4.3.An example of EISCAT observations for a disturbed day of August 05, 1992.Initial median (before polynomial smoothing) profiles calculated over 11:00-12:00 UT time interval are shown.Standard deviations are not experimental, but show the scatter of profiles within the selected 1 h time interval.

Fig. 4 . 4 .
Fig. 4.4.Observed at EISCAT and calculated daytime Ne(h) profiles for quiet (top) and disturbed (bottom) conditions.Profiles for different levels of solar activity (minimum, June 26, 1997; medium, August 02, 1992; and maximum, July 02, 1990, top panel) and different levels of geomagnetic activity are considered (Ap indices for this and previous days are shown next to curves in the bottom panel).Note the variability of Ne(h) profiles under disturbed conditions.

Fig. 4
Fig. 4.5.A comparison of calculated versus MSIS-86 exospheric temperature Tex, shape parameter S, and total neutral density ρ at 300 km for quiet (left-hand panels) and disturbed (right-hand panels) days.Tex is in K, S is in 10 -2 km -1 units, and log ρ is in gcm -3 .Dashed lines indicate a ± 25% band around the model.

Fig. 4 . 7 .
Fig. 4.7.Diurnal variation of calculated Tex, S, O, O2, and N2 to MSIS-86 value ratios.Neutral concentrations are taken at 300 km height.Two upper panels give EISCAT observations of NmF2 and hmF2 for a quiet period of April 1-2, 1992.Solar zenith angle variation is given in the bottom box.

Fig. 4 . 8 .
Fig. 4.8.Observed and calculated at EISCAT (top) and Millstone Hill (bottom) Ne(h) profiles for disturbed April 10, 1990 and quiet reference days.Right-hand panels -calculated and standard model O + /Ne profiles at the two locations.F2-region is seen to be strongly enriched with molecular ions for the disturbed day, while the calculations are close to the models for quiet days.

Fig. 4 . 9 .
Fig.4.9.Calculated and MSIS Tn(h) profiles together with the observed at EISCAT Ti(h) profile for the disturbed day April 10, 1990.Note that MSIS temperatures below 300 km are larger than Ti -the situation is unreal for the conditions in question.

Table 4 .
I. Chemical reactions used in the model.

Table 4 .
II.Average ratios of calculated thermospheric parameters to MSIS-86 values over the sunlit hours.The difference from MSIS values is estimated with the Student criterion.