Can We Intercalibrate Satellite Measurements by Means of Data Assimilation? An Attempt on LEO Satellites

Low Earth Orbit satellites offer extensive data of the radiation belt region, but utilizing these observations is challenging due to potential contamination and difficulty of intercalibration with spacecraft measurements at Highly Elliptic Orbit that can observe all equatorial pitch‐angles. This study introduces a new intercalibration method for satellite measurements of energetic electrons in the radiation belts using a Data assimilation (DA) approach. We demonstrate our technique by intercalibrating the electron flux measurements of the National Oceanic and Atmospheric Administration (NOAA) Polar‐orbiting Operational Environmental Satellites (POES) NOAA‐15,‐16,‐17,‐18,‐19, and MetOp‐02 against Van Allen Probes observations from October 2012 to September 2013. We use a reanalysis of the radiation belts obtained by assimilating Van Allen Probes and Geostationary Operational Environmental Satellites observations into 3‐D Versatile Electron Radiation Belt (VERB‐3D) code simulations via a standard Kalman filter. We compare the reanalysis to the POES data set and estimate the flux ratios at each time, location, and energy. From these ratios, we derive energy and L* dependent recalibration coefficients. To validate our results, we analyze on‐orbit conjunctions between POES and Van Allen Probes. The conjunction recalibration coefficients and the data‐assimilative estimated coefficients show strong agreement, indicating that the differences between POES and Van Allen Probes observations remain within a factor of two. Additionally, the use of DA allows for improved statistics, as the possible comparisons are increased 10‐fold. Data‐assimilative intercalibration of satellite observations is an efficient approach that enables intercalibration of large data sets using short periods of data.

In situ multi-spacecraft measurements are crucial for studying near-Earth radiation.These measurements provide the foundation for validating existing physics-based models of various particle populations, improving our understanding of the underlying physics, and creating more accurate models.Statistical parametrization of the most energetic magnetospheric regions enables the planning of multi-year satellite missions, particularly at MEO and HEO orbits (Friedel et al., 2005).Furthermore, recent studies on Data Assimilation (DA) (Castillo et al., 2021;Cervantes et al., 2020) and assimilative real-time radiation belt forecasting leverage large data sets from multiple spacecraft.However, the quality and reliability of multi-source observations can be affected by several factors.Differences in instrumentation performance or design, lack or degradation of detector shielding, non-standardized instrument calibration (e.g., Cayton & Tuszewski, 2005), and differences in satellite location can all result in significant deviations between measurements from multiple spacecraft.Thus, even observations from similar orbits and magnetospheric regions can vary significantly and require proper intercalibration between the different instruments.
Traditionally, satellite data intercalibrations are performed using satellite conjunctions, which involve comparing real data in magnetic coordinates (e.g., Friedel et al., 2005;Szabó-Roberts et al., 2021;C. Wang et al., 2013) or matching phase space density (PSD) in adiabatic space (e.g., Chen et al., 2005Chen et al., , 2007;;Ni et al., 2011;Zhu et al., 2022).Both approaches require a benchmark instrument (a "gold standard" [GS] as by Friedel et al. (2005)) that provides high-fidelity data and is used to intercalibrate measurements from other instruments.A conjunction between different satellites is defined by imposing strict spatial and temporal criteria on the observations to ensure that physical constraints are met.Then, statistical analysis of the residuals from data comparisons is performed, and scaling factors can be estimated.Although satellite conjunctions have demonstrated reliable results and are an established methodology for satellite data cross-calibration, the strict constraints imposed on the data to make them comparable greatly reduce the number of observations that qualify as a conjunction.This leads to poor statistics and requires large amounts of data.These issues are particularly exacerbated when comparing satellites at very different orbits that observe vastly different magnetospheric regions and particle populations (s.a., LEO vs. HEO, see Figure 1).In such cases, a spacecraft with extensive L-coverage should be used as a reference for intercalibration (Friedel et al., 2005).
To address some of the limitations of data cross-calibration via conjunctions, it would be useful to have an approximation of the state of the entire radiation belts.DA techniques, for example, the Kalman filter (KF) (Kalman, 1960), the Extended Kalman filter (EKF) (Jazwinski, 1970), or the Ensemble Kalman filter (EnKF) (Evensen, 2003), have been utilized in the space weather community since the 2000s to estimate the optimal state of this region using satellite observations and physics-based models (e.g., Bourdarie & Maget, 2012;Drozdov et al., 2023;Godinez & Koller, 2012;Koller et al., 2007;Kondrashov et al., 2011;Naehr & Toffoletto, 2005;Ni et al., 2009;Reeves et al., 2012;Schiller et al., 2012;Shprits et al., 2007Shprits et al., , 2012)).The resulting reconstruction of the system (a time-dependent 3-D PSD volume) is referred to as a data-assimilative reanalysis and represents the state of the radiation belts system that is statistically closest to the "true state."Reanalyses have been used in the past to study the dynamic behavior of the system and to identify missing processes in physics-based models (Cervantes, Shprits, Aseev, & Allison, 2020;Koller et al., 2007;Kondrashov et al., 2007;Schiller et al., 2017).
In this study, we elaborate on an idea proposed by Shprits et al. (2007), and present a new satellite intercalibration method based on the modeling of the outer radiation belt by means of DA.We test our novel intercalibration technique by cross-calibrating six satellites of the NOAA-POES fleet against Van Allen Probes (used here as the reference data set).To do so, a 1 year reanalysis of the radiation belts using Van Allen Probes and GOES data is estimated.By flying the six NOAA-POES satellites through the reanalysis, we can perform on-orbit data comparisons at each POES location, and consequently conduct a statistical analysis of the residuals to estimate the recalibration coefficients.In order to validate our approach, a traditional conjunction study between Van Allen Probes and POES is also carried out.Comparison between the cross-calibration coefficients estimated with both methodologies is presented.
In the next Section, we describe the proposed method.In Section 3, we present the Van Allen Probes and reanalysis data sets.Utilized POES observations and their necessary processing is described in Section 4. Section 5 deals with the POES fly through the data assimilative reanalysis and the statistical analysis of the related on-orbit comparisons.In Section 6, we present the statistical analysis of the comparisons from the conjunction study.General results, final cross-calibration factors and discussion are offered in Section 7, followed by the conclusions and outlook in Section 8.

