GNSS Differential Code Bias Determination Using Rao‐Blackwellized Particle Filtering

The Assimilative Canadian High Arctic Ionospheric Model (A‐CHAIM) is a near‐real‐time data assimilation model of the high latitude ionosphere, incorporating measurements from many instruments, including slant Total Electron Content measurements from ground‐based Global Navigation Satellite System (GNSS) receivers. These measurements have receiver‐specific Differential Code Biases (DCB) which must be resolved to produce an absolute measurement, which are resolved simultaneously with the ionospheric state using Rao‐Blackwellized particle filtering. These DCBs are compared to published values and to DCBs determined using eight different Global Ionospheric Maps (GIM), which show small but consistent systematic differences. The potential cause of these systematic biases is investigated using multiple experimental A‐CHAIM test runs, including the effect of plasmaspheric electron content. By running tests using the GIM‐derived DCBs, it is shown that using A‐CHAIM DCBs produces the lowest overall error, and that using GIM DCBs causes an overestimation of the topside electron density which can exceed 100% when compared to in situ measurements from DMSP.


Introduction
The Assimilative Canadian High Arctic Ionospheric Model (A-CHAIM) is a near-real-time data assimilation model of high latitude ionospheric electron density (Reid et al., 2023).It uses a particle filter technique to assimilate data from ground-based Global Navigation Satellite System (GNSS) receivers, ionosondes, and satellite-borne altimeters on the Jason-3 and Sentinel satellites.All available observations above 45°magnetic latitude within the last 3 hr are assimilated, producing a nowcast and short-term forecast of the high latitude ionosphere.Of the available data sources, ground-based GNSS receivers are by far the most widely distributed and numerous.GNSS receivers provide slant Total Electron Content (sTEC) measurements along line-of-sight from the satellite to the receiver, usually expressed in TEC Units 1 × 10 16 m − 2 (TECU).This measurement is subject to instrumental biases which arise from both the satellite transmitter and receiver hardware, and so to obtain an absolute measurement of ionospheric sTEC these biases must be determined.As A-CHAIM is in continuous operation, and the instrumental biases are subject to change, these biases must be resolved on an ongoing basis.
GNSS satellites broadcast on multiple carrier frequencies in the Ultra High Frequency (UHF) band.These signals are encoded with information about the state of the clock on board the satellite, which when compared to the clock at a receiver allows for the determination of the apparent range between the satellite and receiver.This apparent range, or pseudorange, diverges from the true range due to many factors, including the effects of the ionospheric plasma on the propagation of the signal.As the ionosphere is a dispersive medium, ionospheric group delays and phase advances are dependent on the signal frequency.Using a geometry-free combination of the phase and code observables recorded on each GNSS carrier frequency, where the observables from each frequency are differenced to remove non-dispersive effects, the TEC can (to a first-order approximation) be related to the observables by Equation 1, where A = 40.3,f m is the mth frequency, Δϕ is the difference in the signal carrier phases, DCB rcv and DCB sat are the receiver and satellite Differential Code Biases (DCB) caused by instrumental delays, and W is a phase-leveling term used to correct an integer ambiguity in the phase-derived TEC using the code observables (Reid et al., 2023;Themens et al., 2013).

Existing DCB Estimation Techniques
Several methods exist to resolve the DCBs.From the perspective of an end user, the most straightforward method is to simply use a published estimate, if it exists.The International GNSS Service (IGS) has a network of reference stations, which are used to produce estimates of satellite and receiver biases, satellite clock errors, and orbit determination.The IGS also commissions various worldwide Analysis Centers, which produce their own products using various techniques (Roma-Dollase et al., 2018).These products are widely used, and have proved to be a reliable tool for operational users.As an example, A-CHAIM uses the satellite biases from the Chinese Academy of Sciences (CAS) distributed in the SINEX format (Schaer, 2018), as well as the precise GPS orbits from the IGS.The discrepancies in satellite biases between each Analysis Center are small, on the order of a few tenths of nanoseconds, or less than 1 TECU for the L1/L2 GPS combination (Hernández-Pajares et al., 2009).As the GPS satellites DCBs are common for all receivers worldwide, it is possible to correct any sTEC data using these published satellite DCBs.This is not true for the receiver DCBs, which are unique to each receiver.One clear limitation of published receiver DCBs is that they only exist for those stations included in the IGS data sets.
During August 20th to 10 October 2022, only 76 of the 662 unique stations used in A-CHAIM were also included in the CAS data set.For every other station some other technique would be needed.
Equation 2 is a simplified expression for the sTEC, where the DCBs in Equation 1 have been converted from meters to TECU, and sTEC observed = f 2 1 f 2 2 (Δϕ + W)/ A( f 2 1 − f 2 2 ) .One method to resolve the DCB of an arbitrary receiver is by comparing the data sTEC observed to an ionospheric reference.If the reference value is assumed to be the true ionosphere sTEC true , then subtracting it from the observed TEC provides an estimate of the DCB.
The IGS Analysis Centers also produce maps of vertical Total Electron Content (vTEC), known as Global Ionospheric Maps (GIMs).The GIMs are not independent of the published DCBs, but are instead produced simultaneously as a self-consistent product.The GIMs provide global coverage, and so provide a reference from which any receiver's DCB can be determined (Arikan et al., 2008).These products are distributed by the Crustal Dynamics Data Information System (CDDIS) in the IONosphere Map EXchange format (IONEX) (Noll, 2010;Schaer et al., 2015).An IONEX file contains a series of maps of global vTEC, at fixed time intervals on a geocentic latitude-longitude grid.The IONEX format specifies several techniques to interpolate between times, with the preferred technique being to rotate the maps in local time before performing a spatial interpolation (Schaer et al., 2015).The observed sTEC can then be converted to vTEC using a thin-shell approximation Equation 3, where the ionospheric electron density is assumed to exist entirely in a spherical shell, generally at a fixed altitude h shell .The mapping function from sTEC to vTEC for a thin-shell ionosphere is given in Equation 4where e is the elevation angle of the satellite, R e is the radius of the Earth.The IONEX products from various IGS Analysis Centers used in this analysis are summarized in Table 1 (Feltens, 2007;Ghoddousi-Fard, 2014;Ghoddousi-Fard et al., 2011;Hernández-Pajares et al., 1999, 2009;Iijima et al., 1999;Komjathy et al., 2002;Li et al., 2015;Mannucci et al., 1998;Schaer et al., 2021).These GIM products use several different ionospheric representations internally, with the single thin-shell ionosphere being the most common.Some groups use more complicated representations to capture more of the vertical structure of the ionosphere.The JPL GIM uses three thin shells (Komjathy et al., 2002), and the UPC GIM uses a two-layer voxel technique (Hernández-Pajares et al., 2009).The IGS GIM is produced as a weighted combination of the COD, ESA, JPL, and UPC products.
While each of these products may use different representations of the ionosphere internally, when distributed in IONEX format each GIM is represented as a single thin-shell.vTEC = sTEC ⋅ M(e) (3) To estimate the DCB of a receiver using GIM vTEC, first the point where each sTEC observation intersects the spherical shell is determined, a point known as the Ionospheric Pierce Point (IPP).The vTEC at each of these points in space and time can then be determined by interpolating the GIM.This vTEC can then be converted to an expected sTEC value using the mapping function in Equation 3. The resulting DCB can then be estimated by comparing the observed sTEC to the GIM vTEC as in Equation 5, where w i is a weight for each observation based on the error from the phase-leveling process (Reid et al., 2023).In all comparisons below, this leveling process is performed by considering a full 24 hr of observations, spanning a single daily GIM file.Equation 5 is a modification of the technique specified in Schaer et al. (2015) to include error weighting.
This technique is not limited to the IONEX format.Other groups produce vTEC maps, notably the Madrigal vTEC maps from the Haystack Observatory at MIT (Coster, 2022;Rideout & Coster, 2006).These maps use data from many thousands of receivers, and use a sophisticated bias estimation technique, along with a Chapman-layer vertical parameterization (Vierinen et al., 2016).The same thin-shell bias estimation technique used for the IONEX vTEC maps can also be used with the Madrigal vTEC maps.
It has been shown that there are global biases on the order of a few TECU between different IGS Analysis Centers (Hernández-Pajares et al., 2009;Li et al., 2015;Mazzella Jr, 2012).Validation efforts often use vTEC measurements from the Jason altimeter satellites (Hernández-Pajares et al., 2009, 2020;Li et al., 2015;Roma-Dollase et al., 2018), partial profiles from Incoherent Scatter Radar (Themens et al., 2013), or assuming some fixed vTEC at high latitudes to provide a reference (Stephens et al., 2011;Yizengaw et al., 2008).These validations are always subject to some small ambiguity on the order of a few TECU.Lighting-derived vTEC maps (Lay et al., 2022) could be used as an independent measurement to assess individual GIMs, but due to the limited lightning activity at high latitudes are not directly applicable to the region of interest.With the current literature, it is not clear which GIM product would produce the best result.Single station techniques such as the minimization of standard deviation Ma and Maruyama (2003) are able to provide an estimate without any external reference.The performance of this technique in the high latitude region was assessed in Themens et al. (2013) and was found to produce discrepancies as large as 6 TECU when compared with ISR-derived TEC.The Self-Calibration Of pseudoRange Errors (SCORE) algorithm (Bishop et al., 1996) uses conjunctions between satellites to fix the receiver DCB without an external ionosphere model, assuming a thin-shell ionosphere and sun-fixed variation.The SCORE method has been shown to produce DCBs estimates accurate to within 1 TECU for midlatitude receivers under smooth ionospheric conditions, and when rays which intersect the plasmasphere are excluded (Lunt et al., 1999b).
Another limitation of the IONEX GIMs is that they cannot distinguish between ionospheric TEC and plasmaspheric TEC (Lunt et al., 1999d).Plasmaspheric TEC is expected to contribute up to 50% of sTEC for midlatitude stations during winter night (Lunt et al., 1999a).Anghel et al. (2009); Anghel et al. (2008); Carrano et al. (2009) used Kalman Filter-based techniques to estimate both receiver DCBs and the relative contribution of ionospheric and plasmaspheric electron content.These studies also found that neglecting the plasmaspheric total electron content resulted in an overestimation of ionospheric TEC at midlatitudes.Mazzella (2009) used a variant of the SCORE method, called SCORE for Plasmasphere and Ionosphere (SCORPION), which also found that neglecting the plasmasphere produced an overestimation in midlatitude ionospheric vTEC.This is in agreement with simulations of the SCORE technique conducted with a model plasmasphere in Lunt et al. (1999b), and in experimental results in Lunt et al. (1999cLunt et al. ( , 1999d)).Mazzella (2012) in a direct comparison of SCORE, SCOR-PION, and the UPC, COD, and JPL GIMs, found that plasmaspheric TEC causes a latitudinally dependent error in DCB determination at midlatitudes.In all of these studies, the reported effect on the thin-shell vTEC is relatively small, on the order of 1 or 2 TECU, but consistent.

