Double-thin-shell approach to deriving total electron content from GNSS signals and implications for ionospheric dynamics near the magnetic equator

A new technique was developed to estimate the ionospheric total electron content (TEC) from Global Navigation Satellite System (GNSS) satellite signals. The vertically distributed electron density was parameterized by two thin-shell layers (double-shell approach). The spatiotemporal variation of TEC (strictly speaking, partial electron content) associated with each shell was approximated by the functional fitting of spherical surface harmonics. The major improvements over the conventional single-shell approach were as follows: (1) the precise estimation of TEC was achieved; (2) the estimated TEC was less dependent on the choice of shell heights; and (3) the equatorial anomaly was captured more correctly. Furthermore, higher and lower shells exhibited a different pattern of local time vs latitude variation, providing information on the ionosphere–thermosphere dynamics.


Introduction
Ionosondes have been used as a standard technique of measuring ionospheric characteristics for more than seven decades. Over the last two decades, beacon signals transmitted from the global navigation satellite system (GNSS) satellites have been widely utilized for ionosphere studies. These two methods detect different parameters of the ionosphere. The primary parameter obtained from ionosondes is the electron density (NmF2) at the peak height, and the primary parameter measured by the satellite beacon signals is total electron content (TEC), in which electron density is integrated along the satellite-to-receiver path. Thus, they are often used complementarily for ionosphere and upper-atmosphere studies. NmF2 is directly calculated from the ionospheric critical frequency, foF2, while TEC derivation from beacon signals requires complicated procedure.
TEC is calculated from measurements of the differential propagation delay between two frequencies that are transmitted from GNSS satellites and propagated through the dispersive medium. A basic problem is that the signal delays in the electronics of satellites and receivers are different for the two frequencies, which is often referred to as an inter-frequency instrumental bias. Various algorithms have been proposed to remove instrumental biases using a single receiver (Lanyi and Roth 1988;Arikan et al. 2008;Zhang et al. 2009;Choi et al. 2013;Prasad et al. 2016), a regional or local network of receivers (Otsuka et al. 2002;Ma and Maruyama 2003;Ma et al. 2005;Choi et al. 2011), and a global receiver network (Mannucci et al. 1998;Hernández-Pajares et al. 1999). Regarding the global receiver network, several organizations provide TEC in the form of global ionosphere maps (GIMs) along with instrumental biases (see Hernández-Pajares et al. 2009). However, different approaches based on different assumptions and procedures yield different TECs and biases regardless of whether a single receiver or a network of receivers is used (Brunini et al. 2005;Hernández-Pajares et al. 2009). Thus, there is scope for the TEC derivation to improve.
The primary observable is the differential group delay ( �τ ) including instrumental biases. The ionospheric delay (or slant TEC) is a function of the inclination of the ray path, while the instrumental biases are independent of the inclination. Because of this feature, two variables can be separated. In the separation process, the slant TEC must be converted to the vertical value or vertical TEC at a reference point on the ray path by multiplying by the cosine of the local zenith angle ( χ ). For this purpose, the ionosphere is expressed by a thin shell and χ is represented by the value at the point, where the ray path intercepts the shell (ionosphere piece point (IPP)). Most techniques use the thinshell approximation with a constant shell height. However, there are some problems with the thin-shell approximation: (1) the ionospheric peak height and vertical electron density profile largely vary with the location and time; (2) a different shell height connects the observed slant TEC with a different geographical location of the IPP; and (3) the slant TECs differ at the same IPP but different ray path directions under horizontal gradient of TEC. Thus, there is no theoretical criterion for choosing the shell height. A detailed discussion of these problems is given by Hernández-Pajares et al. (1999).
The ionospheric group delay measured by the modulation code gives an absolute value of slant TEC except that instrumental biases are included. However, the data are subject to noises and the multipath effect. On the other hand, the carrier phase advance gives the same ionospheric effect as the group delay but with a different sign, and with less noises and a reduced multipath problem. Because the origin of each series of continuous phase data is unknown (phase ambiguity), the absolute value is obtained by leastsquares fitting to the code data, which is called code-phase leveling and commonly used to derive TEC precisely. Code-phase leveling, however, is not free of noises or the multipath effect on the code data (Hernández-Pajares et al. 1999;Ciraolo et al. 2007). Furthermore, receiver biases are not necessarily constant during a day, which is another source of the leveling error and, even for a combination of the same satellite and receiver, two different arcs are affected by different leveling errors (Ciraolo et al. 2007) In this paper, to resolve the problems inherent in the conventional techniques, (1) the three-dimensional ionosphere was parameterized by two thin shells and the spatiotemporal distribution of TEC (more strictly, partial electron content) associated with each shell was modeled by spherical harmonics; and (2) the arc biases were defined, which contain instrumental biases, the leveling error, and the cycle-slip correction error. When there were gaps in an arc, that arc was divided into two or more independent arcs. The algorithm is described in "Algorithm" section. The technique is first validated for simulated slant TEC data in "Validation by IRI simulation" section. The application to an actual dataset follows in "Application to actual observations". The discussion and interpretation of the results are given in "Ionospheric dynamics" section. Further application to a dataset under different geophysical conditions is shown in "Application under different geophysical conditions" section. The new technique and its advantage are summarized in "Summary". The leveling-free calculation in which the phase data were directly used without code-phase leveling is described in the appendix.

