-
PDF
- Split View
-
Views
-
Cite
Cite
Dan Aksim, Dmitry Pavlov, Improving the solar wind density model used in processing of spacecraft ranging observations, Monthly Notices of the Royal Astronomical Society, Volume 514, Issue 3, August 2022, Pages 3191–3201, https://doi.org/10.1093/mnras/stac1229
- Share Icon Share
ABSTRACT
Solar wind plasma as a cause of radio signal delay plays an important role in solar and planetary science. Early experiments studying the distribution of electrons near the Sun from spacecraft ranging measurements were designed so that the radio signal was passing close to the Sun. At present, processing of spacecraft tracking observations serves a different goal: precise (at metre level) determination of orbits of planets, most importantly Mars. The solar wind adds a time-varying delay to those observations, which is, in this case, unwanted and must be subtracted prior to putting the data into the planetary solution. Present planetary ephemerides calculate the delay assuming a symmetric stationary power-law model for the solar wind density. The present work, based on a custom variant of the EPM lunar–planetary ephemeris, questions the accuracy and correctness of that assumption and examines alternative models based on in situ solar wind density data provided by OMNI and on the ENLIL numerical model of the solar wind.
1 INTRODUCTION
Two-frequency observations of the Viking spacecraft in 1976 (lander and orbiter) were once used for analysis of radio signal propagation delays due to solar wind plasma (Callahan 1977; Muhleman & Anderson 1981). There were also plans for dual-frequency ranging for the Mars Express spacecraft (Patzold et al. 2004; Remus et al. 2014), but no information is known about the outcome. Apart from those, the spacecraft tracking measurements are made on single frequency; thus, determining the delay in the solar wind directly from the measurements themselves is not possible. While, in principle, other radio observations could help, such as VLBI (Aksim et al. 2019) or observations of scintillations or radio sources (Shishov et al. 2016), in practice there are too few of them, and even fewer, if any, on the dates of spacecraft tracking measurements. Three options remain for accounting for the delay in the solar wind when determining planetary orbits for ephemerides:
Using a simplified model of the solar wind electron density distribution, and then determining its single free parameter (a reference electron density) from the spacecraft tracking observations.
Using in situ electron density measurements taken by spacecraft such as the Advanced Composition Explorer (ACE) and Wind to account for large-scale temporal variations of the said electron density in interplanetary space, while keeping a simplified symmetric model for the electron density distribution at any moment.
Using a more complicated dynamical model of the solar wind, which accounts for corotating interaction regions existing in the solar wind. Such models of the solar wind are often built on top of models of the solar corona, whose parameters are fitted to magnetograms, either daily or mean (averaged over one solar cycle). This approach does not necessarily preclude using the in situ data.
The symmetric models most often include the term of electron density proportional to 1/r2, where r is the distance to the centre of the Sun. This term assumes constant radial velocity of electrons.1 The factor of this term (reference electron density) is a free parameter of the model that is determined from spacecraft ranging observations.
INPOP presumably has two models that are used interchangeably (Verma et al. 2013): (i) a two-parameter model with one term proportional to 1/r2 and the other proportional to 1/r4; and (ii) a two-parameter model with the term proportional to 1/rϵ, where ϵ is the second determined parameter. The parameters are fitted per each solar conjunction of each planet of interest, and also separately for spacecraft being in the fast-wind zone and in the slow-wind zone.
Outside planetary ephemerides, other solar wind electron density models exist; see e.g. Leblanc, Dulk & Bougeret (1998), where 1/r2, 1/r4, and 1/r6 terms are combined. However, the 1/r4 and higher-order terms have a negligible effect on spacecraft tracking signals that pass farther than 15 solar radii (|$15 \, \mathrm{R}_{\odot }$|) from the centre of the Sun. Signals that pass closer than that cannot be used for precise orbit determination anyway because of the temporal and spatial instabilities of the delay in the solar wind in such proximity.
There exist other models with fractional powers of r close to 2; see e.g. Bougeret, King & Schwenn (1984) with r−2.10, Mann et al. (1999) with r−2.16, or Verma et al. (2013) with values of ϵ ranging from 2.10 to 2.50. Such models in general do not possess physical meaning and, as we will now show, do not offer prospects of more accurate representation of spacecraft signal delay.
A major problem in the determination of solar wind parameters from spacecraft ranges (aside from the temporal, latitudinal, and longitudinal variations) is that not only the factor of the 1/r2 (or 1/rϵ) term must be determined from observations, but also the whole set of parameters that form the ephemeris solution. These include parameters of orbits of the planets, the mass of the Sun, the masses of dozens of asteroids, and, most importantly, signal delays (biases) that come from spacecraft transponders going out of calibration and from systematic errors of signal registration in radio observatories (Kuchynka, Folkner & Konopliv 2012).
In fact, just having one factor and one offset to fit for the delay of spacecraft ranges in the solar wind already makes the power of r largely irrelevant. Fig. B1 in Appendix B provides an example where the 1/r2.5 model is fitted to data generated using the 1/r2 model well enough for practical purposes.
We emphasize that there is no power law that can be preferable in a symmetric model of solar wind electron density, either from the position of physical principles or from the position of the processing of signal delays.
1.1 Medium-term density variations
Medium-term variations in the solar wind are ones that have a duration similar to the period of solar rotation (the Sun’s synodic rotation period is 27.2753 d). One obvious problem with the approach of using the 1/r2 model is the assumption that the solar wind is radially symmetric. The asymmetry of the solar wind can be clearly seen, for instance, in the LASCO coronagraph images.2 The fact that spacecraft ranging observations are sensitive enough to be affected by this asymmetry was reported by Kuchynka et al. (2012), who analysed 10-d rolling averages of the observation residuals and noticed the presence of periodic variations with a period of 27 d (see Appendix A for a brief introduction to radio ranging residuals).
The way to account for these periodic variations is to adopt, instead of the 1/r2 model, a three-dimensional numerical model of the solar wind that would simulate (to some level of agreement) the corotating structures of the solar wind based on observations of the Sun.
In this work, we examine the results obtained with the WSA-ENLIL model, which consists of two parts: the ENLIL solar wind model (Odstrcil 2003) and the WSA (Wang–Sheeley–Arge, Wang & Sheeley 1990; Arge & Pizzo 2000; Sheeley 2017) model of the solar corona.
ENLIL is a magnetohydrodynamic (MHD) solar wind model that calculates 3D distributions of the solar wind’s density, velocity, and magnetic field by solving a system of differential MHD equations using a coronal solution given by WSA as a boundary condition at 0.1 au. WSA is a semi-empirical solar corona model that uses line-of-sight magnetic field measurements of the solar photosphere assembled into synoptic magnetograms and simulates the solar wind plasma outflow up to 0.1 au.
The primary purpose of WSA-ENLIL is the prediction of solar wind structures reaching Earth (and, with the help of the CONE model, prediction of geomagnetic storms caused by coronal mass ejections). However, there is no reason that the model cannot be used for other purposes, such as the processing of spacecraft ranging observations.
Fig. 1 shows WSA-ENLIL density distributions in the solar equatorial plane for 2017 January 24 and May 10. The spiral shapes in the map are corotating solar wind structures.