Rationale and Methodology
For lab-calibration procedures, the instrument is exposed to a radioactive source with a well-known spectrum (or signal) and then the measurement is compared to the expected signal.In the case of satellite observations such a procedure is not feasible, because lab recreation of the space conditions is not possible.The problem, however, would be solved, if one could have an approximation of the space environment (the radiation source), in which the non-calibrated spacecraft (NS) is immersed.In this case, having the entire state of the radiation belt system or at least an approximation of it would allow us to easily compare observations, thus avoiding the limitations tied to conjunction cross-calibrations.
Data assimilation techniques enable us to estimate such a state-approximation by blending physics-based models and satellite observations in an optimal way.The information contained in the satellite data will propagate to other areas of the modeling space, giving us a time dependent global reconstruction of the system that is statistically closest to the true state of the system, a so-called reanalysis (RA).Once this reconstruction has been estimated, we can fly satellites/instruments at different orbits through it and compare the real observations (j NS ) with the state-estimate (j RA ) at all locations, energies and equatorial pitch-angles.The idea is to find factors η, such that for each time, location and energy of the instrument it holds: We rename η as R DA , the flux ratio between reanalysis and observations.Note that, R DA may be influenced by a variety of factors, such as geomagnetic activity (or K p ), energy (E), and even location (L*) and equatorial pitch-angle (α eq ).However, the extent to which these factors contribute to R DA can only be assessed through a statistical analysis of all the resulting ratios.
The step-by-step procedure can be summarized as follows: 1. Choose a reference data set to be used as the GS.Ideally, the GS data is pitch-angle resolved, has high energy resolution, provides large L*-coverage and observes the most dynamic regions of the radiation belts (i.e., satellites at HEO would be most suitable here), as this will reflect in the quality of the RA.EnKF), and estimate the RA of the radiation belts for the desired period of time. 5. Convert RA into electron fluxes in observational space.6. Process and constrain NS observations if necessary (e.g., for LEO satellites the use of trapped electron data is greatly important.More details in Section 4.1).7. Fly NS satellite through the RA.This is equivalent to an interpolation of GS-data into the grid of NS satellite.8. Estimate the ratios R DA at each NS-time and -location.9. Perform statistical analysis of R DA in dependence of L*, E, α eq , and K p to determine the most important parameters influencing the ratios R DA .10. Estimate recalibration coefficients and their uncertainties in dependence of parameters found in the previous step.For this use suitable statistical measures depending on the shape of the obtained distributions, for example, statistical mean We validate our approach by presenting a comparison of the recalibration factors obtained through a traditional geomagnetic conjunction study.

Reference Data Set and Reanalysis Data
For this study, we choose the instruments onboard Van Allen Probes as our reference GS data set, and use these observations together with those from GOES 13 and 15 to estimate a data assimilative reanalysis of the radiation belt region for the period of October 2012 to September 2013.A comparison between the POES and Van Allen Probes data sets, and the Van Allen Probes + GOES reanalysis is displayed in Figure 2. Simple visual inspection of the figure clearly shows the need for these data sets to be intercalibrated.An overview of these data sets is given in this section.
In this study, we used MagEIS measurements from probes A and B averaged over 30 min.An example of the Van Allen Probes data set used in this work is presented in Figure 2b for fixed energy (∼1 MeV) and α eq < 15°.
The GOES fleet are a series of meteorological geostationary satellites operated by the U.S. NOAA at nearly geosynchronous orbit (Data Book GOES, 2005).Each GOES spacecraft hosts a Magnetospheric Electron Detectors (MAGED) and two Energetic Proton, Electron, and Alpha Detectors (EPEAD).MAGED consists of nine solid-state-detector telescopes, five in the east-west (equatorial) plane and the other four in the north-south (meridional) plane, measuring differential electron fluxes at energies of: 0.03-0.05,0.05-0.1,0.1-0.2,0.2-0.35,and 0.35-0.6MeV (Hanser, 2011;Rodriguez, 2014a).In addition, the EPEADs measure MeV electron and proton integral flux data in two energy ranges: >0.8 and >2 MeV.To perform the data assimilative reanalysis, we use MAGED and EPEAD pitch-angle resolved electron flux measurements from GOES 13 and 15.The observations are averaged over 30 min.EPEAD integral fluxes and pitch-angles are obtained by averaging the measurements of the East and West telescopes (Rodriguez, 2014b).Integral fluxes as a function of energy are fitted to a power law in order to extend up to 1 MeV energies.We use the 90° pitch-angle differential flux data from MAGED and fit the two integral channels of EPEAD to an exponential function to obtain differential flux at the interpolated energies.