Algorithm
The technique was developed to apply to local receiver networks along a meridian, although it is applicable in single-receiver cases. To adopt the surface harmonics fitting approach, data obtained in a whole day were prepared and the local time of IPP or the angular local time, which covers the range 0-2π , was used as the longitudinal argument ( φ ) of the surface harmonic function. The colatitude argument ( θ ) was the co-angle of the modified dip latitude (modip), which is suitable for the analysis of the equatorial and lowlatitude ionosphere. The differential carrier phase advance ( −�τ ) between the two frequencies of GNSS satellites was calculated from RINEX data files. The consecutive phase dataset for a given receiver and satellite combination was numbered (arc number).

Conversion of measurements
The ionosphere effect on radio wave propagation was assumed to be the summation of the contributions in two thin-shell layers located at two constant heights (doubleshell approximation). The azimuth ( γ ) and elevation ( η ) of the satellite were converted to the IPP parameters, as shown in Fig. 1. The latitude ( i ) and longitude ( ζ i ) of the IPP in each shell and the local satellite zenith angle ( χ i ) ( i = Lo and Hi for the lower and higher shells, respectively) were calculated using the same equation set as that used in the single-thin-shell approximation.
(1) i = sin −1 (sin R cos α i + cos R sin α i cos γ ) (2) ζ i =ζ R + sin −1 sin α i sin γ cos i Here, R and ζ R are the latitude and longitude of the receiver, respectively, R E is the earth's radius, and h si is the shell height. The coordinates of the IPP are further converted to two arguments, colatitude ( θ i ) and angular local time ( φ i ), as follows: where modip ( i , ζ i ) was calculated by the International Geomagnetic Reference Field model (IGRF2012) (Thébault et al. 2015) and T is UT in seconds. In the application of (6), the ionosphere was assumed to be rigid in a reference frame rotating with the sun.

Surface harmonics fitting
Ionosphere pierce points in the thin shell are scattered over a wide range of latitude and longitude (and thus local time) around the ground receiver, and they move with time. Thus, vertical TECs associated with each shell, TEC * VLo and TEC * VHi , at different IPPs and times are represented by the surface harmonics of the colatitude ( θ i ) and local time in terms of angle ( φ i ) as follows: (3) where P m n (cos θ i ) is the associated Legendre polynomial with degree n and order m. The coefficients A nm , B nm , Ã nm , and B nm are unknown parameters, which should be determined from the observations. Once the harmonics fitting is achieved for each shell, the vertical TEC at any given subionospheric location ( θ , φ ) is obtained by adding two functional values.
The maximum degree and order were chosen as N = 9 and M = 7 , respectively, in all runs in this paper. These values were a compromise between the reproducibility of small-scale spatiotemporal variation of TEC and credibility of the results, which depends on the number of IPPs and noises. The available observations are the slant TECs, TEC obs sl , which are calculated from the differential phase advance. Here, the basic problem is The arc bias ( b arc ) is considered to be a constant during an arc of continuous phase data.