Effects of DCB Errors on Electron Density
It would be possible to use a GIM with Equation 5 to fix the receiver DCBs, but relying on an external model for calibration would result in A-CHAIM inheriting the limitations of that model.The potential effects on the assimilation must be considered when incorporating external references.For example, the satellite DCBs provided by all IGS Analysis Centers are constrained to have a zero mean across all satellites.This is a mathematical convenience to fix a free parameter, and imposing this constraint has no impact on the assimilation.If the mean satellite DCB were altered by some value δTECU, the receiver DCBs would appear to change by an equal and opposite amount − δTECU to stay consistent with the observed TEC in Equation 2. If the relative differences between satellites are accurate, the zero mean constraint imposed on the satellite DCBs should not introduce systematic errors.Once the mean value of the satellite DCBs is fixed, the reconstructed ionospheric state is now sensitive to systematic errors in the receiver DCBs.If there were a systematic error in the receiver DCBs then this error would need to be absorbed by the ionospheric model itself.
The potential effects of such an error can be demonstrated using a simple toy model.Figure 1 shows an A-CHAIM ionospheric profile, which is assumed to be the true ionospheric state for this example.In this idealized situation, a GNSS receiver provides a noiseless, vertical TEC measurement subject to some DCB.The parameters which control the shape of the profile are adjusted so that the predicted data matches the observed data.Any error in the receiver DCB must produce errors in the reconstructed electron density.Receiver DCB errors on the order of a few TECU can have surprisingly large effects on the electron density profile.In the first plot of Figure 1, the bottomside thickness H Bot is modified.The impact on the overall profile is dramatic, with small changes in DCB producing clearly unphysical profile shapes.In the second plot, the peak density NmF2 is modified, producing moderate errors as a percentage of electron density.If the bottomside ionosphere is well-constrained by ionosonde measurements, as in A-CHAIM, then the only part of the profile which can absorb any potential error is the topside.This is demonstrated in the third plot of Figure 1.As the topside density is lower, a greater proportional adjustment is needed to maintain consistency with the observed TEC.It is preferable to solve for the ionospheric state and the receiver DCBs simultaneously, as the sensitivity of the ionospheric reconstruction to receiver DCBs is on the same order of magnitude as the differences between DCB products.Other data assimilation models have found improved performance when solving for receiver biases self-consistently (Dear & Mitchell, 2006).

Method
A-CHAIM uses a technique known as Rao-Blackwellized particle filtering (Doucet & Johansen, 2009) to resolve the receiver DCBs.This allows A-CHAIM to efficiently determine the DCBs in real time, by finding the optimal set of DCBs for each ensemble member.A test environment was prepared to perform offline runs of the A-CHAIM system, using data from 20 August 2022 though 10 October 2022.All of the GNSS and altimeter data used in the tests were those collected in near-real-time by the online assimilation system (Reid et al., 2023).In total, 662 unique GNSS stations provided data during the test period.Due to network outages, many ionosonde measurements were unavailable to the online system during this period, so these were added back in for the sake of the tests.
In this study a total of 12 test runs of the A-CHAIM system were performed, as summarized in Table 2. Four of the test runs used the Rao-Blackwellized DCB estimation technique presented in this work.The first two were conducted by assimilating the usual set of instruments available to A-CHAIM in near-real-time (Reid et al., 2023).One of these test runs is identical to the operational A-CHAIM system, which does not include a plasmasphere model.The second, and all subsequent test runs, includes the addition of the Neustrelitz plasmasphere model (NPSM) (Jakowski & Hoque, 2018).This will help determine the influence of plasmaspheric Total Electron Content on the A-CHAIM DCBs.
Figure 1.Simulated effect of DCB errors on reconstructed ionospheric density.Each E-CHAIM profile is altered to maintain a consistent observed vTEC (integrated up to 2,000 km), given an incorrect receiver DCB.The plots show the results for three different parameters, H Bot , NmF2 and H Top respectively.The black profile represents the "true" ionospheric profile.Each profile is drawn at 1 TECU steps, marked with a colored point.The profiles are colored according to the percent error in the reconstructed electron density at that altitude.If the bottomside ionosphere is fixed by an ionosonde measurement of NmF2, the only remaining parameter to adjust is the topside thickness H Top , which causes dramatic changes at high altitudes.
Note.All runs used the full complement of instruments that are normally available to the real-time system.vTEC refers to using the Madrigal vTEC data as an assimilated measurement.
The third and fourth test runs which used the Rao-Blackwellized DCB estimation included an additional data set in the assimilation.These runs assimilated the vTEC values from the Madrigal vTEC maps as measurements, by integrating the full electron density profile for each grid point from 0 km to GPS altitudes at 20,000 km.A-CHAIM uses a pre-filtering technique on a subset of the assimilated data to improve sampling efficiency (Reid et al., 2023), and the Madrigal vTEC measurements were included in this pre-filtering data, along with the ionosondes.This will constrain A-CHAIM to match the Madrigal vTEC values as closely as possible, while leaving the DCBs to adjust accordingly.In testing it was found that incorporating the hundreds of Madrigal vTEC values tended to overwhelm the limited number of ionosondes, causing A-CHAIM to overestimate NmF2.The second vTEC run used an adjusted pre-filtering technique to preferentially modify the A-CHAIM topside thickness H Top to match Madrigal.The resulting DCBs from all four Rao-Blackwellized A-CHAIM runs will then be compared to those determined by leveling to GIMs as in Equation 5.For determining the GIM-derived DCBs, only sTEC observations with an elevation angle greater than 30 o were used.
The final 8 test runs of the A-CHAIM system used GIM-leveled DCBs rather than using the Rao-Blackwellized DCB estimation technique.If any of the GIMs provide an advantage over the self-consistent DCBs, the resulting test run should produce superior representation of the ionospheric electron density.This analysis also provides an opportunity to validate the GIMs produced by the various groups listed in Table 1.Given the example presented in Figure 1, the reconstructed topside electron density should be sensitive to small errors in DCBs.This will provide a quantitative measure of the effects of external DCBs on GNSS sTEC data assimilation.

