Propagating EUV solar flux uncertainty to atmospheric density uncertainty

Available online at www.sciencedirect.com

ScienceDirect Advances in Space Research 63 (2019) 3936–3952 www.elsevier.com/locate/asr

Propagating EUV solar flux uncertainty to atmospheric density uncertainty Fabian Schiemenz a,⇑, Jens Utzmann a, Hakan Kayal b a b

Airbus Defence and Space GmbH, D-88039 Friedrichshafen, Germany Julius-Maximilians-Universita¨t Wu¨rzburg, D-97074 Wu¨rzburg, Germany

Received 13 August 2018; received in revised form 20 February 2019; accepted 28 February 2019 Available online 11 March 2019

Abstract The call for realistic covariances/uncertainties in orbit determination/estimation and propagation is becoming louder. The main reason for unrealistic uncertainties is, in most cases, a disregard of force model uncertainties, i.e. the calculated covariances are based only on the measurement uncertainties. For this reason the covariances are often scaled in practice to avoid misjudgment of orbit knowledge. J. Emmert has recently shown how solar flux uncertainties can be propagated to in-track orbit position errors in the case of long-term orbit propagation. In his method, the density error resulting from solar flux errors is obtained via the underlying atmospheric model. In this paper, a universal analytic approximation of density uncertainty is presented, which is a fast, yet reliable, alternative to the case-based propagation via the atmospheric model. The major benefit of the proposed analytic estimations is that they seamlessly integrate with orbit estimation/propagation. The uncertainty estimate can directly be computed without an additional call to the atmospheric model. Ó 2019 COSPAR. Published by Elsevier Ltd. All rights reserved.

Keywords: Solar EUV uncertainty; Thermospheric density uncertainty; Error propagation; NRLMSISE-00; DTM-2012; DTM-2013

1. Introduction Knowledge of the orbital uncertainty is crucial for many space surveillance applications. The association of tracks, computation of collision probability, anomaly detection and the sensor tasking and scheduling all rely on accurate uncertainty information of the orbit state vector (Poore et al., 2016). As has been determined in a previous study (Schiemenz et al., 2018), many operational systems still rely on simple weighted least squares orbit determination, together with the assumption of Gaussian uncertainties and a construction of the covariance information based solely on measurement uncertainty. In order to obtain more realistic uncertainty information it is required to, at

DOI of original article: 10.1016/j.asr.2019.02.039.

⇑ Corresponding author.

E-mail address: [email protected] (F. Schiemenz). https://doi.org/10.1016/j.asr.2019.02.040 0273-1177/Ó 2019 COSPAR. Published by Elsevier Ltd. All rights reserved.

least partially, account for force model uncertainties. The main contributor of force uncertainty is atmospheric drag. Semi-empirical models, such as NRLMSISE-00, commonly used the F10.7 proxy to describe the extreme ulraviolet (EUV) heating due to its consistency and availability. It has however been recognized that the choice of the proxy/proxies itself has a major significance to the overall density uncertainty (Marcos, 2006; Dudok de Wit and Bruinsma, 2011). Hence, in the development of the JB2008 model, significant focus has been put on the inclusion of heating-processes at other wavelengths. It therefore requires multiple solar proxies (F10.7, S10 (He-2), M10 (MG-2) and Y10 (X-ray flux)) (Bowman et al., 2008). While the incorporation of multiple proxies generally improves model results, Bruinsma and de Wit also mention that this approach has its limit in terms of statistical significance (de Wit and Bruinsma, 2017). Indeed the DTM-2012 model, which still uses F10.7 as the single index

F. Schiemenz et al. / Advances in Space Research 63 (2019) 3936–3952

for the solar forcing, claimed to be the most accurate model at the time of its publication (Bruinsma, 2013). A further improvement was announced with the release of DTM2013 (Bruinsma, 2015). Different strategies have been used in the past to account for various parts of density uncertainty. Some groups started by assuming a certain error process model directly at the density level (e.g. Wilkins and Alfriend, 2000), others analyzed the atmospheric density uncertainty by comparing the outputs of various models (Sagnieres and Sharf, 2017). Emmert et al. addressed the computation of intrack orbit uncertainties due to density uncertainties (Emmert et al., 2014). A more rigorous follow-up paper has been published later in Emmert et al. (2017). The strongest reduction in density uncertainty in recent years has been achieved with the development of HASDM, which is able to remedy the static nature of the models (Storz et al., 2005). The method makes use of drag-derived density data from calibration-satellites with known area to mass ratio in order to debias the underlying atmospheric model. The Combined Space Operations Center (CSpOC) currently uses the technique to determine near-realtime density corrections to JB2008. While HASDM could theoretically be deployed also against other models, it was traditionally developed against the Jacchia models. Furthermore the need of the calibration satellite measurements restricts the application of HASDM to a very limited community that has access to the required data. In this work we do not consider the density uncertainty of the model itself, but of its input parameters. Any uncertainty in the input parameters will also lead to an additional density uncertainty. This is especially true for orbit forecasting, where estimates of future solar indices are required; however also has its justification in cases in which solar proxy measurements are available, due to the accuracy of the proxy measurements themselves and the inherent undersampling (Tapping, 2013). Essentially the propagation of solar flux errors to intrack orbit position errors requires first a mapping from solar flux errors to atmospheric density errors and then from the density errors to orbital position and velocity errors. In this paper we focus on the first step: the propagation of solar flux errors to errors in exospheric temperature and subsequently to neutral particle density. Three popular atmospheric models are compared: NRLMSISE-00 (Picone, 2002), DTM-2012 (Bruinsma, 2013) and DTM2013 (Bruinsma, 2015). The choice of these models is justified by the fact that NRLMSISE-00 is, despite its age, still widely used in many operational applications and both DTM-2012 and DTM-2013 claim to be overall more accurate than JB2008. In contrast to DTM-2012, DTM-2013 uses the F30 index to account for the EUV heating. The paper is structured as follows: after the mapping from absolute errors in the extreme ultra-violet (EUV) flux input to relative errors in the atmospheric density is derived in Sections 2–4, the computation of the corresponding relative density error variance, based on a given absolute EUV

3937

error variance, is presented. The validation and tuning of the derived equations is the topic of Section 5. The paper ends with a conclusion. While the developments of Section 3 are independent of the chosen model, the contents of Section 2 only apply to models that use either F10.7 or F30 as the single proxy for the EUV heating. Note that in this work a small delta (d) denotes absolute errors, a small epsilon (
) relative errors (
q ¼ dq ) and a caret is used for estiq mated quantities that are derived from baseline parameters. 2. Exospheric temperature error from extreme ultra-violet flux input error 2.1. Derivation of sensitivity

@T ex @E

In low Earth orbits (LEO) (height above ground  200– 1500 km) temperature is approximately constant with altitude and approaches an asymptotic boundary known as ‘‘exospheric temperature”, commonly denoted by T ex or T inf (Emmert, 2015). The studied models (NRLMSISE00, DTM-2012 and DTM-2013) all define exospheric temperature via Eq. (1): T ex ¼ T 0 ð1 þ Gð LÞÞ

ð1Þ