Orthogonal network
The surface harmonics fitting (or the determination of coefficients A nm , B nm , Ã nm , and B nm ) and the estimation of b arc were simultaneously carried out by an orthogonal network (Yang and Tseng 1996) with an additional network for bias estimation (Ma et al. 2005), as shown in Fig. 2. The IPP coordinates θ i and φ i are fed to the nodes of the orthogonal function layer (surface harmonics) and each node calculates cos mφ i P m n (cos θ i ) or sin mφ i P m n (cos θ i ) . For convenience of explanation, (7) and (8) are rewritten as where K = M(2N − M + 1) + N + 1 . In (10) and (11), w j and u j correspond to A nm , B nm , Ã nm , and B nm . Fig. 1 Geometry of double-shell approximation of ionosphere. The azimuth ( γ ) and elevation ( η ) angles of the satellite ray path and the differential delay ( �τ ) between two signals are measured. The ray path crosses shells at the height h si (Lo and Hi for the lower and higher shells, respectively). The point of intersection is referred to as the ionosphere pierce point (IPP). RX is the receiver The diurnal variation of TEC is roughly described by the combination of a half cosine part corresponding to daytime and a flat part corresponding to nighttime such as in the Klobuchar model (Klobuchar 1987). A harmonics fitting of local time variation requires higher-order terms to reproduce the transfer between day and night. Furthermore, under the solar minimum condition, the nighttime flat part is close to zero TECu and the harmonics fitting with a lower degree often yields a negative value in nighttime. Thus, for efficient functional fitting with the minimum number of harmonics and to avoid negative values of TEC, a nonlinear activation function σ (softplus function) is introduced in the output layer as follows: where i = Lo and Hi for the lower and higher shells, respectively. The graph of the function σ is shown in Fig. 3. When x is larger than a few TECu, TEC Vi ≃ TEC * Vi and the activation function has no effect on the output. When x is negative, the functional value remains positive and reaches zero at a large negative x value. In the (12) σ (x) = ln 1 + e x (13) TEC Vi =σ TEC * Vi Fig. 2 Orthogonal network for functional fitting. Each node of the orthogonal layer calculates a component of the surface harmonics. The weighted sum of the outputs of the orthogonal layer for each shell is further transferred by the function σ in the output layer. The sum of the two shells gives the vertical TEC. Weights and biases are adjusted to minimize the residual error double-shell parameterization of the ionosphere, only the sum of the contributions from the two shells is compared with the observations. Therefore, even if the sum is positive, TEC associated with one of the two shells can be negative. The softplus activation function in the output layer avoids an unrealistic local minimum in which one of the two shells becomes negative. Thus, (10)-(13) are used instead of (7) and (8).
Functional fitting was performed to minimize the squared residual error of the slant TEC below.
The weights and biases were updated by the gradient descent algorithm as described below and the update was repeated until the residual error converged.
The weights and biases at the step s are updated to those at the step s + 1 with the learning rates µ and μ as In the update process, it is important to choose appropriate learning rates to ensure stability. For a given n, the amplitude of P m n varies greatly with m, because it is not normalized. Instead, µ was scaled by the mean of the square of the spherical surface harmonics over the surface ( S nm ).
The update was carried out in two stages. First, weights and biases were updated after each pattern was presented in a randomized order (online-mode learning) (Haykin 1994). The residual error decreased rapidly in the first ∼ 100 repetitions and then gradually. The error had almost converged by approximately 1000 to 1500 repetitions. However, the error fluctuated with the update step depending on the order of pattern presentation, which was shuffled at each epoch of repetition. Further updates were carried out by batch-mode learning, i.e., averages of w j , u j , and b arc calculated for all patterns were used for updates (Haykin 1994). In all runs in this paper, µ 0 was 0.00001 and 0.001 in the online-and batch-mode learning, respectively, with μ = K µ.
A drawback of the harmonics fitting applied here is the discontinuity of TECs at the start and end of the 1-day period, which were forced to be the same in the fitting with the periodic function, while the ionospheric daily variations are not exactly periodic. To mitigate this problem, the daily dataset was chosen to start immediately before sunrise and to end at the same time on the next day, around which TEC generally takes a minimum and the day-to-day variability is the smallest. (17)

Validation by IRI simulation
When evaluating the vertical TEC estimated from slant TEC measurements, the fitting error is often used as a measure of the correctness of the results. However, it is not a suitable parameter when both biases and TECs are simultaneously determined, because the bias and TEC errors compensate for each other. Thus, to validate the developed technique, we needed to evaluate the bias and vertical TEC errors separately. However, this is not easy from actual observations. The validation was first performed for simulated slant TEC data, and the application to actual observations will be discussed in the next section. To demonstrate the advantage of the doubleshell approach, the same evaluation was conducted for the single-shell approach using a reduced version of the double-shell model. The single-shell model was realized by simply removing the higher-shell part in the diagram shown in Fig. 2.
The International Reference Ionosphere model (IRI2016) (Bilitza et al. 2017) was used to calculate slant TECs. For a realistic simulation, the slant paths were adopted from the actual observation geometry at intervals of 5 min. The observations were carried out using eight receivers composed of a network aligned along the 100 • E meridian near the magnetic equator in Southeast Asia, as summarized in Table 1. GPS data were selected for the period from 22 UT, 16 March to 22 UT, 17 March 2014, starting and ending immediately before sunrise. The total number of arcs was 564. The IPPs at 300 and 600 km are shown in Fig. 4. The electron density was integrated along the ray path up to a height of 2000 km, which is the upper boundary of the IRI model. Arc biases were randomly generated within ±25 TECu and added to the calculated slant TECs. Vertical TECs were directly calculated for use as a reference for evaluating the estimated vertical TECs. A map of the directly calculated TEC (IRI-direct TEC) as a function of local time and latitude is shown in Fig. 5a. The dash-dotted line indicates the magnetic equator. Along with the diurnal variation, the characteristics of the equatorial anomaly with a trough over the magnetic equator and crests away from the magnetic equator are depicted.