Reanalysis Data Using Van Allen Probes and GOES
In this study, we estimate a data assimilative reanalysis of the outer radiation belt for the period of October 2012 until September 2013 following Cervantes et al. (2020).We assimilate the observations of Van Allen Probes (probes A and B), as well as GOES-13 and GOES-15 into the VERB-3D code (Shprits et al., 2009;Subbotin & Shprits, 2009) using a 3D split-operator KF with a timestep of model and assimilation of 1 hr.Assimilation of GOES data allows us to update the outer boundary of the reanalysis using realistic conditions, and thereby improving the overall accuracy of the reconstruction.More details about the implementation of the splitng-operator KF can be found in Shprits et al. (2013) and Shprits, Castillo, et al. (2023).In order to assimilate flux measurements, these need to be converted to PSD in coordinates of phase space (L*, μ, K).To calculate μ, in situ magnetic field measurements from Van Allen Probes are used.For the calculation of K and L*, we use the magnetic field model T89 (Tsyganenko, 1989) and IRBEM-ONERA library (Boscher et al., 2022).Differential fluxes (j) are converted to PSD (f) in units of (c/cm/MeV) 3 following Rossi and Olbert (1970) by f = j/p 2 .
The VERB-3D code computes the numerical solution of the bounce-averaged Fokker-Planck-equation (Shprits et al., 2008;Subbotin et al., 2010) using a fully implicit finite differences method on a high resolution grid with (29 × 101 × 91) points for (L* × E × α eq ), respectively.VERB-simulations include radial, energy and pitch-angle diffusion, as well as losses to the magnetopause.The radial diffusion coefficient is calculated after Brautigam and Albert (2000) in terms of L* and used by the VERB-code for all K p values.The plasmapause position is calculated after Carpenter and Anderson (1992).The bounce-averaged diffusion coefficients for hiss and dayside and nightside chorus waves are computed using the Full Diffusion Code (Shprits & Ni, 2009), and with the parameterizations provided by Orlova et al. (2014), and Orlova and Shprits (2014), respectively.The range of L* reaches values from 1 to 6.6 and for equatorial pitch angles from 0.7° to 89.3°.The energy at the outer radial boundary (L* = 6.6) is defined in the range of 0.01-10 MeV.At the low energy boundary, the energy varies in dependence of the L* value, because electrons are energized during their transport to lower L-shells (e.g., Subbotin & Shprits, 2009), and correspond to μ ≈ 9 MeV/G for electrons at α eq = 90°.The magnetopause position is estimated as the Last Closed Drift Shell calculated with the IRBEM-ONERA library using the T89 model.For further details about the reanalysis, the boundary and initial conditions, we refer the reader to the work by Cervantes et al. (2020).
The resulting assimilated state of the radiation belts is then a time-dependent three-dimensional PSD volume.In order to compare this state to POES measurements, we convert the assimilative reanalysis to differential flux in the coordinates of the observational space (L*, E, α eq ) by f = j/p 2 .A fragment of the electron fluxes from the reanalysis data set used in this study is displayed in Figure 2c for fixed energy and equatorial pitch-angles α eq < 15°.

POES Data Set
Our goal is to test our new intercalibration approach to intercalibrate electron flux data from six satellites of the POES fleet, that is, MetOp2, NOAA−15, 16, 17, 18, 19 (an overview is given in Table 1).In this study, we focus on the observations over the time period of 01 October 2012 until 30 September 2013.
The particle flux data set provided by the POES fleet has gained particular importance due to its large temporal coverage, extensive L*-distribution, and short orbital period.These spacecraft are in Sun-synchronous LEO at about 850 km altitude and have an orbital period of ∼100 min.Since the launch of NOAA-15, the fleet carries the Space Environment Monitor (SEM-2) instrument package (Evans & Greer, 2000), which contains the Medium Energy Proton and Electron Detector (MEPED), and the Total Energy Detector.The SEM-2 MEPED instrument consists of eight particle detector systems: two proton solid-state detector telescopes (each ± 15° wide), two electron solid-state detector telescopes (each ± 15° wide) and four omni-directional (dome) proton detector systems.The electron/proton telescopes are mounted with different orientation in order to observe different particle populations: (a) the 0°-telescope has the central axis of its field of view rotated 9° in the XZ plane pointing away from the local zenith, (b) the 90°-telescope is oriented almost perpendicular to the 0°-telescope with the central axis of its field of view rotated 9° in the YZ plane pointing away from the antiram direction.Original SEM-2 MEPED electron data are reported in three integral electron channels (E1, E2, E3) with a nominal energy range of 0.03-2.5,0.1-2.5, and 0.3-2.5 MeV, respectively (Evans & Greer, 2000;Peck et al., 2015).MEPED count rates (counts/s) are reported in 16 s intervals (Codrescu et al., 1997).
POES observations have been reported to suffer from a number of issues that make their use rather challenging.The rotation angles of the telescopes allow for a clear field of view and for monitoring a mixture of particle populations.Thus, the 0°-telescopes observe mostly particles in the atmospheric loss cone (LC) and only at the geomagnetic equator trapped populations are measured, while the 90°-telescopes monitor trapped particles at high latitudes and L* > 1.4 (Evans & Greer, 2000).Additionally, Rodger, Clilverd, et al. (2010) documented proton contamination of the SEM-2 MEPED electron data, as the detectors respond to protons with energies of up to 2.7 MeV.The amount of contamination varies for each electron energy channel (Yando et al., 2011), but electron data from the 90°-telescopes are of good quality with only 3.5% (on average) to 7% (disturbed times) contamination occurring beyond L = 7. Radiation damage, due to long-term exposure, may also affect the electron detectors, but its impact on the measurements is expected to be rather negligible (Asikainen & Mursula, 2011;Galand & Evans, 2000;McFadden et al., 2007).
In order to address some of the issues mentioned in the previous paragraph, we use the corrected differential electron fluxes estimated by Peck et al. (2015).The authors reduced proton contamination of the MEPED E1 to E3 electron channels.Additionally, using the information about relativistic electrons embedded in the observations of both P6 proton detectors (integral proton channel (P6) with a nominal energy range of 0.03 MeV to > 6.9 MeV), the authors produced a virtual fourth electron channel (E4) with energies between 0.3 MeV and −2.5 MeV, centered at ∼0.612 MeV (Green, 2013).The count rates estimated for the E1-E4 electron energy channels were then used to calculate continuous spectra over the energy range from 0.025 to 10 MeV (total of 27 energy channels).Peck-corrected MEPED data set also contains error estimates accounting for measurement errors and for errors in the fitting of the spectral distributions.An example of the electron fluxes measured by MEPED onboard NOAA-16 used in this study are displayed in Figure 2a for ∼1 MeV energy and equatorial pitch-angles α eq < 15°.

