Imaging of fast moving electron-density structures in the polar cap

The imaging of fast-moving electron-density structures in the polar cap presents a unique set of challenges that are not encountered in other ionospheric imaging problems. GPS observations of total electron content in the po- lar cap are sparse compared to other regions in the Northern Hemisphere. Furthermore, the slow relative motion of the satellites across the sky complicates the problem since the velocity of the plasma can be large in compar- ison and traditional approaches could result in image blurring. This paper presents a Kalman-filter based method that incorporates a forward projection of the solution based on a model plasma drift velocity field. This is the first time that the plasma motion, rather than just integrations of electron density, has been used in an ionospheric imaging algorithm. The motion is derived from the Weimer model of the electric field. It is shown that this novel approach to the implementation of a Kalman filter provides a detailed view of the polar cap ionosphere under severe storm conditions. A case study is given for the October 2003 Halloween storm where verification is provided by incoherent scatter radars.


Introduction
Imaging the electron density in the Earth's ionosphere presents a particular challenge at high latitudes.This is for two reasons: firstly the electron density can be moving at a large enough velocity to cause 'blurring' (Watermann et al., 2002) and secondly the ground-based datasets result in a more sparse geometrical coverage than at lower latitudes.This results in two requirements for an imaging algorithm beyond those needed for imaging the mid-latitude ionosphere.Firstly, the algorithm must work on short time-windows of data or it must accom-modate movement within the imaging.Secondly it also needs to compensate for missing data in a self-consistent manner.While much of the early literature on ionospheric imaging deals with the incomplete and biased geometry of the measurements (e.g., Yeh and Raymund, 1993), this paper introduces a new algorithm that addresses the motion of the plasma and the sparsity of data in the high-latitude regions.
A number of authors have applied Kalman filters to accommodate time changes in ionospheric imaging Wilson et al. (1995) represented the ionosphere as a shell at a fixed altitude and Juan et al. (1997) applied a Kalman filter to full 3D imaging.Howe et al. (1998) showed simulation results that combined and extended ideas from the vertical-slice ionospheric tomography work of Fremouw et al. (1992) and the two-dimensional horizontal mapping of GPS TEC Wilson et al. (1995) into a full three-dimensional algorithm.Bust et al. (2001Bust et al. ( , 2004) ) proposed that meteorological forecasting concepts could be applied to ionospheric specification through variational analysis and formed an exponential decay term to incorporate time evolution into their GPS imaging.More recently several other publications (e.g., Hajj et al., 2004;Schunk et al., 2004) have proposed other forms of Kalman filters to specify the time dependent character of the ionosphere in data assimilation.
The Kalman filter method itself provides a statistically optimal solution to a linear inverse problem where estimates of the covariance in both data and a priori models are known.A significant problem with the approach is that such information is not readily available for the ionosphere and prior estimates of the ionospheric electron density from models are often systematically biased.In many cases a direct inversion solution to the inverse problem using a mapping to a set of orthogonal basis functions can be used (e.g., Howe et al., 1998).The use of orthogonal functions to specify the vertical profile of electron density in GPS ionospheric imaging successfully is demonstrated in Mitchell and Spencer (2003).Spencer et al. (2004) showed that a reformulation of the Kalman filter approach using orthogonal profiles with electron-density gradients rather than absolute values can be used successfully in imaging the electron density over the mainland U.S.A.
This paper introduces a new a priori constraint for ionospheric imaging; the incorporation of a model of the plasma motion into a Kalmanfilter based ionospheric imaging technique.This is necessary for imaging fast-moving polar cap structures such as those found under certain storm conditions where the interplanetary magnetic field is southward.These structures are forced by the electric and magnetic fields in the polar cap to move at a velocity given by the cross product of the two fields.Using a model of the electric potential in the polar cap, provided by Weimer (1995), an estimate of this drift velocity was used to forward project the estimate of the ionosphere at each point in time.This approach avoids the need for a prior electron density model for the polar cap since the 'model' is constructed from a projection of the prior solution.
The following section provides a detailed description of the theory of the method.Results for the storm of 30 October 2003 are presented followed by a discussion of these results and a comparison with data obtained from two inco-herent scatter radars: the European Incoherent Scatter Radar (EISCAT, 69.6°N, 19.2°E) and the EISCAT Svalbard Radar (78.2°N, 15.8°E).