Single-shell results
In single-shell approaches, the shell height, where IPPs are located is usually fixed slightly above the F-peak. However, there is no solid basis for this height, and moreover, the height of the F-peak varies greatly, especially in the equatorial region. In the case of the IRI model, the peak height varied from 280 to 455 km under the current conditions in the area shown in Fig. 4.
For the above reason, three different shell heights, 350, 450, and 550 km, were examined. The reconstructed maps of TEC as a function of local time and latitude are shown in Fig. 5b-d. Figure 5e-g shows the reconstruction errors (difference between the reconstructed and IRI-direct TEC maps). Two types of systematic error were revealed: (1) TEC was underestimated when the shell height was small, and TEC increased with increasing shell height; and (2) the ratio of TEC at the equatorial anomaly crests to that at the trough was small irrespective of the shell height. In other words, the equatorial ionization anomaly was flattened by the reconstruction. Figure 6a-c show the bias errors for the three shell heights. The horizontal axis shows the randomly generated arc biases and the vertical axis shows the estimated arc biases. When the shell height was chosen as 350 km, the biases were overestimated and the TECs were underestimated, as shown in Fig. 5e. In contrast, the biases were underestimated for a shell height of 550 km and the TECs were overestimated, as shown in Fig. 5g. Thus, the bias and overall TEC errors compensated for each other. When the shell height was taken as 450 km, the bias error was widely scattered around the diagonal line corresponding to the flattening of the equatorial anomaly. Figure 7 shows the residual error of fitting, the root mean square error (RMSE) of the arc bias, and the RMSE of the slant TEC for different shell heights. The shell height giving the minimum fitting error was ∼ 400 km, but the

Double-shell results
In the double-shell approach, it is rational to choose one shell below the F peak and the other shell above the F peak. In addition, the TECs associated with each shell should not differ greatly. Thus, the selection of heights is limited compared with that in the case of the singleshell approach. The three cases examined were the combinations of 250/600, 300/600, and 300/700 km for the lower/higher shells. The results are shown in Fig. 8. The left panels show the vertical TEC maps for each  shell height combination and the right panels show the errors. For comparison, the same color code as shown in Fig. 5e-g is used in the error maps. At a glance, it can be seen that the three vertical TEC maps are almost identical, and the errors were less than ±1 TECu irrespective of the shell heights. The equatorial anomaly was correctly reproduced. Arc bias errors are shown in the lower three panels of Fig. 6 for the above three combinations of shell heights. The biases were also correctly estimated for all combinations and the scatter was much smaller than that for the single-shell model in the upper panels.

Application to actual observations
The advantage of the double-shell approach has been validated by simulation using the dataset generated by the IRI model. The technique was then applied to the actually observed slant TEC. The dataset has already been described, for which the ray path geometry was used to simulate the slant TEC by the IRI model. The maximum geomagnetic three-hourly K p index during 16-17 March 2014 was 1 + and the solar activity was moderate (solar radio flux, F 10.7 , at 20 UT on 17 March was 136.4 sfu (= 10 −22 Wm −2 Hz −1 )). An application to the observations under different conditions will be presented later.
To maintain consistency with the IRI simulations, the arc biases were initialized to zero (initial guess) and slant TECs with minimal biases were obtained using  Ma and Maruyama (2003). These slant TECs are not very accurate at low latitudes, where the TEC gradient is large. Prior to applying the code-phase leveling, cycle slips were corrected by the procedure described by Horvath and Crozier (2007). When there were gaps in an arc, that arc was divided into two or more independent arcs after the code-phase leveling was done, which allowed to correct the cycle-slip correction error. The accumulated error of the cycle-slip correction, the code-phase leveling, and the instrumental bias estimation were corrected as a unified residual arc bias. The leveling-free calculation in which phase data are directly used is described in the appendix. Differently from the order in the previous section, the double-shell results are first described, then a possible error in the single-shell results is discussed.

Double-shell results
If the two shell heights are suitably chosen, the TEC associated with the lower shell (referred to as the lower-shell TEC or ls-TEC hereafter) decreases and that associated with the higher shell (referred to as the higher-shell TEC or hs-TEC hereafter) increases when the ionosphere is raised by the E×B drift. On the other hand, downward plasma flux along the magnetic field line due to the gravitational force and neutral drag may decrease hs-TEC and increase ls-TEC. Thus, the TEC maps for the lower and higher shells should behave complementarily for several features.
We postulate that ls-TEC and hs-TEC should be similar amplitudes, and the two combinations of shell heights were taken as 300/600 and 280/700 km for h sLo /h sHi . Figure 9a-c show vertical TEC, ls-TEC, and hs-TEC for shell heights of 300/600 km. Figure 9d-f show the same results for shell heights of 280/700 km. The two results of vertical TEC were almost identical and the maxima of ls-TEC and hs-TEC were similar for both combinations.
Notable features of the TEC maps in Fig 9a and d were the high equatorial anomaly crests around 12-16 LT and the deep equatorial hole-like trough centered at 21 LT. The patterns of the ls-TEC and hs-TEC maps were different from each other and also from the TEC maps for either height combination. The equatorial anomaly was prominent in ls-TEC, but hs-TEC exhibited only small crests at around 13 LT. The evening equatorial hole observed in ls-TEC was wider and deeper than that in vertical TEC, but nothing was observed in hs-TEC. In the development stage of the hole in ls-TEC, hs-TEC reached a secondary peak at 18.5 LT. Such distinctive features of ls-TEC and hs-TEC are discussed later from the standpoint of ionospheric dynamics.

