A framework for estimating spherical vector fields using localized basis functions and its application to SuperDARN data processing

Abstract

A technique for estimating a plasma drift velocity distribution in the ionosphere is presented. This technique is based on a framework for representing a global vector field on a sphere by using a set of localized basis functions which is newly derived as a variant of the spherical elementary current system (SECS). A vector field on a sphere can be divided into its divergence-free (DF) component and curl-free (CF) component. The DF and CF components can then be represented by weighted sums of the DF and CF vector-valued basis functions, respectively. While the SECS basis functions have a singular point, the new basis functions do not diverge over a sphere. This property of the new basis function allows us to achieve robust prediction of the drift velocity at any point in the ionosphere. Assuming that the ionospheric plasma drift velocity has no divergence, its distribution can be represented by a weighted sum of the DF basis functions. The proposed technique estimates the ionospheric plasma drift velocity distribution from the SuperDARN data by using the DF basis functions. Since there are some wide gaps in the spatial coverage of the SuperDARN, an empirical convection model is combined with the framework based on the new basis functions. It is demonstrated that the proposed technique is useful for the estimation and modeling of the ionospheric plasma velocity distribution.

Introduction

The global ionospheric convection pattern is a fundamental property of the ionosphere and the magnetosphere, and many studies have conducted modeling of the global convection pattern. One notable empirical model was provided by Heppner and Maynard (1987), which improves on the earlier model of Heppner (1977) by using electric field measurements from Dynamic Explorer 2 (DE2). Weimer (2001) presented an improved electric field model based on the DE2 data as a function of solar wind parameters and geomagnetic conditions. Ruohoniemi and Greenwald (1996) derived empirical convection maps based on Goose Bay HF radar observations. However, since the ionospheric convection field varies greatly according to geomagnetic disturbances, these empirical models are not always guaranteed to well represent variations of the convection pattern.

In order to monitor the variations of the ionospheric convection pattern, global observation networks have been constructed over the last couple of decades. The Super Dual Auroral Radar Network (SuperDARN) (e.g., Greenwald et al. 1995; Nishitani et al. 2019) is one such global network, consisting of HF radars covering high latitudes, and recently it covers even the middle latitudes. Each radar of this network observes the backscatter from field-aligned irregularities of electron density in the ionosphere and obtains the line-of-sight component of the electron drift motion, which follows the $\varvec{E}\times \varvec{B}$ drift in the F-region ionosphere. The radar beam typically scans in 16 directions over a 52° sector every 1 or 2 min. The SuperDARN provides global information on the ionospheric plasma drift velocity distribution which is associated with the electric field distribution in the ionosphere. However, there are some wide gaps in the spatial coverage of the SuperDARN. In addition, each radar gives only the line-of-sight component of the drift velocity and data are frequently missing. Thus, it is not necessarily easy to retrieve a global convection pattern from the SuperDARN data.

There have been proposed several techniques for obtaining the variations of the global convection pattern. Richmond and Kamide (1988) proposed a method to estimate the variations of the global ionospheric potential distribution by combining various observations including HF radar observations and geomagnetic observations. Although their method is not necessarily for analyzing the SuperDARN data, it is applicable to obtain the ionospheric potential from the SuperDARN data. Ruohoniemi and Baker (1998) presented a method for obtaining the global convection pattern, which is specialized to use the SuperDARN data. The method for SuperDARN data was further improved by employing a set of empirical convection models based on data from the polar cap through the mid-latitudes (Thomas and Shepherd 2018). These methods basically use spherical harmonic functions as basis functions to represent the global potential distribution. However, since each spherical harmonic function represents some global pattern, spurious structures can be produced in the data gap region as a result of unevenly distributed data. In addition, it would be difficult to consider relatively small-scale local structures which could be resolved by the SuperDARN. Cousins et al. (2013b) used empirical orthogonal functions based on spherical harmonic functions (Matsuo et al. 2002; Cousins et al. 2013a), and Gjerloev et al. (2018) employed a similar method to obtain the global convection map. Although these approaches would improve an estimate in the data gap region by using empirical orthogonal functions, they may still suffer from the problems due to the use of global basis functions.

Amm (1997) proposed another approach based on a different set of basis functions called spherical elementary current systems (SECS), which was originally introduced for representing ionospheric electric currents. Each SECS basis function is symmetric with respect to its source point, and its value decreases as the distance from its source point increases. The problems of spurious structures in the data gap due to the use of spherical harmonic functions could be avoided by using the SECS basis functions. Amm et al. (2010) demonstrated SECS basis functions are useful for obtaining a regional velocity field from the SuperDARN data. Reistad et al. (2019a, b) performed a statistical analysis of ionospheric convection sources using SECS functions. However, an original SECS basis function has a singular point at which the value of the basis function diverges. This problem is unfavorable when using the analysis results as a reference to compare with other observations such as spacecraft observations. For example, if the footprint of a spacecraft comes to the neighborhood of a singular point by chance, the data from spacecraft observation cannot be compared with an SECS analysis result. Thus, in many applications of the SECS functions, it was essential to appropriately treat the singularities as discussed by Vanhamäki and Juusola (2020).

The present study introduces an alternative set of localized basis functions as a variant of SECS basis functions. The basis functions proposed in this study do not have singular points in order to enable us to obtain reliable estimates and predictions for ionospheric properties. We then propose a technique for estimating the global drift velocity distribution from the SuperDARN data based on the newly introduced basis functions. The use of localized basis functions would make it difficult to effectively predict the drift velocity in the data gap region. However, another existing model such as an empirical convection pattern model can be incorporated into a model with localized functions in order to fill in the data gap. There exists some other localized approaches. Fiori et al. (2010) employed localized spherical cap harmonic functions for analyzing the SuperDARN data. Bristow et al. (2016) proposed another localized approach which discretizes space into cells and assumes local divergence-free for each cell. On the other hand, our approach which is an extension of the SECS method is easy to apply and it would allow a flexible modeling of the ionosphere. Although a basic idea on the proposed technique is briefly described in the paper by Seki et al. (2018), the present paper provides a detailed explanation of the idea of the set of basis functions and their application to estimating the ionospheric drift velocity distribution from the SuperDARN data.

Spherical elementary current systems

Amm (1997) proposed sets of SECS basis functions for representing a vector field on a sphere. There are two sets of SECS basis functions. One comprises curl-free (CF) basis functions:

$$\begin{aligned} \varvec{v}^{\text{cf}}(\varvec{r}, \varvec{r}_i)=\varvec{e}_{\theta , i}\frac{I^{\text{cf}}}{4\pi R} \cot \frac{|\Delta \theta |}{2}, \end{aligned}$$

(1)

and the other comprises divergence-free (DF) basis functions:

$$\begin{aligned} \varvec{v}^{\text{df}}(\varvec{r}, \varvec{r}_i)=\varvec{e}_{\phi , i}\frac{I^{\text{df}}}{4\pi R} \cot \frac{|\Delta \theta |}{2}, \end{aligned}$$

(2)

where $\varvec{r}$ and $\varvec{r}_i$ are vectors pointing to positions on a sphere with radius R (i.e., $|\varvec{r}|=|\varvec{r}_i|=R$), and $\Delta \theta$ is the angle between the vectors $\varvec{r}$ and $\varvec{r}_i$, that is,

$$\begin{aligned} \Delta \theta = \arccos \frac{\varvec{r}\cdot \varvec{r}_i}{R^2}, \end{aligned}$$

(3)

and the vectors $\varvec{e}_{\phi }$ and $\varvec{e}_{\theta }$ are defined as follows:

$$\begin{aligned} \varvec{e}_{\phi , i}= & {} \frac{\varvec{r}_i\times \varvec{r}}{R^2\sin (\Delta \theta )}, \end{aligned}$$

(4)

$$\begin{aligned} \varvec{e}_{\theta , i}= & {} \frac{\varvec{e}_{\phi , i}\times \varvec{r}}{R}. \end{aligned}$$

(5)

The vector $\varvec{e}_{\phi }$ and $\varvec{e}_{\theta }$ correspond to unit vectors in the azimuthal and polar-angle directions of the spherical coordinate system where the pole is in the direction of the vector $\varvec{r}_i$. We hereinafter omit the constant factors in Eqs. (1) and (2) because they would not affect the following discussions, that is,

$$\begin{aligned} \varvec{v}^{\text{cf}}(\varvec{r}, \varvec{r}_i)= & {} \varvec{e}_{\theta , i} \cot \frac{|\Delta \theta |}{2}, \end{aligned}$$

(6)

$$\begin{aligned} \varvec{v}^{\text{df}}(\varvec{r}, \varvec{r}_i)= & {} \varvec{e}_{\phi , i} \cot \frac{|\Delta \theta |}{2}. \end{aligned}$$

(7)

Amm and Viljanen (1999) proposed an expansion of a vector field on a sphere in a series of SECS basis functions:

$$\begin{aligned} \varvec{V}(\varvec{r})=\sum _{i=1}^n \left[ w_i^{\text{cf}}\varvec{v}^{\text{cf}}(\varvec{r}, \varvec{r}_{i}) + w_i^{\text{df}}\varvec{v}^{\text{df}}(\varvec{r}, \varvec{r}_{i})\right] , \end{aligned}$$

(8)

where i is an index of a basis function, n is the number of the basis functions used in this expansion, and $\varvec{r}_i$ is the position of the center of the ith basis function. If basis functions located at various positions are distributed over the sphere, then various spherical vector fields can be represented by tuning the weights $\{w_i^{cf}\}$ and $\{w_i^{\text{df}}\}$. Equation (8) thus provides a flexible representation of a vector field on a sphere. However, this model is problematic when it is used for practical data analysis. The SECS basis functions $\varvec{v}^{\text{cf}}(\varvec{r}, \varvec{r}_i)$ and $\varvec{v}^{\text{df}}(\varvec{r}, \varvec{r}_i)$ diverge to infinity when approaching $\varvec{r}_i$. Therefore, $\varvec{V}(\varvec{r})$ in Eq. (8) would diverge to $\infty$ or $-\infty$ at any $\varvec{r}_{i}$. This singularity would deteriorate the robustness of the analysis. The prediction using Eq. (8) would thus be unreliable in the vicinity of $\varvec{r}_{i}.$

Extension of spherical elementary current systems

A vector field on a sphere can be decomposed into CF and DF components. The basic idea of the SECS expansion is to represent each of the vector field components by a linear combination of CF and DF localized basis functions. The CF SECS basis function was chosen so that its divergence is constant on the sphere except at $\varvec{r}=\varvec{r}_i$, and the DF SECS basis function was chosen so that its curl is constant except at $\varvec{r}=\varvec{r}_i$. As described above, however, the SECS basis functions diverge to infinity at $\varvec{r}=\varvec{r}_i$. This section explores a set of localized basis functions which satisfy the CF or DF condition and take only finite values over the sphere.

Before obtaining alternative set of basis functions, we consider scalar functions $\Psi ^{\text{cf}}(\varvec{r} )$ and $\Psi ^{\text{df}}(\varvec{r} )$ defined on a sphere. From these scalar functions, the following two kinds of spherical vector fields can be derived:

$$\begin{aligned} \varvec{V}^{\text{cf}} (\varvec{r})= & {} -\nabla _\Omega \,\Psi ^{\text{cf}}, \end{aligned}$$

(9)

$$\begin{aligned} \varvec{V}^{\text{df}} (\varvec{r})= & {} -\varvec{e}_r\times \nabla _\Omega \,\Psi ^{\text{df}}, \end{aligned}$$

(10)

where $\nabla _\Omega \,$ denotes the gradient operator on the sphere and $\varvec{e}_r$ denotes the radial unit vector at $\varvec{r}$. The vector field $\varvec{V}^{\text{cf}}$ satisfies the CF condition $\nabla _\Omega \,\times \varvec{V}^{\text{cf}}=0$, where $\Psi ^{cf}$ gives its scalar potential. The vector field $\varvec{V}^{\text{df}}$ satisfies the DF condition $\nabla _\Omega \,\cdot \varvec{V}^{\text{df}}=0$, where $\Psi ^{\text{df}}$ gives its stream function. We approximate each scalar function by a series of radial basis functions as follows:

$$\begin{aligned} \Psi ^{\text{cf}}(\varvec{r})= & {} \sum _{i=1}^n w_i^{\text{cf}}\psi ^{\text{cf}}(\varvec{r}, \varvec{r}_i) , \end{aligned}$$

(11)