Normalized solar wind electron number density calculated by WSA-ENLIL with daily-updated GONGz magnetograms. The data are plotted in the helio-equatorial plane with a latitude of 0° in the HEEQ coordinate system (see Thompson 2006).
1.2 Long-term density variations
Along with the previously mentioned 27-d periodic variations in the solar wind’s density, there are also long-term variations that happen on a larger scale, i.e. on a scale of years. These long-term variations can be seen by analysing in situ density measurements applying a rolling average filter with a window size large enough to smooth out any medium-term data variations.
Fig. 2 shows 365-d rolling average of the solar wind electron density at 1 au calculated twice: (i) from the OMNI data set (King & Papitashvili 2005; Papitashvili, Bilitza & King 2014) using values of proton density and the alpha–proton ratio and assuming charge neutrality (see Section 3.2), and (ii) from the ENLIL simulations (see Section 3.3). The plot reveals that the average electron density in the long term can deviate up to 20 per cent from its mean value, so the assumption of some ‘nominal’ electron density for all time would be unjustified in the processing of spacecraft radio ranges.

Solar wind electron density at 1 au (365-d rolling average), years 2006–2017.
Our suggestion is that these long-term variations must be taken into account regardless of the particular solar wind model used. Specifically, if a power-law model like N0/r2 is used, one constant N0 multiplier would result in constant density at 1 au, which would not be in agreement with in situ observations. Therefore, N0 must be made proportional to the in situ density smoothed in a way that would eliminate any medium-term disturbances and variations, such as 27-d periodic ones resulting from the Sun’s rotation.
The WSA-ENLIL average density curve shown in Fig. 2 was calculated with fixed values of dfast and vfast, i.e. with one constant value of momentum flux. The figure shows that when momentum flux is fixed, WSA-ENLIL’s average electron density at 1 au (i) is not constant and (ii) does not match the OMNI numbers.
Quite obviously, with average density being non-constant, simply setting dfast proportional to electron density derived from OMNI would do very little to bring the two curves into agreement. Instead, knowing that variations in density are defined by variations in WSA velocity, we presume that the first step should be to use values of vfast based on the WSA velocity distribution vr(ϑ, φ), possibly choosing the distribution’s average as the parameter’s value. Doing so would produce solutions with a near-constant average density at 1 au, and, at that stage, dfast may be set to values proportional to the OMNI density. Such a set-up should, in theory, produce solutions that match the observational in situ data.
However, since ENLIL simulations require time to complete and solutions may be scaled post-factum, the simplest way – used in this work – to bring ENLIL data into agreement with OMNI is to remove the density trend shown in Fig. 2 directly from the solution Ne(t, r, ϑ, φ) and to apply to the data the same procedure as that used for the N0 parameter of the N0/r2 model.
For a detailed description of the in situ data used, see Section 3.2.
2 RADIO SIGNAL PROPAGATION IN THE SOLAR WIND
The term τdist in the above expression corresponds to the time delay due to the distance between |$\mathbf {r}_a$| and |$\mathbf {r}_b$|, i.e. the time that it would take for the signal to travel the same path in vacuum, and τcor is the additional dispersive time delay due to the lower propagation velocity in the solar wind.
In the case of a power-law electron density model Ne(r) ∝ rα, the integral in equation (4) may be expressed in a closed form (see Aksim et al. 2019). In the general case, specifically if the electron density is defined by a grid of numerical values, which is the case with ENLIL and other numerical models, the integral has to be evaluated numerically. To calculate the integral numerically, we used Simpson’s rule (section 5.1 in Kahaner, Moler & Nash 1989) with a step size of 1/160 au, equal to the radial grid spacing of the simulation data used (see Section 3.3), which allows us to exclude possible errors originating from the grid.
3 DATA
Despite this work being focused on Mars spacecraft observations, each of the considered solar wind models was tested as part of its respective planetary solution. Each solution was fitted not only to Mars spacecraft observations, but to the full set of observations usually required to build ephemerides, including not only spacecraft ranging, but also optical observations and other kinds of observations. For details, we refer the reader to Pitjeva (2013), Pitjeva & Pitjev (2014).
3.1 Spacecraft ranges
The Mars spacecraft ranges that were crucial to this work were taken from the webpage of the Solar System Dynamics (SSD) group at NASA JPL,4 which includes the Mars Global Surveyor (MGS, 1999–2006), Odyssey (2002–2013), and Mars Reconnaissance Orbiter (MRO, 2006–2013). The extended set of ranges for MRO and Odyssey up to the end of 2017 was kindly provided by William Folkner. The provided ranges are normal points (NPs) made from raw observations. A normal point represents the distance from an Earth radio observatory to the centre of Mars, corrected for the delay in ionosphere and troposphere.
Following the decision by Kuchynka et al. (2012), pre-2005 Odyssey and MRO ranges acquired by Deep Space Network stations 26 and 43 were removed from analysis, as well as ranges acquired by station 54 before 2008 October.
MGS, Odyssey, and MRO ranges that pass closer than 30|$\, \mathrm{R}_{\odot }$| to the Sun suffer too much from solar wind instabilities and have been removed from analysis. Ranges that pass farther than 30|$\, \mathrm{R}_{\odot }$|, but closer than 60|$\, \mathrm{R}_{\odot }$|, have been treated separately (but equally with other ranges).
Mars Express (MEX) and Venus Express (VEX) ranges (Morley & Budnik 2009) that are also sensitive to the delay in the solar wind were downloaded from the Geoazur website.5 They too are published as observatory–Mars distances corrected for ionosphere and troposphere, but otherwise they are ‘raw’ observations, so normal points had to be made from them before processing. VEX ranges between 2010 August 23 and October 31 were excluded from analysis due to poor residuals, the reason believed to be bad orbit determination due to trajectory control manoeuvres being more frequent than usual;6 VEX ranges between 2011 October 2 and December 5 were also excluded for the same reason. Due to generally big noise in VEX residuals, the decision was made to exclude VEX ranges that pass closer than |$90\, \mathrm{R}_{\odot }$| to the Sun, so that they do not interfere with less noisy observations in the comparison of the solar wind models. The MEX ranges were treated similarly to Odyssey and MRO, with cut-off values of |$30$| and |$60\, \mathrm{R}_{\odot }$|.
MESSENGER (2011–2014, Park et al. 2017) and Cassini (2004–2014) spacecraft ranges (normal points) were also taken from the aforementioned SSD webpage; however, the published values have already been reduced for the delay in the solar wind (calculated by the JPL DE solar wind model), so they did not participate directly in the analysis of solar wind models. However, these observations are very important for Solar system ephemerides as a whole; in particular, they reduce the uncertainty of the mass of the Sun.
In this work, different models of solar wind electron density are compared in their application to building planetary ephemerides. One of the models, WSA-ENLIL, depends (in its WSA part) on the GONG magnetograms that did not exist before 2006 September (see Section 3.3), so the comparison is done for 2006–2017. The spacecraft observations that touch this time span are summarized in Table 1. The reception (Rx) and transmission (Tx) frequencies are listed so as to emphasize that the frequency is important for calculation of the delay of radio signals in plasma; see equation (3).
Spacecraft observations used in the comparison of solar wind models. NPs = normal points, Rx = reception, Tx = transmission, σ = a priori error.
Spacecraft . | Time span . | # of NPs . | σ, one-way . | Frequency . | ||
---|---|---|---|---|---|---|
. | . | |$\gt 60 \, \mathrm{R}_{\odot }$| . | |$30\!-\!60 \, \mathrm{R}_{\odot }$| . | . | Rx . | Tx . |
Odyssey | 2002–2017 | 7988 | 1784 | 0.5 m | 7155 MHz | 8407 MHz |
MRO | 2006–2017 | 1924 | 561 | 0.5 m | 7183 MHz | 8439 MHz |
MEX | 2005–2015 | 2888 | 718 | 1.5 m | 7100 MHz | 8400 MHz |
VEX | 2006–2013 | 1294a | – | 3 m | 7166 MHz | 8419 MHz |
Spacecraft . | Time span . | # of NPs . | σ, one-way . | Frequency . | ||
---|---|---|---|---|---|---|
. | . | |$\gt 60 \, \mathrm{R}_{\odot }$| . | |$30\!-\!60 \, \mathrm{R}_{\odot }$| . | . | Rx . | Tx . |
Odyssey | 2002–2017 | 7988 | 1784 | 0.5 m | 7155 MHz | 8407 MHz |
MRO | 2006–2017 | 1924 | 561 | 0.5 m | 7183 MHz | 8439 MHz |
MEX | 2005–2015 | 2888 | 718 | 1.5 m | 7100 MHz | 8400 MHz |
VEX | 2006–2013 | 1294a | – | 3 m | 7166 MHz | 8419 MHz |
Note. aThe Sun–signal distance for all VEX normal points is greater than |$90 \, \mathrm{R}_\odot$|.
Spacecraft observations used in the comparison of solar wind models. NPs = normal points, Rx = reception, Tx = transmission, σ = a priori error.
Spacecraft . | Time span . | # of NPs . | σ, one-way . | Frequency . | ||
---|---|---|---|---|---|---|
. | . | |$\gt 60 \, \mathrm{R}_{\odot }$| . | |$30\!-\!60 \, \mathrm{R}_{\odot }$| . | . | Rx . | Tx . |
Odyssey | 2002–2017 | 7988 | 1784 | 0.5 m | 7155 MHz | 8407 MHz |
MRO | 2006–2017 | 1924 | 561 | 0.5 m | 7183 MHz | 8439 MHz |
MEX | 2005–2015 | 2888 | 718 | 1.5 m | 7100 MHz | 8400 MHz |
VEX | 2006–2013 | 1294a | – | 3 m | 7166 MHz | 8419 MHz |
Spacecraft . | Time span . | # of NPs . | σ, one-way . | Frequency . | ||
---|---|---|---|---|---|---|
. | . | |$\gt 60 \, \mathrm{R}_{\odot }$| . | |$30\!-\!60 \, \mathrm{R}_{\odot }$| . | . | Rx . | Tx . |
Odyssey | 2002–2017 | 7988 | 1784 | 0.5 m | 7155 MHz | 8407 MHz |
MRO | 2006–2017 | 1924 | 561 | 0.5 m | 7183 MHz | 8439 MHz |
MEX | 2005–2015 | 2888 | 718 | 1.5 m | 7100 MHz | 8400 MHz |
VEX | 2006–2013 | 1294a | – | 3 m | 7166 MHz | 8419 MHz |
Note. aThe Sun–signal distance for all VEX normal points is greater than |$90 \, \mathrm{R}_\odot$|.
3.2 In situ solar wind measurements
The low-resolution (daily) version of NASA/GSFC’s OMNI data set7 was used as the source of in situ data in this work. The data set is combined8 from in situ data measured by multiple spacecraft. The ones that serve as sources of the particle density data for 2006–2017 are ACE (Stone et al. 1998) and Wind (Wilson et al. 2021).
OMNI provides the proton density Np and alpha–proton ratio Na/Np. To calculate electron density Ne, we used the plasma quasi-neutrality condition Ne = Np + Na/2.
The daily density data provided by OMNI are volatile because of the nature of the solar wind structures; see the ‘raw’ electron density data derived from OMNI in Fig. 3. Since we are interested in the (average) electron density along the signal path, our analysis cannot benefit from the information about the daily variations at the end of the signal path (near Earth). On the other hand, we assume that a sufficiently smoothed electron density curve will represent long-term variations of electron density of the solar wind in the whole inner Solar system, including the Earth–Mars line. The Savitzky–Golay filter (Savitzky & Golay 1964) was chosen as one of the most widely used filters for time series. The parameters of the filter were chosen in an ad hoc manual procedure: a window of 511 d and a polynomial order of 5. The filtered data are shown on Fig. 3.