Theory
Dual-frequency GPS observations provide an estimate of the integrated total electron content along a path from satellite to receiver.Representing the ionosphere as a set of three-dimensional voxels, the line integral from each of the receivers is approximated as a summation.Figure 1 shows a typical example of the ray path intersections with a shell at 400 km altitude over a period of ten minutes from the 61 receivers used in this study.The horizontal size of the voxels of 2°i n latitude and longitude is chosen heuristically on the basis of the spatial density of information available as shown in fig. 1.In order to ensure the sizes of the voxels are approximately uniform across the polar cap the grid is defined over the equator and then rotated so as to be centered over the geomagnetic pole.Menke (1989) gives an excellent discussion of resolution issues in tomographic imaging applied to geophysical systems.For inversions using ground-based GPS receivers the lack of horizontal ray paths through the ionosphere severely limits the resolution of the technique in the radial, hence the reformula- tion of the problem using an orthogonal basis set to represent the vertical profile.A value of 40 km for the radial size of the voxels was chosen since it is sufficient to ensure accurate numerical estimates of the ray path integrals.
The problem of reconstructing the electron density field from the ray path integral observations may be represented as a linear inverse problem where, H defines the per voxel path integral contributions for each ray, x the unknown electron densities and z the observed line integrals.A single update of the Kalman filter from time t-1 to time t is expressed mathematically as (2.6) The state transition, or forward, model matrix A defines how the prior solution is projected into the future.P represents the error covariance of x and Q the error covariance of the state transition model.Given the matrix H and its associated covariance R, the Kalman gain, K, is estimated enabling the solution for x and P at time t to be obtained.Variances in the observations, R, are estimated from the variance of the inter-frequencybias (IFB) corrected differential pseudo-range from the differential-phase observations.The determination of the IFBs and their effect on the TEC measurement is addressed in Dear and Mitchell (2007).The error covariance matrices P and Q are band limited sparse matrices with a non-zero term defined over 3 neighboring pixels in latitude and longitude.The error covariance in Q is defined by a truncated Gaussian function with the variance set to 5 TECU (in accordance with Dear and Mitchell, 2007) and a half-width (t) of 200 km.The correlation distance of 200 km is chosen so as to be the same order as the spa-

Hx z =
tial density of information available as shown in fig. 1.There is no literature we can find specifically relating to the issue of correlation distance in the polar cap and obviously, the local correlation will be spatially and temporally variable.
The covariance in the radial direction is set to zero since the stability of the algorithm in the radial direction is provided using an orthogonal mapping described below.
Typically in the application of the Kalman filter approach to the ionosphere the forward projection of the state is carried out using a prior model estimate of the three-dimensional electron-density field.Such an approach was considered unsuitable for imaging the polar cap where electron-density models are unable to predict the distribution of localized fast moving structures. he method developed here involves using a projection of the prior solution based on an estimate of the velocity of the plasma due to the E×B drift.Weimer (1995) provides a model of the electric field, E, in the polar cap as a function of date, time, the solar magnetic field components By, Bz and solar wind velocity.The latter three parameters can be obtained from ACE satellite data enabling the changing electric field in the polar cap to be modeled as a function of time.Combined with an estimate of the Earth's magnetic field, B, the velocity of the plasma is given by . (2.7) For the current approach the Earth's magnetic field over the polar cap was assumed vertical with a magnitude, B, of 50000 nT.The plasma drift velocity will then be in the local horizontal only.The forward projection A i,j at point r for a time step between updates of the filter of ∆t is defined as .(2.8) Thus the forward projection at a point is given by a Gaussian convolution of neighboring voxels in the local horizontal at an offset distance given by an estimate of the distance the plasma moves between filter updates.No convolution is applied in the radial direction.The initial starting estimate of the electron density field used in the reconstructions is the IRI 95 model (Bilitza, 1997).
Representation of the ionosphere as a discrete set of voxels is computationally expensive considering the high spatial resolution required to image structures over the polar cap.Consequently the state vector, x, is mapped to a set of empirical ortho-normal functions, EOFs, in the radial direction obtained from a set of Chapman functions with peak heights ranging from 250 km to 550 km.From these model profiles a set of two EOF's was derived using singular value decomposition.An ortho-normal matrix, M, is constructed to enable the mapping from a space of discrete voxels in three-dimensions to a space of discrete EOFs with independent EOF's for each latitude and longitude in the grid.Following the approach outlined in Spencer et al. (2004), the matrices required to describe the filter state in the domain of the basis functions, denoted using a dash, are given by (2.9) (2.13) After updating the filter the voxel electron densities are then recovered from the mapped solution using . (2.13)