$$\begin{aligned} \Psi ^{df}(\varvec{r})= & {} \sum _{i=1}^n w_i^{\text{df}}\psi ^{\text{df}}(\varvec{r}, \varvec{r}_i). \end{aligned}$$

(12)

Each of these equations can be regarded as a radial basis function network (e.g., Bishop, 2006). The radial basis functions $\psi ^{\text{cf}}$ and $\psi ^{\text{df}}$ yield the following localized vector-valued basis functions satisfying CF and DF conditions, respectively:

$$\begin{aligned} \varvec{v}^{\text{cf}}(\varvec{r}, \varvec{r}_i)= & {} -\nabla _\Omega \,\psi ^{\text{cf}}(\varvec{r}, \varvec{r}_i), \end{aligned}$$

(13)

$$\begin{aligned} \varvec{v}^{\text{df}}(\varvec{r}, \varvec{r}_i)= & {} -\varvec{e}_r\times \nabla _\Omega \,\psi ^{\text{df}}(\varvec{r}, \varvec{r}_i). \end{aligned}$$

(14)

Equtions (9) and (10) can thus be rewritten in the forms

$$\begin{aligned} \varvec{V}^{\text{cf}}(\varvec{r})= & {} \sum _{i=1}^n w_i^{\text{cf}}\varvec{v}^{\text{cf}}(\varvec{r}, \varvec{r}_i), \end{aligned}$$

(15)

$$\begin{aligned} \varvec{V}^{\text{df}}(\varvec{r})= & {} \sum _{i=1}^n w_i^{\text{df}}\varvec{v}^{\text{df}}(\varvec{r}, \varvec{r}_i), \end{aligned}$$

(16)

and thus a vector field $\varvec{V}(\varvec{r})$ can be decomposed as follows:

$$\begin{aligned} \varvec{V}(\varvec{r})=\sum _{i=1}^n w_i^{\text{cf}}\varvec{v}^{\text{cf}}(\varvec{r}, \varvec{r}_i) + \sum _{i=1}^n w_i^{\text{df}}\varvec{v}^{\text{df}}(\varvec{r}, \varvec{r}_i). \end{aligned}$$

(17)

Equtions (13) and (14) provide sufficient conditions of CF and DF basis functions, respectively. If we choose

$$\begin{aligned} \psi ^{cf}= \psi ^{\text{df}}= 2\log \left| \sin \frac{\Delta \theta }{2}\right| , \end{aligned}$$

(18)

then we obtain the SECS CF and DF basis functions in Eqs. (6) and (7).

A different choice of $\psi ^{\text{cf}}$ and $\psi ^{\text{df}}$ provides a different set of basis functions. One practical choice would be a spherical Gaussian function as follows:

$$\begin{aligned} \psi ^{\text{cf}}(\varvec{r}, \varvec{r}_i)= \psi ^{\text{df}}(\varvec{r}, \varvec{r}_i) = \exp \left[ \eta \left( \cos \Delta \theta -1\right) \right] = \exp \left[ \eta \left( \frac{\varvec{r}\cdot \varvec{r}_i}{R^2}-1\right) \right] , \end{aligned}$$

(19)

where $\eta$ is a scale parameter which gives the spread of the spherical Gaussian function. As $\eta$ is smaller, the spread becomes larger. When Eq. (19) is used, Eqs. (11) and (12) are regarded as Gaussian radial basis function networks bounded on a sphere for representing $\Psi ^{\text{cf}}$ and $\Psi ^{\text{df}}$. In principle, the Gaussian radial basis function network can represent any scalar function in a finite region with arbitrary accuracy by tuning the scale parameter and the number of basis functions (Park and Sandberg 1991, 1993). According to Eqs. (13) and (14), the spherical Gaussian function in Eq. (19) yields the following CF and DF basis functions:

$$\begin{aligned} \varvec{v}^{\text{cf}}(\varvec{r}, \varvec{r}_i)= & {} \varvec{e}_{\theta }\left( \eta \sin \Delta \theta \right) \exp \left[ \eta \left( \cos \Delta \theta -1\right) \right] = \varvec{e}_{\theta }\left( \eta \sin \Delta \theta \right) \exp \left[ \eta \left( \frac{\varvec{r}\cdot \varvec{r}_i}{R^2}-1\right) \right] , \end{aligned}$$

(20)

$$\begin{aligned} \varvec{v}^{\text{df}}(\varvec{r}, \varvec{r}_i)= & {} \left( \eta \varvec{r}_i\times \varvec{r}\right) \exp \left[ \eta \left( \cos \Delta \theta -1\right) \right] = \left( \eta \varvec{r}_i\times \varvec{r}\right) \exp \left[ \eta \left( \frac{\varvec{r}\cdot \varvec{r}_i}{R^2}-1\right) \right] . \end{aligned}$$

(21)

Figure 1 shows the shape of basis stream functions and the profile of the azimuthal component of DF basis functions. In the left panel of Fig. 1, the dotted line indicates the profile of a basis stream function for an original SECS function in Eq. (18), and the solid line indicates that for a proposed basis function in Eq. (19), where $\eta = 131.4$. In the right panel, the dotted line indicates the profile of the azimuthal component of a DF SECS basis function in Eq. (7) and the solid line indicates a proposed DF basis function in Eq. (21). While the SECS functions diverge as $\Delta \theta$ is approaching 0, the proposed basis functions take finite values at any point. Accordingly, an expansion in the form of Eq. (17) would take finite values at any point on the sphere as far as all the weights $\{w_i^{\text{cf}}\}$ and $\{w_i^{\text{df}}\}$ are finite.

Estimation of ionospheric plasma velocity field

As an application of the newly introduced basis functions in Eqs. (20) and (21), we estimate a temporal evolution of the ionospheric plasma velocity field from the SuperDARN data. The plasma velocity field in the ionosphere can be reasonably assumed to have no divergence. We can thus consider a stream function of the velocity field at the time $t_k$ in the form of Eq. (12):

$$\begin{aligned} \Psi _k(\varvec{r})=\sum _{i=1}^n w_{k,i}\psi (\varvec{r}, \varvec{r}_i), \end{aligned}$$

(22)

where a spherical Gaussian function in Eq. (19) is chosen as the stream function $\psi$ here. Since we consider a DF situation where the CF component can be ignored, the superscript df is omitted. The velocity field at $t_k$ can be derived from the stream function in Eq. (22) as follows:

$$\begin{aligned} \varvec{V}_k(\varvec{r})=\sum _{i=1}^n w_{k,i}\varvec{v}(\varvec{r}, \varvec{r}_i), \end{aligned}$$

(23)

where $\varvec{v}(\varvec{r}, \varvec{r}_i)$ is the DF basis function in Eq. (21). We placed 2500 basis functions over the region of the northern hemisphere at higher geomagnetic latitudes than 40°. The poles of the basis functions $\{\varvec{r}_i\}$ are distributed uniformly but at random locations over this region. Since the basis functions are distributed fairly uniformly, Eq. (23) corresponds to a Monte Carlo approximation of the following convolution integral (Heaton et al. 2014)

$$\begin{aligned} \varvec{V}_k(\varvec{r}) =\int\limits_S w(\varvec{r}') \varvec{v}(\varvec{r}, \varvec{r}')d\varvec{r}', \end{aligned}$$

(24)

where S indicates the region above $40^{\circ }$ in geomagnetic latitude. A spherical Gaussian function has a scale parameter $\eta$ which gives its spread. As $\eta$ is smaller, the spread becomes larger. If the spread is taken to be too small, each localized basis function could not fill the gaps with its neighbors and erroneous artifacts would be produced when the expansion of Eq. (23). Meanwhile, if the spread is too small, the spatial resolution of Eq. (23) becomes poor. Here we set $\eta$ at 131.4, which is determined so that the exponent of Eq. (19) becomes half of the peak value at $\Delta \theta =5^{\circ }$. This means the typical spread of a localized basis function in Eqs. (20) and (21) is $5\pi R/180$ in geodesic distance on the sphere. This spatial spread is large enough to enable us to represent a spatially smooth solution with 2500 basis functions distributed in the region of interest.

We want to obtain such weights $\{w_{k,1}, \ldots \}$ that $\varvec{V}_k(\varvec{r})$ will fit to the observations from the SuperDARN. We denote the jth record of the line-of-sight plasma velocity as $y_{k,j}$ and the location of that record as $\varvec{r}_{k,j}$. A prediction for $y_{k,j}$ based on Eq. (23) is given as

$$\begin{aligned} y_{k,j}^{\text{pred}}=\varvec{e}_{k,j}\cdot \varvec{V}_k({\varvec{r}_{k,j}})= \sum _{i=1}^n w_{k,i}\varvec{e}_{k,j}\cdot \varvec{v}({\varvec{r}_{k,j}}, \varvec{r}_i), \end{aligned}$$

(25)

where $\varvec{e}_{k,j}$ denotes a unit vector in the line-of-sight direction of a radar beam providing a velocity measurement at the point $\varvec{r}_{k,j}$. Denoting the observation noise as $\varepsilon _{k,j}$, the observed line-of-sight velocity $y_{k,j}$ can be written as

$$\begin{aligned} y_{k,j}=\sum _{i=1}^n w_{k,i}\varvec{e}_{k,j}\cdot \varvec{v}(\varvec{r}_{k,j}, \varvec{r}_i) +\varepsilon _{k,j}. \end{aligned}$$

(26)

Combining all the records of line-of-sight plasma drift velocity at the time $t_k$ into one vector $\varvec{y}_k$, Eq. (26) for all the records can be written in the following form:

$$\begin{aligned} \varvec{y}_k=\mathsf {H}_k\varvec{w}_k+\varvec{\varepsilon }_k, \end{aligned}$$

(27)

where $\varvec{w}_k$ is a vector consisting of all the weights contained in Eq. (23), $\varvec{\varepsilon }_k$ is a vector consisting of $\varepsilon _{k,j}$ for all $y_{k,j}$, and $\mathsf {H}_k$ is a matrix of which each element $H_{k,ji}$ is defined as

$$\begin{aligned} H_{k,ji}=\varvec{e}_{k,j}\cdot \varvec{v}(\varvec{r}_{k,j}, \varvec{r}_i) . \end{aligned}$$

(28)

If the estimate of $\varvec{w}_k$ is obtained, then we can basically predict the plasma drift velocity at an arbitrary point $\varvec{r}$ by using Eq. (23). However, there are some wide gaps in the spatial coverage of the SuperDARN, and the radars could miss velocity measurements even in the fields of view. It is difficult to appropriately determine the weights for the entire region in the hemisphere. In order to resolve the problems of observation gaps and missing data, we use an empirical model of the ionospheric electric potential as additional information. We decompose the weight $\varvec{w}_{k}$ into the background component $\varvec{\zeta }_{k}$ and the disturbance component $\varvec{\beta }_k$ as follows:

$$\begin{aligned} \varvec{w}_{k}=\varvec{\zeta }_{k}+\varvec{\beta }_{k}. \end{aligned}$$

(29)

The background component $\varvec{\zeta }_k$ can be obtained by assuming that the plasma drift velocity is equal to the $\varvec{E}\times \varvec{B}$ velocity which is given by the empirical model of the electric potential distribution. In this study, the electric potential model by Weimer (2001) is used to obtain $\varvec{\zeta }_k$. The OMNI hourly solar wind data were used as the input of the Weimer model. The background potential distribution at each minute is then obtained by linear interpolation of the hourly distribution of the Weimer model.

If $\varvec{\zeta }_k$ is given, Eq. (27) can be rewritten as

$$\begin{aligned} \varvec{y}_k=\mathsf {H}_k\varvec{\zeta }_k+\mathsf {H}_k\varvec{\beta }_k+\varvec{\varepsilon }_k. \end{aligned}$$

(30)

We then estimate the disturbance component $\varvec{\beta }_k$ by using a Bayesian method, which is based on the posterior distribution given the observations $\varvec{y}_k$. The posterior distribution $p(\varvec{w}_k|\varvec{y}_k)$ is obtained according to Bayes’ theorem:

$$\begin{aligned} p(\varvec{\beta }_k|\varvec{y}_k)=\frac{p(\varvec{y}_k|\varvec{\beta }_k)\, p(\varvec{\beta }_k)}{p(\varvec{y}_k)} . \end{aligned}$$

(31)