Processing of POES Observations
For a proper comparison of the Van Allen Probes and POES data sets some considerations need to be taken into account, and consequently further processing and/or constraining of the observations has to be performed.(Asikainen & Mursula, 2011;Lam et al., 2010).

Table 1 National Oceanic and Atmospheric Administration Polar-Orbiting Operational Environmental Satellites Satellites Used in This Study and Their Characteristics
All POES data are processed with the IRBEM-ONERA library using the magnetic field model (T89) (Tsyganenko, 1989).We first constrain the POES data to observations at equatorial pitch-angles α eq ≥ 6° because the smallest pitch-angle channel of MagEIS can detect α eq ∼ 6° based on the center point.Only time intervals of quiet to low geomagnetic activity are used (i.e., times when K p ≤ 4 − ) to reduce possible inaccuracies of the magnetic field model.Additionally, we restrict the L*-range to values between 3 and 6.6 R E , as we want to focus on observations of the outer radiation belt.Figure 3 presents the L* and α eq -distributions of the raw (Figures 3a and 3c) and the constrained (histograms b and d) data sets.The final overlap of the distributions for the constrained data suggests that comparison of Van Allen Probes and POES observations for the studied time period is only feasible for L* = 3-5 R E and α eq = 6°-12°.
As previously mentioned, the POES-fleet observes a mixture of electron populations, therefore we only use measurements from the 90°-telescopes.Since these observations are very close to the LC, we need to isolate the measured populations and remove drift-and bounce LC (DLC and BLC, respectively) measurements from our data sets.The purpose of this step is twofold: (a) DLC and BLC observations from POES cannot be compared to Van Allen Probes measurements because Van Allen Probes does not resolve the LC and (b) the use of only trapped particles allows us to rely on Liouville's theorem to map PSD at the geomagnetic equator.
The approach used to isolate POES populations used in this work is similar to the one presented by Shprits, Michaelis, et al. (2023), and is described in the next paragraphs.Measurements of the MEPED detector for each energy channel are reported as the total counts per second estimated over 8 consecutive integration periods of 2s.Due to the wide angle of aperture of the detector and the integration time for the measurement, a large range of electrons with local pitch-angles between α loc ± 15° can enter the detector, so that the measurement of the central angle may be biased.For this reason, using the local pitch-angle from the central-angle measurement α c = α loc , we estimate the other two possible edge values for the local pitch-angle at satellite position (assuming a symmetric detector opening), that is, α min = α loc − 15° and α max = α loc + 15°.Using the conservation of the first adiabatic invariant (μ), we can calculate the corresponding magnetic field intensity at the mirror point for each of these pitch-angle values, (i.e., B c , B min , B max , respectively) using IRBEM-ONERA library.For the characterization, we only use the minimum of the three values (here notated as B M = min(B c , B min , B max )), thereby imposing the strongest assumption to ensure that measurements labeled as trapped are accurate.However, an unambiguous characterization of the observed electron populations is rather impossible.The intensity of the Earth's magnetic field at 100 km altitude (B foot ) is estimated using the IGRF-12 model (Thébault et al., 2015).
We then determine if a particle precipitates into the atmosphere or not, as follows: • The BLC is defined as the range of pitch-angles at satellite location with mirror points below the atmosphere in either hemisphere.These particles will precipitate into the atmosphere within one bounce period.For each measurement, we find the minimum B foot value between both hemispheres and compare this value to B M .It holds: if B foot ≤ B M , the particle bounces below the atmosphere and will be lost, therefore the measurement is labeled as BLC.• The DLC is defined as the range of α loc at fixed drift-shell, that reach altitudes lower than ∼100 km at the South Atlantic Anomaly (SAA) and will therefore precipitate into the atmosphere within one drift period.We estimate the L-shell (McIlwain value) for each POES measurement using IGRF.We then find the minimum B foot for the given L-shell along constant longitude (longitude of satellite location).This is the magnetic field intensity at the SAA (B SAA ) and we compare it to B M .It holds: if B SAA ≤ B M , the particle drifts below 100 km at the SAA and it will be lost, therefore the measurement is labeled as DLC.
• If the measurement is not labeled as BLC nor as DLC, it will be labeled as TRAPPED.Only these data are used for the present work.
The obtained geographical distributions of the electron populations agree well with those obtained by Rodger, Carson, et al. (2010) (see Figure 4).Only trapped data are used for the comparison with Van Allen Probes + GOES-reanalysis, Van Allen Probes observations, and for the respective estimation of recalibration coefficients.