Particle Filtering
A-CHAIM is vertically parameterized as a semi-Epstein layer, whose shape is controlled by a set of harmonic expansions of several key ionospheric parameters (Reid et al., 2023).These parameters include the peak density of the F2 layer, NmF2, the altitude of the F2 peak hmF2, as well as thickness parameters H Top and H Bot which control the shape of the topside and bottomside ionosphere, respectively.The electron density at any point ⃗ r is therefore a nonlinear function N e (x s , ⃗ r) of these harmonic coefficients x s , which can be described as a vector in a state space X s .The subscript s is included to indicate these elements are part of the ionospheric state.Given some set of ionospheric observations y n at time t n , A-CHAIM estimates an optimal set of coefficients x s,n , and so can produce a model of ionospheric electron density.To model the time evolution of the ionosphere, a sequence of states x s,1:n = {x s,1 , x s,2 , …, x s,n− 1 , x s,n } must be determined.This requires evaluating p(x s,1:n |y 1:n ), the probability distribution of a sequence of states x s,1:n conditioned on the observations y 1:n .

p( x s,1:n |y
Evaluating Equation 6is a highly nonlinear inverse problem, and so in A-CHAIM a particle filter technique is used (Reid et al., 2023).A particle filter is a Monte Carlo technique which uses an ensemble of sample points, or particles, X i s ∈ X s with associated statistical weights W i to approximate a distribution π(x s ) on X S (Doucet & Johansen, 2009).
The particles X i s are sampled from an importance distribution q(x s,1:n |y 1:n ).Choosing an optimal importance distribution is a critical part of particle filter design, but the precise form of q(x s,1:n |y 1:n ) is not important in this context.
The forecast model f(x s,n |x s,n− 1 ), gives the probability of transitioning from a state x s,n− 1 to a state x s,n .This, with the likelihood function p(y 1:n |x s,1:n ) allows A-CHAIM to constantly update the weights of the particles w i with Equation 10.

Space Weather
10.1029/2023SW003611 After normalizing the weights w n (x s,1:n ) the sum in Equation 7takes the following form: Equation 12 provides an empirical approximation to the full probability distribution Equation 6, allowing for the estimation of statistical properties of the ionospheric state.While the unnormalized weights w i n are calculated using the particle filter, only the normalized weights W i n are used directly in the approximate solution p( x s,1:n |y 1:n ) .

Rao-Blackwellized Particle Filtering
As a Monte Carlo technique, the performance of a particle filter is directly dependent on the number of particles in the ensemble.As the number of dimensions of the state space X s increases, the number of particles required to sample the space appropriately increases dramatically, posing significant challenges for data assimilation models (Reid et al., 2023;van Leeuwen et al., 2019).This presents a problem when attempting to fix the DCBs in A-CHAIM.x s,1:n ∈ X s are the components of the state space which control the ionospheric density in A-CHAIM as above.By including the DCBs as parameters to be determined, a new set of numbers x b,1:n ∈ X b is added to the state.The subscript s indicates values belonging to the ionospheric state and subscript b corresponds to those of the DCBs.The new state space is the product of these two sets of parameters Equation 13.
The probability distribution p(x 1:n |y 1:n ) can therefore be rewritten as Equation 14.
While it would be possible to approximate Equation 14 with a particle filter, each GNSS receiver added to the A-CHAIM data set adds a new DCB, increasing the number of dimensions of the state space.The resultant undersampling can produce the undesirable situation where including additional data to the assimilation produces a worse outcome.The total number of particles useable by the system is limited by the computational resources available, and so it is not practical to compensate by naively adding more particles.An efficient solution to this problem is Rao-Blackwellized particle filtering, which enables an analytical solution to the DCBs (Doucet & Johansen, 2009).The conditional probabilities in Equation 14 can be expressed as: Space Weather

10.1029/2023SW003611
By re-arranging and combining Equation 15with Equation 14, the problem can be restated in the form Equation 16.
The first term on the right hand side of Equation 16 is simply Equation 6 for the ionospheric state alone.When expanded fully, the new empirical distribution becomes Equation 17.
The unnormalized weights w i n can now be calculated using Equations 18 and 19.
Recall that the X i s,k in Equation 19 are the X i k of the original particle filter Equation 10.The new weights are therefore identical to the original weights, with the exception of the factors p( X i b,k |y k ,X i s,k ) .If this expression can be evaluated analytically, the additional dimensions of state space added by the DCBs do not need to be sampled by the particle filter.

DCB Modeling
In A-CHAIM it is assumed the DCB behaves as a Gaussian random walk with an average step size ̅̅̅̅̅ ̅ Q b √ = 0.05 TECU over a 5 min assimilation step.With Q b as a diagonal square matrix of size N rcv , the number of DCBs, the forecast model for the DCBs is expressed as Equation 20.
The forward model operator H n (x) allows the reproduction of a set of measurements for some x, which takes the form Equation 21 for sTEC observations.Here 1 ≤ k ≤ m is an index over all m observations in y n .
While the integrals through the modeled ionosphere Equation 21 are nonlinear with respect to the ionospheric state x s , H n (x) is linear with respect to the DCBs.Determining p( X i b,k |y k ,X i s,k ) for a fixed X i s is therefore a simple linear Gaussian problem, which can be solved with a Kalman filter.At time t 0 , or whenever a new DCB is added to the state, an initial guess X i b,0 is found by leveling to the current state of the ionosphere as in Equation 4, and an initial error covariance P b,0 = (2 TECU) 2 .As shown in Section 3.1, this is consistent with the actual errors produced by the initial estimation procedure, at least for an assimilation run which has already been initialized.The predicted bias xb,n and bias covariance P n are as follows, where Q b,n is a diagonal process noise covariance chosen to keep P b,n ≥ (0.05 TECU) 2 .This is considerably larger than the equivalent Q b,n = (0.001TECU) 2 used in Carrano et al. (2009).That study was conducted under tightly controlled conditions with well-known receivers, and a larger uncertainty is used in A-CHAIM to account for older or unknown hardware.For a random walk with a step size of 0.05 TECU every 5 min, the standard deviation at the end of a full day is σ day = (0.05 TECU) ̅̅̅̅̅̅̅ ̅ 288 √ ≈ 0.85 TECU.A-CHAIM expects 67% of receiver DCBs that have converged to a steady-state to stay within 0.85 TECU of their starting value over a single day.

Space Weather
The DCBs are simply added to the observations, so the measurement matrix H i n has a similarly simple form.
The Kalman gain K i n can therefore be evaluated with the observation error covariance R i n .This allows the calculation of the optimal estimator Xi b,n and posterior covariance Pi n .
If all P i 0 are initialized with the same values for each particle, these equations become simpler still.Neither the Kalman gain K i n , nor any of its constituent matrices have any dependence on y n or X n .As a result, K i n and Pi n will always be identical for every particle, and need only be calculated once per assimilation step.In the following expression the i superscript is dropped to reflect this feature.It is then straightforward to calculate Equation 28.
As A-CHAIM now has access to the optimal estimator Xi b,n , it would be inefficient to use any other choice of DCB.By setting X i b,n = Xi b,n , Equation 28 simplifies still further to Equation 29.As Pn is identical for all particles, so too is p( Recall from Equation 12 that only the normalized weights are used to calculate statistical moments.As Equation 29 is a constant, identical for all particles, Equation 19 can therefore be simplified to Equation 30.This is the same expression as Equation 10, the expression for the particle weights before the DCBs were added to the state space.By using the Rao-Blackwellized particle filter, choosing P 0 to be identical for every particle, and taking the optimal estimator Xi b,n , the DCBs can be factored out of the problem entirely.
It should be stressed that this complete factorization of the DCBs out of Equation 30 is not a typical result of Rao-Blackwellized particle filters, but occurs here due to the specific nature of the problem, and by somewhat careful choice of initial conditions for P 0 .This nearly trivial technique adds very little computational cost to the assimilation.An additional benefit is that any uncertainty in the DCBs does not contribute to the observation errors of the sTEC measurements.If the DCBs are fixed through some external model, every observation from that receiver has a covariant measurement error term from the uncertainty in the DCB correction.By moving the DCBs into the state, this difficult covariance problem disappears, as all errors associated with the DCBs are now contained in P n .Without Rao-Blackwellized particle filtering this benefit would be offset by the increased complexity of the increased size of the state space.With the technique above, both the DCB estimation problem and the correlated errors they create are simply factored out, for essentially no cost and few assumptions about the underlying behavior of the DCBs.

Results
There are several ways to assess the validity of the DCB estimation technique used in A-CHAIM.As A-CHAIM has already demonstrated success in reconstructing ionospheric electron density, it can be inferred that the DCBs are at least consistent with the modeled ionospheric state (Reid et al., 2023).This does not provide evidence that the modeled DCBs are consistent with real receiver behavior.The values of the A-CHAIM DCBs determined through the Rao-Blackwellized particle filter technique can be compared to those determined through other means, such as by using the IGS GIMs.For each GIM in Table 1, the daily DCB was determined for each receiver using the near-real-time data as in Equation 5.The difference between A-CHAIM DCBs and GIM DCBs, for all receivers, is shown in Figure 2. In Figure 2 and hereafter, the quantity Δ DCB is defined as Equation 31.If the A-CHAIM DCBs and the GIM DCBs were in agreement, it is expected that the mean Δ DCB would be near zero.
It might be expected, given the known biases between GIMs (Hernández-Pajares et al., 2009), that the A-CHAIM DCBs would show the best agreement with a GIM in the middle of the range.As this is not the case, this suggests some systematic effect causing the A-CHAIM DCBs to diverge from all GIM DCBs.The cause of this discrepancy must be determined, for which there are only four possibilities.The first is that the A-CHAIM DCBs are entirely inconsistent with expected DCB behavior, which is addressed in Section 3.1.The second possible source of this discrepancy is that A-CHAIM produces an ionospheric state which is systematically inconsistent with reality, and therefore produces DCBs that are inconsistent with the true values.As A-CHAIM has demonstrated success in ionospheric specification (Reid et al., 2023), the potential role of the plasmasphere is examined in Section 3.2.This discrepancy could also be a result of systematic underestimation of TEC, which is examined in Section 3.3.The third possible source of a systematic error is that the GIM-leveling procedure used to produce the comparison data set is biased.This possibility is examined in Section 3.4.The final possibility is the inverse of the second, that the GIM products are systematically inconsistent with the ionospheric state, and so produce DCBs which are inconsistent with reality.This will be examined in Section 3.5.

DCB Convergence Time and Stability
Unlike many other DCB estimation techniques, A-CHAIM does not require the DCB to be static over a full day.The A-CHAIM DCBs are free to migrate with every 5-min assimilation step.One advantage of this technique is that it allows for intraday variability (Coster et al., 2013;Zhang et al., 2019), but small-scale perturbations in the ionosphere could conceivably contaminate DCBs causing rapid fluctuations or other unrealistic behavior.
Another advantage of the A-CHAIM DCB estimation technique is that it is not necessary to keep an entire day of GNSS data to be able to produce a DCB estimate.Of course, if the convergence time for the DCBs is greater than a day, then this technique would introduce needless error.
Ideally, the DCBs should converge rapidly to some value, and then vary slowly thereafter.When a new GNSS station is added to the A-CHAIM data set, x 0 = DCB rcv is first estimated by leveling the vTEC to the current state of the A-CHAIM ionosphere, analogous to the procedure in Equation 5.As a full day of data is not available, only the sTEC for the current 5-min assimilation window is used.The error covariance is then initialized as P 0 = (2 TECU) 2 .When the model is operating normally, the initial estimate for the DCB is able to use the information provided by the A-CHAIM ionospheric state to produce a reasonable first estimate.This is not true when the entire model is initialized, as A-CHAIM will be identical to the background model E-CHAIM (Themens et al., 2017(Themens et al., , 2018(Themens et al., , 2019)), and could potentially have large errors in TEC at midlatitudes.The time history of each A-CHAIM DCB during the test period is plotted in Figure 3.The first plot of Figure 3 includes all receivers present when the model was initialized on August 20th.It is immediately apparent that a large population of DCBs were initialized with errors on the order of 10 TECU, suggesting that P 0 = (10 TECU) 2 might be more appropriate.The second plot of Figure 3 shows those stations which were added from August 21st onward, which do not exhibit the initialization errors.Most of the DCBs in this population were initialized within 2 TECU of their final values, justifying the choice of P 0 = (2 TECU) 2 , at least for stations added after the start of an assimilation run.In practice, the choice of P 0 controls how quickly the DCBs are allowed to converge, and choosing a slightly lower than optimal convergence rate will only affect the first few hours of an assimilation run.In operational use a complete reset is a rare event, and so the choice of P 0 = (2 TECU) 2 over some more complicated model has minimal impact on timescales beyond a day.
Figure 4 shows the time evolution of several example stations.Two test runs of A-CHAIM are plotted, one with and one without the Neustrelitz plasmasphere model (NPSM) (Jakowski & Hoque, 2018).These values are also plotted in Figure 4, along with the published DCBs from CAS in SINEX format.The first plot of Figure 4 shows the high latitude station IQAL, with an initial DCB close to the true value.The second plot shows the midlatitude station GODE, initialized over 10 TECU away from the true DCB, but converging rapidly.The large error covariance at initialization P o allows the DCBs to change rapidly early on.As the DCBs become more certain, and P n becomes small, the DCBs are less malleable.Large errors in the estimated DCB which occur during initialization are quickly corrected, but afterward the DCBs are expected to change slowly.This expectation is not always valid, as any changes to the receiver hardware at a ground station can have a dramatic effect on the DCB.
After an antenna or cable is swapped, for the purpose of DCB estimation the receiver should be considered an entirely new entity.With hundreds of GNSS receivers across a dozen networks such hardware changes are not rare, but difficult to detect reliably.Of the 662 receivers in this study, 15 had a DCB change of greater than 10 TECU.One such example, MAR6, is shown in the third plot of Figure 4.After the DCB change, the error in the MAR6 DCB was comparable to the initial error of the GODE DCB.Unlike GODE, which converged in a matter of hours, the MAR6 DCB took several days to converge to the new value.Hardware changes are also the principal mechanism which produce the extreme outliers in Figure 3, as the median DCB for a station with a hardware change is not representative of the DCB before or after.

Effects of Plasmaspheric Total Electron Content
If A-CHAIM were systematically underestimating TEC, this would create an overestimation of the DCBs in order to keep the assimilation self-consistent.One potential region of electron content that is neglected in A-CHAIM is the plasmasphere.A-CHAIM has a maximum altitude of 2,000 km.In normal A-CHAIM operation, the electron density above this altitude is assumed to be negligible.The plasmasphere is constrained to mid-and low latitudes, and so the vertically integrated electron density above 2,000 km is small.Any GPS observation where the ray path leaves the assimilation region through the southern boundary are excluded in A-CHAIM, and as a result most southward-pointing rays with significant plasmaspheric TEC are not assimilated.However, as shown in Figure 5, there are rays which terminate on the upper boundary of A-CHAIM which pass through much of the plasmasphere.
To quantify the impact of plasmaspheric plasma on bias estimation, test runs of A-CHAIM were conducted with and without the Neustrelitz plasmasphere model (NPSM) (Jakowski & Hoque, 2018).The inclusion of a plasmasphere had very little effect on the DCB estimation, as summarized in Figure 6.The effect of the plasmasphere was a bias of 0.04 TECU, which is in good agreement with the values in Stephens et al. (2011) which used upwardlooking sTEC measurements from the COSMIC satellites.This comparison shows that the effect of plasmaspheric TEC on A-CHAIM is at least an order of magnitude smaller than the difference between the A-CHAIM DCBs and the GIM-leveled DCBs.
Figure 7 shows the sTEC predicted by the NPSM model above 2,000 km altitude.Each plot corresponds to a hypothetical receiver at a different geographic latitude.The lack of measurements at high elevations is a result of the 55°inclination of the GPS satellites.The black curve shows the lower limits of the A-CHAIM assimilation region, at 45°magnetic latitude and 2,000 km altitude.Any ray outside of this region is not assimilated by A-CHAIM, which excludes nearly all plasmaspheric TEC at all latitudes.This essentially reproduces the situation in the simulation study of Lunt et al. (1999b), where it is shown that excluding low-latitude rays which pass through the plasmasphere is needed for accurate DCB determination.

Effects of Madrigal vTEC Measurements
An additional two A-CHAIM test runs were conducted using the Rao-Blackwellized DCB estimation procedure.
In Figure 2

Self-Consistency of GIM-Leveled DCBs
The GIM-leveling process is not without error, and some of this apparent underestimation could be an artifact of this technique.There is also some ambiguity in how best to interpolate the GIMs (Schaer et al., 2015).To ensure that the leveling technique used in this study is self-consistent, the GIM-leveled DCBs can be compared to the official IGS DCBs published in SINEX format.Several of the GIM products used in this analysis are distributed with receiver DCBs.As A-CHAIM uses the Chinese Academy of Sciences (CAS) satellite DCBs, the SINEX DCBs and the GIM-leveled DCBs should be self-consistent.
The CAS GIM-leveled DCBs are compared to the CAS SINEX DCBs in Figure 8.There is a small bias of − 0.7 TECU, and a standard deviation of 1.7 TECU, meaning the CAS GIM-leveled DCBs tend to slightly underestimate the CAS SINEX values.The variance is almost entirely attributable to noise.If it is assumed that most DCBs are static, then the standard deviation of the DCBs for each station is a reasonable proxy for the error.The CAS SINEX DCBs tend to be somewhat noisy, with the average station having an standard deviation of 1.1 TECU.The CAS GIM-leveled DCBs are slightly noisier, with a standard deviation of 1.3 TECU.Added in quadrature, this gives an expected error of 1.7 TECU for the difference summarized in Figure 8, exactly what is observed.
It is not obvious what the source of this − 0.7 TECU bias is, although from an operational perspective the precise mechanism is unimportant.It is noteworthy that there is a non-zero mean, and that reconstructing the DCBs from a GIM product can introduce a bias.If all GIM products in this study gained a similar bias from the leveling process, then a small amount of the offset between A-CHAIM DCBs and the GIM-leveled DCBs could be due to this leveling error.This is insufficient to explain all of the observed differences, and suggests that the sTEC processing technique used in A-CHAIM is not biased by more than 1 TECU from an IGS standard.

Effect of GIM-Leveled DCBs on A-CHAIM
Rather than using the Rao-Blackwellized DCB estimation procedure, test runs of A-CHAIM can be conducted using the GIM-leveled DCBs found for each day during the August 20th through 10 October 2022 test period.The details of each test run are summarized in Table 2.If the GIM-leveled DCBs are more accurate, they should produce an improvement in reconstructed electron density.For this assimilation experiment a test run of A-CHAIM was conducted for each GIM in Table 1, using the GIM-leveled DCBs rather than solving for them with the Rao-Blackwellization method.Each of these new runs also included the NPSM plasmasphere, to isolate the effect of the imposed DCBs.

Space Weather
10.1029/2023SW003611 shows the same, but from the A-CHAIM (+vTEC) run which assimilated the Madrigal vTEC measurements.The third plot shows the observed sTEC corrected by the CAS IONEX-derived DCB, and the vTEC from the A-CHAIM test run using these DCBs.The vTEC from the CAS IONEX GIM is also plotted, as a black dashed line.The fourth plot shows the same, but using the JPL IONEX GIM.
The A-CHAIM vTEC in each of the plots in Figure 9 changes.In the first plot, the A-CHAIM vTEC is slightly, but consistently, lower than the Madrigal vTEC.Once corrected by the A-CHAIM DCBs, the resulting sTEC occasionally dips below zero during period of very low ionospheric density.This is to be expected, as the phaseleveling procedure used to convert GNSS observables to sTEC is subject to error, typically on the order of 1.5 TECU to 3 TECU per continuous arc (Reid et al., 2023).In the second plot the A-CHAIM vTEC overlaps the Madrigal vTEC, as this run assimilated these vTEC values as measurements.Each of these runs used the Rao-Blackwellized DCB estimation procedure, and the A-CHAIM state and DCBs were produced as a selfconsistent set.The sTEC obs − DCB shown in the third and fourth plots are the data assimilated by A-CHAIM, as the DCBs were fixed with the external GIM.To stay self-consistent, the A-CHAIM state must change to match the observed sTEC, and therefore shifts to match the IONEX vTEC.
The electron density produced by each of the A-CHAIM test runs can be compared to measurements.Three kinds of measurements are used in this analysis, autoscaled ionosonde NmF1 and NmF2, and in situ electron density measurements from the Defense Meteorological Satellite Program (DMSP).The ionosonde measurements were assimilated by each of the test runs, and so each run should show good agreement.The in situ electron density is not assimilated, and provides an independent reference.The results for each run are summarized in Table 3.
Most test runs showed comparable performance when compared to NmF1, with the exception of the +vTEC run.This run significantly overestimated NmF1, with an RMSE more than twice as large as any other run.The UPC run also showed a slight overestimation of NmF1 when compared to ESA and COD, runs with a similar ΔDCB.
For the other GIM-derived runs, there is a slight trend for runs with more negative ΔDCB to underestimate NmF1.
When examining NmF2, the +vTEC run is again an outlier, and to a lesser extent the UPC run.Otherwise, there is no obvious trend in mean error of NmF2 with DCB.The mean error is small compared to the Root Mean Square Error (RMSE) of any run other than +vTEC, so any bias is a small contributor to the overall error.All runs have a small RMSE when compared to typical values of NmF2, which are on the order of ×10 11 to ×10 12 m − 3 .The differences between runs is modest, as expected, as these NmF2 observations were assimilated into the model.A stronger trend is present in the DMSP in situ data.As the DCBs decrease, the mean error increases, indicating a tendency to overestimate the topside electron density.Both A-CHAIM runs tend to slightly underestimate the topside density, with the Madrigal GIM run having almost zero bias.The two runs which incorporate the Madrigal vTEC measurements have the smallest RMSE overall.For the other runs using GIM DCBs, the RMSE increases steadily with the mean error, as this bias is a major component to the error in the reconstruction.For in situ measurements, the A-CHAIM runs which fit self-consistent DCBs have the smallest RMSE, along with the run using DCBs derived from the Madrigal vTEC maps.
For a more detailed view, Figure 10 shows the mean errors for each test run, for both the assimilated NmF2 measurements and the DMSP in situ electron density.Rather than averaging over all available data, the observations and DCBs were binned into groups by latitude and longitude before being averaged.This highlights the effects receiver DCBs have on their local area.The left plot shows the assimilated NmF1 measurements, with the slight negative trend with decreasing ΔDCB visible, along with the outliers +vTEC and UPC.In NmF2, most runs show no bias with ΔDCB, other than the outliers +vTEC and UPC.
The third plot of Figure 10 shows the error in DMSP in situ electron density.A very clear linear relationship exists between the DCBs and the in situ error.After limiting consideration to only those regions with GNSS receivers, the A-CHAIM runs without Madrigal vTEC data show no overall bias.The Madrigal DCBs now show a slight overestimation, and all other GIMs showing significant overestimation.Again, the +vTEC and UPC runs are outliers, showing better topside performance than runs with similar ΔDCBs, at the cost of decreased NmF1 and NmF2 performance.This is evidence of a bifurcation in A-CHAIM, with one attractor preferring to compensate

Space Weather
10.1029/2023SW003611 for underestimated DCBs by increasing NmF2, and the other by increasing the topside thickness.This behavior emerged naturally in the case of UPC, but was selected specifically by choosing to force the topside of A-CHAIM to match the Madrigal vTEC in +vTEC HTop.
The differences between each of the test runs can easily be seen when the results are mapped as in Figure 11.Each test run is shown with three maps in a row, showing the DCB offset from the A-CHAIM (NPSM) DCBs, the mean percent error in NmF2, and the mean percent error in DMSP in situ electron density.The effects on NmF2 are predictably small, although when expressed as percent error all runs tend to overestimate NmF2, a departure from the mixed behavior seen in Figure 10.This is consistent with a small positive bias in NmF2, which would result in a strong overestimation as a percent of NmF2 at night when densities are small, but be negligible compared to other sources of error during the day.As expected, the UPC run overestimates NmF2 to a greater extent than the others, particularly at high latitudes.+vTEC also shows significant overestimation in NmF2, but at the lower latitude ionosondes.
The variations in DMSP in situ electron density errors are much more dramatic, with the JPL and EMR products saturating the color scale.All GIM runs show their worst overestimation in the lower latitude band over North America, which is well populated by ionosondes and has the greatest density of GNSS receivers.The best performance appears to be in the circle around the polar cap, though this is not due to any effect of the GNSS receivers.This region is not readily observed by GNSS sTEC due to orbital geometry, and has enhanced densities which are not well captured by E-CHAIM.As such, A-CHAIM tends to strongly underestimate this region, as seen in Reid et al. (2023).Any regional overestimation of topside thickness will coincidentally act to correct this error.Both A-CHAIM runs tend to slightly underestimate the topside in the Russian sector, and in the polar cap.
As the Madrigal GIM run slightly overestimated topside thickness, it had a better mean error of − 3.8 × 10 8 m − 3 than the A-CHAIM run at − 1.9 × 10 9 m − 3 , albeit with a worse RMSE of 9.2 × 10 9 m − 3 compared to 8.9 × 10 9 m − 3 (Table 3).This underestimation in the polar cap is also why the two runs which assimilated Madrigal vTEC had the best DMSP RMSE overall, but this apparent improvement disappeared when the results were binned into regions where receivers were actually present.In the regions where rays pass through the ionosphere, the run with the lowest overall RMSE +vTEC HTop overestimates the topside electron density.The UPC run here maintains its outlier status, with relatively good topside performance outside of the midlatitude American sector.The regions where UPC NmF2 performance is worst tend to be where in situ measurements were best.
When Madrigal vTEC measurements are included in the assimilation, the resulting DCBs are in good agreement with the GIM-derived DCBs from Madrigal.When Madrigal vTEC is assimilated with no other constraints, as in the +vTEC test run, the resulting errors in NmF1 and NmF2 were significant.When the topside thickness was adjusted in +vTEC HTop, the resulting errors were very similar to the run which used the Madrigal-derived DCBs.This can be clearly seen in Figure 10.The +vTEC HTop run overestimated midlatitude topside density

Space Weather
10.1029/2023SW003611 to a greater extent than the Madrigal GIM run, but also produced slightly more negative ΔDCBs.In the areas where measurements are present, including the Madrigal vTEC measurements did not improve overall performance.To match Madrigal vTEC, either midlatitude NmF1 and NmF2 were overestimated, or midlatitude topside electron density was overestimated.

Discussion
In this study, the A-CHAIM run which includes the NPSM plasmasphere performs slightly better than the run without the NPSM in the topside, and slightly worse in NmF2.It is unclear if this is due to chance, as particle filters are a Monte Carlo technique and do not produce perfectly repeatable results.The operational A-CHAIM model does not currently include a plasmasphere model, and more study is needed to determine if one needs to be included.The relative contribution of the plasmasphere to TEC will vary seasonally, and with solar and geomagnetic activity, and so a long-term validation is required.As shown in Figure 1, a small error in DCB can result in a modest error in NmF2, or a significant error in topside electron density if NmF2 is fixed.In the idealized example, a DCB error of − 5TECU results in a nearly 100% overestimation in electron density at 800 km altitude.As NmF2, hmF2, NmF1 and hmF1 measurements are assimilated in A-CHAIM, the bottomside profile should be well constrained in regions where there are many ionosondes.The in situ measurements from the DMSP satellites, which orbit at approximately 800 km, allow us to recreate Figure 1 at a large scale in Figure 10.All of the color scales are identical to those in Figure 1, to allow for straightforward comparison.By choosing the A-CHAIM (NPSM) DCBs as the reference value, they are acting as the true DCBs for this analysis.As can clearly be seen in Figure 10, when a DCB is biased relative to the A-CHAIM DCBs, it causes an error in the electron density.This error in the reconstruction is consistent with the expected result if that DCB were biased relative to the truth by the same amount.An offset from the A-CHAIM DCBs is indistinguishable from an error.
The ΔDCBs for each GIM product show spatial structure, and not a simple flat offset as might be expected.If a feature appears in every GIM DCB map, it may be due to some odd behavior of A-CHAIM, or it may be due to some shared assumption in the GIMs, such as the thin-shell approximation.While the GIM-leveled DCBs used every available sTEC measurement, A-CHAIM can only use those rays which stay fully inside the assimilation region, meaning receivers near the lower boundary cannot use southward rays.Several of the ΔDCB maps, COD, CAS, and IGS show a strong negative band at the southern limits of the assimilation region.However this band is absent in Europe, being limited to the American sector, and does not appear in other products such as UPC.In fact, some products show a positive band at the lower limits in Europe, including Madrigal, UPC, and EMR, while others such as JPL or ESA show no such banding.This band structure is therefore not likely to be a feature of the A-CHAIM DCBs.
Figure 12 shows the ΔDCB of each GIM, binned by latitude.This curve is then decomposed into the mean ΔDCB, and the residual variation ΔDCB − ΔDCB.A-CHAIM (NPSM) is not compared to itself, and is omitted from the figure.All GIM-derived DCBs show a downward trend below 45°in latitude.This effect is smallest in EMR, UPC, and Madrigal, with a drop ∼1 TECU at the lowest latitudes, whereas the other GIM show a more pronounced drop of ∼2 − 3 TECU.ESA, COD, CAS, IGS, and JPL all show nearly identical curves.This curve is similar both in shape and in magnitude to that found in Figure 11 of Mazzella (2012), indicating that this latitudinal effect is likely due to aliasing of plasmaspheric TEC into the midlatitude ionosphere in the GIMs.As seen in Figure 7, stations at approximately 45°latitude and below begin to have a significant plasmaspheric contribution to sTEC.Notably, the two GIMs which do not use a thin-shell approximation, Madrigal and UPC, show the least evidence of plasmaspheric influence.Madrigal uses two orders of magnitude more GNSS receivers than any other GIM product in this study, and has a altitude-dependent error term, (Vierinen et al., 2016) which may help limit the influence of low-elevation equatorward rays with significant plasmaspheric TEC.The good latitudinal agreement of the EMR GIM, when compared to other single thin-shell GIMs, may be attributable to the comparatively large number of high-latitude GNSS stations used in that product (Ghoddousi-Fard, 2014).The southward rays of these additional high-latitude receivers may provide a greater constraint to the midlatitude vTEC.The IGS GIM, as a weighted average of the COD, ESA, JPL, and UPC products, appears to have inherited the strong plasmaspheric effects found in three of the four constituents.Of all of the IGS Analysis Centers in this study, UPC has consistently been an outlier, which may be attributable to its unique use of voxels.It has the best latitudinal agreement with the A-CHAIM DCBs, comparable to those of the Madrigal-derived DCBs.
For the GIM runs, as the ΔDCBs became more negative the runs tended to underestimate NmF1.This is counterintuitive, as the system was simultaneously boosting the topside electron density to match the underestimated DCBs.A-CHAIM was hollowing out the bottomside electron density profile to insert it in to the topside, which would have the effect of raising the apparent shell height of the ionosphere.Figure 13 shows an electron density profile measured with the Millstone Hill ISR.The ISR is co-located with the MHJ45 Digisonde, one of the instruments assimilated into A-CHAIM, and so should have a very well constrained bottomside ionosphere.This ISR profile had a close conjunction with the DMSP F17 satellite, which passed approximately 5°t o the east.The A-CHAIM profile for each test run is also plotted.When examining the topside, at the altitude of DMSP, the differences between each test run are immediately apparent.The four test runs which used the Rao-Blackwellized DCB estimation method, and the test run which used the Madrigal-derived DCBs, are all clustered together, and in excellent agreement with the density measured in situ.All other runs, those derived from IONEX GIMs, significantly overestimate the topside density.This overestimation extends down even to the region

Space Weather
10.1029/2023SW003611 immediately above the F2 peak.Below the F2 peak, the runs which overestimated the topside are now underestimating the bottomside, as both the bottomside thickness is reduced and hmF2 is increased.Each run underestimated the E-region density, but the test runs which overestimated the topside had a more severe underestimation of the E-region.The test runs which best matched the DMSP in situ density also show the best agreement with the ISR profile.

Conclusion
The Rao-Blackwellized particle filter technique used in A-CHAIM is able to produce stable estimates of receiver DCBs.After a period of initial convergence, lasting less than a day, these DCBs generally evolve slowly to follow the real drifts in the receiver.When hardware changes cause a large change in DCB, A-CHAIM is able to adjust within a few days.This is an area for potential improvement, by re-initializing the DCB for a receiver if some error threshold is detected.By running tests both with and without the NPSM plasmasphere model, the effect of the plasmasphere on DCBs is estimated to be on the order of 0.05 TECU at high latitudes, in agreement with previous estimates (Stephens et al., 2011).One additional advantage of this technique is that it only requires concurrent access to the sTEC of a single assimilation step, rather than to a full day of observations.This makes it suitable for real-time applications.
The DCBs produced by A-CHAIM show systematic biases relative to the DCBs estimated using GIMs, resulting in differences ranging from − 1 TECU to − 7 TECU.Some portion of this bias may be attributable to the GIMleveling process itself, as the CAS SINEX DCBs showed a bias of − 0.7 TECU relative to the CAS GIM-leveled DCBs.This does show that the A-CHAIM data processing pipeline produces sTEC data which is able to match the CAS DCBs using the CAS GIM to within 1 TECU, and therefore phase-leveling is not a source of significant bias.When constrained to match the vTEC of the Madrigal vTEC maps, A-CHAIM is able to reproduce the Madrigalderived DCBs within 0.4 TECU.When the latitudinal differences between the A-CHAIM DCBs and the GIMderived DCBs were examined, the results were in excellent agreement with previous studies of plasmaspheric effects on midlatitude vTEC estimation (Anghel et al., 2009;Carrano et al., 2009;Lunt et al., 1999b;Mazzella, 2012).This effect was most pronounced in GIMs which use a thin-shell representation of the ionosphere.The remaining differences between GIMs were comparable with the relative biases between GIMs seen in previous studies (Hernández-Pajares et al., 2009;Li et al., 2015;Mazzella, 2012).DCBs generated using the Madrigal vTEC maps produced the closest agreement, often within 1 TECU of the A-CHAIM DCBs outside of the auroral oval.
Eight additional test runs of A-CHAIM were performed, each using the DCBs from a different GIM product, to assess the impact these would have on the assimilation.All test runs had good performance when compared to the assimilated autoscaled NmF2, which indicated the assimilation was able to stay self-consistent.The two test

Space Weather
10.1029/2023SW003611 runs of A-CHAIM with Rao-Blackwellized DCB estimation, with and without the NPSM, had the smallest RMSE of any runs when compared to in situ electron density measurements from the DMSP satellites.In test runs which used GIM-derived DCBs, the performance of the model in reconstructing the topside electron density deteriorated dramatically.In the regions with GNSS receivers, any divergence from the A-CHAIM DCBs resulted in overestimation of topside electron density.The resulting overestimation of topside electron density is consistent with the GIM-derived DCBs being underestimates of the true receiver DCBs.This suggests that the GIMs are overestimating high-latitude electron content, both through plasmaspheric effects at midlatitudes and overall.The typical shell height of 450 km specified by IONEX maps is also a likely source of error, as these test runs forced A-CHAIM to underestimate the bottomside electron density even while overestimating topside density.
Any mapping from sTEC to vTEC, even those with a vertical parameterization of the ionospheric profile, neglects the effects of gradients in the ionosphere (Vierinen et al., 2016).vTEC is a weighted spatial average of the ionosphere along the line of sight of the GNSS receiver.If strong or persistent gradients are present, this spatial average is not necessarily representative of any particular point along line of sight.This is a fundamental limitation of vTEC, and a particular challenge at high latitudes where orbital geometry requires all rays to be southward.Lower latitudes will have greater electron density, and this will be averaged in to the high latitude vTEC.It is inevitable that high latitude vTEC tends to be overestimated.
The DCB estimation technique presented here is not subject to the limitations of vTEC.A-CHAIM assimilates sTEC directly by linearly integrating through a 3D ionosphere, and therefore does not create the spatial averaging of vTEC-based techniques.The A-CHAIM ionospheric profile is updated continuously by measurements, including independent data sources like ionosondes which help constrain the F2-layer peak, as well as the bottomside thickness (Reid et al., 2023).By comparing with the topside electron density measured in situ, the entire electron density profile is well characterized.This includes the plasmasphere, through the inclusion of the NPSM.
If the error in the A-CHAIM DCBs is much greater than 1 TECU, then A-CHAIM would need to be significantly, systematically, and globally underestimating some part of the ionospheric electron density which is not already measured.
Satellite Data Processing and Distribution (2020).Altimeter data from the Sentinel-3 satellites are provided by the ESA Copernicus Data Space Ecosystem (European Space Agency, 2023).Real time solar proton flux provided by the NOAA SWPC (National Oceanic and Atmospheric Administration (NOAA) & National Aeronautics and Space Administration (NASA), 2022).Auroral electron energy and flux measurements provided by JHU/APL (Paxton, 2022).In situ measurements from the DMSP missions are provided by National Centers for Environmental Information (NCEI) (2022a).Millstone Hill ISR data provided by Phil Erickson, MIT/Haystack Observatory (Erickson, 2022).Categorical color palettes from the coloropt library (Tsitsulin, 2021).

Figure 2 .
Figure2.A comparison of DCB estimation techniques during the August 20th through 10 October 2022 test period.The DCBs from the A-CHAIM run with NPSM plasmasphere were subtracted from other estimates, for all unique stations and all times.The edges of each box represent the 25th and 75th percentiles, with the outer limits at the 5th and 95th percentiles.

Figure 3 .
Figure 3. Superimposed epoch analysis of the DCB convergence time in A-CHAIM.The first 7 days of DCB estimates for each receiver are superimposed, based on when that receiver first appeared in the data set.The median value of the DCBs for each receiver from day 2 through day 14 has been subtracted, showing the overall convergence envelope.

Figure 4 .
Figure 4. Time series of GNSS DCBs for three stations, estimated by various techniques during the August 20th through 10 October 2022 test period.The top figure shows the DCBs for the IQAL station, which shows typical behavior for a high-latitude station.The second figure shows GODE, which was initialized with a large error of nearly 15 TECU, which converged in less than a day.The third figure shows MAR6, which had a large change in DCB around September 2nd.MAR6 did not provide data to A-CHAIM from September 2nd through September 8th.

Figure 5 .Figure 6 .
Figure 5. A-CHAIM ionospheric and NPSM plasmaspheric electron density along a slice of fixed longitude at 60°E.The edges of the A-CHAIM domain are indicated with a solid black line.The solid white line shows the possible positions of a GPS satellite, as determined by orbital geometry.The dashed white lines show the line-of-sight of a GNSS receiver at the extreme southern boundary of A-CHAIM, and at 0°elevation at the geographic North pole.
the run labeled A-CHAIM (+vTEC) assimilated Madrigal vTEC as a normal measurement, and the run A-CHAIM (+vTEC HTop) attempted to preferentially adjust the A-CHAIM topside thickness H Top to match the Madrigal data.It is clear from Figure2that the assimilation of the Madrigal vTEC as a measurement caused a change in the Rao-Blackwellized DCBs.The mean ΔDCB for A-CHAIM (+vTEC) was − 1.59 TECU, and for A-CHAIM (+vTEC HTop) the mean ΔDCB was − 1.51 TECU.These are in excellent agreement with each other, and very similar to the ΔDCB of − 1.18 TECU from the Madrigal-derived DCBs.These sets of DCBs were produced through entirely different means, one set by performing the standard GIM-leveling procedure with a thin-shell ionosphere, and the other two within the A-CHAIM model.By constraining the A-CHAIM profile to match the Madrigal vTEC, the resulting Rao-Blackwellized DCBs are within 0.4 TECU of the Madrigal-derived DCBs.

Figure 7 .
Figure 7. Sky maps of plasmaspheric sTEC predicted by the NPSM model from 2,000 to 22,000 km.The black region shows the projection of the A-CHAIM boundary, only rays within the encircled region are assimilated into A-CHAIM.The gray circle shows the rays above 30°elevation.

Figure 9
Figure9shows the observed sTEC at the ALBH GNSS receiver from 1 September 2022, through 3 September 2022.The top plot shows the sTEC corrected by the DCBs from A-CHAIM, along with the vTEC produced by vertically integrating A-CHAIM at the station location.The vTEC from Madrigal is also shown.The second plot

Figure 9 .
Figure9.Observed sTEC and modeled vTEC at the ALBH GNSS receiver for a period of 3 days from 1 September 2022 through 3 September 2022.Each plot shows the vTEC from the corresponding A-CHAIM test run, and the IONEX GIM used to produce the DCBs (if applicable).All vTEC values are for the single point at the station coordinates.The Madrigal vTEC is shown on each plot for reference.

Figure 10 .
Figure10.Effects of GNSS Receiver DCBs on reconstructed electron density.Ionosonde NmF1 and NmF2 was the same data assimilated into each A-CHAIM run.DMSP uses combined in situ measurements from F-16, F-17, and F-18.Each observation type and DCB offset were binned by 2°latitude and 4°longitude and averaged over the August 20th to 10 October 2022 test period.Each small circle corresponds to a bin with at least one electron density measurement and GNSS receiver.The mean value for each run is indicated with a large outlined circle.

Figure 11 .
Figure 11.Global comparisons of each of the A-CHAIM test runs, averaged over the whole test period from August 20th through October 10th.The left map for each run shows a marker for each GNSS receiver, colored according to the offset of the DCB used from the A-CHAIM DCB.The center map shows each assimilated ionosonde, colored with the mean percent error in reconstructed NmF2.The right map shows the mean percent error in DMSP in situ electron density measurements.All colors are on the same scale as in Figure 1.

Figure 12 .
Figure 12.Differences between A-CHAIM DCBs and GIM-derived DCBs as a function of latitude, binned by 1°of geographic latitude.Error bars represent the standard deviation of all binned DCBs.The dotted line shows the mean ΔDCB of all receivers.The solid line shows the binned ΔDCB with the mean subtracted, to highlight the latitudinal variations.

Figure 13 .
Figure 13.Ionospheric electron density profile measured by the Millstone Hill ISR at 11:26:02 UT on 29 August 2022.This profile had a near conjunction with the DMSP F17 satellite, being displaced by ∼5°in longitude.The A-CHAIM profile for each of the test runs in this study are plotted, showing the effects of the different DCB estimation techniques on both the bottomside and topside thickness.

Table 1
Vertical TEC Products Used to Test A-CHAIMs DCB Determination Nrcv indicates the average number of GNSS stations reported in each daily file over the 20 August 2022 through 10 October 2022 test period.All files were obtained from CDDIS in the IONEX format, except for the Madrigal vTEC maps, which were obtained from the CEDAR Madrigal database in HDF5 format.

Table 3
Summary of Different DCB Estimation Techniques, and Their Effects on Reconstructed Electron Density When Used inA-CHAIM DCBs referenced to the A-CHAIM (NPSM) run, averaged over all receivers.All DCB values are in TECU, all other values in ×10 10 m − 3 .