In the following, a Gaussian distribution with mean $\varvec{\mu }$ and covariance matrix $\mathsf {P}$ is denoted as ${\mathcal {N}}(\varvec{\mu }, \mathsf {P})$. If $\varvec{\varepsilon }_k$ in Eq. (30) obeys a Gaussian distribution ${\mathcal {N}}(\varvec{0}, \mathsf {R}_k)$, $p(\varvec{y}_k|\varvec{\beta }_k)$ becomes the Gaussian distribution with mean $\mathsf {H}_k\varvec{\zeta }_k+\mathsf {H}_k\varvec{\beta }_k$ and variance $\mathsf {R}_k$, that is

$$\begin{aligned} p(\varvec{y}_k|\varvec{\beta }_k)= {\mathcal {N}}(\mathsf {H}_k\varvec{\zeta }_k+\mathsf {H}_k\varvec{\beta }_k , \mathsf {R}_k). \end{aligned}$$

(32)

In addition, we assume that the prior distribution $p(\varvec{\beta }_k)$ is Gaussian as follows:

$$\begin{aligned} p(\varvec{\beta }_k)={\mathcal {N}}(\varvec{\beta }_{b,k}, \mathsf {P}_{b,k}). \end{aligned}$$

(33)

The posterior distribution $p(\varvec{\beta }_k|\varvec{y}_k)$ then becomes a Gaussian distribution, where its mean vector $\bar{\varvec{\beta }}_{k}$ and covariance matrix $\bar{\mathsf {P}}_k$ are

$$\begin{aligned} \bar{\varvec{\beta }}_{k}= & {} \varvec{\beta }_{b,k}+(\mathsf {P}_{b,k}^{-1}+\mathsf {H}_k^T\mathsf {R}_k^{-1}\mathsf {H}_k)^{-1} \mathsf {H}_k^T\mathsf {R}_k^{-1} (\varvec{y}_k-\mathsf {H}_k\varvec{\zeta }_{k}-\mathsf {H}_k\varvec{\beta }_{b,k}) , \end{aligned}$$

(34)

$$\begin{aligned} \bar{\mathsf {P}}_{k}= & {} (\mathsf {P}_{b,k}^{-1}+\mathsf {H}_k^T\mathsf {R}_k^{-1}\mathsf {H}_k)^{-1}. \end{aligned}$$

(35)

The mean of the posterior in Eq. (34) can be regarded as an estimate of the disturbance, and it provides an estimate of the weight $\bar{\varvec{w}}_k$ as follows:

$$\begin{aligned} \bar{\varvec{w}}_k=\varvec{\zeta }+\bar{\varvec{\beta }}_k. \end{aligned}$$

(36)

Once $\bar{\varvec{w}}_k$ is obtained, the stream function can be estimated by Eq. (22), and the plasma velocity field $\varvec{V}_k(\varvec{r})$ can be estimated by Eq. (23). The covariance matrix of the posterior in Eq. (35) represents the uncertainty of the estimate of $\bar{\varvec{\beta }}_k$. The uncertainty of $\Psi _k(\varvec{r})$ can be evaluated by introducing the following vector:

$$\begin{aligned} \varvec{a}(\varvec{r})=\left( \begin{array}{c} \psi (\varvec{r}, \varvec{r}_1) \\ \psi (\varvec{r}, \varvec{r}_2) \\ \vdots \\ \psi (\varvec{r}, \varvec{r}_n) \end{array} \right). \end{aligned}$$

(37)

Using the vector $\varvec{a}$, the variance of $\Psi _k(\varvec{r})$ becomes

$$\begin{aligned} \sigma _{\Psi _k(\varvec{r})}^2=\varvec{a}(\varvec{r})^T\bar{\mathsf {P}}_k\varvec{a}(\varvec{r}). \end{aligned}$$

(38)

This can be used as a proxy of the uncertainty of $\Psi _k(\varvec{r})$. The uncertainty of each component of the plasma velocity $\varvec{V}_k(\varvec{r})$ can also be evaluated in a similar manner.

The prior distribution of the disturbance $p(\varvec{\beta }_{k})$ can be determined in various ways. One effective way is to determine $p(\varvec{\beta }_{k})$ based on the observations in the previous time intervals. This means that the conditional distribution $p(\varvec{\beta }_{k}|\varvec{y}_{1:k-1})$ is used as a prior distribution, where $\varvec{y}_{1:k-1}$ denotes the sequence of the observations $\{\varvec{y}_1, \varvec{y}_2, \ldots , \varvec{y}_{k-1}\}$. We then obtain

$$\begin{aligned} p(\varvec{\beta }_k|\varvec{y}_{1:k}) =\frac{p(\varvec{y}_k|\varvec{\beta }_k)\, p(\varvec{\beta }_k|\varvec{y}_{1:k-1})}{p(\varvec{y}_k|\varvec{y}_{1:k-1})} . \end{aligned}$$

(39)

The conditional distribution $p(\varvec{\beta }_{k}|\varvec{y}_{1:k-1})$ can be obtained from $p(\varvec{\beta }_{k-1}|\varvec{y}_{1:k-1})$, which corresponds to the distribution of Eq. (39) for the previous time interval. We assume the distribution $p(\varvec{\beta }_{k-1}|\varvec{y}_{1:k-1})$ is Gaussian as follows:

$$\begin{aligned} p(\varvec{\beta }_{k-1}|\varvec{y}_{1:k-1})= {\mathcal {N}}(\varvec{\beta }_{k-1|k-1}, \mathsf {P}_{k-1|k-1}) . \end{aligned}$$

(40)

In addition, the probability distribution of $\varvec{\beta }_{k}$ given the disturbance for the previous time interval, $\varvec{\beta }_{k-1}$, is assumed to be a Gaussian distribution

$$\begin{aligned} p(\varvec{\beta }_k|\varvec{\beta }_{k-1})= {\mathcal {N}}(\alpha \varvec{\beta }_{k-1}, \mathsf {Q}) . \end{aligned}$$

(41)

We then obtain a predictive distribution of $\varvec{\beta }_k$ as a Gaussian distribution:

$$\begin{aligned} p(\varvec{\beta }_{k}|\varvec{y}_{1:k-1})= {\mathcal {N}}(\varvec{\beta }_{k|k-1}, \mathsf {P}_{k|k-1}) , \end{aligned}$$

(42)

where the mean vector $\varvec{\beta }_{k|k-1}$ and covariance matrix $\mathsf {P}_{k|k-1}$ are given as