OMNI data of electron density near Earth, years 2006–2017. The filter used is Savitzky–Golay with a window of 511 d and a polynomial order of 5.
Alternatively to OMNI, in situ electron density data are also provided by the ESA/NASA Solar and Heliospheric Observatory (SOHO, Hovestadt et al. 1995). These data9 are not part of the OMNI combination. The results obtained with the SOHO data (not shown in the text) turned out to be somewhat worse than those obtained with OMNI.
3.3 Numerical solar wind density maps
WSA-ENLIL data were produced by request via the ‘runs on request’ service on the Community Coordinated Modeling Center (CCMC) website.10 The whole data set covers dates between 2006 September and 2017 December and consists of 12 separate runs, each one covering one year.
The grid size in the WSA-ENLIL runs was set to 256 × 10 × 90 (r × ϑ × ϕ) with grid boundaries 0.1–1.7 au, −20–+20° latitude, 0–360° longitude (in the HEEQ coordinate system). The grid spacing, therefore, is 1/160 au in the radial direction and 4° in the latitudinal and longitudinal directions. The output frequency for the simulations was set to once every 12 h.
The solar surface magnetic field data used for WSA coronal solutions consisted of daily updated synoptic maps produced by the Global Oscillation Network Group (GONG, Hill 2018).
It must be noted that there are two kinds of GONG magnetograms, called GONGb and GONGz. They were supposed to be identical with the only difference being that GONGz has the non-solar magnetic-field bias removed. However, it was found in late 2015 that GONGz magnetograms suffered from oversubtraction in the polar regions. The correction of the oversubtraction gradually went into effect starting from 2016 August and reaching all GONG stations in 2017 April. Unfortunately, it is not possible to retroactively correct the past GONGz data (see section 5.4 in Mays et al. 2020). So the numerical simulations in 2017 and later are expected to be more accurate when GONGz is used; before 2017, the expected improvement of GONGz is not reached and GONGb is more accurate.
4 ANALYSIS
As said above, full planetary ephemeris solutions were obtained for each of the considered solar wind models. All solutions were based on EPM and the era-8 software (Pavlov & Skripnichenko 2015).
4.1 Dynamical model
EPM has a single dynamical model of the Solar system that includes all planets, the Moon, the Sun, asteroids (see Section 4.2), a discrete model of the Kuiper belt (160 uniformly distributed points of equal mass), and 30 individual trans-Neptunian objects (Pitjeva & Pitjev 2018b). 16 bodies (the Sun, the planets, Pluto, Ceres, Pallas, Vesta, Iris, Bamberga) obey Einstein–Infeld–Hoffmann equations of motion. Other bodies, for the sake of performance, interact with those 16 bodies with only Newtonian forces, and do not interact with each other. It has been checked that taking into account the full set of relativistic interactions affects the Earth–Mars distance by a few centimetres, which is negligible with the present observations.
Apart from point-mass interactions, the model includes additional accelerations from solar oblateness and the Lense–Thirring effect. Earth also gets ‘point-mass figure’ accelerations that come from the Sun, Venus, Mars, Jupiter, and the Moon (Pavlov, Williams & Suvorkin 2016).
The planetary part of the EPM model has over 600 parameters that are determined simultaneously in the least-squares fit (see Appendix A). Here are the parameters that are important for the analysis of radio ranges from Earth to Mars orbiters:
orbital elements for all planets at epoch,
solar oblateness factor,
masses of 279 asteroids,
distance-scale parameter, historically known as correction to the value of au,
single factor of electron density, which applies to the delay in solar wind plasma,
biases to compensate for calibration errors or clock offsets on Earth or on spacecraft: 15 per-station biases for MGS/Odyssey, 14 per-station biases for MRO, single transponder delay for MGS, single transponder delay for Mars Express.11
The distance-scale parameter represents the correction to GM⊙, the mass of the Sun. It is numerically harder to determine the mass of the Sun directly because of correlations with semimajor axes of planets. Determination of GM⊙ (and, in specific solutions, its time derivative) allows one to study the physical properties of the Sun and to test general relativity (Pitjeva et al. 2021).
It is important to process different kinds of observations simultaneously to get realistic values of parameters and error estimates. For example, MESSENGER observations help greatly to determine the solar oblateness factor, because Mercury is close to the Sun. The solar oblateness factor affects the orbits of Earth and Mars. To determine the distance-scaling factor, one has to take ranging data for MESSENGER, MEX, the Mars orbiters, and Cassini. Earth’s orbit, whose precision is obviously important for processing of the Earth–Mars ranges, depends on the Moon’s orbit around Earth, which is determined from lunar laser ranging observations. Optical observations are important for determination of the orbits of Uranus and Neptune, which in turn affect Mars’s orbit on long time spans.
4.2 Perturbing asteroids
An important part of the model, especially in regard to the Mars spacecraft observations, are the point masses that represent the asteroids in the Main Asteroid Belt. In EPM2017, 301 largest asteroids were present in the dynamical model; the remaining part of the Main Asteroid Belt was modelled as a discrete rotating annulus consisting of 180 uniformly distributed points of equal mass (Pitjeva & Pavlov 2017; Pitjeva & Pitjev 2018a). Masses of 16 individual asteroids were fixed to the values determined from observations of spacecraft or natural satellites orbiting them; 30 masses were determined as part of the ephemeris solution (i.e. by the perturbations that they inflict on orbits of inner planets); other asteroids’ masses were determined as mean densities of three taxonomic classes (C, S, M).
In the updated ephemeris dynamical model used in this work, the number of individual asteroids is 279. The source list of asteroids was merged (with obvious removal of duplicates) with 343 asteroids of the DE430 model (Folkner et al. 2014) and 287 asteroids from table A.1 of Kuchynka et al. (2010). The former is believed to be the list of most perturbing asteroids, while the latter is believed to be the list of most ‘non-ring-like-acting’ asteroids; i.e. their cumulative effect on inner planets cannot be modelled by a uniform ring. The resulting list contained 379 asteroids. Then, 100 asteroids whose masses were determined negative (though always within uncertainty) were excluded from the model. Such negative masses are an artefact of the least-squares method and they come from the fact that there are not enough observational data to properly determine them.
None of the 279 masses were fixed in planetary solutions in this work; rather, all of them were determined along with other planetary parameters by a weighted least-squares method enhanced with Tikhonov regularization (Kuchynka & Folkner 2013; Kan & Pavlov 2022). The a priori masses and uncertainties used in the regularization were taken from different sources. The masses of 17 asteroids are known with small uncertainties, thanks to observations by spacecraft or their natural satellites (in comparison with EPM2017, the binary asteroid Euphrosyne was added to the list). For 77 asteroids, the a priori values were those determined from observed deflections of orbits of other asteroids during close approach. For the remaining 185 asteroids, the estimates were taken based on diameters observed in the infrared (Tedesco et al. 2004; Masiero et al. 2011) and densities assigned to the three taxonomic classes.
4.3 Solutions
Three solutions were made, with identical sets of observations and parameters, differing only in the solar-wind delay model for the time span of 2006–2017.
In Solution I, the delays were calculated in accordance with the simplest model of electron density, symmetric and stationary, Ne = CN0/r2, where N0 (the reference electron density at r = 1 au) was made equal to the 2006–2017 time-average of OMNI data, which was found to be 5.97 cm−3. The formal error of this average value was calculated to be 0.07 cm−3 from the formal errors of the proton density and alpha/proton ratio, σ(Np) and σ(Na/Np). The single electron density factor C (which works as the factor of delay; see the τcor term in equation 4) was fitted to the observations. Ideally, C would be equal to 1; the purpose of C is to account for uncertainty in the observed electron density, and also for the fact that the assumption of stationary electron density in 2006–2017 is not very realistic.
In Solution II, the delays were calculated in accordance with a symmetric but non-stationary electron density model, Ne = CN1(t)/r2, where the function N1(t) is is equal to the smoothed electron density values derived from OMNI (see Section 3.2 and Fig. 3). The single factor C was fitted to observations.
In Solution III, the delays were calculated in accordance with the WSA-ENLIL model of electron density, scaled so that the (365-d rolling average of) electron density near Earth is equal to the smoothed OMNI values. Similarly to Solutions I and II, a single factor C of the electron density was fitted to observations.
The values of the fitted electron density factor C in the three solutions are listed in Table 2. While its reference value is 1 in all solutions, there is a difference in how it is fitted. In Solution I, C was not constrained, since the assumption of stationary electron density in 2006–2017 is already not very realistic. In Solutions II and III, the Tikhonov regularization was used to effectively constrain C so that its a priori standard deviation is 1.17 per cent. The value was taken equal to the relative standard deviation of the electron density in the OMNI data (see the description of Solution I above).
Weighted root-mean-square (WRMS) statistics and values of important parameters in the three solutions. The listed uncertainties are 3σ.
. | Solution I . | Solution II . | Solution III . |
---|---|---|---|
. | 1/r2 . | 1/r2 + OMNI . | WSA-ENLIL . |
WRMS | |||
Odyssey | 0.564 m | 0.563 m | 0.563 m |
– conj.a | 1.16 m | 1.15 m | 1.15 m |
MRO | 0.619 m | 0.620 m | 0.619 m |
1.20 m | 1.19 m | 1.19 m | |
MEX | 2.46 m | 2.46 m | 2.46 m |
1.96 m | 1.96 m | 1.95 m | |
VEX | 6.24 m | 6.04 m | 6.04 m |
Density | 1.31 ± 0.81 | 1.037 ± 0.03 | 1.037 ± 0.03 |
factor C | |||
|$\Delta \mathrm{GM}_{\odot }^b ,\ \mathrm{km}^3\, \mathrm{s}^{-2}$| | 0.972 ± 0.526 | 1.457 ± 0.463 | 1.470 ± 0.463 |
Variance of | 1.1127 | 1.1124 | 1.1124 |
unit weight | |||
(reduced χ2) |
. | Solution I . | Solution II . | Solution III . |
---|---|---|---|
. | 1/r2 . | 1/r2 + OMNI . | WSA-ENLIL . |
WRMS | |||
Odyssey | 0.564 m | 0.563 m | 0.563 m |
– conj.a | 1.16 m | 1.15 m | 1.15 m |
MRO | 0.619 m | 0.620 m | 0.619 m |
1.20 m | 1.19 m | 1.19 m | |
MEX | 2.46 m | 2.46 m | 2.46 m |
1.96 m | 1.96 m | 1.95 m | |
VEX | 6.24 m | 6.04 m | 6.04 m |
Density | 1.31 ± 0.81 | 1.037 ± 0.03 | 1.037 ± 0.03 |
factor C | |||
|$\Delta \mathrm{GM}_{\odot }^b ,\ \mathrm{km}^3\, \mathrm{s}^{-2}$| | 0.972 ± 0.526 | 1.457 ± 0.463 | 1.470 ± 0.463 |
Variance of | 1.1127 | 1.1124 | 1.1124 |
unit weight | |||
(reduced χ2) |
Notes. a‘conj.’ (conjunction) observations are those where the closest distance from the Earth–Mars signal to the Sun is between 30 and 60 solar radii.
b|$\Delta \mathrm{GM}_{\odot } = \mathrm{GM}_{\odot } - 132\,712\,440\,041\ \mathrm{km}^3\, \mathrm{s}^{-2}$|.
Weighted root-mean-square (WRMS) statistics and values of important parameters in the three solutions. The listed uncertainties are 3σ.
. | Solution I . | Solution II . | Solution III . |
---|---|---|---|
. | 1/r2 . | 1/r2 + OMNI . | WSA-ENLIL . |
WRMS | |||
Odyssey | 0.564 m | 0.563 m | 0.563 m |
– conj.a | 1.16 m | 1.15 m | 1.15 m |
MRO | 0.619 m | 0.620 m | 0.619 m |
1.20 m | 1.19 m | 1.19 m | |
MEX | 2.46 m | 2.46 m | 2.46 m |
1.96 m | 1.96 m | 1.95 m | |
VEX | 6.24 m | 6.04 m | 6.04 m |
Density | 1.31 ± 0.81 | 1.037 ± 0.03 | 1.037 ± 0.03 |
factor C | |||
|$\Delta \mathrm{GM}_{\odot }^b ,\ \mathrm{km}^3\, \mathrm{s}^{-2}$| | 0.972 ± 0.526 | 1.457 ± 0.463 | 1.470 ± 0.463 |
Variance of | 1.1127 | 1.1124 | 1.1124 |
unit weight | |||
(reduced χ2) |
. | Solution I . | Solution II . | Solution III . |
---|---|---|---|
. | 1/r2 . | 1/r2 + OMNI . | WSA-ENLIL . |
WRMS | |||
Odyssey | 0.564 m | 0.563 m | 0.563 m |
– conj.a | 1.16 m | 1.15 m | 1.15 m |
MRO | 0.619 m | 0.620 m | 0.619 m |
1.20 m | 1.19 m | 1.19 m | |
MEX | 2.46 m | 2.46 m | 2.46 m |
1.96 m | 1.96 m | 1.95 m | |
VEX | 6.24 m | 6.04 m | 6.04 m |
Density | 1.31 ± 0.81 | 1.037 ± 0.03 | 1.037 ± 0.03 |
factor C | |||
|$\Delta \mathrm{GM}_{\odot }^b ,\ \mathrm{km}^3\, \mathrm{s}^{-2}$| | 0.972 ± 0.526 | 1.457 ± 0.463 | 1.470 ± 0.463 |
Variance of | 1.1127 | 1.1124 | 1.1124 |
unit weight | |||
(reduced χ2) |
Notes. a‘conj.’ (conjunction) observations are those where the closest distance from the Earth–Mars signal to the Sun is between 30 and 60 solar radii.
b|$\Delta \mathrm{GM}_{\odot } = \mathrm{GM}_{\odot } - 132\,712\,440\,041\ \mathrm{km}^3\, \mathrm{s}^{-2}$|.
Pre-2006 spacecraft ranges, which include MGS and parts of Odyssey and MEX, were processed with the symmetric stationary solar wind model, with different factors fitted for different periods of time between Sun–Mars conjunctions. One factor was fitted for 2002 January 1–2003 December 31, another for 2004 January 1–2006 September 29.
A different factor was fitted for pre-2002 MGS ranges. The same factor was used for pre-1990 radar observations of Mars, Mercury, and Venus, although they are hardly sensitive to the delay in the solar wind because of low accuracy.
5 RESULTS
Fig. 4 shows the residual (post-fit) differences of the Odyssey (blue points) and MRO (red points) ranges plus the values of the model delay in the solar wind. For better visibility, only ranges from mid-2013 to mid-2017 are shown. The model delay in the solar wind is shown separately with green curves, so that one can see how the model delay curve is fitted to observations.