POES Fly Through Across the Reanalysis
In this and the following sections, we present the formal tests and results of our intercalibration approach on the NOAA-16 satellite data set.The results obtained for the other satellite missions mentioned in Table 1 are summarized in Supporting Information S1.
Since the reanalysis represents the "optimal state" of the outer radiation belt (i.e., the closest to the true state) at all times and locations, we can fly each POES satellite through this global reconstruction.A spacecraft fly through across the data assimilative reanalysis is equivalent to an interpolation of the assimilated electron fluxes onto the spatial/temporal-grid of the POES fluxes.For the fly through, POES data are binned into 1h time bins (i.e., the time step of the reanalysis) and the (L*, E, α eq )-nodes in the VERB-grid closest to the satellite measurement are labeled.To obtain the flux value of the reanalysis at the satellite location, we perform three 1D interpolations using piece-wise cubic splines.We interpolate electron fluxes over 1D intervals enclosing the measured POESdata point and at least five RA grid nodes around it.Since the VERB-code only models diffusion of energetic particles trapped in the radiation belts without convection, we focus on radiation belt energies from ∼0.2-1 MeV (i.e., energy channels 10 to 17 of the Peck-corrected data).
We then extract the corresponding flux values of the reanalysis (j RA ) at POES location (L*), energy (E) and pitch-angle (α eq ), and compare them with the actual flux values measured by the LEO satellites (j POES ) at same location, energy and pitch-angle.Figure 5 shows the 2D-histogram of L* and α eq values, at which fly through data are available.We find a total of 23,664 data points available for comparison in the ranges of L* = 3.2-4.4and α eq = 6°-12°.
Since we now have two flux values at same location, we can estimate the flux ratios (R DA ) between the reanalysis fluxes (j RA ) and the measured fluxes (j POES ) for each time-bin (reanalysis time (t RA )), satellite location (L*), energy channel (E) and equatorial pitch-angle (α eq ), as follows: (,  * , , ) = (,  * , , ) (,  * , , ) (1) We analyze the distributions of R DA in dependence of E, α eq , L*, and K p , in order to determine the influence of each of these parameters on the flux ratios.The histograms of R DA in dependence of the energy channel are presented in Figure 6.The distributions show slightly skewed bell shapes with clear peaks.The spread and skewness of the distributions appears to be larger for E ≤ 0.5 MeV.We estimate the median of R DA over time for each energy channel E i (red line), that is, and use the Median Absolute Deviation (MAD) (green lines) to estimate the median variation of the residuals around the median of the distribution.For skewed distributions the MAD is more robust than the standard deviation, because it is more resilient to outliers, and it is defined as the median of the absolute deviations from the median of the data, as follows (Rousseeuw & Croux, 1993): (2) The median of R DA for energies <0.7 MeV remains close to 1 (note that the x-axis is log 10 (R DA )), but at higher energies it shows a clear increase up to values of ∼2 for E = 0.973 MeV.The lower MAD values constantly fall around 0.8-0.9, but noticeably increase above 1 for E = 0.779 MeV and E = 0.973 MeV.For most energy channels, the upper bounds of the MAD oscillate around 2-3, reaching highest values (>4) at E = 0.779 and E = 0.973 MeV.These features suggest a strong dependence of the R DA on the energy channel.We further study the dependence of R DA on α eq for each energy channel, as shown in the 2D-histograms in Figure 7I.The red dashed line represents the median and the magenta dashed lines are the MAD of the distributions (note that the y-axis is log 10 (R DA )).
Here, the skewness and spread of the distributions also appear to decrease with increasing energy.Clusters in the data can be well seen for all energy channels at least up to α eq = 9°, with highest sample density around α eq = 6°-7°.The median of the distributions seems to decrease with increasing value of α eq in a non-linear way at all energies.For E < 0.3 MeV, the median of R DA moves from values close to ∼1 at α eq = 6° down to ∼0.6 at α eq = 11°.Furthermore, for E > 300 and E < 0.7 MeV, the median of R DA also peaks around 1 at α eq = 6°, but it reaches down to ∼0.2 at α eq = 11°.Higher energy channels show larger values for the median of R DA with the maximum being >2 at α eq = 6° and the minimum falling close to 1 at α eq = 11°.For all the energy channels, the upper limit of the MAD remains around 0.3 above the median, while the lower bound decreases rapidly with increasing value of α eq , so that it cannot be estimated for α eq = 11° in most of the cases.
Similar trends in the skewness and spread are observed in Figure 7II, which displays the 2D-histograms R DA versus L* for each energy channel.These distributions also show clear data bulks between L* = 3.2-4.2with peaks at L* = 3.6-4.0for all energies.The median curves of R DA present inverse parabolic behavior that seems to flatten at E = 0.973 MeV.The median reaches its minimum at L* = 3 and increases within one order of magnitude until it finds its maximum at L* = 4 and then begins to decrease at L* = 4.2.The median at L* = 4 oscillates close to 1 for E < 0.6 MeV, but increases its value above 2 at higher energies.The trends in the MAD are similar to those seen in Figure 7I, which is expected due to the inverse proportionality of L* and α eq .
Finally, we analyze the variation of R DA with respect to the geomagnetic activity index K p (see Figure 7III).The same trends in the skewness and spread with regard to the energies observed before, are also seen here.However, in this case the spread of the distributions appears to be less than one order of magnitude.The histograms show clear bulks of samples between K p = 0-3.Unlike the previous cases, the median of R DA does not show much variation and oscillates around 1 for all K p values and E < 0.7 MeV.At higher energies, the median curve also increases its values slightly showing a small peak at K p ∼ 0.3, but remaining rather constant otherwise.The MAD shows larger uncertainties in the upper limits around the median, but remains within 0.4 of the median values.
The spread of the MAD also decreases noticeably with increasing energy.
The analysis of R DA presented in this section suggests a strong dependence on the energy channel, L*-location and α eq .In contrast, the value of K p shows a rather small, if not negligible, influence on the flux-ratios.Before we look  deeper into these parameters and their influence on R DA , we check in the next section if a traditional conjunction approach delivers similar insights into the behavior of the flux-ratios.