Single-shell results
In the single-shell runs, the same shell heights as those in the IRI simulation were chosen. The TEC maps are shown in the left panels of Fig. 10 for shell heights of 350, 450, and 550 km. Although the correct values of vertical TEC are unknown, unlike in the IRI simulations, the double-shell results are presumed to be close to the correct TEC in analogy with the IRI simulations. Then, the error of the single-shell results was evaluated with reference to the double-shell results. The right panels of Fig. 10 are TEC, the difference between the left panels and the double-shell TEC with shell heights of 300/600 km. When the shell height was taken as 350 km, TEC was negative everywhere, as shown in Fig. 10d, and the largest errors appeared near the equatorial anomaly crests and in a wide range of latitudes in the evening hours. When the shell height was 450 km, the daytime equatorial and nighttime TECs were close to the reference. However, TEC remained negative near the anomaly crests with the largest error around the evening hours, as shown in Fig. 10e. When the shell height was increased further to 550 km, TEC increased further. At around 20 LT, the error was still negative near the anomaly crests and positive around the equator.
The error of the single-shell model applied to the actual data was generally similar to that in the IRI simulations, but it was larger and more complex than in the simulations. The error in Fig. 10e with a shell height of 450 km ranged from −8 to 4 TECu, while the error in the IRI simulation with the same shell height ranged from −4 to 2 TECu, as shown in Fig. 5f. The largest error in the actual observations was concentrated around the hour of the evening enhancement in the right panels of Fig. 10.

Residual biases in preprocessed slant TECs
Although it was not the purpose of this paper to validate any other techniques of TEC estimation, it is worthwhile to examine the residual arc biases in the preprocessed slant TECs. The biases contained errors originated in a limit of fixed single-shell height (this study), the intra-day variations of instrumental biases (Ciraolo et al. 2007), inadequate application of the simple method of bias estimation at low latitudes (Ma and Maruyama 2003), as well as the cycle-slip correction and the code-phase leveling. The residual arc biases are shown in Fig. 11. The error was large and negative for CMU0, UDON, and NKSW, which were close to the equatorial anomaly crest. The error scattered mostly between −25 and 15 TECu, and the random biases ( ±25 TECu) added to slant TECs in the IRI simulations in "Validation by IRI simulation" section were on a level with this.

Ionospheric dynamics
In the previous sections, it was shown that the doubleshell model was more suitable for capturing the vertical TEC variations correctly than the single-shell model, especially in the equatorial anomaly region. Furthermore, ls-TEC and hs-TECs varied in a different way with the local time and latitude, which can be ascribed to the ionospheric dynamics. In other words, the different behaviors of TECs associated with each shell are expected to provide information on the ionospheric dynamics such as the E×B drift and the field-aligned flux induced by the gravitational force and neutral drag. From this perspective, the results in Fig. 9 for the individual shells were interpreted.

Complementary features between two shells
The diurnal variations of vertical TEC, ls-TEC, and hs-TEC away from the equator ( 12 • N magnetic latitude, 15.9 • N geographic latitude as shown by the red dashed line in Fig. 9a-c) and at the magnetic equator ( 7.6 • N geographic latitude as shown by the dash-dotted line in Fig. 9a-c) are shown in Fig. 12 for the shell height combination of 300/600 km.
TEC first started to increase in the lower shell at around sunrise and then in the higher shell with a delay of 1-2 h. The increase in ls-TEC is an immediate response to the start of photoionization at 150-200 km. However, the increase in hs-TEC is due to the refilling of the topside ionosphere by the upward plasma flux along the magnetic field line, which caused the delayed response. The upward E×B drift might also have contributed to the increase in hs-TEC.
A notable change occurred at the off-equator latitude when ls-TEC reached a peak slightly before 11 LT. After that time, ls-TEC and hs-TEC varied complementarily or out of phase with each other for two cycles until 22 LT. A similar variation started 1.5 h earlier at the magnetic equator and continued until 23 LT with the exception of the period of 18.5-20 LT, during which ls-TEC markedly decreased. These complicated variations were the manifestation of the ionospheric dynamics, and the difference in TEC between the two shells was further interpreted by focusing on events from hour to hour.