Odyssey and MRO range residuals without correction for the delay in the solar wind (blue and red points, respectively), their 5-d rolling averages (black curve), and model values (green curve) of delays in the solar wind fitted to observations via the factor C. The top plot corresponds to Solution I (CN0/r2 model), the middle plot corresponds to Solution II (CN1(t)/r2 model), and the bottom plot corresponds to Solution III (WSA-ENLIL model). The time span is from mid-2013 to mid-2017. The vertical grid lines are placed 27.2753 d apart so as to assist in examining the periodic variations due to the Sun’s rotation.
The black curves show the time-averaged spacecraft range residuals. They have visible medium-term variations with period close to 27 d (the period changes as Earth and Mars move on their orbits; the vertical lines on Fig. 4 are separated by the solar rotation period). The amplitude of variations increases near Sun–Mars conjunctions, which is expected because, during the conjunctions, Earth–Mars signals pass close to the Sun, and at small distances from the Sun the electron density is higher than at large distances. Observations with signal passing closer than |$30\, \mathrm{R}_{\odot }$| to the Sun are not shown.
On the bottom plot of Fig. 4 it can be seen that the WSA-ENLIL solar-wind delay model also has 27-d variations, as expected. In some cases there is a good match between peaks in the model and in the residuals, but it often goes ‘out of sync’ or the amplitudes become very different. One explanation for this is the unfortunate fact that all the GONG magnetograms are obtained from Earth – and the daily-updated map of the magnetic field that is fed to the ENLIL’s MHD equations is correct only on the visible side of the Sun. On the non-visible side, it is a prediction, often not a very accurate one.
Before 2017, we can also attribute the differences between WSA-ENLIL model data and real spacecraft delays to the GONGb magnetograms containing non-solar magnetic-field bias (see Section 3.3). For 2017, when GONGz magnetograms were used, we support the first explanation with anecdotal evidence. A good match between the WSA-ENLIL model and spacecraft data can be clearly seen in the bottom plot of Fig. 4 for the first three months of 2017, while in the following three months this is no longer the case. Fig. 1 shows two example maps for the two quarters of 2017. One can see that on the first map the corotating structures that intersect the Earth–Mars line originate from the visible half of the Sun’s surface; while on the second map the Earth–Mars line has changed in the way that it is also intersected by a corotating structure whose origin is on the non-visible half. (If the pictures were animated, the structures would rotate anticlockwise in time.) We admit that this cannot be accepted as strong evidence and hope for a future possibility to analyse post-2017 spacecraft data together with WSA-ENLIL simulations based on properly fixed GONGz magnetograms.
The Mars Express data are shown in Fig. 5. One can see, apart from the expected peaks with a period of approximately 27 d, long-term variations that far exceed the peaks. The cause of the variations certainly cannot be an inaccuracy of the EPM’s dynamical model of the motion of Mars, reductions, or the solar-wind delay model, because there are no variations of that magnitude in the Odyssey and MRO residuals. It is also unlikely that the variations come from a radio observatory on Earth: the same Deep Space Network stations observe the Mars Express, Odyssey, and MRO spacecraft. The most reasonable explanation would be that there are problems in orbit determination for the spacecraft. In any case, the comparison of solar-wind delay models against the MEX observations currently makes little scientific sense, at least with the ranges that pass farther than 30 |$\, \mathrm{R}_{\odot }$| from the Sun.