The function Gð LÞ represents spherical harmonic variations in the parameter L, which is a placeholder for a vector of the environmental parameters (latitude, local solar time, solar flux, and geomagnetic activity) (Bruinsma, 2015). T 0 may be considered as the global and temporal average, prior to adaption by Gð LÞ. The strongest driver of the exospheric temperature is the solar energy input. Consequently the sensitivity of the exospheric temperature with respect to the solar flux is of major interest in order to derive an estimation of how solar flux errors propagate to exospheric temperature errors. Today’s atmospheric models typically contain two parameters for the solar flux input: a daily flux input of the previous day and an 81-day average of the flux (which corresponds to an average condition of three solar rotations). The DTM atmospheric models use the average of the past 81 days, whereas NRLMSISE-00 makes use of a centered 81-day average input to the model. The reason for the two parameters is a separation of faster and slower variations in the solar flux input. Both types have different impacts on the thermosphere. The slow variations are represented by the running 81-day average and fast variations via the difference between the daily values and the 81-day average (de Wit and Bruinsma, 2017). Note that the term ”solar flux” in this paper generally refers to the daily flux input argument (denoted E) and not the 81-day-average1 (denoted E), as this is the parameter with the largest error in orbit propagation/estimation (the average flux parameter is more stable).

1 Prior 81-day average for the DTM models and central 81-day average for NRLMSISE-00

3938

F. Schiemenz et al. / Advances in Space Research 63 (2019) 3936–3952

Our first goal is to derive a description of how the exospheric temperature varies for a given solar flux input. We are therefore interested in a description of dT ex ð dEÞ. From Eq. (1) it follows that dT ex T 0 ðGðL2 Þ  GðL1 ÞÞ ¼ dE dE

ð2Þ

In order to proceed with Eq. (2), it is necessary to extract the dependency of the spherical harmonic variations with respect to the solar flux. This is done by analyzing the sensitivity of Gð LÞ in NRLMSISE-00 and DTM-2012/2013 with respect to E. Neglecting all harmonic variations which are not in the solar flux, E, the average contribution of the solar flux to GðL2 Þ  GðL1 Þ may be given in the following form (the asterisk denotes that harmonic variations other than due to the flux input are neglected):

2 dGI ¼GðL2 Þ  GðL1 Þ ¼ c1 E2  E þ c2 E2  E h

2 i  c1 E 1  E þ c2 E 1  E

ð3Þ

h
2
2 i ¼c1 dE þ c2 E2  E  E1  E

ð4Þ

c1 and c2 are model-specific coefficients and dE ¼ E2  E1 . c2 is typically two to three orders of magnitude (OOM) smaller than c1 , which results in a mostly linear dependency between exospheric temperature and the solar flux. Due to its multiplication with the squared flux differences however, the contribution of the second-order term has a strong dependence on dE2 and therefore quickly gains importance in case of increasing absolute solar flux errors. Assuming that the absolute flux error in the input parameter, dE, is small (dE/10 sfu), we may drop the smaller c2 -term and thereby obtain a first-order description of dGI ð dEÞ: dGI ffi c1 dE

ð5Þ

The performance of the first-order approximation, i.e. the effect of the omitted second-order term on the subsequent estimation of the relative density error, is further discussed in Section 5.1 when the developed equations are validated. Inserting Eq. (5) into Eq. (2) yields:

dT ex ðGðL2 Þ  GðL1 ÞÞ c1 dE ffi T0 ¼ c1 T 0 ffi T0 ð6Þ dE dE dE Eq. (6) represents a first order approximation of the sensitivity of the exospheric temperature with respect to the daily solar flux input to the atmospheric model. It has been obtained by assuming that the error in the absolute solar flux input (dE) is sufficiently small to drop the secondorder term in Eq. (4). 2.2. Model-specific parameters and sensitivity In this section we introduce the values of c1 and c2 for each model, which allows to evaluate the sensitivity of the exospheric temperature to absolute errors in the daily solar flux. 2.2.1. NRLMSISE-00 Using the definition of exospheric temperature in NRLMSISE-00 results in h

 2
2 i dGI ¼ c01 þ c01 c3 E  150 dE þ c2 E2  E  E1  E ð7Þ Comparing Eq. (7)  with Eq. (4) we observe that
c1 ¼ c01 þ c01 c3 E  150 . The coefficient c1 hence contains a dependency on the average solar flux input2 E. Solving Eq. (6) for the change in the exospheric temperature results in dT ex ffi T 0 c1 dE ¼ s1 dE

ð8Þ

Inserting the definition of c1 results in 
 dT ex ffi T 0 c01 þ T 0 c01 c3 E  150 dE

ð9Þ

Next, define s01 ¼ T 0 c01 ; s2 ¼ T 0 c2 ; s3 ¼ T 0 c01 c3 , which yields 
 dT ex ffi s01 þ s3 E  150 dE ¼ s1 dE ð10Þ The s1 -parameter is computed as the product of the modelcoefficient c1 and the baseline temperature T 0 . It represents the first order-dependency of the change in the exospheric temperature, given an absolute error in the daily solar flux input. A study of the NRLMSISE-00 model coefficients reveals the following parametrization:

Table 1 NRLMSISE-00: model-coefficients required for the computation of dT ex . T0 1027.318

c01

c2

1:65686  10

3

c3 6

5:53887  10

3

2:46423  10

s01

s3

1:7021

4:1944203  103

s1


 1:7021  4:1944203  103 E  150

Table 2 DTM-2013: model-coefficients required for the computation of dT ex . T0

c1

c2

c3

s1

s2

1029.77

2:75746  103

5:68252  106

0:0

2:83955

0:00585169

2

NRLMSISE-00 uses the central 81-day average flux.

F. Schiemenz et al. / Advances in Space Research 63 (2019) 3936–3952

3939

Table 3 DTM-2012: model-coefficients required for the computation of dT ex . T0 1029.77

c1

c2 3

2:51565  10

6

21:4824  10

Emmert et al. mention an approximate dependency of 1.7 between dT ex and dE for NRLMSISE-00 (Emmert et
al., 2014).  This value is obtained by dropping the term s3 E  150 in Eq. (10). This is motivated by the fact that
 s01  s3 E  150 for medium solar activity (120 sfu/E /180 sfu), however it is an unnecessary approximation if E is known. Also note that c2  c1 , as mentioned earlier, which is why the first-order approximation of the change in exospheric temperature is appropriate, as long as the solar flux errors range in the order of a few solar flux units. 2.2.2. DTM-2013 The coefficients c1 and c2 may also be found from the model coefficients file of DTM-2013. The DTM-models have no dependency of c1 on the average flux level. DTM-2013 uses the same model-implementation as DTM-2012, however the model-coefficients file is different. Also the F30 flux is the primary proxy of the model to estimate the EUV energy input. As F30 is subsequently scaled to F10.7 units via (de Wit and Bruinsma, 2017): F 30I ffi 1:554  F 30  1:6

ð11Þ

we can derive a ratio when using F30 scaled to F10.7 units and another for the unscaled F30 flux value. Using F 30I (i.e. considering the situation after scaling F30), we find from the model coefficients file: which when introduced to Eq. (10), results in dT ex ffi c1 T 0 ¼ s1  2:840  2:8 dE

ð12Þ