Development of equatorial anomaly
The equatorial anomaly started developing at around 10 LT, as shown in Fig. 9a. At the same time, ls-TEC at the equator started decreasing and formed a bite-out around noon in Figs. 9b and 12b. On the other hand, hs-TEC at the equator continued to increase and reached a peak slightly after the lower-shell bite-out in Figs. 9c and 12b. Figure 13a shows the latitudinal variation of the TECs at 11 LT (dashed lines) and 12 LT (solid lines). During the intervening hour, hs-TEC increased and ls-TEC decreased at all latitudes. The chain of TEC variations is ascribed to the intense upward E×B drift, consistent with the development of the equatorial anomaly.
After 12.5 LT, ls-TEC started to increase and hs-TEC started to decrease in Fig. 12. This is interpreted as a weakening of the eastward electric field and the dominant field-aligned plasma diffusion toward lower latitudes in the topside ionosphere. While the downward plasma flux maintained the large crest TEC in Fig. 9a, the crest latitude receded toward the equator as a consequence of the weakening of the electric

Evening enhancement
A transition occurred at 17 LT when hs-TEC was increasing again and ls-TEC was decreasing in Fig. 12. The secondary peak of hs-TEC appeared at 18-19 LT in Fig. 9c, which was earlier and more intense at latitudes closer to the magnetic equator. Slightly delayed from the hs-TEC peak, a large depression appeared in ls-TEC in Fig. 9b. This is clearly due to the effect of the evening enhancement of the upward E×B drift (Fejer et al. 2008). The latitudinal variations of the TECs at 17 and 18 LT are shown in Fig. 13b. Unlike in Fig. 13a, the decrease in ls-TEC was significant and the twoshell TEC never increased, which is due to the absence of further photoionization. The center of the equatorial depression or the hole of ls-TEC was at 20 LT in Figs. 9b and 12b, approximately one hour earlier than that of the center of the two-shell TEC depression in Figs. 9a and 12b. The reversal of E×B drift was estimated to occur at 20 LT. Figure 13c and d exhibit the reversal of the magnitudes of ls-TEC and hs-TEC in intervening hours caused by the large downward E×B drift. Associated with the reversal of the electric field, the equatorial anomaly crests receded toward the equator after 20 LT in Fig. 9a.

Hemispheric asymmetry
Another interesting feature observed in Fig. 9 is the north-south asymmetry of TEC. The asymmetry was prominent at the anomaly crests in the ls-TEC maps of Fig. 9b and e, in which the northern anomaly crest was higher than the southern anomaly crest during daytime and in the evening hours. The asymmetry was reversed at 21 LT and continued until midnight. In the meridional plots, ls-TEC was less than hs-TEC near the magnetic equator in daytime and evening hours, as shown in Fig. 13a-c. The magnitudes of ls-TEC and hs-TEC reversed at off-equatorial latitudes in the northern hemisphere. In Fig. 13a, during the development phase of the equatorial anomaly, hs-TEC increased and ls-TEC decreased at all latitudes in the hour from 11 to 12 LT, which has already been ascribed to the upward movement of the F peak due to the vertical E×B drift. The increase in hs-TEC and the decrease in ls-TEC were more significant at off-equatorial latitudes in the southern hemisphere than in the northern hemisphere. In other words, the F peak in the southern hemisphere rose more rapidly than that in the northern hemisphere, which is ascribed to the effect of northward field-aligned plasma flow.
In Fig. 13b, the equatorial anomaly almost disappeared in the two-shell TEC. The curve of hs-TEC was convex upward and that of ls-TEC was convex downward, resulting in arch-shaped iso-electron density contours across the magnetic equator. The decrease in hs-TEC and the increase in ls-TEC toward off-equatorial latitudes were large in the northern hemisphere. Plasma pileup in the lower shell and a depression in the higher shell occurred in the northern hemisphere, which are ascribed to the downward (northward) fieldaligned plasma flow.
At 20 LT, the lower shell became almost empty around the equator, and the equatorial anomaly developed greatly in the two-shell TEC as shown in Fig. 13c. The plasma pileup persisted in the lower shell in the northern hemisphere, but hs-TEC was symmetric about the equator. The equatorward neutral wind might be dominant in the northern hemisphere. Two hours later, the relative magnitudes of ls-TEC and hs-TECs was reversed, as shown in Fig. 13d, indicating the reversal of the E×B drift. The large reduction in TEC in the northern hemisphere is striking. For comparison, the two-shell TEC at 20 LT is reproduced in Fig. 13d by the purple dotted line. At the same time, the northsouth asymmetry of ls-TEC was reversed. hs-TEC also exhibited a greater reduction in the northern hemisphere, which suggests that the reduction in TEC was due to chemical recombination under a large downward E×B drift. The reduction in TEC in the southern hemisphere was not significant in spite of the large downward E×B drift. The north-south asymmetry of TEC reduction is ascribed to the ion drag of northward neutral wind, such that the downward field-aligned plasma flow enhanced the descent of the F layer due to the downward E×B drift in the northern hemisphere and, in contrast, the upward field-aligned plasma flow resisted the descent of the F layer due to the downward E×B drift in the southern hemisphere.
From the examination of the meridional variation, the hemispherical asymmetry of ls-TEC and hs-TEC was ascribed to the field-aligned plasma flow combined with the vertical E×B drift and chemical recombination. The plasma flow was due to the ion drag of the thermospheric neutral wind, which was directed northward in most of the hours until ∼ 24 LT.