MEX range residuals without correction for the delay in the solar wind (red points), their 5-d rolling averages (black curve), and model values (green curve) of solar-wind delays fitted to observations via the factor C. The plot corresponds to Solution III (WSA-ENLIL model). The time span is from mid-2013 to mid-2017. The vertical grid lines are placed 27.2753 d apart so as to assist in examining the periodic variations due to the Sun’s rotation.
The statistics of the post-fit results and values of important parameters are given in Table 2. Solution II (time-dependent symmetric model) is slightly better than Solution I (stationary symmetric model) overall by the reduced χ2 metric and particularly better near conjunctions. Solution III (WSA-ENLIL-based solar-wind delay model) is no better than Solution II, most probably because of the aforementioned problems with magnetograms.
In Solutions II and III, the electron density delay factor C was constrained by Tikhonov regularization to be close to the reference value of 1. An additional experiment (not shown) was made where the Tikhonov regularization was not applied to C in Solution II. The resulting value of C turned out to be 1.35 ± 0.81 (as opposed to 1.037 ± 0.03 in the ordinary Solution II). This confirms that proper usage of absolute in situ measurements of electron density really helps to bring the solution, including the GM⊙ estimate, closer to reality. Also, the constraint on the factor helps to reduce the uncertainty of GM⊙ by ∼12 per cent (the uncertainty is 0.526 in Solution I and 0.463 in Solution II).
Another important improvement of Solution II over previous EPM solutions is that only a single parameter related to the solar wind is fitted in the solution, as compared to more than 50 in EPM2017 (see B and |$\dot{B}$| in equation 1 per year per planet) or 10 (B per conjunction from 2002 to 2018 and single B prior to 2002) in a later EPM version that was not released (Pavlov 2019; Pitjeva et al. 2021).
6 CONCLUSIONS AND FUTURE WORK
Conclusions:
Interplanetary spacecraft ranging observations are sensitive to signal delay due to the solar wind and, in particular, to its medium-term variations caused by the rotation of the Sun and its long-term variations caused by changes in the Sun’s activity.
The ENLIL model of the solar wind, built on daily magnetograms done by the GONG network, can serve as a basis for a numerical solar-wind delay model. Such a model, though, can have or not have a good match to the MRO and Odyssey ranges. Spacecraft ranges that pass through the parts of the electron density map that are effectively built from visible disc magnetograms have (with otherwise equal conditions) better residuals than the ranges that pass through the parts of the map that are built from older magnetograms rotated to the present date.
Mathematical models of ephemerides can benefit from taking into account a solar-wind delay model that has long-term variability determined from the OMNI data set of in situ observations.
Constraining the (usually unconstrained) electron density factor in a planetary ephemeris solution to a reasonable range allows one to reduce the error in the determined gravitational parameter of the Sun by ∼12 per cent.
Using the WSA-ENLIL-based solar-wind delay model to account for medium-term variations can improve the residuals in individual cases, but does not improve the solution overall.
Interplanetary spacecraft ranging observations are sensitive enough to electron density so that they are used in this work to perform comparative analysis and validation of solar wind models. They can be used for validation of future solar wind models; and potentially they can be used as a source of data for solar wind models.
ACKNOWLEDGEMENTS
The authors would like to thank their colleagues from the IAA RAS: Elena Pitjeva for continuous support of this work and a lot of helpful advice regarding the planetary ephemeris solutions, and Margarita Kan for implementing the Tikhonov regularization method in the era-8 software and putting together the list of a priori masses of asteroids.
Planetary observations, including spacecraft ranges that are essential for this work, were collected at the NASA SSD webpage supported by William Folkner and at the Geoazur website supported by Agnès Fienga.
The authors acknowledge use of NASA/GSFC’s Space Physics Data Facility’s OMNIWeb service, and OMNI data, and are grateful to the OMNIWeb team.
Solar wind simulation results have been provided by the Community Coordinated Modeling Center at Goddard Space Flight Center through their public ‘runs on request’ system.
The numerical heliospheric code ENLIL has been developed by Dusan Odstrcil who is currently at the George Mason University.
The WSA model, which gives the boundary conditions for the ENLIL model, was developed by Yi-Ming Wang and Neil R. Sheeley Jr at the Naval Research Laboratory and later improved by Nick Arge, who worked at the National Oceanic and Atmospheric Administration in Boulder at the time and now works at the Goddard Space Flight Center.
The authors are indebted to Dusan Odstrcil for dedicated numerical experiments and detailed explanations of inner workings of the ENLIL model, to Leila Mays for organizing the numerical simulations at the CCMC, and to Maria Kuznetsova for leading the CCMC itself.
DATA AVAILABILITY
The Mars Global Surveyor (MGS), Mars Reconnaissance Orbiter (MRO), and Odyssey spacecraft ranging data underlying this article are available at the webpage of the Solar System Dynamics (SSD) group at NASA JPL.12
The extended set of ranges for MRO and Odyssey up to the end of 2017 was kindly provided by William Folkner of the SSD group by permission and may be shared on reasonable request to the SSD team of NASA JPL.
Mars Express (MEX) and Venus Express (VEX) ranges were downloaded from the Geoazur website.13
The WSA-ENLIL data were produced by request via the ‘runs on request’ service at the Community Coordinated Modeling Center (CCMC) website.14 The data are available for download through the ‘runs on request’ service;15 the names of the runs are:
Dan_Aksim_122319_SH_2 (2006)
Dan_Aksim_122319_SH_3 (2007)
Dan_Aksim_122319_SH_4 (2008)
Dan_Aksim_122319_SH_5 (2009)
Dan_Aksim_122319_SH_6 (2010)
Dan_Aksim_122319_SH_7 (2011)
Dan_Aksim_122319_SH_8 (2012)
Dan_Aksim_122319_SH_9 (2013)
Dan_Aksim_120920_SH_1 (2014)
Dan_Aksim_122319_SH_11 (2015)
Dan_Aksim_122319_SH_12 (2016)
Dan_Aksim_112119_SH_3 (2017)
The OMNI in situ data are available publicly through the OMNIWeb website.16 Specifically, the file ‘omni_01_av.dat’, containing daily averages of solar wind parameters, comes from the low-resolution version17 of the OMNI data set.
Footnotes
The connection between density distribution Ne(r) and radial velocity vr(r) is expressed by the equation r2vr(r)Ne(r) = const, which follows from the continuity equation (Aksim et al. 2019).
Conforming to the tradition of measuring time delays in units of length, by ‘time delay τ’ we mean not time delay specifically, but rather optical path-length, which is measured in metres.
This information was obtained in personal communication from Elena Pitjeva, who in turn obtained it from Trevor Morley.
Methods used for the combination are outlined at http://omniweb.gsfc.nasa.gov/html/omni_min_data.html
It is not possible, nor necessary, to determine a transponder delay for MRO, because it is absorbed by the per-station biases in the solution. For the same reason, of the MGS/Odyssey pair, only the MGS transponder delay is determined.
REFERENCES
APPENDIX A: RADIO RANGING RESIDUALS AND DETERMINATION OF PARAMETERS
A radio signal travels twice the distance between the observatory on Earth and the spacecraft. The time delay is being measured. The ranging residual is the difference between the ‘observed’ delay (normal point; see Section 3.1) and the model (‘computed’) value of the delay. The process of calculation of ‘observed minus computed’ value of the delay is called reduction.
The precision of the normal point is sub-metre; for the model value to match that precision, its formula must account for the speed of light, relativistic signal delay caused by the Sun and other massive bodies, the proper time of the observer, delay in the troposphere and ionosphere, and delay in the solar wind. All the mentioned delays are to be calculated twice – for the up and down legs of the signal. A detailed description can be found e.g. in Moyer (2003). For each normal point, two equations are solved iteratively: one to find the time of retransmission of the signal by the spacecraft; and the other to find the time of reception of the retransmitted signal.
Observations must also be corrected for the spacecraft transponder delay (the short moment of time between reception of the signal by the spacecraft and transmission of the return signal). There is always some delay; it is measured before launch, but varies due to space conditions as the spacecraft makes its way to the planet. The transponder delay thus itself must be determined from radio ranges.
In addition to transponder delays, there are biases that come from miscalibrations on radio observatories. Following the decision by Kuchynka et al. (2012), two sets of biases were determined for Deep Space Network (DSN) stations: one for the MGS and Odyssey spacecraft, another for the MRO spacecraft.
One has to transform the position of the station in the terrestrial reference frame to the geocentric celestial reference frame (GCRS) at times t1 and t3. The standard IAU2000/2006 precession–nutation model (Wallace & Capitaine 2006) and IERS EOP series (Bizouard et al. 2011) are used for this. Also, displacements of the station due to solid-Earth tides and pole tides are to be calculated according to the IERS Conventions 2010 (Petit & Luzum 2010).
Δgrav is the delay of the signal due to space–time curvature; it was calculated in this work with a theoretical approximation (Kopeikin 1990); delays from the Sun, Earth, the Moon, Jupiter, and Saturn were added up.
Δatm is the delay of the signal in the ionosphere and troposphere, usually calculated according to the respective IERS models (Petit & Luzum 2010). In this work, the given observations (normal points) have already been corrected for the atmospheric delay.
Δsol is the delay of the signal in the solar wind, the subject of this work.
While the ‘observed’ values are known and do not change, the ‘computed’ values depend on model parameters (see Section 4.1). The model parameters are to be fitted so that the differences between the observed and computed, i.e. residuals, are minimized. The values of residuals obtained after the fitting (minimization) are called post-fitting residuals.
In reality, the ‘true’ thermal noise of measurement is hard to separate from other systematic errors unaccounted for in the model. The main source of those errors is the solar wind plasma itself; there can also be errors from incompleteness of the set of perturbing asteroids (see Section 4.2) or from wrong estimates of asteroid masses or orbits. Finally, systematic errors (other than the aforementioned constant offsets) can come from hardware on the spacecraft or at the radio observatory. A frequent approach to overcome this problem, which is also used in this work, is to choose the formal errors to be somewhat smaller than a posteriori estimates of the standard deviation of noise. These estimates come from post-fit residuals (see Section 5). The decision to make the formal errors smaller than the post-fit residuals was made in order to account for systematic errors in processing otherwise unaccounted for.
The second term of equation (A3) comes from the Tikhonov regularization modification of least-squares. μj is the a priori value of the parameter βj, while εj is the prior uncertainty of that value (taking εj = ∞ will mean effectively no prior value).
APPENDIX B: Difficulties of determination of solar wind model parameters simultaneously with other parameters of the planetary ephemeris solution
Apart from the solar wind plasma delay, observations of radio signal delays always contain unknowns that are not directly related to the solar wind, but have a similar effect on the ranges. It is important not to be overly confident in e.g. the determined value of the power ϵ in the 1/ rϵ model of solar wind density. Fig. B1 provides a synthetic example with ϵ = 2.5. Let us assume that the ‘real’ solar wind behaves according to the 1/r2 model, but we fit the 1/r2.5 model to observations. That is, we must fit the electron density factor, but also an offset (bias) that is inherent to spacecraft observations. It turns out that the 1/ r2.5 model can be made virtually equivalent to the 1/r2 one as long as a factor and a bias are fitted, and signal paths that pass closer than 15 solar radii to the Sun are excluded from the analysis due to their bad quality, as they normally are in ephemeris solutions (in fact, the barrier is often 30 solar radii or more).