The value of c2 is very close to that of NRLMSISE-00, hence we expect similar and small deviations from the ratio. The value of the latter however is significantly larger than for NRLMSISE-00 (2.8 vs. 1.7), which indicates that DTM-2013 models a stronger response of the exospheric temperature to the solar energy input. If one is interested in an approximate ratio for the unscaled F30 flux, the conversion factor needs to be multiplied to Eq. (12), as dEI ¼ 1:554  dE. Hence we have: dT ex ffi 2:75746  103  1029:77  1:554  4:4 dEI

ð13Þ

2.2.3. DTM-2012 From the model-coefficients file we obtain: Using the coefficients of Table 3 in Eq. (6), results in: dT ex ffi c1 T 0 ¼ 2:51565  103  1029:77 ¼ s1  2:6 dE

ð14Þ

Two things strike the eye when comparing the DTM-2012 coefficients with those of NRLMSISE-00. First, the ratio is

c3

s1

s2

0:0

2:59054

0:0221219

similar to DTM-2013 and second, c2 is roughly four times larger than for NRLMSISE-00 and DTM-2013. Due to this relative magnitude of c2 , a widely linear relationship cannot be expected and large deviations from the ratio of 2.6 are to be anticipated when using the first-order approximation for larger dE. 2.3. Verification of dT ex ð dEÞ To verify the findings of Sections 2.1 and 2.2 the three models under consideration need to be called with global and temporal average conditions enabled. This allows to disable the sensitivity of the computed atmospheric density with respect to all parameters other than the solar flux. To achieve these conditions, one option is to make use of the switches in NRLMSISE-00 and to change the modelcoefficients file of DTM-2012 and DTM-2013 to zero the dependency of all environment parameters other than the solar flux in L (see also Section 6). Another option is to run a Monte-Carlo simulation with all dependencies enabled. By covering the full theoretical range of the parameters in L in discrete steps and averaging, the results ultimately approach those of the first method. In this section we present the output of the Monte-Carlo iterations, as it also shows the variability in dT ex ð dEÞ for the changes of the other parameters in L. Figs. 1a, b and 2 show the results for an 81-day average solar flux value (central for NRLMSISE-00 and previous for the DTM models) of 183.5 sfu. 2.3.1. NRLMSISE-00 Recall that for NRLMSISE-00 we expect an average
 slope of 1:7021–4:1944203  103 E  150 , which equates to 1.562 for E ¼ 183:5 sfu (red line in Fig. 1a). After averaging all Monte-Carlo iterations the ratio dTdEex was found to be 1.566, which closely matches the expected value of 1.562. Fig. 1a also allows to see the effect of the neglected second-order term in Eq. (4). While the expected dT ex ð dEÞ matches the observed dT ex ð dEÞ closely for small solar flux errors (jdEj 6 10 sfu), the simulated distribution ‘‘bends away” from the expected linear relationship as the magnitude of the absolute daily solar flux error dE increases, with the maximum deviation occurring in the tails of the solar flux axis. 2.3.2. DTM-2013 For DTM-2013 we expect an average slope of 2.8 after scaling F30 to F10.7 units. Judging from Fig. 1b the situation is very similar to NRLMSISE-00. For small dE the change in exospheric temperature is accurately estimated.

3940

F. Schiemenz et al. / Advances in Space Research 63 (2019) 3936–3952

Fig. 1. Performance of the first-order description of dT ex ðdEÞ vs. the Monte-Carlo iterations for NRLMSISE-00 and DTM-2013. The blue bars indicate the variation caused by the environmental parameters L for a certain value of dE. The red line depicts the expected result according to Eq. (10). (For interpretation of the references to colour in this figure legend, the reader is referred to the web version of this article.)

Fig. 2. DTM-2012: Exospheric temperature vs. solar flux for the considered model parameter ranges.

With increasing magnitude of dE however the justification for dropping the second-order term vanishes, resulting in the Monte-Carlo iterations to bend away from the expected line. The relative deviation is similar to that of NRLMSISE-00, as the neglected c2 -parameter is of similar magnitude. 2.3.3. DTM-2012 The expected dT ex ð dEÞ for DTM-2012 is 2.6. Recall, that the neglected c2 -parameter was comparatively large in case of the DTM-2012 coefficients with respect to NRLMSISE00 and DTM-2013. Consequently we expect a pronounced deviation from the linear trend with growing dE. As expected, the sensitivity-plot (Fig. 2) shows a pronounced bending of the Monte-Carlo iterations away from the expected slope. The intensity of this behavior is related to the magnitude of the c2 coefficient. Consequently the domain in which the first-order approximation of dT ex ð dEÞ applies, is limited to a smaller region than for NRLMSISE-00 and DTM-2013. With increasing magni-

tude of dE the delta in the exospheric temperature even decreases again, which indicates that the second ordereffect dominates the first order trend. 3. Relative density error from exospheric temperature error In this section we will derive the sensitivity of atmospheric density to exospheric temperature and use this sensitivity to derive an expression of how the relative density error relates to an absolute error in exospheric temperature. The derivation is valid for all models assuming hydrostatic balance, i.e. that the pressure gradient and gravitational forces cancel. Consider the column of infinitesimal height dz shown in Fig. 3. Hydrostatic balance is expressed as: dP ¼ qg dz

ð15Þ

where g is the local gravitational acceleration, z represents geometric height and P is atmospheric pressure which

F. Schiemenz et al. / Advances in Space Research 63 (2019) 3936–3952

@P @P @P dq þ dT þ dM @q @T @M kT qk qkT ¼ dq þ dT  2 dM M M M

3941

dP ¼

ð22Þ

Inserting Eq. (22) into (19) we obtain

Fig. 3. Hydrostatic Balance: pressure gradient and gravity cancel.

decreases with altitude. It makes sense to transform Eq. (15) into a form in which g does not decrease with height. This is achieved by transitioning to geopotential, which is defined by: Z z /ð zÞ ¼ gdh ) d/ ¼ gdz ð16Þ 0

Substituting Eq. (16) into Eq. (15) yields: dP ¼ qd/

ð17Þ

Geopotential height is the geopotential divided by a reference gravity: fð z Þ ¼

/ ð zÞ ) d/ ¼ g0 df g0

ð18Þ

Inserting Eq. (18) into (17) results in the desired form of the hydrostatic balance equation: dP ¼ qg0 df

ð19Þ

All models considered use z0 ¼ 120 km as reference altitude. At this height 2 m m g0 ¼ 9:807 s2  ð6371=ð6371 þ 120ÞÞ ffi 9:45 s2 . Next we need to use the ideal gas law in order to replace the dependency on pressure with one on temperature: P ¼ nkT ¼ q

kT M

ð20Þ

n ¼ Mq is the mean particle density, M the mean particle mass, k the Boltzman constant and T denotes temperature. Eq. (20) represents an expression for P, whereas Eq. (19) is written in terms of dP. Hence we need to determine dP from Eq. (20). Also we must not forget that the atmosphere consists of a mixture of different gases, each with different q; T and M. Consequently atmospheric pressure varies with density, temperature and mean mass, where the latter is computed according to: Pnspecies s¼1 ns M s M¼ P ð21Þ nspecies s¼1 ns To determine dP from P we need to apply the total differential, which results in