Conjunction Study Between Van Allen Probes and NOAA-16
In this section, we analyze the behavior of flux-ratios obtained from a geomagnetic conjunction study performed between the NOAA-16 and Van Allen Probes (A and B) satellites.In this case, we choose Van Allen Probes observations to be the "gold standard," which we use as a reference to carry out on-orbit comparisons with NOAA-16 measurements in geomagnetic space (Friedel et al., 2005).For a pair of (Van Allen Probes, NOAA-16) observations to be considered a conjunction, the following conditions should be met: (a) The location of both satellites must be within ±0.1 L*, (b) ideally the observed electrons have the same equatorial pitch-angles: ±0.5°α eq , (c) the energy of the measurements has a maximum deviation of ±10%: E VAP = E POES ± 10%, (d) the conjunction must occur within a time frame of Δt = ±1 hr, and (e) the conjunction occurs during low to moderate levels of geomagnetic activity: K p ≤ 4 − .
Figure 8 presents the 2D-histogram of L* and α eq values, at which the geomagnetic conjunctions are found.We have a total of 1,129 conjunctions between Van Allen Probe-A and NOAA-16 (Figure 8a) and, 1,131 conjunctions between Van Allen Probe-B and NOAA-16 (Figure 8b), in the ranges of L* = 3.6-4.4and α eq = 6°-8°.Bins with the largest number of data points are centered around L* = 3.8 and α eq = 8.5°.
Since we now have comparable pairs of (Van Allen Probes, NOAA-16) observations, we can perform flux-comparisons at same satellite location and estimate the flux ratios (here notated as R Conj ) between Van Allen Probes measured fluxes (j VAP ) and POES measured fluxes (j POES ) for each time-bin (Van Allen Probes time (t VAP )), satellite location (L*), energy channel (E) and equatorial pitch-angle (α eq ), as follows: (3) Similar to the previous section, we analyze the statistical dependence of R Conj on E, L*, α eq , and K p . Figure 9 shows the histograms of R Conj per energy channel.Since the distributions are rather irregular and show large spread, we estimate their peak as the median of R Conj over time for each energy (i.e., Q 2 (R Conj (E i )) = median(R Conj (E i ))) (indicated by the red bar); and their deviation through the MAD (green lines) is estimated by: For energies <0.7 MeV, the value of the median remains close to 1, but it increases for higher energies reaching a maximum at E = 0.973 MeV.While the upper bound of the MAD seems to stick constantly close to the median for all energies, the lower bound becomes too small for several energies and cannot, therefore, be displayed in log 10 scale.However, a clear peak in sample density is observed at L* = 3.6-4.For E < 0.7 MeV, the value of the median of R Conj seems to remain constantly around 1 or increases with increasing L* value, showing a peak at the L* = 4.0 bin and then decreasing again.MAD values for the upper bound remain around 2 units above the median, but the lower limit becomes too small for L* < 3.8 in most energy channels.
Additionally, Figure 10III shows the 2D-histograms of R Conj in dependence of the geomagnetic index K p .While a clear peak in sample density can be observed at K p = 0, the distributions show large spread and for K p > 1 no clear peak can be seen.The median value at the bulk of the samples is very close to 1 for all energy channels.However, the curve of the median oscillates in rather random way at higher K p values, so no clear trend can be observed.While the upper bound of the MAD closely follows the median value, the lower MAD limit becomes too small for the log-scale.

Results and Discussion
Taking into account the statistical analyses presented in Sections 5 and 6, here we compare the median values of R DA and R Conj (denoted by Q 2 (R DA ) and Q 2 (R Conj ), respectively), in dependence of E, L*, α eq , and K p (i.e., the red lines in the previous histograms).We discuss our findings and estimate final intercalibration coefficients for NOAA-16.The uncertainties of both data sets are quite large to the upper limits of the median.Lower bound uncertainties remain small, but they do increase for Q 2 (R DA ) at E > 0.6 MeV.
The upper limit uncertainty for Q 2 (R DA ) remains around a factor of ∼2 for E < 0.7 MeV, but increases up to a factor of ∼3.5 for higher energy channels.
The upper bound uncertainties of Q 2 (R Conj ) are generally larger than those of Q 2 (R DA ), but remain within a factor of ∼2-2.5.
We further study the behavior of the median values of R DA and R Conj with respect to L*, α eq , and K p for each energy channel.The median values of R DA and R Conj in dependence of α eq and energy channel are presented in Figures 12c and 12d, respectively.R DA -median curves clearly resemble the trends observed in Figure 12a.In general, Q 2 (R DA ) increases with decreasing value of α eq for fixed energy.For E < 0.6 MeV, most Q 2 (R DA ) values are below 1, indicating that POES measurements tend to be larger than the reanalysis.MEPED fluxes at α eq = 6°-7° appear to be very close to the reanalysis fluxes below 0.7 MeV.The largest difference between the data assimilative output and the POES measurements is observed above E = 0.7 MeV.Trends of Q 2 (R Conj ) in Figure 12d coincide well with those in Figure 12b.For E < 0.4 MeV, Van Allen Probes measurements at α eq = 7° are higher than MEPED fluxes, but at α eq = 6°, the opposite is the case.For E > 0.6 MeV and at α eq = 6°-7°, MEPED fluxes underestimate Van Allen Probes observations.The curves of Q 2 (R DA ) and Q 2 (R Conj ) in dependence of K p and energy channel are displayed in Figures 12e and 12f, respectively).For most energy channels, Q 2 (R DA )-curves are equal or close to 1.At E > 0.6 MeV, we observe an increase in the median value, suggesting that POES underestimates Van Allen Probes fluxes at these energies.Q 2 (R Conj )-values move close to 1 only for E < 0.6 MeV.At E > 0.6 MeV, an increase in Q 2 (R Conj )-values up to a factor of 2 is well observed.
With increasing K p -value the statistical significance of the K p -bins is strongly reduced (i.e., points per bin ≤10), which resembles in the irregular behavior of the curves.Therefore, we only plot the results for K p ≤ 1.