Earth–Mars signal delays with the 1/r2 and 1/r2.5 density models. The multiplier (0.37) and the additive constant (1.16 m) for the 1/r2.5 curve are chosen so that the curve fits the 1/r2 curve on the 57974–58119 interval, which corresponds to the Earth–Mars signal passing farther than |$15\, \mathrm{R}_\odot$| from the Sun.
Another example concerns the fitted electron density factor itself. Even if the power of r in the model is correct, it is important not to be overly confident in the determined factor. In a planetary ephemeris solution, GM⊙ is an unknown parameter that has to be determined from observations. In EPM, GM⊙ is not routinely determined directly; a distance-scale parameter is determined instead. In any case, an incorrect value of GM⊙ adds an error to the computed range, proportional to the range itself. The delay of radio signals in the solar wind is largest near Sun–Mars conjunctions; it so happens that the Earth–Mars distance follows the same pattern. That, together with the unknown bias that also has to be fitted, makes the two sources of error not perfectly separable, as illustrated by a synthetic example in Fig. B2. Let us assume that the blue curve is the ‘real’ delay in the solar wind. If no reduction for solar-wind delay was made in the residuals, the residuals will lie on the blue curve. But if, instead of taking into account the solar-wind delay, we decide to fit the distance-scale factor and bias, we will end up subtracting the red curve from the residuals. The blue curve is not very far off from the blue curve. This leads to the conclusion that an incorrectly determined value of the electron density factor can be partially absorbed by a small scale factor (and hence a small change in the mass of the Sun) and small offset.

Red curve (change of Earth–Mars range due to scale) fitted to the blue curve (model values of the delay in the solar wind), with a factor of (1 + 9.67 × 10−12) and an offset of 0.97 m. The two curves are correlated at 87 per cent.