Day-to-day variability
The double-shell approach to TEC estimation has been shown to provide information on ionospheric dynamics. Forces acting the ionosphere are the thermospheric neutral wind and the zonal electric field, which vary from day to day as well as with the season. Figure 14 shows TEC maps for 7 consecutive days centered on the day analyzed in the previous sections. The solar activity during these days was moderate from F 10.7 = 136 to 151 sfu and the maximum K p index was 3. In the two-shell TEC maps of the left column, a similar pattern of the daytime equatorial anomaly was repeated. The amplitude and time of the anomaly crests exhibited day-to-day variability, reflecting the variability of the zonal electric field. An equatorial hole of ls-TEC at around 20 LT in the middle column and the preceding enhancement in hs-TEC in the right column, which are the feature of the evening enhancement of the zonal electric field, were also repeated every day. There observed a notable feature in the north-south asymmetry of shell-TECs. At 10-11 LT (vertical orange dotted line in the spaces between ls-TEC maps), the off-equatorial ls-TEC was higher in the northern hemisphere than that in the southern hemisphere until 18 March 2014 and it reversed on 19 and 20 March. Contrary, at 13 LT (vertical purple dotted line in the spaces between hs-TEC maps), the off-equatorial hs-TEC was higher in the southern hemisphere than that in the northern hemisphere before In analogy with the previous analysis, the shell-TEC maps suggest that the thermospheric wind pattern was changed around these days.

Application under different geophysical conditions
In the previous section, the ionospheric dynamics were inferred from the double-shell results under the conditions of the equinox and a moderate solar activity. The vertical ionospheric drift induced by the zonal electric field and the field-aligned plasma flow caused by the thermospheric neutral wind vary with the season and the solar activity. To ascertain the validity of the interpretation, the double-shell model was applied to the observations in the June solstice and low solar flux conditions from 22 UT, 21 June to 22 UT, 22 June 2014. F 10.7 on 22 June was 94.2 sfu and the maximum K p index was 2.
The results are shown in Fig. 15. The maximum TEC near the equatorial anomaly crests was ∼ 50 TECu and greatly decreased from the March results shown in Fig. 14, reflecting the low solar activity. A notable difference between the TEC maps for the two periods was an absence of the evening hole in two-shell and ls-TECs as well as an absence of the enhancement of hs-TEC preceding the hole. In the empirical model of Scherliess and Fejer (1999), the evening enhancement of the equatorial vertical drift is prominent in the equinoxes (March-April, September-October) in both low and high solar flux conditions. On the other hand, the evening enhancement disappears in the June solstice (May-August) in low solar flux conditions ( F 10.7 = 90 ). The disappearance of the feature indicating the evening enhancement in Fig. 15 was consistent with the Scherliess and Fejer model. Figure 16a shows the latitudinal variation of twoshell, ls-, and hs-TECs at 12 LT. ls-TEC in the southern hemisphere was higher than that in the northern hemisphere, while hs-TEC in the southern hemisphere was lower than that in the northern hemisphere. For  Fig. 16c. The hemispheric asymmetry in the June solstice was opposite to that in March. Similarly, at 22 LT, hs-and ls-TECs in the northern hemisphere were higher than those in the southern hemisphere, as shown in Fig. 16b, which were also opposite to the March results in Fig. 16d. In analogy with the discussion in "Hemispheric asymmetry", the field aligned plasma flow was reversed in the two seasons. We should note that the magnetic equator at these longitudes is at ∼ 7.6 • N and the sun is magnetically in the southern hemisphere in March.
Finally, double-and single-shell models were compared. Figure 17a-c are the single shell TECs for the shell heights of 350, 450, and 550 km and Fig. 17d-f are the difference ( TEC ) between the single-and doubleshell TECs. The general tendency of lower (higher) TEC for the smaller (larger) shell height was the same as the March results in Fig. 10. A notable difference between the June solstice and March was the absence of a large error in the evening hours in June, which was ascribed to the absence of the evening enhancement of the vertical drift in the June solstice. The largest negative error was near the southern anomaly crest in daytime from 12 to 16 LT in the June solstice, as shown in Fig. 17d-f. On the other hand, the daytime largest negative error was near the northern anomaly crest in March, as shown in Fig. 10d-f. In both cases, the layer height was depressed in the areas of a large negative error in the single-shell results. As already discussed, the depression of layer height was due to the poleward neutral wind. Through the analyses of different seasons and solar activities, it was generally concluded that the large errors in the single-shell model were caused by a large change in the layer height due to the zonal electric field and the thermospheric meridional winds.