Results
Figures 2 and 3 show comparisons between the vertical electron density obtained from two incoherent scatter radars and those obtained from the inversions.The vertical axis ranges from an altitude of 100 km to 800 km and the horizontal axis shows universal time for 30 October 2003.In both cases it can be seen that extreme enhancements in electron density occur between 18 UT on 30 October and 00 UT on 31 October.These enhancements are also associated with large peak heights often in excess of 500 km.It is interesting to note the good agreement in the large-scale features, for example the pre-noon enhancement and the successive passing of enhancements during the evening.Some smaller scale features observed in the ISRs, at  the sub-grid level, are not expected to be imaged by our method.The sequence of images in fig.4a-h show both the vertical total electron content and volume visualizations of the evolution of the elec-tron-density field over a period from 18 UT on 30 October to 00 UT on 31 October.The dramatic increase in the height of the plasma over the US mainland can be clearly seen to the left of the images.This uplifted plasma is then convected around the polar cap towards Northern Europe.Figure 5 illustrates the solution obtained at 22 UT from a separate run of the filter using no input observations and is shown for comparison with fig.4g.The lack of significant cross-polar flow and peak height enhancement clearly illustrates that the results presented are primarily data driven and not the product of prior assumptions in the model.

Discussion and future work
The October 2003 storm has warranted much attention in the literature and was featured in a special issue of Geophysical Research Letters (for example see Mitchell et al., 2005).In this pa-   per we did not aim to address the science; the storm is being used as an extreme example of plasma convection to demonstrate a technique to image the polar cap ionosphere.Nevertheless, our results do provide further evidence that convecting plasma appearing over Europe (Stolle et al., 2006) and causing GPS scintillations (as reported by Mitchell et al., 2005) originated from the storm enhanced density observed in the American sector.This supports results by Foster et al. (2005), who have shown using ISRs and other observations that another significant storm on November 20 2003 exhibits a flow of plasma from the U.S.A. mid-latitudes across the polar cap into the European sector.This preliminary result for the new algorithm presented here is very promising but it remains to be tested on a number of different cases.Particularly important is the fact the method is able to image the extreme heights obtained by the plasma under storm conditions as seen by ISR.This result showing the peak height uplift is demonstrated to be data (not model) driven by the auxiliary run of the imaging technique here using only the initiation and convection models.
There is no reason to limit the application of this new algorithm to storm conditions -in fact the technique will help to image the polar cap for all the recent solar maximum even in cases of slower convection.Improvements to the results are envisaged by the incorporation of velocity vectors from the SuperDARN radars (Greenwald et al., 1995) or from an empirical model of the electric field (Ridley et al., 2000) into the algorithm.It is intended that the new technique presented here will serve as a useful complement to other ionospheric instruments, allowing a three-dimensional visualization of the large-scale plasma over the entire Northern Hemisphere under all ionospheric conditions.

Fig. 1 .
Fig. 1.Ray path intersections with a 400 km altitude shell for a typical inversion.Data color coding represents slant TEC in TEC units.

Fig. 2 .
Fig. 2. Electron density as a function of height and universal time from the EISCAT radar (69°N, 19°E), above, and inversion, below.

Fig. 3 .
Fig. 3. Electron density as a function of height and universal time from the ESR radar (78°N, 16°E), above, and inversion, below.

Fig. 5 .
Fig. 5. Volume visualization of electron density at 22 UT obtained from a separate run of the filter using no input data.Iso-contours at 5E11m -3 , transparent, and 10E11m -3 , opaque, are shown (for comparison with fig.4g).

Fig
Fig. 4e-h.Volume visualizations of electron density from the evening of October 30 2003.Iso-contours at 5E11m -3 , transparent, and 10E11m -3 , opaque, are shown.The images are presented at two hourly intervals starting from 18 UT.