Automated identification and tracking of polar-cap plasma patches at solar minimum

A method of automatically identifying and tracking polar-cap plasma patches, utilising data inversion and feature-tracking methods, is presented. A well-established and widely used 4-D ionospheric imaging algorithm, the Multi-Instrument Data Assimilation System (MIDAS), inverts slant total electron content (TEC) data from groundbased Global Navigation Satellite System (GNSS) receivers to produce images of the free electron distribution in the polar-cap ionosphere. These are integrated to form vertical TEC maps. A flexible feature-tracking algorithm, TRACK, previously used extensively in meteorological storm-tracking studies is used to identify and track maxima in the resulting 2-D data fields. Various criteria are used to discriminate between genuine patches and “false-positive” maxima such as the continuously moving day-side maximum, which results from the Earth’s rotation rather than plasma motion. Results for a 12-month period at solar minimum, when extensive validation data are available, are presented. The method identifies 71 separate structures consistent with patch motion during this time. The limitations of solar minimum and the consequent small number of patches make climatological inferences difficult, but the feasibility of the method for patches larger than approximately 500 km in scale is demonstrated and a larger study incorporating other parts of the solar cycle is warranted. Possible further optimisation of discrimination criteria, particularly regarding the definition of a patch in terms of its plasma concentration enhancement over the surrounding background, may improve results.