kT qk qkT dq dq þ dT  2 dM ¼ qg0 df () M M q M g0 M dT dM ¼ df  þ kT T M Realizing that as d lnðqÞ ¼ 

dx x

¼d

R

dx x

ð23Þ

¼ d lnð xÞ, we may write Eq. (23)

  Mg0 T df  d ln M kT

ð24Þ

Integrating formula (24) with respect to geopotential height from f0 to f yields Z T ð fÞ M ð fÞ g 0 f M 0 df þ ln  lnðqðfÞÞ ¼ lnðqðf0 ÞÞ  ln T ð f0 Þ M ð f 0 Þ k f0 T ð25Þ Eq. (25) describes the profile of a fully mixed atmosphere. This is true up to  100 km altitude. Between 100 km and 120 km height the mixture of the gases transitions into an atmosphere that becomes dominated by molecular diffusion (Emmert, 2015). In the LEO-altitudes of interest (> 200 km) molecular diffusion and gravity cause the constituent species to separate according to their mass and follow individual hydrostatic balance constraints. Hence the concept of ‘‘diffusive equilibrium” is typically applied to Eq. (25) when working with the upper atmosphere. Diffusive equilibrium means to add a thermal diffusion coefficient to Eq. (25) and to evaluate the equation separately for each species (each is in an individual hydrostatic equilibrium). Written in terms of individual species (index i) and considering altitudes higher than 200 km, Eq. (25) becomes: Z
 T ð fÞ M i g 0 f 1 0 lnðqi Þ ¼ ln qi;0  ð1 þ ai Þ ln  0 df |fflfflfflffl{zfflfflfflffl} T ð f0 Þ k T f ð Þ f0 thermaldiffusion

ð26Þ Note that, since Eq. (26) is per species, we have M ¼ M i ¼ const: and ln MMððff0ÞÞ ¼ ln 1 ¼ 0, where i denotes the ith species. The property we are actually interested in is the sensitivity of the atmospheric density to changes in the exospheric temperature, as this sensitivity allows to relate an (absolute) error in exospheric temperature to a (relative) error in atmospheric density. To do so, we start with a description of neutral atmospheric density in LEO altitudes, based on the individual species in diffusive equilibrium:

3942

F. Schiemenz et al. / Advances in Space Research 63 (2019) 3936–3952

X

nspecies

q¼

d lnðqÞ dT ex

qi ()

i

nspecies X d ln expðlnðqi ÞÞ ¼ dT ex i

! ð27Þ

Introducing Eq. (26) to Eq. (27), we obtain: nspecies X
 d d T ð fÞ lnðqÞ ¼ ln exp ln qi;0  ð1 þ ai Þ ln dT ex dT ex T ð f0 Þ i

 Z f M i g0 1 0  ð28Þ 0 df k f0 T ð f Þ 0

Realizing that dxd lnð f ð xÞÞ ¼ ff ððxxÞÞ ; f0 ¼ fðz0 Þ ¼ 0 and enforcing the constraint of the constant exospheric temperature (T ðfÞ ¼ T ex ¼ const:), we can transform Eq. (28) into hP h
 ii nspecies M i g0 T ex d exp ln q ð Þln  f  1 þ a i i i;0 T0 kT ex dT ex d h
 i lnðqÞ ¼ Pn M i g0 species T ex dT ex exp ln q  ð1 þ ai Þln  f i

i;0

T0

kT ex

ð29Þ which, after evaluation of the differential and simplification, results in Pnspecies d lnðqÞ ¼ dT ex

i

qi;0 ð1þai Þ

T0

h i h i

ð1þai Þ M i g0 exp  MkTi gex0 f þ T ex fexp  MkTi gex0 f kT 2ex h i qi;0 ð1þai Þ T exp  MkTi gex0 f ð1þai Þ ex

ð1þai Þ

ð1 þ ai Þ T exT ex Pnspecies i

T0

ð30Þ

We shall next assume that all species have the same thermal diffusion coefficient. While this may seem crude, it is close to the parametrization of nowadays models, which commonly employ a value different from 0.0 only for the light species helium and hydrogen (with the exception of NRLMSISE-00 which also uses a non-zero thermal diffusion coefficient for argon). The consequence of this assumption is the introduction of a small error in the mapping of dT ex to the relative density error in altitude regions which are dominated by helium and hydrogen (z > 600–800 km, depending on solar condition) due to a different weighting of the thermal component in Eq. (30) for H and He. Consequently we set ai ¼ a and by splitting the sum obtain:

Pnspecies d lnðqÞ ffi dT ex

i

 ð 1 þ aÞ

qi;0 ð1þaÞ T0

d  ð 1 þ aÞ g 0 f lnðqÞ ffi þ 2 dT ex T ex kT ex h i Pnspecies M i g0 f M n exp  i i;0 i kT ex h i Pn M i g0 f species n exp  i;0 i kT ex

ð32Þ

Observing

Eq. (32) closely, we realize that h i M i g0 f ni ð zÞ ¼ ni;0 exp  kT ex is the barometric formula per species and hence describes the particle density at height z. Consequently Pnspecies d  ð 1 þ aÞ g 0 f M i ni ð z Þ i lnðqÞ ffi þ 2 P ð33Þ nspecies dT ex T ex ni ð z Þ kT ex i Comparing Eq. (33) with Eq. (21) it becomes evident that the last factor of the second summand is nothing else than the mean mass computed over all species at altitude z. The sensitivity of log density to exospheric temperature is hence approximated via d g M ð zÞ ð 1 þ aÞ lnðqÞ ffi 0 2 f  dT ex T ex kT ex

ð34Þ

The geopotential height in Eq. (34) can be calculated via: h iz R z GM 1 0 GM  0 dz R þz z0 ðR þz0 Þ2 z0 fð z Þ ¼ ¼ GM g0 2 ðR þz Þ

0

R þ z 0 ¼ ð z  z0 Þ R þ z

ð35Þ

where R denotes the radius of Earth. Inserting Eq. (35) into formula (34) we obtain d g M ð zÞ R þ z0 ð 1 þ aÞ lnðqÞ ffi 0 2 ð z  z0 Þ  dT ex T ex R þ z kT ex

ð36Þ

Setting a ¼ 0, as discussed above, results in a description of the logarithmic density sensitivity: d g M ð zÞ R þ z0 1 lnðqÞ ffi 0 2 ð z  z0 Þ  dT ex T ex kT ex R þ z

ð37Þ

The final step is to obtain an expression of the atmospheric density error, given an exospheric temperature error. To first order the change of a function with respect to an error in its input is given by the derivative with respect to the

h i P h i n q M i g0 ð1þaÞ M i g0 exp  MkTi gex0 f þ i species ði;01þaÞ T  f exp  f 2 ex kT ex kT ex T h 0 i Pnspecies qi;0 ð1þaÞ M i g0 exp  kT ex f ð1þaÞ T ex i

ð1þaÞ

T ex T ex

ð31Þ

T0

Assuming the error in the mean mass due to an error exospheric temperature to be negligible, Eq. (31) can be further simplified to

parameter of interest times the error in that argument. Applied to the log density we obtain:

F. Schiemenz et al. / Advances in Space Research 63 (2019) 3936–3952

d lnðqÞ ffi

d lnðqÞdT ex dT ex

ð38Þ

Inserting Eq. (37) into Eq. (38) yields: d lnðqÞ ffi

! g 0 M ð z Þ R þ z 0 1 ð z  z0 Þ  dT ex T ex kT 2ex R þ z

ð39Þ

For functions described in terms of the natural logarithm of its argument, the absolute error in the logarithm is equal to the relative error in the argument, since: d lnðqÞ ffi

d ln q dq dq ¼ ¼
q dq q

ð40Þ

Introducing Eq. (40) to (39) we obtain an expression which relates the relative density error in percent to an absolute error in exospheric temperature in Kelvin: ! g0 M ð zÞ R þ z0 1 dT ex
q ffi ð z  z0 Þ  ð41Þ T ex kT 2ex R þ z The relative density error increases with height z and mean molecular mass M ð zÞ, which however decreases with altitude. The profile also depends on exospheric temperature. Fig. 4 contains plots of
q vs. altitude using NRLMSISE00 configured with global and temporal average conditions, different solar flux-levels (min, mean, max) and dT ex ¼ 10 K. Eq. (41) requires a thermospheric density model to compute T ex and M ð zÞ, however no additional model calls are necessary, as T ex and M ð zÞ may be taken from the model call required for the orbit propagation/estimation. It is approximately valid for all models assuming hydrostatic balance and diffusive equilibrium. An alternative to Eq. (41) is to use an atmospheric density model to first compute a baseline atmospheric density for a given exospheric temperature and altitude and then to observe the relative density errors for a set of different absolute errors in exospheric temperature (”model-based propagation”).

3943

4. Relative density error and variances from solar irradiance error Combining Eqs. (8) and (41) we can formulate an expression of the relative density error (in percent) with respect to an absolute solar irradiance error (in sfu). The first-order relationship reads ! g 0 M ð z Þ R þ z 0 1
q ffi s1 ð z  z0 Þ  dE ð42Þ T ex kT 2ex R þ z where the s1 coefficient is to be taken from Tables 1, 2 or 3, depending on the model used. The variance associated with Eq. (42) is obtained easily via the relationship VarðaX Þ ¼ a2 Varð X Þ. Considering dE as the random variable, we obtain the following estimate of the associated variance: !2
 1 2 g 0 M ð z Þ R þ z 0 Var
q ffi s1 ð z  z0 Þ  Varð dEÞ T ex kT 2ex R þ z ð43Þ Eq. (43) relates the variance of the solar flux error to the variance in the relative density error. Since we have a linear dependency of both quantities, the relative density error will follow the same stochastic process as the absolute solar flux error. The modeling of the latter depends on the usecase and is treated to a large extend in (Emmert et al., 2017). In case of simulations, daily solar flux measurement data is generally available. The inter-day variations of the flux to this parameter are best modeled via white noise, whereas for long-term propagation and operational usecases only propagated flux data is available. In these cases the solar flux error is best modeled as a Brownian motion process (Emmert et al., 2017). We shall consider both cases in the next section. 5. Validation The validation of Eqs. (42) and (43) is performed against the atmospheric model-based error propagation. Considering dE to follow a white-noise process and a Brownian motion process, we compare the relative errors and the associated variances. Note that for the validation against the models all variations, except for the flux dependency, have been disabled. 5.1. Limits of analytic approximation

Fig. 4. Relative density error vs altitude according Eq. (41) for different solar conditions.

The NRLMSISE-00 atmospheric model implements logical switches that allow to en-/disable variations (symmetrical annual, semidiurnal, . . .). For this model switches 2–15 were disabled, which results in a sensitivity of the model only to the solar flux inputs and altitude. After setting up the absolute solar EUV error process according to a Brownian motion or white noise process model, we computed the baseline case, which represents our estimate of

3944

F. Schiemenz et al. / Advances in Space Research 63 (2019) 3936–3952

Fig. 5. Comparison of model-derived relative density errors and uncertainty with analytic approximations.

Fig. 6. Second order description of dT ex ðdEÞ (blue) vs. first order approximations (green and red) for global average conditions and z ¼ 500 km, as well as E ¼ 200 sfu. The second order equation perfectly matches dT ex as no approximation is involved. (For interpretation of the references to colour in this figure legend, the reader is referred to the web version of this article.)

the following parameters: M b ; T exb and qb . The true values of these quantities are computed using the atmospheric model and a true daily flux argument of Et ¼ Eb þ dE. The relative errors are always computed relative to the baseline quantity (
X ¼ dX ) in order to obtain unbiased relXb ative errors.3 As an initial test we performed 500 Monte-Carlo iterations using E ¼ 80 sfu, Eb ¼ 110 sfu and z ¼ 500 km as baseline parameters. We then compared

q with ^
q (according Eq. (42)), as well as r
q with
 r ^
q (according Eq. (43)). The results of this initial test are shown in Fig. 5. The relative solar EUV irradiance

3

If the relative errors were computed against the true quantity, this would bias the relative errors. Consider e.g. X b ¼ 5 and dX ¼ 1. Then dX 1 1 1 X ¼ 6 and  4, whereas errors relative to the baseline result in 5.

errors are depicted in the first subplot. The process model is a random walk model with a standard deviation of 11 sfu after 7 days of simulation. The second subplot shows the model-derived relative density errors b (qq ). In both cases the red line shows the 3r deviation qb of the 500 member ensembles. The third subplot depicts the relative density errors according to the analytic approximation (Eq. (42) using the baseline parameters). Also the estimated 3r deviation according Eq. (43) is shown in blue. For a comparison of the predicted and the actual uncertainty, the blue curve of subplot three is also added to subplot two. Two main points strike the eye: (see Fig. 6).  The uncertainty predicted by Eq. (43) matches the uncertainty of the predicted relative density errors according Eq. (42) (subplot three).

F. Schiemenz et al. / Advances in Space Research 63 (2019) 3936–3952

 The uncertainty according Eq. (43) slightly overestimates the actual uncertainty by a constant factor of ^
q ¼ 1:11 (subplot two).
q In order to further analyze the constant overestimation ^
of this case, we next compared
qq using Eq. (41) to determine the contribution to the overestimation of the respective parts of Eq. (42) (the mapping from the absolute EUV flux error to an absolute exospheric temperature error and the mapping from the absolute exospheric temperature error to the relative density error). We found that ^
also Eq. (41) leads to
qq ¼ 1:11. This indicates that the overestimation is purely due to the equation describing the relative density error from an absolute error in exospheric temperature, meaning that the first-order mapping from dE to dT ex is sufficient in this case. In a second test using the parameters z ¼ 500 km, E ¼ 120 sfu and Eb ¼ 90 sfu, ^
we found
qq ¼ 0:87 using Eq. (42), which shows an underestimation by 13%, when both error mappings are combined. Using Eq. (41) with dT ex obtained from the model, reveals that the mapping from the absolute exospheric temperature ^
error to the relative density error itself leads to
qq ¼ 1:038. Hence, in this case the dominant contribution of the underestimation is the first-order approximation of dT ex ð dEÞ according to Eq. (8). These initial tests show that the analytic expressions both over- and underestimate the model-propagated uncertainties. To get a more complete view, the performance was evaluated for different altitude/flux-combination, resulting ^
in a grid of t :¼
qq -values. t > 1:0 indicates overestimation and t < 1:0 underestimation. In each case the daily solar

