Modeling the plasmasphere to topside ionosphere scale height ratio

A new model of plasmasphere to topside ionosphere scale heights ratio is developed, based on topside electron density (Ne) profiles deduced from the International Satellites for Ionospheric Studies (ISIS)-1 satellite measurements. The model is able to improve operational algorithms for space weather predictions. The topside ionospheric and plasmaspheric scale heights are determined by the lowest and largest gradients of measured profiles, respectively, converted in dh/dlnNe units. The new model depends on four parameters: the month of the year (M), the local time (LT), the geomagnetic latitude (glat), and the ln(O) density (zO) at the O-H ion transition height. It is designed to replace the old one-dimensional model of the ratio in the TaD (TSMassisted Digisonde) profiler. The parameters M, LT, and glat are approximated by trigonometric basis functions, while zO is described by a polynomial. A series of models were produced with different number of coefficients (number of terms) of the basis functions. Comparison between models revealed that those with larger number of coefficients can produce unrealistic extremes of the model curves due to the non-uniform sampling of data along the axes. Further considered is the simplest model approximating M, LT, and glat by simple 24 sinusoidal functions and linearly depending on zO. The model description and its 54 coefficients are given in Appendix 1 and can be used by other users for reconstruction of plasmasphere density profiles. The main variation of the ratio along geomagnetic latitude at fixed values of the other model parameters is illustrated in a series of plots.