Discussion
The comparisons presented in the previous sections clearly show how the data-assimilative method is able to compare more data points (Figure 5) than the conjunction study (Figure 8), thereby consistently improving the statistics for the intercalibration.This is because the reanalysis provides a global reconstruction of the entire space of the radiation belts, allowing us to compare most of the real observations at all satellite locations.In Figures 6  and 7, an increase in spread and skewness of the R DA distributions below E < 0.5 MeV is well observed.This is not the case for R Conj (Figures 9 and 10).The reason for this may lay in the physics used by the VERB-3D code, which as a diffusion model is more suitable to model energetic particles.The comparison in Figure 11 clearly shows the potential of our data-assimilative intercalibration approach.Differences between Q 2 (R DA ) and Q 2 (R Conj ) may be related to the very different statistics of both data sets.All conjunction statistics contain less data points than the statistics of the data-assimilative method.Another possibility is a bias coming from the way the on-orbit comparisons are estimated.By just comparing the observations in space and time, we neglect the dependence of the instrument's response on the hardness of the real energy spectrum.For instance, if due to a loss process, low energy particles are removed from the environment, the net energy of the spectrum will increase.The observed dependence of Q 2 (R DA ) and Q 2 (R Conj ) on L* and α eq (Figures 12a-12d, respectively) further supports this hypothesis.Such a dependence was also reported by Peck et al. (2015), in comparison with the data set from the Detection of Electro-Magnetic Emissions Transmitted from Earthquake Regions satellite Instrument for Detecting Particles (Sauvaud et al., 2006).Since the original energy channels of the POES measurements were derived as integral fluxes over broad ranges of energy, the effect of spectrum hardening could be particularly high on the effective energy of the POES data set.Restriction of the K p values to ≤4 − may help reduce the effect of hardening, however, the large width of the real energy channels, the large field of view of the detector and possible remaining contamination can cause the observations to be dominated by higher energy particles.
Values of Q 2 (R DA ) are similar for all L* and α eq .This is not the case for Q 2 (R Conj ), where values in dependence of L* show higher maxima than those in dependence of α eq .The proximity of POES pitch-angle measurements to the LC may be the reason for this result.In the VERB-code the LC is modeled for a dipole field using a pitch-angle dependent exponential decay, which does not account for the energy dependent losses at 100 km altitude (Marshall & Bortnik, 2018).Additionally, classification of observations as trapped contains unavoidable inaccuracies.On the other hand, Van Allen Probes observations in the smallest pitch-angle channel are also very close to the LC, such that measurements from these channels probably contain LC particles, even though the central angle of the instrument may be outside of the LC.Optimal performance of the KF requires highly accurate modeling of the system as well as accurate knowledge of the covariances.Unfortunately, for dynamic systems that cannot be reproduced in a laboratory and for which only sparse observations are available (e.g., the radiation belts), such accuracy levels cannot be achieved and errors may be present in the reanalysis.Since we interpolate the reanalysis onto the POES grid, errors or inaccuracies contained in the reanalysis will also appear in the fly through data and become noticeable after comparison with the real data set.Therefore, the use of a data-assimilative intercalibration approach also enables us to learn about possible improvements in the physics of our model.In general, Q 2 (R DA ) values in L* and α eq are lower than those of Q 2 (R Conj ).This is potentially an indication of inaccuracies in the latitudinal dependencies of the used diffusions coefficient of the VERB-3D code, which determine the shape of the pitch-angle distribution.In the future, more advanced diffusion coefficients such as Drozdov et al. (2017), D. Wang et al. (2019), andSaikin et al. (2022) may deliver better agreement.
The analysis on the K p dependence of Q 2 (R DA ) and Q 2 (R Conj ) presents large increases in both curves for E > 0.6 MeV.Since the last integral channel (E4) of the original SEM-2 data is centered at about 0.612 MeV, we find that this is an indication of a possible bias in the Peck-corrected differential fluxes, perhaps related to the spectral fit.While this data product delivers large amounts of observations and the possibility to work with higher energies, the broad width of the energy channels of the original POES data set may impose some limitations to extensions of the observations to higher energies.
Our results show that the highest dependence of Q 2 (R) is on energy, L* and α eq .Since in Figure 12 the inverse relation between L* and α eq is easily observed, for the purpose of this study, we present final recalibration coefficients (values of Q 2 (R DA )) only in dependence of energy and L* in tabular form (see Table 2).

Conclusions
In the present study, we have shown the potential of a data-assimilative satellite intercalibration approach.The proposed method was tested and validated using measurements of energetic electrons in the radiation belt region   , and Van Allen Probes.Using our intercalibration approach, we are able to considerably improve the statistics of on-orbit data comparisons.Satellite intercalibration via data assimilative fly through requires therefore shorter periods of data than comparisons through conjunctions.Our comparative analysis clearly shows that due to very few conjunctions, flux-ratios may be falsely estimated, while using data-assimilative intercalibration shows that Peck-corrected POES data are already in good agreement with Van Allen Probes observations below E ≈ 0.6 MeV (i.e.R DA ≈ 1), and can be used to reconstruct the global state of the radiation belts.For higher energy channels the data sets are within a factor of 2, so that intercalibration is required, as shown by both methods in this study.The recalibration factors estimated with our data-assimilative method are consistent with the results from the conjunction study.
The results of this study are encouraging as large satellite data sets can be efficiently and automatically intercalibrated with this technique.In particular, the intercalibration of POES against Van Allen Probes allows us to use the recalibrated POES data set for DA studies using the assumption that both data sets have similar uncertainties.
In the future, we plan to extend the pitch-angle distribution of the Peck-corrected POES data sets using Smirnov et al. (2022) approach and perform global reconstruction of the radiation belts using our recalibrated data set.We also want to perform a similar analysis using original uncorrected SEM-2 integral fluxes, including lower ring current energies.In this study, we have excluded such a comparison since it would only concern one energy channel for radiation belt energies.Additionally, we plan to perform assimilation of the intercalibrated POES data set using our full-3D EnKF (Castillo et al., 2021) to explore how much the quantification of the model error and model bias can be improved by introducing these parameters into an augmented state vector.We also look forward to using this intercalibration method with other satellite fleets providing large data sets, such as GPS.

Figure 3 .
Figure 3. Data distributions: L* and equatorial pitch-angle (α eq ) observed by Van Allen Probes (a, b) and Peck-corrected data of NOAA-16 for 01 October 2012 until 30 September 2013.(a, c) L* and α eq -distributions of raw data, respectively.(b, d) L* and α eq -distributions of constrained data sets for intercalibration.