3945

first step a second-order expression of the absolute exospheric temperature error with respect to the absolute solar flux error has been derived. The following derivations represent improvements with respect to the first-order approximations given in Section 2.1. Reconsider Eq. (4). In Section 2.1 the second-order term of dGI has been dropped. Once dE however becomes large, a higher-order relationship of dT ex ð dEÞ is required in order to compute an accurate change in exospheric temperature. Instead of dropping the c2 -term, we may continue Eq. (4) using E2 ¼ E1 þ dE, which leads to: h
i 2 dGI ¼c1 dE þ c2 ðE1 þ dEÞ  E  E21 þ 2E1 E  E2 ð44Þ
 ¼c1 dE þ c2 dE2 þ 2c2 dE E1  E

ð45Þ

In case of small disturbances (E1  E), the last term may be dropped. Since however we are seeking to fully optimize the computation of dT ex , we shall continue with Eq. (45). We obtain the corresponding change in exospheric temperature via: dT ex dGI ¼ T0 dE dE


 c1 dE þ c2 dE2 þ 2c2 dE E1  E () dT ex ¼ T0 dE 
 ¼ s1 þ 2s2 E1  E dE þ s2 dE2

ð46Þ

Consequently the relative density error may be reformulated as ð47Þ

flux argument, E, was chosen to vary according to the solar flux error process model about E, i.e. Eb ¼ E, which resembles periods of low disturbances. This choice of E is in favor of the first-order description of dT ex ð dEÞ, as the dependence of dT ex on E is part of the second-order term, which has been dropped so far. It hence represents the optimum performance achievable via a first-order description of dT ex ð dEÞ. The obtained performance using t ¼ 0:941; tmin ¼ 0:551; tmax ¼ NRLMSISE-00 was: 1:692; rt ¼ 0:256, showing that, on average, the results are matching well, however larger over- and underestimations are possible for certain Eb ¼ E; z-combinations. 5.2. Improved analytic approximations 5.2.1. Second order description of dT ex ð dEÞ As the performance of Eqs. (42) and (43) was not yet satisfactory, we sought for possible improvements. As a

where dE and dE2 follow a Gaussian distribution. Next,

 0 ð z  z0 Þ define s12 :¼ s1 þ 2s2 E1  E and f ¼ g0kTM2ðzÞ RR

þz þz ex

 T1ex . Then the variance of the second-order approximation may be computed via

 2 Var
q ¼ ðs12 fÞ Varð dEÞ þ 2s12 s2 f2 Cov dE; dE2
 þ ðs2 fÞ2 Var dE2 ð48Þ Define the random variables X :¼ dE ¼ r0 Z þ l and Y :¼ X 2 ¼ dE2 , where Z follows a standard normal distribution. Using the definition of covariance and the moment-generating function of the normal distribution we find: 
 2 Cov X ; X 2 ¼ Cov r0 Z þ l; ðr0 Z þ lÞ
 ¼ r03 Cov Z; Z 2 þ 2lr02 ð49Þ

3946

F. Schiemenz et al. / Advances in Space Research 63 (2019) 3936–3952


     Since Cov Z; Z 2 ¼ E Z 3  E½Z E Z 2 and E½Z ¼ 0, as  3 well as E Z ¼ 0, the covariance vanishes, i.e.
 Cov dE; dE2 ¼ 0. Consequently

 2 2 Var
q ¼ ðs12 fÞ Varð dEÞ þ ðs2 fÞ Var dE2 ð50Þ h
 i
    2 2 2 Since Var Y 2 ¼ E Y 2  E½Y ¼ E X 2  E X2 ¼   2 E X 4  Varð dEÞ we may insert the fourth moment of
 the normal distribution to find Var Y 2 ¼ 3r4  r4 ¼ 2r4 , where r is the standard deviation of X. Hence
 2 ð51Þ Var dE2 ¼ 2Varð dEÞ Substituting Eq. (51) into Eq. (50) and inserting the definition of s12 and f, we find that

estimations. The tuning-function will be optimized for non-disturbed conditions (Eb ¼ E), since this is the most general scenario. Nevertheless, also tuning for offsets between Eb and E is possible. Operational scenarios could then switch between tuning-functions depending on the current condition. The computation of t ¼ ^
q =
q over an altitude/flux-grid allows to evaluate the performance of the analytic approximations and at the same time to compensate for the discrepancy. Based on such a grid, tuningfactors may either be derived via interpolation or surfacefitting. While interpolation is potentially more accurate, surface fitting is faster and yields a generic functional description. Also it is the purpose of the tuning-function to compensate for general trends and not to match specific conditions. Using MATLAB, the computed tuning-factor

ð52Þ

where s1 and s2 are to be taken from the model coefficients (Tables 1, 2 or 3). ^
Using Eqs. (47) and (52) the average
qq of the LEO-grid statistics remains equal. This is not surprising for two reasons: first, for the test Eb was chosen to match E and second, dE is comparatively small with a standard-deviation of 11 sfu after seven days of propagation. If, however we introduce an offset between Eb and E, the second order description of dT ex is able to show its superior performance over the first order approximation. For Eb ¼ E þ 50 sfu, for example, the grid performance is: t2nd ¼ 0:84; t1st ¼ 1:24; rt2nd ¼ 0:19; rt1st ¼ 0:26; tmin2nd ¼0:54;tmin1 ¼0:66; st tmax2nd ¼ 1:66; tmax1st ¼ 1:73. The results of Eq. (47) are identical to those obtained with Eq. (41), since the second order description of dT ex is free of approximations. Despite the improvement of dT ex via a second order expression, Eqs. (47) and (52) still may result in overestimations for high solar flux levels and low altitudes, as well as underestimations for medium altitudes and low flux levels, as well as high flux levels and high altitudes. From a mathematical perspective, further enhancements would require improving the equation describing
q ðdT ex Þ. Another option, however, is to make use of the fact that the ratio of the analytically-derived relative density errors and the modelpropagated ones is constant for each parameter combination. This finding motivates the introduction of a tuningfunction that is able to further optimize the analytically computed relative density errors. 5.2.2. Altitude- and flux-dependent tuning factor As has been noted before, the discrepancy of the analytic relative flux error estimation to the model-derived errors is constant with time for each altitude and flux combination. In the following we hence develop a tuningfunction that compensates for the major over- and under-

grid was loaded and fifth-order polynomial surface-fits were computed. The grid points for NRLMSISE-00, as well as the results of the surface fits, are shown in Fig. 7. The obtained surface-equation is given in formula (53) (fit-accuracy: R2 ¼ 0:9801). Due to the multiplicative nature of the tuning-factor, its consideration in Eqs. (47) and (52) is straightforward. The performance with active tuning for a variety of different conditions, is shown in Fig. 8. After tuning, the performance over the grid is t ¼ 0:991; tmin ¼ 0:831; tmax ¼ 1:255; rt ¼ 0:046.
 t z;E  E ¼ 1:925  0:001588z þ 0:002612E  8:272  106 z2 þ3:25  105 zE  1:25  105 E2 þ 1:079  108 z3 þ1:029  108 z2 E  1:258  107 zE2 þ 5:122  108 E3 þ5:793  1012 z4  1:541  1010 z3 E þ 5:648  1010 z2 E2 8:038  1010 zE3 þ 8:011  1010 E4  7:754  1015 z5 þ1:17  1013 z4 E  4:658  1013 z3 E2 þ 7:964  1013 z2 E3 4:04  1013 zE4  7:987  1013 E5 ð53Þ