Introduction
Plasmasphere thermal plasma density distribution is an important factor in transionospheric radio propagation, satellite communications, Global Navigation Satellite System (GNSS) applications, and many other human activities in space.It has been a subject of continuous study and modeling.In addition to earlier observations of the topside ionosphere and lower plasmasphere by means of rockets (Gringauz et al. 1960) and whistlers (Carpenter & Smith 1964), many direct measurements and models have contributed to the present knowledge of the dynamic nature of the plasmasphere.A first representative description of plasmaspheric density was given by Bent et al. (1972) based on whistler measurement.They provided graphs of exponential height distribution of density with scale height depending on the critical frequency foF2, the geomagnetic latitude, and the solar flux index (F10.7).Rawer et al. (1978Rawer et al. ( , 1981) incorporated the Bent model (Bent et al. 1972) into the International Reference Ionosphere (IRI) model, transforming the graphic presentation of density distribution into analytical expressions.Numerous verifications of the IRI model have revealed that its topside and plasmasphere model provided unrealistically high values of plasma density, especially at high solar activity, although substantial improvements have been made during the years.Bilitza (2004) made an important improvement in the IRI topside and plasmasphere model by introducing a correction factor based on ISIS-1 and 2 and Alouette-1 and 2 topside sounder data.This correction factor varies with height and modified dip latitude and brings the IRI topside and plasmasphere model in good agreement with the data.The new corrected model is included in the IRI2007 version (Bilitza & Reinisch 2008).Gallagher et al. (1988) developed an empirical model of plasmasphere low energy plasma based on measurements from the Retarding Ion Mass Spectrometer (RIMS) experiment onboard the Dynamic Explorer-1 satellite.This model provides the electron density as a function of day of the year, universal time (UT), height, and geographic coordinates.The Gallagher et al. (1988) plasmasphere model was used in the Parameterized Ionospheric Model (PIM) as a plasmasphere extension of the global ionospheric model (Daniell et al. 1995).For user-specified geophysical conditions and spatial coordinates, PIM produces electron density profiles (EDPs) between 90 and 25,000 km altitude.Gallagher et al. (2000) developed the Global Core Plasma Model (GCPM) which provides empirically derived core plasma density and ion composition (H + , He + , and O + ) as a function of geomagnetic and solar conditions throughout the inner magnetosphere.The GCPM model uses the Gallagher et al. (1988) plasmasphere model for high altitudes and merges with the IRI model at 500-600 km altitude.
A new empirical model of plasmasphere density distribution has been developed recently, based on the Radio Plasma Instrument (RPI) onboard the Imager for Magnetopauseto-Aurora Global Exploration (IMAGE) satellite (Reinisch et al. 2001;Huang et al. 2004) which provides N e profiles along magnetic field lines (Ozhogin et al. 2014).To connect the plasmasphere model to the topside ionosphere, Reinisch et al. (2007) introduced a modified Chapman function, named Vary-Chap, with a continuously varying scale height H(h).These authors modeled the ratio of H(h) to the Digisondederived scale height at the F2 peak (H m ) as a function of season, latitude, and local time.
In this paper we present a new model of plasmasphere to topside ionosphere scale heights ratio, based on topside electron density (N e ) profiles deduced from the International Satellites for Ionospheric Studies (ISIS)-1 satellite measurements.This model is a part of the Topside Sounder Modelassisted Digisonde (TaD) Profiler software providing the electron density profile from the F layer peak to GNSS orbit heights.The TaD software is running operationally in the European Digital Upper Atmosphere Server, DIAS (Belehaki et al. 2006(Belehaki et al. , 2007)), to provide the maps of the electron density over preselected heights in the topside ionosphere and plasmasphere, and maps of the Total Electron Content (TEC) over Europe (http://dias.space.noa.gr).The new model proposed here describes the scale height ratio as a function of four parameters and aims to replace the old one-dimensional model ratio currently used in the TaD Profiler.The chronology of developing the family of models to which the TaD profiler belongs is described below in more detail.

TaD background
A series of models have been developed by utilizing data from topside sounders aboard Alouette and ISIS satellites.The TaD model belongs to this group of models.Its core component, the Topside Sounder Model (TSM) originally proposed by Marinov et al. (2004) and Kutiev et al. (2006), reproduced the topside electron density scale height (H T ) and the O + -H + (upper) transition height (h T ), based on 172,622 measured N e profiles by topside sounders onboard Alouette-1a, -1b, -1c, and -2 and ISIS-1 and -2 satellites (Bilitza et al. 2003).The H T and h T models were represented by five-dimensional polynomials expressing both quantities as functions of month of the year, geomagnetic latitude, local time, solar flux (F10.7), and geomagnetic index (K p ). Later, Kutiev & Marinov (2007) developed a model of the ratio R T = H T /h T , based on ratios of individual profiles and combined the three models into a single model named Topside Sounder Model (TSM).The key element in the modeling approach is the assumption that the lowest gradient on the topside N e profile represents the altitude gradient of the O + density profile.Assuming an exponential distribution of the O + density, the transition height h T is found at the height where the upward extrapolated O + density becomes half of the measured N e density.Kutiev et al. (2006) offered a method for obtaining the shape of the vertical plasma distribution in the topside ionosphere and plasmasphere by introducing the TSM parameters H T and h T into three well-known formulas describing the vertical plasma distribution: a-Chapman, sech-squared, and exponential (see Stankov et al. 2003).The analytical formulas containing H T and h T , named the Topside Sounder Model Profiler (TSMP), provide the shape of the vertical plasma distribution.As was recently shown by Verhulst & Stankov (2014, 2015), during non-disturbed conditions the exponential and a-Chapman formulas offer better approximations of the topside F region electron density among other analytical profilers.The a-Chapman formula was used to describe the O + distribution, while the H + distribution was presented by an exponential formula with scale height H H+ = 16H T (Kutiev et al. 2009a).To obtain the density distribution, the F2 layer peak density (NmF2) and height hmF2 are required to be specified at its lower boundary.Kutiev et al. (2009a) used the TSM ratio R T = H T /h T to link the TSMP to the Digisonde topside profiling technique.The combined profiler, named TaD (TSM-assisted Digisonde) profiler, used the F2 layer peak density NmF2, peak height hmF2, and the scale height H m , all provided by the Digisonde software, to reconstruct the electron density profile (EDP) to 20,000 km height.
The Digisonde software (Reinisch & Huang 2001) extends the measured bottomside profile in the topside by assuming an a-Chapman shape with a scale height H m obtained from the bottomside profile at hmF2.The scale height H m is a neutral scale height, so well above hmF2 the Digisonde a-Chapman EDP decreases with a gradient corresponding to a doubled H m .Kutiev et al. (2009a) found that the doubled H m is systematically lower than the TSM scale height H T by a factor of ~1.25.In order to adjust TSMP to Digisonde measurements, Kutiev et al. (2009b) multiplied H m by a correction factor k = 2.5 and obtained the corresponding transition height h T through the ratio h T = kH m /R T , with the ratio R T taken from TSM.Having obtained the scale height H T = kH m and transition height h T , the TaD profiler extends the measured bottomside EDP to the topside F region and plasmasphere.It should be noted that while TSMP describes the shape of the topside EDP only, TaD provides EDP from hmF2 up to GPS orbit heights over Digisonde locations.
The correction factor k, which converts H m into H T , is statistically obtained by comparing Digisonde data with the corresponding TSM product.Although the correction factor improves the statistical reliability of TaD profiles, it cannot guarantee reliable reconstruction of EDP in all individual cases.To further improve the TaD accuracy, Kutiev et al. (2012) developed a method which compared the EDP integral (bottomside + topside profiles) with GNSS-derived total electron content (TEC GNSS) obtained at the location of the ionosonde.In this method the profile parameters k, H m , h T , and H p are iteratively varied until the difference between the integral and the measured TEC becomes negligible.It was found that H T changes contribute more than 90% to the total change of the integral and therefore in the final method only H T was used in the adjustment.

Database
The ISIS-1 satellite had a perigee of 3600 km and the topside sounder onboard the satellite was the only sounder systematically providing electron density profiles (EDP) in the plasmasphere.In order to derive a reliable scale height value from the H + part of EDP, we selected those profiles which are well extended into the plasmasphere (their uppermost values had altitudes above 2500 km).The selected database contains 14,641 measured EDPs. Figure 1 shows the distribution of measured EDPs in geographic coordinates (left) and local time (on the right).Clearly, the global distribution is not uniform, instead it clusters at Asian equatorial latitudes and in the American sector.The local time coverage is acceptable in equatorial and polar latitudes but poor at mid-latitudes.
The database established for this development is formed by the topside scale height H T , the transition height h T , the ratio R T = h T /H T , the plasmasphere scale height H p , and the ratio R H = H p /H T , extracted from each individual measured profile.
The method of extracting the topside parameters H T , h T , and R T was described by Marinov et al. (2004), Kutiev et al. (2006), and Kutiev & Marinov (2007).The topside (O + density) scale height H T was defined at the minimum gradient dh/d(ln O + ) of the topside part of the profile plotted in height/lnNe scale.Gradients were calculated at the height of each data point of the measured profile and those of them with values within 30% from the minimum gradient were averaged and then converted to the topside scale height H T .The 30% allowance is set to account for the increase of the scale height due to plasma temperature (T p ) increase in the range.It should be noted that the scale height defined here is a measure of the decay of the vertical ion density distribution with altitude and it differs from the theoretical scale height (kT p /mg) of plasma distribution along magnetic field lines at diffusive equilibrium.The data points, whose gradients were used to determine H T , were then approximated by a regression line to represent the O + distribution above hmF2.The transition height h T was found by extending the regression line upwards to the height where it reaches a density value of one-half from that of the measured profile.The scale height H T , transition height h T , and their ratio R T were separately modeled by five-dimensional polynomials (Kutiev & Marinov 2007), unified in a single Topside Sounder Model (TSM).Kutiev et al. (2009b) used a similar approach to define the plasmasphere scale height H p .They identified H p by the maximum gradient of the measured profile lying above h T and also averaged the gradients at several data points around the height of maximum gradient to obtain the H p value.In the present paper we obtain H p in a similar way, using this time a dynamic filter to select the gradients around the maximum one.At the first step the filter gathers data points with gradients 1% below the maximum one.The procedure repeatedly increases the band of allowance stepwise by 1% until the number of data points inside the band reaches 6.This is usually achieved in 3-4 steps.Then H p is calculated by taking a linear regression over the gradients in the allowance band.
Figure 2 shows an ISIS-1 measured N e profile in altitude/ lnNe scales and the decomposed O + and H + distributions [in cm À3 ].The gray line represents the height h T at which the extrapolated upward O + density becomes one-half of N e .The H + distribution is independently drawn downward from the top of the N e profile and it does not necessarily meet the O + profile at the transition height h T .When this is the case, the H + density at h T is predominantly lower than that of O + .Belehaki et al. (2012) have studied these cases and suggested that this is due to the presence of helium ions He + .They modeled the density ratio g = n(H + )/n(O + ) at the transition height as a function of geomagnetic latitude, local time, and O + density.Then the He + density distribution was presented as an a-Chapman shape centered at h T with peak density n(He + ) = (1 À g)n(O + ) and scale height four times that of H T .In some applications of TaD profiler, as, for example, the adjustment of TaD integral to TEC GNSS measurements, the He + contribution is considered negligible and g is taken equal to 1.
The complete expression of the TaD profiler is given by:  where N O+ is the O + density, g is the ratio of H + to O + density, H p is the H + (or plasmasphere) scale height, and the scale height of He + is taken equal to 4H T in accordance with the He atomic mass.The first term on the right represents the O + density and the second and third terms represent H + and He + density distributions as functions of altitude.Therefore the electron density profile N e is a sum of O + , H + , and He + partial distributions.The increment of exponents contains the absolute value of the height difference |h À h T | in order to force the corresponding densities to decrease with decreasing height below h T , as it really happens.The model assumes constant scale heights with increasing altitudes, although in reality the scale heights increase with altitude, due to the increasing plasma temperature.This controversy is solved by calculating the scale heights as regression lines over selected height ranges of the measured profiles.The scale heights so obtained contain the changes due to temperature increase and can be regarded as representative of the whole respective ion distributions.
Figure 3 shows all 14,641 H T and H p values contained in the database plotted as a function of geomagnetic latitude.H T (green dots) increases toward the poles with an additional increase within ±10°around geomagnetic equator (in all figures in the paper H T is shown as HT due to graphic software limitations).H p (red dots, scaled on the right) exhibits a large increase in low and equatorial latitudes.In analogy with the ratio R T = H T /h T , here we model the ratio R p = H p /H T instead of H p itself in order to make the profile parameters dependent on H T .This is especially important when the profile parameters are varied to adjust its height integral with the independently measured TEC.

Modeling R p = H p /H T ratio
The old model (Kutiev et al. 2009b) is one-dimensional, depending solely on the geomagnetic latitude glat (H p / H T = 9cos 2 (glat) + 4).Here we develop a model of the ratio R p = H p /H T as a function of four parameters: month of year (M), local time (LT), geomagnetic latitude (glat), and ln(O+) density at transition height h T (denoted as zO).Along with the external parameters M, LT, and glat, on which R p strongly depends, the internal to the profile parameter zO is also involved to account for variations induced by ionospheric disturbances.The ratio R p calculated from all individual H T and H p pairs in the database is approximated by trigonometric and polynomial basis functions with varying number of terms.For variations along M, LT, and glat axes, the base functions are trigonometric (Fourier expansion terms), while along the zO axis the basis function is polynomial.Mathematical description of the model is given in Appendix 1.We varied the number of coefficients (or terms) of each parameter and obtained a set of models with their individual model errors.As an example, eight models with their number of coefficients and the model errors (absolute and relative standard deviations) are given in Table 1A of Appendix 1, along with the errors of the old model.Figure 4 compares the selected eight models along the geomagnetic latitude with the fixed values of zO = 10 at noon time of months M = 0.75 (January), 3.75 (April), 6.75 (July), and 9.75 (October).It is obvious that the basis functions with larger number of coefficients allow more extremes along the respective axes (the number of extremes is equal to the number of coefficients minus one), which usually ripples the model curves.The large maximum around 30°glat in January is a clear argument to question the use of higher order models.It is seen from Table 1A that the model error decreases with increasing number of coefficients and the 8th model performs the best.The wave-like structure of larger terms models, however, is probably an effect of the nonuniform distribution of the data along the model axes.The fitting procedure forces the models to get closer to the data clusters and can produce as a result unrealistic extremes.For further use, we favor the simplest model (3, 3, 3, 2) which has two extremes of each trigonometric function and linear dependence on zO.The model has an acceptable number of coefficients (54) which allow an easier reproduction from other users.Coefficients are given in Table 2A of Appendix 1.
Figure 5 is a summary plot of all R p values extracted from the database (blue dots) and of the model (3, 3, 3, 2) predictions (red dots), so that each data value corresponds a model value calculated for the same conditions.In the left panel, the model values (further on the numbers of selected model (3, 3, 3, 2) will be omitted) fit the data at all latitudes.The old model, which is shown in the right panel, provides a good approximation to the average data in the latitude bins (see Kutiev et al. 2009b), but significantly overestimates the ratio of scale heights at higher latitudes.
Figure 6 provides another view of the model performance.The data rows are sorted by the increasing values of the model ratios.In the left plot, data is sorted by the increasing values of the old model, while in the right plot the sorting is made by the increasing values of the new model.This is the reason why the model curves have a slope of monotonically increasing angle, although they are composed by a large number of dots.The X axes order data (blue dots) according to the model values from 1 to 14,641.As mentioned above, the old model overestimates the lower data values, but is close to the average data when the scatter increases.The total standard deviation of the old model is 0.55.The new model approximates data much    better (standard deviation is 0.40).The models with higher number of coefficients fit data with better accuracy as seen from Table 1A.

The R p behavior
The main variations of the R p ratio for several fixed values of the model parameters are presented in Figures 7-9.All curves are plotted versus geomagnetic latitude along which R p variations are strongest.Figure 7 contains four plots showing R p at four local times LT = 0.75 (night),6.75 (morning),12.75 (noon),and 18.75 (evening).Each of the panels contains four curves representing months: M = 0.75 (January), 3.75 (April), 6.75 (July), and 9.75 (October), as zO parameter is fixed to 10.The symbols are shown on the top left of each panel.
All curves have a maximum near the equator and a minimum around the poles.Seasonal variations are significant, especially at noon.The maximum of most of the ratios exceeds 20, which is higher than the maximum of the old model.Individual ratios may exceed 60, as seen from Figure 5, which means that at certain condition the plasmaspheric density is almost constant with altitude.
Next Figure 8 compares the ratios at selected four local times, as each panel represents the conditions at the four selected months.As expected, the daytime ratio is always smaller than the nighttime ratio and its maximum is shifted to the winter hemisphere.Figure 9 illustrates the dependence on the O + density at the transition height.Each panel shows the basic variations of R p along geomagnetic latitude for 3 values of zO: 9, 10, and 11, for noon and night at the months of January and July.The O + density at h T affects the ratio mostly during the day.

Discussion
The empirical model of the ratio R p is based on the individual ratios of plasmasphere and topside ionosphere scale heights.The model is free from any theoretical assumptions and its accuracy depends only on how well the data are approximated with the chosen basis functions.The fact that the data do not  1A, the overall model error of 40% is considerably lower than the error of 55% calculated for the old model, which is onedimensional.Actually, the R p data exhibit a slight double-crest distribution around the geomagnetic equator, as seen in Figure 5, which is an effect of the visible increase of H T at the equator.The models with 4 and even more coefficients usually predict this structure, but they are also able to produce additional maxima and minima at other latitudes which cannot be controlled.After extensive inspection of the models behavior, we finally decided to rely on the simplest model (3, 3, 3, 2) which behaves steadily along all parameter axes.
The model concerns the ratio H p /H T instead of H p itself mainly because it is designed to replace the old model of R p in the TaD profiler and to be applied in operational space weather products in the DIAS system and in the European Ionosonde Service of the European Space Agency Space Situational Awareness Programme (Belehaki et al. 2015).A model of H p solely, in our view, is not so useful if the lower boundary of plasmasphere (namely h T ) and the density there are not known.The main contributor to the total electron content (TEC) is the F region of ionosphere where H T plays a key role.The adjustment procedure (Kutiev et al. 2012), incorporated into the TaD software, adjusts the integral of the TaD profile to the measured GNSS TEC at the same location by varying the value of H T .The transition height h T can also be determined by H T through the ratio R T = H T /h T .In this sense, H T is a driver of the whole TaD profile and determining the plasmasphere scale height H p again through H T is a reasonable approach.Using Digisonde bottomside profiles, the TaD profiler reconstructs accurately the N e profiles up to the GNSS heights.The model coefficients, provided in Appendix 1, could be helpful for users aiming to reconstruct N e profiles in the ionosphere and plasmasphere.Kutiev et al. (2009b) have shown that in the equatorial and low latitudes the ratio R p had a shape similar to that of magnetic field lines and can be well approximated with cos 2 (glat).The present model has generally the same form, but the cosine function is scaled by a factor that depends on four parameters (in the old model this factor is a constant, equal to 9).The R p behavior is physically meaningful, showing that inside the inner plasmasphere (L < 3) the electron density is higher and often its height distribution is almost constant with altitude.At higher geomagnetic latitudes (above ±50°) the plasmaspheric density decreases rapidly with altitude due to the lower scale height.

Summary
The new empirical model of the plasmaspheric to the ionospheric scale heights ratio is based on topside sounder profiles measured onboard the ISIS-1 satellite.It approximates the ratios of independently extracted H + and O + scale heights from each individual profile as a function of geomagnetic latitude, month of the year, local time, and ln(O + ) density at the transition height.The first three model parameters are approximated by trigonometric base functions, while the last parameter is approximated by a polynomial function.By fitting the data with various numbers of terms of the base functions, a set of models were produced providing curves with different number of extremes.After extensive inspection of the different models, the simpler model, representing the data by a simple sinusoid along each axis, was favored for further analysis.The new model is designed to replace the old one-dimensional model of the ratio currently used in the TaD profiler.The overall model error of the new model is 40%, significantly smaller compared to the 55% of the old one.

Fig. 2 .
Fig. 2. Deriving H T and H p .Red dots: ISIS-1 measured N e profile.Dashed lines represent O + and H + distributions decomposed from the measured profile and the green line is their sum.Horizontal gray line shows the height h T at which the O + density becomes equal to N e /2.Density is expressed in cm À3 .

Fig. 3 .
Fig.3.Topside scale height H T (green dots) and H + scale height H p (red diamonds) versus geomagnetic latitude.H p is scaled on the right.In all figures H T is marked as HT.

Fig. 5 .
Fig. 5.The R p values extracted from the database (blue dots) and corresponding R p models along geomagnetic latitude.With the red dots, we provide the model (3, 3, 3, 2) R p values (to the left) and the old R p model, 9cos 2 (glat) + 4 (to the right).

Fig. 4 .
Fig. 4. The ratio R p calculated by the eight models shown in the Appendix 1 as a function of geomagnetic latitude.Plots represent the models at the four selected months, zO = 10 and LT = 12.75 (noon).

Fig. 7 .
Fig. 7. R p ratio versus geomagnetic latitude for 4 months of the year.The panels represent seasonal variations for morning, noon, evening, and night conditions.

Fig. 8 .
Fig. 8. R p ratio versus geomagnetic latitude for four local times.The panels represent the local time curves for four different months of the year.

Fig. 9 .
Fig. 9. R p ratio versus geomagnetic latitude for 3 values of zO.The panels represent zO curves for noon and night in January and July.