Figure 4 .
Figure 4. Global distribution of electron populations in the radiation belts as observed by the 90°-telescopes Medium Energy Proton and Electron Detector onboard NOAA-16, as of Peck-corrected Space Environment Monitor data.DLC = Drift loss cone and BLC = Bounce loss cone.

Figure 5 .
Figure 5. Fly through data: 2D-histogram of L* versus α eq covered by the fly through of NOAA-16 for the period of October 2012 until September 2013.A total of 23,664 data points are available, color-coded is the number of data points per bin.

Figure 6 .
Figure 6.Distribution of R DA in energy for NOAA-16: Histograms of R DA (in log 10 scale) versus number of samples for each energy channel (each R DA unit is divided into 10 bins).The median is indicated by the red lines, while the Median Absolute Deviation is given by the green lines.

Figure 7 .
Figure 7. 2D-Distributions of R DA for NOAA-16: (I) 2D-Histograms of R DA (in log 10 scale) versus α eq for each energy channel (plotted in 1°-bins and R DA -bins of 1.4 width).(II) 2D-Histograms of R DA (in log 10 scale) versus L* for each energy channel (plotted in L*-bins with 0.25R E width).(III) 2D-Histograms of R DA (in log 10 scale) versus K p for each energy channel (plotted in K p -bins of 0.33 width).Color-coded are the number of samples.The median is indicated by the red dashed lines, and the Median Absolute Deviation is given by the magenta dashed lines.

Figure 8 .
Figure 8. Conjunction data: 2D-histogram of L* versus α eq , at which geomagnetic conjunctions between NOAA-16 and (a) Van Allen Probe-A; and (b) Van Allen Probe-B are available for the period of October 2012 until September 2013.The total number of conjunctions is displayed in the lower left part of each plot, color-coded is the number of data points per bin.

Figure
Figure 10I presents the 2D-histograms of R Conj in dependence of α eq .Although the distributions show high spread and nonuniform behavior, a clear peak can be seen between α eq = 8°-9°.The median of R Conj (red dashed line) appears to remain constant around a value of 1 for energies below 0.7 MeV, showing a decrease in value at the 6° bin.At higher energies the value of the median increases, as also observed in the previous figure.The MAD bounds (magenta dashed lines) indicate higher deviation to the upper values of R Conj .Figure 10II displays the 2D-histograms of R Conj in dependence of L*.The distributions are again rather irregular and show large spread.

Figure 9 .
Figure 9. Distribution of R Conj in energy for NOAA-16: Histograms of R Conj (in log 10 scale) versus number of samples for each energy channel (each R Conj unit is divided into 10 bins).The median is indicated by the red lines, while the Median Absolute Deviation is given by the magenta lines.

Figure 10 .
Figure 10.2D-Distributions of R Conj for NOAA-16: (I) 2D-Histograms of R Conj (in log 10 scale) versus α eq for each energy channel (plotted in 1°-bins and R Conj -bins of 1,4 width).(II) 2D-Histograms of R Conj (in log 10 scale) versus L* for each energy channel (plotted in L*-bins with 0.25R E width).(III) 2D-Histograms of R Conj (in log 10 scale) versus K p for each energy channel (plotted in K p -bins of 0.33 width).Color-coded are the number of samples.The median is indicated by the red dashed lines, and the Median Absolute Deviation is given by the magenta dashed lines.
Figures 12a and 12b show the median of R DA and R Conj (respectively) in terms of L* and energy.Q 2 (R DA ) curves are smooth and present similar trends as those seen in Figure 11 for all L*-bins.The values of Q 2 (R DA ) increase with increasing L*-value for fixed energy, but remain between ∼0.5 and ∼1.5 at 0.2 MeV, and reach ∼1.2-2.5 at 0.973 MeV.At L* ≤ 3.8 and for E < 0.6 MeV, MEPED slightly underestimates the reanalysis fluxes.For E ≥ 0.8 MeV, this underestimation is seen in all L*-bins and also maximum values of Q 2 (R DA ) are observed here.Below L* = 3.6 and for E < 0.6 MeV MEPED consistently overestimates the reanalysis fluxes.The curves of R Conj -median values (Figure 12b) are less smooth than those of Q 2 (R DA ), and no clear trends are observed.For most L* and energy values, MEPED underestimates Van Allen Probes fluxes, only at L* = 4 below 0.5 MeV mild overestimation or agreement are observed.Highest Q 2 (R Conj ) values are at E > 0.6 MeV for most L*-values.

Figure 11 .
Figure 11.Values of Q 2 (R DA ) and Q 2 (R Conj ) versus Energy.Plotted in linear scale are the median values of R DA (pink diamonds) and R Conj (red diamonds) estimated for NOAA-16 in dependence of the energy channel.Error bars are estimated from the corresponding Median Absolute Deviation values and displayed for Q 2 (R DA ) in pink color and for Q 2 (R Conj ) in blue.

Figure 12 .
Figure 12.Values of Q 2 (R DA ) and Q 2 (R Conj ) versus L*-bin, α eq and K p .Curves of Q 2 (R DA ) in dependence of the energy channel for NOAA-16, color-coded are the curves (a) for each L*-bin, (c) for each α eq -bin, (e) for each K p -bin.Curves of Q 2 (R Conj ) in dependence of the energy channel for NOAA-16, color-coded are the curves (b) for each L*-bin, (d) for each α eq -bin, (e) for each K p -bin.The Y-axes in all plots is in linear scale.
Columns are satellite name, altitude, inclination angle, local time of the ascending node (LTAN), and the intervals of the data used in this study Note.The coefficients are given in terms of energy and L*.

Table 2
Recalibration Coefficients for NOAA-16: Final Intercalibration Coefficients (Q 2 (R DA )) Estimated for NOAA-16 Using Our New Data Assimilation Approach from POES satellites (NOAA