Fig. 7. Fifth order surface-fit of the tuning-factor grid (grey dots). The fit is very accurate and only deviates notably for very low fluxes at high altitudes.

F. Schiemenz et al. / Advances in Space Research 63 (2019) 3936–3952

3947

Fig. 8. Comparison of NRLMSISE-00 derived relative density errors (second subplot) and the proposed analytic estimations according Eqs. (47), (52) and (53) (third subplot). The second and fourth plot show the performance for a white-noise solar flux error process model. Case (a) depicts quiet conditions at 300 km altitude. At 500 km altitude, case (b) represents low solar activity with medium 81-day average flux levels. Case (c) shows high solar activity with a medium 81-day average at 700 km altitude. Case (d) finally depicts high, yet non-disturbed solar activity at the upper LEO region. In all cases the blue curve in the second subplot (a copy of the blue curve in the third subplot) matches the red one, i.e. the tuned second-order equations match the model-propagated relative density errors. (For interpretation of the references to colour in this figure legend, the reader is referred to the web version of this article.)

The min and max values are obtained for single parameterpairs for which the surface-fit has problems to capture the grid dynamics. The small standard deviation however shows that basically over the entire grid a close match with the model-propagated equations is obtained. The tuning factor function was derived under the assumption that E  E, as this is the most generic condition. Obviously, if this condition is not met, the performance is degraded. Nevertheless it is also possible to compute multiple tuning-functions for different relationships between Eb and E which can then be activated depending on the given inputs Eb and E. 6. Performance using DTM In contrast to NRLMSISE-00, the DTM models do not include switches to quickly simulate global and temporal average conditions. The only option is hence to modify the model-coefficients file in such a way that all variations other than those caused by the solar flux inputs add zero contribution to Gð LÞ. For this purpose the components T_Lat, T_kp, T_SLat, T_SASLat, T_NSLat, T_SAN-

SLat, T_DiAn, T_SDiAn, T_TDi, T_AMg, T_Lon and T_dPhas of the model parameters tt, h, he, ox, az2, o2, az, t0 and tp have to be set to zero. Also the flux coefficients T_flux(7:12) have to be set to zero, since they model cross-couplings between latitude and solar flux. 6.1. DTM-2012 We first present the results obtained using DTM-2012 for its similar domains of under- and overestimation with respect to NRLMSISE-00. On average, the untuned approximation of the relative density error matches the model-derived relative density error well (t ¼ 0:993; tmin ¼ 0:56; tmax ¼ 1:63; rt ¼ 0:28). Also for DTM-2012 the analytic approximations both over- and underestimate the model-derived densities (c.f. Fig. 9c and d). Overestimation tends to occur for low altitude and increasing flux-levels, whereas underestimation occurs for both medium altitude and low flux levels, as well as high altitudes and high flux levels. This is the same behavior as has been observed for NRLMSISE-00. Since the second order description of dT ex ð dEÞ is mathematically

3948

F. Schiemenz et al. / Advances in Space Research 63 (2019) 3936–3952

Fig. 9. Performance using DTM-2012. (a) Demonstrates the accurate second-order description of dT ex ðdEÞ. The tuning-factor grid and surface-fits are shown in (b). As is the case for NRLMSISE-00, the analytic expressions, on average, match the model-propagated relative density errors, however both over- and underestimation is possible ((c) and (d)).

identical to Eq. (41), the only contributor to over-/ underestimation is Eq. (41) itself. The second order description of dT ex ð dEÞ is very important for DTM-2012 due to the comparatively large c2 -component with respect to the other models, which causes a strong non-linear behavior of dT ex once dE becomes large (see Fig. 9a). To compensate for the maximum deviations of the approximated and model-derived relative density errors, a tuning-function was computed in the same way as for NRLMSISE-00. The fifth-order polynomial surface fit (Fig. 9b) achieved an accuracy of R2 ¼ 0:9731. The quality of the surface-fit is similar to NRLMSISE-00 (similar R2 values). The performance of the analytic expressions using DTM-2012 may be found in Appendix A for the same cases as given in Fig. 8. It generally matches the actual uncertainty well, even under conditions producing ”artifacts”. The grid-statistics after tuning is: t ¼ 0:996; tmin ¼ 0:836; tmax ¼ 1:458; rt ¼ 0:059. tmax represents a single peak for E ¼ 67 sfu and z ¼ 1200 km and is

due to the degraded surface-fit for this parameter combination. 6.2. DTM-2013 For DTM-2013 we show the results after the solar flux conversion to F10.7 units, i.e. the same flux range as for the other two models is covered. The same function as for DTM-2012 could be used to activate flux-sensitive global and temporal average behavior. In contrast to DTM-2012 and NRLMSISE-00, DTM-2013 underestimates the relative density error on average (t ¼ 0:82; tmin ¼ 0:53; tmax ¼ 1:06; rt ¼ 0:144). Nevertheless it shows the smallest standard deviation (only half of the other two models) and does not have the tendency of overestimation at low altitudes and high flux-levels. The strongest underestimation occurs, as is the case for the other two models, at medium altitudes and low flux values, as well as for high altitudes and high flux values. Due to the compar-

F. Schiemenz et al. / Advances in Space Research 63 (2019) 3936–3952

3949

Fig. 10. DTM-2013: Second-order description of dT ex ðdEÞ and fifth order surface fit of the tuning-factor grid.

atively small magnitude of c2 , the second-order description of dT ex may also be replaced with the first order approximation in case of small to medium disturbances (see Fig. 10a for the different representations of dT ex ð dEÞ). Among the models considered, DTM-2013 shows the highest dynamics in the tuning-factors. The fifth-order surface fit achieves R2 ¼ 0:9121 and is depicted in Fig. 10. With tuning the analytic expressions are able to accurately represent the relative density uncertainty obtained by the modelbased propagation. The results of the tuned analytic expressions for the same cases as given in Fig. 8 may be found in Appendix B. They agree perfectly in three out of four cases and slightly overestimate the modelpropagated uncertainties in the remaining one. The statistics after tuning are: t ¼ 0:995; tmin ¼ 0:836; tmax ¼ 1:447; rt ¼ 0:059. tmax is obtained again for E ¼ 67 sfu and z ¼ 1200 km and is due to the degraded surface-fit for this parameter combination, which is also visible in the front-left of Fig. 10b.

7. Conclusion The sensitivity of exospheric temperature to changes in the daily solar flux input has been derived for three major semi-empirical thermospheric models. These estimations have been used to derive analytic expressions of the resulting relative density error, which were compared to purely model-propagated relative density errors for the cases of a white noise and a random walk absolute solar flux error process model. The analytic equations describing the relative density error and its uncertainty were found to, on global and temporal average scales, agree well with the modelpropagated relative density errors. Nevertheless the possibility of over- and underestimation was identified for certain flux/altitude combinations. While the development of