$$\begin{aligned} \varvec{\beta }_{k|k-1}= & {} \alpha \varvec{\beta }_{k-1|k-1}, \end{aligned}$$

(43)

$$\begin{aligned} \mathsf {P}_{k|k-1}= & {} \alpha ^2\mathsf {P}_{k-1|k-1}+\mathsf {Q}. \end{aligned}$$

(44)

Similarly to Eqs. (34) and (35), the posterior mean $\varvec{\beta }_{k|k}$ becomes

$$\begin{aligned} \varvec{\beta }_{k|k}= \varvec{\beta }_{k|k-1}+(\mathsf {P}_{k|k-1}^{-1}+\mathsf {H}_k^T\mathsf {R}_k^{-1}\mathsf {H}_k)^{-1} \mathsf {H}_k^T\mathsf {R}_k^{-1} (\varvec{y}_k-\mathsf {H}_k\varvec{\zeta }_{k}-\mathsf {H}_k\varvec{\beta }_{k|k-1}) , \end{aligned}$$

(45)

and the posterior covariance matrix $\mathsf {P}_{k|k}$ is

$$\begin{aligned} \mathsf {P}_{k|k}= (\mathsf {P}_{k|k-1}^{-1}+\mathsf {H}_k^T\mathsf {R}_k^{-1}\mathsf {H}_k)^{-1} . \end{aligned}$$

(46)

Applying Eqs. (43)–(46) recursively, $\varvec{\beta }_k$ for each time interval can be estimated. This is the Kalman filter algorithm for estimating $\varvec{\beta }_k$ based on the sequence of observations $\varvec{y}_1, \ldots , \varvec{y}_k$. The estimate of $\varvec{\beta }_k$ provides the estimate of $\varvec{w}_k$ according to Eq. (29). The posterior distribution of $\varvec{w}_k$ is thus given as

$$\begin{aligned} p(\varvec{w}_k|\varvec{y}_{1:k})={\mathcal {N}}(\varvec{w}_{k|k}, \mathsf {P}_{k|k}) , \end{aligned}$$

(47)

where $\varvec{w}_{k|k}=\varvec{\zeta }_k+\varvec{\beta }_{k|k}$. The stream function is estimated according to Eq. (22):

$$\begin{aligned} \varvec{\Psi }_{k|k}(\varvec{r})=\sum _{i=1}^n w_{k|k,i}\psi (\varvec{r}, \varvec{r}_i), \end{aligned}$$

(48)

where $w_{k|k,i}$ denotes the ith element of the vector $\varvec{w}_{k|k}$. The plasma drift velocity distribution is also estimated according to Eq. (23):

$$\begin{aligned} \varvec{V}_{k|k}(\varvec{r})=\sum _{i=1}^n w_{k|k,i}\varvec{v}(\varvec{r}, \varvec{r}_i). \end{aligned}$$

(49)

As done in Eq. (38), the variance of the $\Psi _{k}$ given $\varvec{y}_{1:k}$ can be obtained as follows:

$$\begin{aligned} \sigma _{\Psi _{k|k}(\varvec{r})}^2=\varvec{a}(\varvec{r})^T\bar{\mathsf {P}}_{k|k}\varvec{a}(\varvec{r}), \end{aligned}$$

(50)

and the variance of $\varvec{V}_k$ given $\varvec{y}_{1:k}$ can also be obtained using a similar equation.

In the beginning of the Kalman filter algorithm, $\varvec{\beta }_{0|0}$ and $\mathsf {P}_{0|0}$ were set as follows:

$$\begin{aligned} \varvec{\beta }_{0|0}&=\varvec{0},&\mathsf {P}_{0|0}&=400\mathsf {Q}. \end{aligned}$$

(51)

Although the initial variance $\mathsf {P}_{0|0}$ is very large, it becomes a reasonable value after several steps. In order to estimate a sequence of $\varvec{\beta }_k$, the parameter $\alpha$ in Eq. (41) and the covariance matrices $\mathsf {Q}$ and $\mathsf {R}_k$ must also be given in advance. In order to ensure spatial smoothness among the elements of $\varvec{\beta }_k$, the matrix $\mathsf {Q}$ is set as follows:

$$\begin{aligned} \mathsf {Q}=\sigma _Q^2 \left( \begin{array}{cccc} C_Q(\varvec{r}_1, \varvec{r}_1) &{} C_Q(\varvec{r}_1, \varvec{r}_2) &{} \cdots &{} C_Q(\varvec{r}_1, \varvec{r}_n) \\ C_Q(\varvec{r}_2, \varvec{r}_1) &{} C_Q(\varvec{r}_2, \varvec{r}_2) &{} &{} C_Q(\varvec{r}_2, \varvec{r}_n) \\ \vdots &{} &{} \ddots &{} \vdots \\ C_Q(\varvec{r}_n, \varvec{r}_1) &{} C_Q(\varvec{r}_n, \varvec{r}_2) &{} \cdots &{} C_Q(\varvec{r}_n, \varvec{r}_n) \end{array} \right), \end{aligned}$$

(52)

where $C_Q(\varvec{r}_i, \varvec{r}_j)$ is a correlation function between $\varvec{r}_i$ and $\varvec{r}_j$. We assume $C_Q(\varvec{r}_i, \varvec{r}_j)$ is a spherical Gaussian function:

$$\begin{aligned} C_Q(\varvec{r}_i,\varvec{r}_j)= \rho (\varvec{r}_i,\varvec{r}_j) \exp \left[ \kappa \left( \frac{\varvec{r}_i\cdot \varvec{r}_j}{R^2}-1\right) \right] , \end{aligned}$$

(53)

where we assume $\kappa =14.7$ so that the exponent becomes half of the peak at $\Delta \theta =15^{\circ }$. The function $\rho$ in Eq. (53) forces the stream function $\Psi _k(\varvec{r})$ at the boundary ($40^{\circ }$ geomagnetic latitude) to be near zero. We define $\rho$ as follows:

$$\begin{aligned} \rho (\varvec{r}_i,\varvec{r}_j) = \frac{(\sin \lambda _i-\sin 40^{\circ })(\sin \lambda _j-\sin 40^{\circ })}{(1-\sin 40^{\circ })^2}. \end{aligned}$$

(54)

The matrix $\mathsf {R}_k$ is set as follows:

$$\begin{aligned} \mathsf {R}_k=\sigma _R^2 \left( \begin{array}{cccc} \delta _{b_1b_1}C_R(g_1, g_1)&{} \delta _{b_1b_2}C_R(g_1, g_2) &{} \cdots &{} \delta _{b_1b_m}C_R(g_1, g_m) \\ \delta _{b_2b_1}C_R(g_2, g_1)&{} \delta _{b_2b_2}C_R(g_2, g_2) &{} &{} \delta _{b_2b_m}C_R(g_2, g_m) \\ \vdots &{} &{} \ddots &{} \vdots \\ \delta _{b_mb_1}C_R(g_m, g_1)&{} \delta _{b_mb_2}C_R(g_m, g_2) &{} \cdots &{} \delta _{b_mb_m}C_R(g_m, g_m) \end{array} \right), \end{aligned}$$

(55)

where m is the dimension of $\varvec{y}_k$, $\delta _{ij}$ is the Kronecker delta:

$$\begin{aligned} \delta _{ij}= {\left\{ \begin{array}{ll} 1 \quad (i=j), \\ 0 \quad (i\ne j), \end{array}\right. } \end{aligned}$$

(56)

$b_i$ and $g_i$ denote the beam number and the range gate of the ith element of the observation $\varvec{y}$. The beam number indicates the direction of a radar beam and the range gate is associated with the distance from a radar site. We assume $C_R(g_i, g_j)$ can be written in the following form:

$$\begin{aligned} C_R(g_i,g_j)=\exp \left[ -\frac{(g_i-g_j)^2}{2}\right] . \end{aligned}$$

(57)

Equation (55) means the correlations within the same beam are considered while the correlations between different beams are ignored. The parameters $\sigma _Q$, $\sigma _R$, and $\alpha$ are determined by maximizing the marginal likelihood:

$$\begin{aligned} p(\varvec{y}_{1:K}|\sigma _Q, \sigma _R, \alpha )=\prod _k\int p(\varvec{y}_k|\varvec{\beta }_k)\, p(\varvec{\beta }_k|\sigma _Q, \sigma _R, \alpha )\, d\varvec{\beta }_{k} \end{aligned}$$

(58)

(e.g., Morris 1983; Casella 1985). As a result, we chose $\sigma _Q=100$, $\sigma _R=400$, and $\alpha =0.9$.

In order to evaluate how the proposed basis functions works, we performed fitting to the drift velocity distribution obtained from an empirical model by Weimer (2001). The red line in Fig. 2 shows the latitudinal profile of the azimuthal velocity at 3h MLT which was reconstructed by the proposed functions. For reference, the reconstruction with the original SECS functions is overplotted with the blue line, and the drift velocity profile directly derived from the Weimer’s model is overplotted with the green line. In this fitting, we did not consider the temporal evolution; that is, the estimation was performed with Eq. (34), where $\mathsf {P}_{b,k}$ was taken to be as large as $\mathsf {P}_{0|0}$ in Eq. (51). While the Weimer’s model gives a smooth latitudinal profile, the original SECS functions produce a little serrated profile due to their singularities. This problem would be inconvenient when compared with other observations such as spacecraft observations. In comparison with high-resolution data such as total electron content (TEC) data from Global Navigation Satellite System (GNSS) satellites, it would be almost impossible to expect all the data points are apart from a center of any basis functions. The singularities could be eliminated by modifying the original SECS functions in the neighborhood of their source points as discussed by Vanhamäki and Juusola (2020). The proposed basis functions provide another practical way which enables us to obtain a reasonable estimate as smooth as the Weimer’s model.

Application

We estimated the temporal evolution of the ionospheric plasma drift velocity distribution during the magnetic storm on March 17, 2015 by using the method described in the previous section. In advance of applying the above method, the records of line-of-sight velocity of which the absolute values were less than $100\,\mathrm {m/s}$ were excluded because they might be ground scatter echoes. The records of the absolute values that were larger than $2000\,\mathrm {m/s}$ were also excluded as outliers. We used the data from 16 radar sites (Adak Island East, Adak Island West, Christmas Valley East, Christmas Valley West, Fort Hays East, Fort Hays West, Hankasalmi, Hokkaido East, Hokkaido West, Inuvik, Kapuskasing, King Salmon, Prince George, Pykkvibaer, Rankin Inlet, Saskatoon) in operation during this event.

Figure 3 shows the estimated plasma drift velocity distribution and some related quantities at 08:30 UT, when the data coverage on the nightside was relatively good. The upper-left panel indicates the stream function estimated according to Eq. (48), and the upper-right panel indicates the uncertainty of the stream function value as estimated by the variance of the posterior distribution of the stream function, which can be obtained using Eq. (50). The lower-left panel shows the estimated plasma drift velocity distribution with white arrows. The uncertainty of the plasma drift velocity is superposed with a color code in this panel. The lower-right panel shows the number of the available SuperDARN data for past 10 min for each bin, where white indicates no data points being available for that bin. As shown in the lower-right panel, the SuperDARN data points were not evenly distributed, but rather most of the data were obtained on the duskside or the nightside at this time. Accordingly, the uncertainty of the estimate tends to be small on the duskside and the nightside, while it tends to be large on the dawnside and the dayside as indicated in the upper-right panel and the lower-left panels. The uncertainty near the outer boundary is also small because we assume $\Psi _k(\varvec{r})$ at the boundary to be zero.

Figure 4 shows geomagnetic indices (the SYM-H and AE indices) and the southward component of the interplanetary magnetic field observed by the ACE spacecraft from 00:00 to 12:00 UT for reference. At 08:30 UT, a magnetic storm was in progress, and the empirical model by Weimer (2001) predicted that the two-cell convection structure was well developed, which is shown in the upper-left panel of Fig. 3. In order to determine the impact of the SuperDARN data, the drift velocity distribution is decomposed into two components: the background Weimer model component $\zeta _k$ and the disturbance component $\beta _k$. As described above, the disturbance $\beta _k$ is estimated so as to fit the SuperDARN data. Thus, the disturbance component reflects the impact of the SuperDARN data. The left and right panels in Fig. 5 show the background and the disturbance components, respectively. In the polar cap region, the disturbance component shows a flow toward the dayside, which is opposite to the normal two-cell convection as seen in the empirical model. This would suggest that the flow speed in the dayside polar cap region was overestimated by the empirical model and that it was estimated to be smaller by considering the SuperDARN data. On the nightside, the disturbance component was small, which means the empirical model predicted the flow distribution with high accuracy.

Fig. 1

The profile of the basis stream function with respect to $\Delta \theta$ for the original SECS and the proposed basis function (left) and the profile of the azimuthal component of the divergence-free SECS basis function and the proposed divergence-free basis function (right). The scale parameter $\eta$ for the proposed basis functions is set at 131.4. In each panel, the dotted line indicates the original SECS and the solid line indicates the proposed basis function

Fig. 2

The latitudinal profile of the azimuthal drift velocity at 3h MLT obtained by fitting to the Weimer’s empirical model. The red line indicates the reconstruction with the proposed basis functions, the blue line indicates the reconstruction with the original SECS functions, and the green line indicates the Weimer’s model

Fig. 3

The result of the estimation at 08:30 UT on March 17, 2015 for the northern hemisphere. The upper-left panel shows the estimated stream function, the upper-right panel shows the uncertainty of the stream function, the lower-left panel shows the estimated plasma drift velocity distribution (arrows) and its uncertainty (color code), and the lower-right panel shows the cumulative number of data points for each bin for the past 10 minutes. In each panel, the center is the north pole and the outer boundary is the line of $40^{\circ }$ geomagnetic latitude. The unit of the stream function is $R_I\cdot \mathrm {m}/\mathrm {s},$ where $R_I$ is the radius of the ionosphere

Fig. 4

Overview of geomagnetic conditions from 00:00 to 12:00 UT on March 17, 2015. The upper panel shows the SYM-H index, the middle panel shows the AU and AL indices, and the bottom panel shows the z component of the interplanetary magnetic field in GSM coordinates measured by the ACE spacecraft. The interplanetary magnetic field is plotted by assuming the time delay from ACE to the Earth was 45 min

Figure 6 shows the temporal evolution of the estimated plasma drift velocity distribution from 06:00 to 08:30 UT, and Fig. 7 shows the corresponding evolution of the disturbance component. The estimation result shows the ionospheric convection was weak at 06:00 UT. At 06:30 UT, the ionospheric convection was estimated to be enhanced, which was in accord with the southward turning of the solar wind magnetic field around 06:00 UT as shown in Fig. 4. The disturbance component at 06:30 UT was in the opposite sense to the estimated plasma velocity shown in Fig. 7 especially in the polar cap region, which suggests the empirical model overestimated the plasma velocity there. The disturbance component gradually decayed until 07:30 UT on the dayside. In the post-midnight region, the westward disturbance component is maintained from $50^{\circ }$ to $60^{\circ }$ geomagnetic latitude. In the pre-midnight region, the disturbance component is oriented northward. At 08:30 UT, the disturbance component in the sunlit polar cap was visible again. Figure 8 shows a histogram of the deviation of the observed line-of-sight velocity from the estimate obtained with the proposed method with a red line. The deviation of the observation from the prediction with the background Weimer’s model is also shown with blue shade. Large residuals still remain after the fitting probably because we assumed temporal smoothness as in Eq. (41) and spatial smoothness as in Eq. (52) to obtain robust results. However, the residuals certainly became smaller than the prediction with the Weimer’s model. The global observation such as the SuperDARN provides vital information for improving a prediction with an empirical model for each event. The proposed technique is helpful for making use of the global measurement in order to obtain a reliable prediction.

Discussion and concluding remarks

We have proposed a technique for estimating the spatial distribution of ionospheric plasma drift velocity. In this technique, the ionospheric plasma velocity distribution is represented by a sum of the background component and the perturbation component. While the background component is obtained from an empirical convection model, the perturbation component is written by a weighted sum of DF localized basis functions. In modeling of the global ionospheric drift velocity distribution, global basis functions such as spherical harmonic functions were usually used (e.g., Ruohoniemi and Baker 1998; Thomas and Shepherd 2018). However, such global basis functions may cause spurious artifacts in the data gap region when using the SuperDARN data. The use of localized basis functions would be an effective way to avoid such artifacts.

The basis functions introduced in this study can be regarded as a variant of the SECS proposed by Amm (1997). However, the original SECS basis functions have a singular point. This problem would be serious especially when analyzing measurement data with high spatial resolution such as the SuperDARN data. If high-resolution data are used, it would be difficult to place all singular points apart from any observation points. This means that some observation points would be in the vicinity of singular points. In such a situation, analysis with the SECS functions would be strongly constrained by the data obtained in the vicinity of singularities, and representation of large-scale features could be poor.

Fig. 5

The background Weimer model (left panel) and the disturbance component of the estimation (right panel) at the same time as Fig. 3

Fig. 6

The sequence of estimated plasma drift velocity distribution (arrows) and its uncertainty (color code) from 6:00 UT to 8:30 UT on March 17, 2015

Fig. 7

The sequence of the estimated perturbation component for plasma drift velocity distribution from 6:00 UT to 8:30 UT on March 17, 2015

Fig. 8

Histogram of the deviation of the observed line-of-sight velocity from the prediction with the proposed method (red line) and the deviation of the observation from the prediction with the background Weimer’s model (blue shade).

In contrast, the basis functions employed in this paper have been derived from spherical Gaussian functions, and they do not have a singular point. The proposed framework would enable us to obtain a natural smooth solution in analysis of high spatial resolution data. Moreover, it enables us to obtain a reliable estimate of the plasma drift velocity at an arbitrary point of the ionosphere. This allows us to use analysis results as a reference for comparison with other observations such as spacecraft observations. Nowadays a variety of ionospheric and geomagnetic measurements attain high spatial resolution. The proposed framework would also be helpful for comparing such data with high spatial resolution.

This study focused on the estimation of the global plasma velocity in the ionosphere. However, the basis functions introduced in this study can be applied for a variety of purposes as the original SECS functions are used for many studies. For example, since the proposed framework employs local basis functions, we can apply this framework for estimating a regional pattern of the plasma velocity distribution as done by Hori et al. (2018). Modeling of the ionospheric electric current would also be a potentially useful application. The original work by Amm (1997) proposed the SECS functions for representing ionospheric currents and analyzing ground magnetic disturbances due to the ionospheric currents. Many studies have employed the SECS functions for this purpose. Such analyses of ionospheric currents could be improved with the proposed smooth basis functions.