Summary
The thin-shell approximation of the ionosphere is a widely used technique to estimate ionospheric TEC from the trans-ionospheric radio waves from GNSS satellites. In most cases, single-shell approaches are taken, assuming a fixed shell height. Ionosphere pierce points of the slant ray path vary depending on the shell height. In spite of this, there is no clear criterion for choosing this height, and furthermore, the ionospheric height changes widely with the time and place. As a result, the vertical TEC and instrumental biases, which are simultaneously determined, depend on the shell height. To overcome this difficulty, a new technique was developed by taking a double-thin-shell approach.
The vertically distributed electron density was parameterized by two thin-shell layers. The spatiotemporal variation of TEC (strictly speaking, partial electron content) associated with each shell was approximated by the functional fitting of spherical surface harmonics. For this application, dataset were prepared for 24 hours by the receiver network along the meridian at the southeast Asian longitude across the magnetic equator.
To evaluate the performance of the double-shell model, a single-shell model, which was a reduced version of the double-shell model, was also examined. The technique was first validated for a simulated slant TEC using the IRI electron density profile and artificially added random biases. The vertical TEC estimated from the slant TEC was compared with the known values of vertical TEC directly calculated by IRI (IRIdirect TEC). The single-shell model yielded different TECs and biases depending on the shell height. When the shell height was small, TEC was generally low and vice versa. Regardless of the shell height, the equatorial anomaly was flattened compared with the IRI-direct TEC. On the other hand, the double-shell approach successfully reproduced the IRI-direct TEC and randomly added arc biases. Also, the estimated vertical TEC and arc biases were less dependent on the choice of shell heights.
In the application to real data obtained by the receiver network, the general characteristics of the difference between the single-and double-shell models were similar to those in the application to the IRI simulations. The double-shell results provided not only accurate TEC values but also comprehensive information on the ionospheric dynamics, such as the vertical ionospheric drift induced by the zonal electric field and the field-aligned plasma flow induced by the thermospheric neutral drag, that was not extracted by single-shell approaches.
The double-shell algorithm developed in this paper was easily extended to the three-shell model in mathematics. The three-shell model yielded a smaller fitting error than that of the double-shell model. However, the improvement was not so distinct as the improvement brought by the double-shell model. Moreover, an overfitting problem has arisen because of the increase in the degrees of freedom. Careful examinations of each shell-TEC are required and we must choose the model parameters more carefully, which will be a future issue.
The iteration converged in about a minute for the number of receivers and the degree and order of the harmonics in this study (see the appendix more). The Fortran source code used in this study is available on request.

Appendix: Leveling-free calculation
In the main text, slant TECs were preprocessed using the code-phase leveling and the coarse instrumental biases estimated by the method described in Sect. 4 of Ma and Maruyama (2003) (method A). An example of the preprocessed slant TEC is shown in Fig. 18a for SOKA. The discontinuity of arcs at 00 UT was due to the separate data processing for each UT day. The code-phase leveling, however, is not necessarily required to prepare input slant TECs in the arcbias-based model. In an alternative method, the slant TEC arcs were adjusted to zero at the start of each arc (method B), as shown in Fig. 18b for the same receiver. Method B was successful and iteration converged quickly in the single-shell model. But iteration converged slowly in the double-shell model. A two-step method or combined use of single-and double-shell models accelerated the convergence. Approximate arc biases were obtained first by the single-shell model run applied to slant TECs prepared by method B. In the second step, improved arc biases were obtained by the double-shell model run starting from the initial guess of arc biases obtained by step one (method C). The efficiency and performance of three methods are compared below.
Vertical TEC maps obtained by the different methods were compared instead of the fitting error of slant TEC. Figure 19a shows the TEC map obtained by method A with a sufficiently large number of iterations (10000 in the online mode and 5000 in the batch mode), which was used as a reference for evaluating the convergence. The CPU time used for this run was 295 s when the program code was executed on Intel Zeon (2.8 GHz) processor. Figure 19b shows the TEC map obtained by method A with 2000 iterations in both online and batch modes. The TEC error was less than 0.5 TECu in most of the map area, as shown in Fig. 19g, and the iteration almost converged. The CPU time for this run was 79 s. Figure 19c shows the TEC map obtained by method B with the same iterations as Fig. 19b. The equatorial anomaly crest was low in the northern hemisphere and high in the southern hemisphere compared with the reference. The difference is shown in Fig. 19h. The systematic error in north-south asymmetry of TEC was fatal for a discussion on the ionospheric dynamics based on TEC maps like that in "Hemispheric asymmetry". Figure 19d and e are the TEC maps after 5000 and 8000 iterations in the online mode, respectively. Figure19i and j are the corresponding TEC error. There slightly remained a negative error ( < −1 TECu) at the northern end and a positive error ( > 1 TECu) at the southern end of the map, even after 8000 iterations. Figure 19f shows the TEC map by method C after 2000 iterations in both online and batch modes. The CPU time for this run was 90 s, including 10-s step-one calculation by the single-shell model. The error shown in Fig. 19k was on a level with Fig. 19g.