an exact second order description of dT ex ð dEÞ was able to accurately predict changes in exospheric temperature for a variety of solar flux conditions, the over- and underestimation was determined to be caused by the derived relationship for
q ðdT ex Þ. For each parameter combination we found the ratio between the analytic estimates and the model-propagated relative density errors to be constant. Therefore an altitude- and flux-dependent compensationfactor was derived and shown to accurately tune the analytic equations. We stress that the modeled solar-flux errors were of typical magnitude for seven days of propagation. While this time-scale is common to orbit forecasting and collision avoidance scenarios, the relevant timescales for orbit estimation, possibly even simulations of past data with actual space-weather measurements, are commonly of the order of less than a day. Consequently small differences between the model-propagated uncertainties and the analytically-derived ones will not have the time to notably diverge. Therefore also the absolute accuracy of the analytic expressions is better for small propagation-times. This is true, independently of the model used. Among the studied ones NRLMSISE-00 and DTM-2012 showed very similar performance. DTM-2013 had the smallest standard deviation of the tuning-factors over the LEO-grid, however, on average, showed an underestimation of the density error. After tuning, the analytic approximations agree well with the model-propagated errors for a wide variety of typical LEO conditions. All in all, the presented analytic expressions allow to accurately estimate the resulting relative atmospheric density error/uncertainty due to absolute solar flux errors/ uncertainties. In a future work we plan to integrate the presented equations with the propagation of relative density errors to absolute in-track orbital position errors, as presented in Emmert et al. (2017). It is hoped that this approach will notably improve covariance realism.

3950

F. Schiemenz et al. / Advances in Space Research 63 (2019) 3936–3952

Appendix A. Tuned performance of analytic uncertainties (DTM-2012) Fig. A.11.

Fig. A.11. Comparison of DTM-2012 derived relative density errors and the proposed analytic estimations according Eqs. (42) and (43) including the tuning-function. The second and fourth plot show the performance for a white-noise solar flux error process model. Except for the case of 700 km altitude, all other cases are showing an excellent agreement between the model-propagated and the analytically computed uncertainties (the blue and red curves in the second subplots agree). The third case simulates high solar activity with a medium 81-day average level and is worthy of a separate discussion. Both the analytically calculated and the model-propagated errors and uncertainties show a strong shift towards the negative. The reason for this is that no Gaussian temperature errors arise from the random-walk Gaussian solar flux errors. This happens both in the model and in the analytic calculation. The reason is the comparatively large c2 -parameter of DTM-2012. This causes positive flux deltas to no longer result in positive dT ex . For this reason, the distribution of the T ex -error (and consequently also of the relative atmospheric density errors) experiences the same negative shift. This shift results in the 3r uncertainty to under-represent the actual uncertainty, since it is calculated from the non-Gaussian 500 members. In the analytic case, the 3r uncertainty is calculated via Eq. (52). Since it is based on the same shift in dT ex , it has to match the red curve in subplot three. Nevertheless an actual 3r uncertainty cannot be identified anymore due to the shift. Comparing the blue and red curves in subplot two, actually the analytically computed uncertainty represents the maximum deviations into the negative better than the red curve. Hence, even though it disagrees with the red curve, it still is an accurate measure of uncertainty. (For interpretation of the references to colour in this figure legend, the reader is referred to the web version of this article.)

F. Schiemenz et al. / Advances in Space Research 63 (2019) 3936–3952

3951

Appendix B. Tuned performance of analytic uncertainties (DTM-2013) Fig. B.12.

Fig. B.12. Comparison of DTM-2013 derived relative density errors and the proposed analytic estimations according Eqs. (42) and (43) including the tuning-function. The second and fourth plot show the performance for a white-noise solar flux error process model. Despite the slight overestimation of uncertainty for E ¼ 80 sfu, z ¼ 700 km, the overall accuracy of the tuned equations is noteworthy (the blue and the red lines in the second subplots of the figures match). (For interpretation of the references to colour in this figure legend, the reader is referred to the web version of this article.)

3952

F. Schiemenz et al. / Advances in Space Research 63 (2019) 3936–3952

References Bowman, B., Tobiska, W.K., Marcos, F., Huang, C., Lin, C., Burke, W., 2008. A new empirical thermospheric density model JB2008 using new solar and geomagnetic indices. In: AIAA/AAS Astrodynamics Specialist Conference and Exhibit, p. 6438. Bruinsma, S., 2013. The semi-empirical thermosphere model DTM2012. 6th European Conference on Space Debris, vol. 723. Bruinsma, S., 2015. The DTM-2013 thermosphere model. J. Space Weather Space Clim. 5, A1. de Wit, T.D., Bruinsma, S., 2017. The 30 cm radio flux as a solar proxy for thermosphere density modelling. J. Space Weather Space Clim. 7, A9. Dudok de Wit, T., Bruinsma, S., 2011. Determination of the most pertinent EUV proxy for use in thermosphere modeling. Geophys. Res. Lett. 38 (19). Emmert, J., 2015. Thermospheric mass density: a review. Adv. Space Res. 56 (5), 773–824. Emmert, J., Byers, J., Warren, H., Segerman, A., 2014. Propagation of forecast errors from the sun to leo trajectories: how does drag uncertainty affect conjunction frequency? In: AMOS Conference. Emmert, J., Warren, H., Segerman, A., Byers, J., Picone, J., 2017. Propagation of atmospheric density errors to satellite orbits. Adv. Space Res. 59 (1), 147–165.

Marcos, F., 2006. New satellite drag modeling capabilities. In: 44th AIAA Aerospace Sciences Meeting and Exhibit, p. 470. Picone, J.M. et al., 2002. NRLMSISE-00 empirical model of the atmosphere: statistical comparisons and scientific issues. J. Geophys. Res.: Space Phys. 107 (A12). Poore, A.B., Aristoff, J.M., Horwood, J.T., Armellin, R., Cerven, W.T., Cheng, Y., Cox, C.M., Erwin, R.S., Frisbee, J.H., Hejduk, M.D., et al., 2016. Covariance and Uncertainty Realism in Space Surveillance and Tracking. Tech. rep., Numerica Corporation. Sagnieres, L., Sharf, I., 2017. Uncertainty characterization of atmospheric density models for orbit prediction of space debris. In: 7th European Conference on Space Debris. Schiemenz, F., Utzmann, J., Kayal, H., 2019. Survey of the operational state of the art in conjunction analysis. CEAS Space Journal. https:// doi.org/10.1007/s12567-019-00242-2. Storz, M.F., Bowman, B.R., Branson, M.J.I., Casali, S.J., Tobiska, W.K., 2005. High accuracy satellite drag model (HASDM). Adv. Space Res. 36 (12), 2497–2505. Tapping, K., 2013. The 10.7 cm solar radio flux (F10.7). Space Weather 11 (7), 394–406. Wilkins, M., Alfriend, K., 2000. Characterizing orbit uncertainty due to atmospheric uncertainty. In: Astrodynamics Specialist Conference, p. 3931.