Introduction
Polar-cap plasma patches are enhancements of the ionospheric F layer, most often but not solely defined as having a maximum electron concentration at least twice that of the surrounding background. In some past studies an absolute difference between background and patch maximum electron concentration has been used Heelis, 1995, 1998;Dandekar, 2002;Hunsucker, 2002). An additional criterion of a minimum gradient in electron concentration at the edge of the patch has also been used (Coley and Heelis, 1995). Patches have been observed more frequently in the Northern Hemisphere winter than in summer regardless of whether an absolute or relative definition of a patch is adopted (Coley and Heelis, 1995;Dandekar, 2002). Patches have been observed on scales from as small as 50 km up to 1000s of kilometres Hunsucker, 2002). Patches are at least approximately circular in the locally horizontal dimensions and have a lifetime of the order of hours (Hunsucker, 2002). Evidence has been put forward in support of a weak relationship between the Kp index of planetary magnetic field disturbance and the frequency of formation of patches; see Dandekar (2002) and references therein.
In the polar-cap ionosphere, plasma drifts according to the formula V = E × B/B 2 , where V is the velocity of the plasma, E is the electric field strength and B is the magnetic field strength. This is referred to as an E × B drift (Chen, 1984). It is driven by the Earth's magnetic field and an electric field mapped from the magnetosphere. Once patches form, they move as dictated by this E × B drift. In general, patch motion can be described as approximately antisunward (Hunsucker, 2002). The velocity field formed by 198 R. Burston et al.: Automated identification and tracking of polar-cap plasma patches the E × B drift is generally known as a circulation pattern and it varies according to ambient conditions, being strongly affected by the variable inter-planetary magnetic field (IMF). Diagrams of idealised circulation patterns under various conditions of the B z and B y components of the IMF can be found in Hunsucker (2002), along with further details of patch dynamics and references to the primary literature.
Tracking of patches is of practical as well as scientific interest. Patches are associated with the phenomenon of ionospheric scintillation, which is potentially disruptive to satellite-to-ground radio links. This in turn can negatively affect the performance of Global Navigation Satellite Systems (GNSS) used by numerous industries in polar regions, including safety-critical ones such as civil aviation (Skone et al., 2001). Knowledge of the presence, locations and scale of patch activity would be greatly beneficial to these industries for operational planning purposes.
In past studies patches have been tracked using data from coherent scatter radars and incoherent scatter radars (ISR) as well as ionosondes; however the patches themselves have been identified by humans (Clausen et al., 2012;Grant et al., 1995;Mori et al., 2012;Oksavik et al., 2006;Pedersen et al., 1998). An algorithm for detecting patches automatically from Dynamics Explorer 2 satellite data has also been developed (Coley and Heelis, 1995). Since the orbital period of Dynamics Explorer 2 was 98 min and the relevant measurements were in situ, the algorithm could not follow the history of any given patch in detail but rather detected any patches that the satellite passed through. Dynamics Explorer 2 reentered the lower atmosphere in 1983.
An algorithm for automated detection and tracking of patches in data from a single all-sky camera (ASC) using a cross-convolution technique was developed by Hosokawa et al. (2009). The limitations of ASC instruments mean that the algorithm can only be used when the sun and moon are below the horizon and there is no cloud in the field of view. The field of view of any single ASC will not cover the whole of the polar cap. A false-positive detection rate using this method of 5-10 % was also reported.
A previous study  traces the flow of plasma in the polar cap using the output from an ionospheric imaging algorithm, Ionospheric Data Assimilation Three-Dimensional (IDA3D). This algorithm uses data assimilation to produce images of the electron concentration of the ionosphere. This method can give results for the whole polar cap for any time period when sufficient input data are available. In practice this means at least all times since the period studied in the paper (December 2001). Subsequently, Foster and Evans (2009) demonstrated a method of determining a velocity field from any time series of vertical total electron content (vTEC) fields by applying block-based motion estimation. The technique was developed for studying tongues of ionisation, but not other phenomena. Another approach, using optical flow constrained by an ionospheric model, generates a velocity field for the whole image area from a time series of vTEC images (Benton and Mitchell, 2012). The methods discussed in this paragraph all attempt to determine plasma motion in general, not the motion of plasma patches as discrete entities.
The study reported in this paper demonstrates the feasibility of a fully automated technique for identifying and tracking the motion of polar-cap plasma patches. The approach is based on ionospheric imaging combined with feature tracking in spherical data fields. The patches are identified from the output of an algorithm, the Multi-Instrument Data Assimilation System (MIDAS). MIDAS is a tomographic imaging algorithm, with a ten-year history at the time of writing (Mitchell and Spencer, 2003), that has been applied to many ionospheric studies, for example Alfonsi et al. (2011), Muella et al. (2011, Rose et al. (2009Rose et al. ( , 2011 and Yin and Mitchell (2011). It solves an under-determined four-dimensional inverse problem in order to compute a gridbased image of the electron concentration of the ionosphere in a user-selected region. In principle any data that can be put in the form of an in situ electron concentration measurement or a slant total electron content (TEC) measurement with known time and position or path respectively can be assimilated. GNSS signals are ideal for this purpose as reliable observations such as those provided by the International GNSS Service (IGS) are readily available. Since the geometry of the problem and the limitations of data collection ensure that the problem is under-determined, additional constraints are required in order to obtain a unique solution in any given circumstance (Austen et al., 1988). The extant MI-DAS literature shows that the algorithm is capable of identifying structures of approximately 500 km and greater in scale including plasma patches (Burston, 2012;Yin and Mitchell, 2011;Yin et al., 2008Yin et al., , 2009. The TRACK algorithm, which has a nearly 20-year history, identifies and tracks features in diverse spherical data fields (Hodges, 1994(Hodges, , 1995(Hodges, , 1999. It is used to find and follow the patches, in a completely automated fashion. This is different from the approach taken in Bust and Crowley (2007), which back-traces the motion of plasma packets that form patches rather than automatically identifying patches in the ionospheric images and tracking them. This back-trajectory method is often used in atmospheric sciences for following air parcels to examine the origins of water vapour or pollution (Methven et al., 2001). As with Bust and Crowley (2007), the new method described here allows the whole polar cap to be treated, and the only limit on the time periods that can be studied is a sufficiency of input data. Since MIDAS and IDA3D both principally use GNSS ground receiver data, it is expected that the limit for MIDAS is approximately the same as that for IDA3D.
The TRACK algorithm has been formulated to be applicable to a wide range of phenomena but has mainly been used in the atmospheric sciences to track tropical (Bengtsson et al., 2007) and extra-tropical cyclones (Hoskins and Hodges, 2002). Because of the substantial extant literature on applying both MIDAS and TRACK, neither can be considered new; the innovation in this paper is their combined use, which has not been reported in the peer-reviewed literature before.
A proof-of-concept feasibility study covering the Halloween Storm period (late October 2003) was presented at the Institute of Navigation (ION) GNSS 2010 (Burston et al., 2010). It shows that care must be taken to eliminate falsepositive results and was used to define the selection criteria adopted here. These criteria are described in the Method section.
A 12-month period, March 2007-February 2008, inclusive, is studied. The period was chosen to coincide with a series of measurements taken by the European Incoherent Scatter (EISCAT) ISR at Longyearbyen, Svalbard, which were used to validate the MIDAS results.
The tracking technique proves successful, identifying 71 separate patches, in a 12-month period (March 2007 to February 2008, inclusive) at solar minimum but cannot identify patches of less than approximately 500 km in diameter because of the resolution limit of MIDAS. The results show few patches in the period March through June 2007 inclusive and peaks of activity in July and October 2007 and February 2008. This does not conform to the usually observed seasonal bias during the first half of the 12 months under study. However, 71 tracks are probably too few to draw conclusions regarding seasonal occurrence. A clear bias toward patch initiation over North America as compared with Europe and Asia is observed. There is no correlation with Kp, but the index never rises to 6 in the period studied.
This technique offers the possibility of improved plasma patch climatology and a method, independent of coherent scatter radar, of observing patch flow patterns. It also offers a possible measure of ∇n/n , a critical parameter in growth rate equations for instabilities proposed as the primary cause of scintillation associated with patches. A much larger database of patches, from across all parts of the solar cycle, is required in order to perform all these analyses.

Method
MIDAS exists in two main variations, as described in Mitchell and Spencer (2003) and Spencer and Mitchell (2007), henceforth referred to as the 2003 version and the 2007 version. Despite being specifically developed for imaging the polar-cap ionosphere, the 2007 version is not used here. This is because the 2007 version uses a Kalman filter to balance the final solution of the algorithm for any particular epoch between the result of the data inversion and a projection from the previous solution calculated by using the Weimer model. The Weimer model is a statistical model of the polar-cap E × B drift pattern (Eriksson et al., 2002;Weimer, 1995Weimer, , 1996Weimer, , 2001aWeimer, , b, 2005aWinglee et al., 1997). The use of such a solution would make inferences about the motion of patches difficult in that it would be hard to ascertain how much influence the Weimer model had on them. The fact that the solution from the Weimer model is weighted in the ultimate solution by a Kalman filter further complicates the matter as the dependence would vary with every epoch. The version of MIDAS used in this study is the 2003 version with the variation that constraints in the locally horizontal dimensions are imposed by a convolution smoothing technique rather than fitting to spherical harmonics. Because the 2003 version makes no underlying assumptions about the nature of the polar-cap circulation pattern, the results presented can be taken as inferred purely from the input data.
The grid was defined as having divisions every 4 • of latitude and every 4 • of longitude. The grid was defined at the Equator as being bounded at −50 and +50 • latitude and −50 and +50 • longitude then rotated so as to be centred on the geomagnetic north pole. The vertical extent of the grid was from 100 to 980 km in altitude with a division every 40 km. Data from IGS permanent GPS ground receiver stations within the grid, operating during the relevant period, were used as input to the MIDAS algorithm. The ground area of the grid and the locations of every IGS station used are shown in Fig. 1. MIDAS was run for the 12 months 1 March 2007-29 February 2008, inclusive. The time step (epoch) length for the MIDAS run was 10 min. The algorithm solves a four-dimensional problem using an adjustable number of epochs on either side of the target epoch: in this case, seven epochs before and seven epochs after the target epoch, totalling 15 epochs. Note that temporal resolution is determined by the individual epoch length, not by the total time length. This is because resolution is determined by voxel volume, not by overall grid volume. In the 4-D problem voxels have a 10 min length in the time dimension and the overall grid length in the time dimension is 15 epochs × 10 min = 150 min. This means that the input data used in the solution for each time step are not 10 minutes' worth but 150 minutes' worth. Figure 2a shows the input rays for one arbitrarily selected epoch. Figure 2b shows all the rays used to solve for that epoch.
Observations of the highest value of electron concentration in the F2 layer of the ionosphere (NmF2) as observed by the EISCAT ISR at Longyearbyen, Svalbard (78.15 • geographic latitude, 16.02 • geographic longitude), were used to check that the MIDAS reconstruction was reasonably accurate. The data are from the International Polar Year programme, which was run using the 42 m dish. Technical specifications for this ISR installation can be found in Rottger et al. (1995) and Williams (1995). This data set was chosen because it offered the most complete and continuous temporal coverage available from any ISR for the 12-month period studied, whilst operating on an unchanging set of parameters. The highest value of electron concentration in the column of voxels over Svalbard was extracted from the MIDAS grid for each time step and compared with the ISR observation. If more than one ISR observation was available in a given time step, 32 the mean was taken and used for the comparison. Figure 3 shows the availability of the ISR data during this period.
Although the MIDAS algorithm produces a 3-D grid of electron concentration values as output, and the TRACK algorithm can be applied at different levels in a 3-D data field and detect structures in the vertical dimension, it is simpler to establish the technique's effectiveness using 2-D data fields. Hence the 3-D MIDAS grids for each epoch were vertically integrated to give a 2-D field of vTEC values using the formula vTEC = nds, where n is the electron concentration and s is the path length.
The TRACK algorithm works by first identifying features of interest; this can be done in a number of ways, but here the following approach is used. First the individual vTEC fields are spectrally filtered to suppress noise and to remove the large-scale background, which makes identification easier. The spectral filtering uses a 2-D discrete cosine transform (Denis et al., 2002) to transform to spectral space where the spectral coefficients are truncated to provide the smoothing and the low-order coefficients are set to zero to remove the large-scale background before transforming back to real space on the original grid. The spectral band chosen is representative of a spatial resolution of 500-5000 km. The TRACK algorithm processes the filtered data sequentially by first applying a threshold to the data, a value of 0.5 is used, and the grid points above this threshold are grouped into objects using a connected component labelling algorithm (Hodges, 1994) that searches for grid point adjacencies to label points with the same label that belongs to the same object. Each object is then searched for local extrema, in this case maxima, a complex object can have more than one extremum. The tracking is performed by first initialising a set of tracks by linking the identified maxima using a nearestneighbour search, subject to a maximum displacement of 8 • geodesic degrees. For a 10 min time step this is equivalent to a speed of 1.5 km s −1 . The tracks are then refined by minimising a cost function for track smoothness. The cost function is formulated in spherical geometry (Hodges, 1995) in terms of changes in direction and speed, determined over three consecutive time steps, and with weights used to combine these measures of direction and speed change. The total cost function is the sum of these, over all tracks. Full details of the algorithm can be found in Hodges (1994Hodges ( , 1995Hodges ( , 1999. In order to avoid false-positive results, the tracks are postprocessed to retain only those that last longer than 6 time steps (1 h); also tracks within 5 • of the domain boundary are removed. To further isolate the tracks of interest they are referenced to, the vTEC field converted to a percentage of background and those tracks that do not satisfy the condition that the intensity exceeds twice the background for at least 4 consecutive track points are discarded. Finally the track point locations are transformed back to true longitudes and latitudes, and only those tracks that pass within 10 • of the north magnetic pole are selected, taken as (84.2 • N, 124.9 • W), which is suitable for the 2007/2008 period used. Ideally a higher spatial and temporal resolution of the vTEC would be beneficial, but the data production is limited by the availability of the GNSS rays. The tracks can then be used to compute a range of statistics including the spatial distribution, which is the main focus here. This is done using spherical kernel estimators as described in Hodges (1996). Illustrations of each stage of the TRACK process for identifying plasma patches from vTEC data can be seen in Burston et al. (2010).
Polar-cap plasma patches are usually defined as having a maximum plasma concentration of at least twice that of the surrounding background. This definition has been adapted for the case of a 2-D field of vTEC values simply to be a maximum in the vTEC field at least twice that of the background. In turn this requires a definition of the background value. This was calculated individually for each epoch in the following manner: first each value in the field was rounded to the nearest TEC unit (TECu), where 1 TECu = 10 16 electrons m −2 . Next the modal average value is determined. This value is then adopted as the background vTEC value.

Results
There are 52 704 epochs in the MIDAS run. This means that direct graphical comparison of the ISR NmF2 values with the MIDAS values is difficult to interpret. Hence an alternative method of comparison, utilising the mathematical method of  convolution, has been adopted. Let x(n) be the series of values of NmF2 recorded by the ISR and y(n) a series of corresponding MIDAS NmF2 values. Then the curve is the convolution of the ISR NmF2 series with itself. This may then be compared with the curve the cross-convolution of the ISR NmF2 series with the MI-DAS NmF2 series. In Eqs.
(1) and (2), is the epoch number, n is the total number of epochs, x is an ISR observation of NmF2 for a given epoch and y is a MIDAS value of NmF2 for a given epoch. Only epochs where there is an ISR value are included in these series. The closer the crossconvolution result, Z, is to the self-convolution result, X, the more closely the MIDAS values of NmF2 match those from the ISR. The self-convolution of the ISR results is compared with the cross-convolution of the ISR values with the MIDAS values in Fig. 4. It can be seen that, for the most part, the ISR-MIDAS curve is higher than the ISR-ISR curve. This implies that MIDAS is giving larger values for NmF2 than are observed by the ISR for a considerable number of epochs. Detailed inspection shows that MIDAS, in this study, consistently slightly overestimates the nighttime minimum in NmF2, thus explaining the convolution results. This slight bias is unlikely to affect the patch detection and tracking results.
The values of the polar-cap background vTEC calculated in the manner described above ranges from 1 to 21 TECu. The distribution of values is shown in Fig. 5a. Figure 5b shows the distribution of background vTEC values as a percentage of the total number of epochs and the cumulative percentage of epochs. A total of 99.4 % of epochs have a background value of 10 TECu or below; 99.9 % of epochs have a background value of 12 TECu or below.  There are very few patches detected in the first four months. Whilst there are some patches detected during November-January, the major peaks of activity are in July, October and February. This is illustrated in Fig. 6, which shows the distribution of patch occurrences, as identified by TRACK, for the 12 months under study.
An arbitrarily selected sequence of MIDAS images showing a patch from 12 May 2007 is shown in Fig. 7. The red marker indicates the location of the patch maximum as identified by TRACK. This indicates that the correspondence between patches as reconstructed by MIDAS and as identified by TRACK is good.
The Kp index never rises to 6 during the study period. The mean value is 1.6. Figure 8a shows the distribution of epochs of a given Kp value for the entire 12-month period. The distribution of Kp values at times when a patch is first identified by TRACK is shown in Fig. 8b. The mean value is 2.2. Figure 8c shows the percentage of epochs with a given Kp value during which a patch is first identified by TRACK. The value, R = 0.1, shown in panel c, is the linear correlation coefficient for the distributions shown in panels b and c. The low and high values shown are the corresponding 95 % confidence limits for R. A similar analysis was conducted for the Auroral Electrojet (AE) index. Figure 9a shows the distribution of hourly epochs of a given AE index value for the entire 12-month period. The distribution of AE index values at times when a patch is first identified by TRACK is shown in Fig. 9b and 9c. Panel c shows the percentage of epochs with a given AE index value during which a patch is first identified by TRACK.  confidence limits for R. A quite good anti-correlation exists. This is discussed in the Conclusions.
The spatial concentration of all identified tracks is shown in Fig. 10. The unit area used is a spherical cap with a 5 • arc angle. Results are for the whole 12-month study period. Red markers indicate the locations where maxima are first identified by TRACK. The blue marker shows the location of the geomagnetic pole. There is a bias toward North American latitudes in these locations. The magnetic local time (MLT) of first identification of patches shows a bias against formation of patches in the 06:00 to 12:00 quadrant (Fig. 11).
There were insufficient tracks to meaningfully analyse relationships between patch formation and motion and the directions of IMF components. There are four combinations arising from B z and B y independently being either positive or negative, each with its own ideal convection pattern. The convection pattern that obtained cannot be derived solely from the trajectory of an individual patch, however. A large number of trajectories would be required for each set of B z /B y conditions so that a statistical analysis could be conducted. There proved to be too few available in each category from this study alone.

Conclusions
The MIDAS-TRACK method of automating detection and tracking of polar-cap patches has proved successful. Small patches (less than approximately 500 km in diameter) will not be detected by TRACK because of the limited resolution of the MIDAS grid. Patches are also not identified until they show a maximum plasma concentration twice that of the surrounding background. These factors mean that the time of first identification of a patch by TRACK is very likely to be significantly later than the actual time of inception of the patch.
The method used to determine a polar-cap background vTEC value gives realistic results except in a very small proportion of the 52 704 epochs: 0.01 % have a value greater than 12 TECu. The highest value is 21 TECu. It is probable that some of these higher values are inaccurate because the modal average figure is representing the region surrounding the polar cap rather than the polar cap itself. This is most likely to occur both in the summer when vTEC values will be at their highest and at the magnetically quietest times when the polar-cap ionosphere will be at its geographically smallest area. Fewer pixels of the vTEC field will then be contained within the polar cap, and hence the modal average figure is more likely to occur outside the polar cap. Given the extreme rarity of these aberrant values and the small number of tracks detected, the results from TRACK are unlikely to be significantly affected by them.
The seasonal distribution of patch occurrences shown in Fig. 6 exhibits a very low level of activity for the first four months and peaks of occurrence during July, October and February. This does not conform with previous studies such as Coley and Heelis (1995) and Dandekar (2002), but there are probably too few patches in total to draw any inferences. Figure 8 shows no linear correlation between the Kp value and the first identification of patches. This lack of correlation between Kp and patch occurrence is counter to the findings of Dandekar (2002) and other studies cited therein. The Dandekar (2002) study covers a total of 48 months in four 12-month segments each at a different part of the solar cycle. Here, results for a 12-month period at solar minimum, where Kp < 6 throughout, are shown, with only 71 patches in total. The lack of results for Kp ≥ 6 and the small number of patches detected mean that the two studies are not directly comparable. Only by similarly sampling a whole solar cycle could direct comparison between results from the method presented here and the results of Dandekar (2002) be made.
The quite good anti-correlation between AE index value and the first identification of patches has broad confidence limits (see Fig. 9), and it is probably unwise to draw inferences from it for the same reasons as given for the Kp index. Figure 10 shows a bias in the results towards patch formation in the longitudes of North America. Although the number of IGS stations contributing MIDAS input data is larger in northern North America than in northern Europe or northern Asia, Fig. 3 shows the ray coverage over the entire polar cap to be weakest directly over the Arctic Ocean. Since TRACK detects patches over the Arctic Ocean, where the least input data exist, it is concluded that this is a real bias and not an artefact of the data distribution. The offset between geographic and geomagnetic poles does not explain this bias as the distribution is still biased toward whichever pole is considered to be the centre of the pattern. Asymmetry of the polar-cap circulation pattern induced by co-rotation of the ionosphere is also not the cause of this observation as the pattern completes one rotation every 24 h and 366 such complete rotations are included in the analysis. This means that the asymmetry has no longitudinal bias. The MLT of first identification of patches conforms to the results of Coley and Heelis (1995), using Dynamics Explorer 2 data, but does not shed light on this bias. The phenomenon remains unexplained. Study over a greater total time period including all parts of the solar cycle and both poles is required in order to establish whether the phenomenon remains in a much larger data set.
Natural extensions of this work include application of the technique to other parts of the solar cycle and to the Southern Hemisphere. The consequent increased number of patch tracks will improve the climatological analysis. An investigation of results using different criteria for defining a patch must also be conducted -specifically, use of an absolute rise above the background threshold, e.g. 5 TECu, rather than the relative one (twice the background) used here.