Journal of Geodesy https://doi.org/10.1007/s00190-018-1165-8

ORIGINAL ARTICLE

Towards thermospheric density estimation from SLR observations of LEO satellites: a case study with ANDE-Pollux satellite Francesca Panzetta1 · Mathis Bloßfeld1 · Eren Erdogan1 · Sergei Rudenko1 · Michael Schmidt1 · Horst Müller1 Received: 14 December 2016 / Accepted: 29 May 2018 © Springer-Verlag GmbH Germany, part of Springer Nature 2018

Abstract The present contribution investigates the possibility to obtain thermospheric neutral density estimates using satellite laser ranging (SLR) observations of low Earth orbiters (LEOs). This approach is based on the analysis of the satellite atmospheric drag, driven by the fact that the drag force is the largest non-gravitational perturbation acting on LEOs. Due to the uncertainty of current thermospheric models, it is the main error source in the LEO orbit determination process. Moreover, the drag is physically related to the thermospheric density distribution, the interaction of the satellite surface with the surrounding thermosphere and thermospheric winds. For this investigation, a spherical satellite called “Atmospheric Neutral Density Experiment-Pollux” (ANDE-P) developed by the Naval Research Laboratory (USA) is adopted as a case study. The satellite flew at the very low altitude of about 350 km. The most important perturbing acceleration at this altitude, the atmospheric drag, is easier to model for a spherical satellite like ANDE-P than for a satellite of complex geometry. A precise orbit determination of ANDE-P was performed with the DGFI Orbit and Geodetic parameter estimation Software (DOGS) over a period of 49 days (16 August 2009 until 3 October 2009) using SLR observations to this satellite and different thermospheric models. In total, we tested four thermospheric models, namely CIRA86, NRLMSISE00, JB2008 and DTM2013. Correspondingly, scale factors of these reference models are estimated with a 6-h resolution. The results confirm that the estimation of force model parameters from SLR measurements of the ANDE-P satellite is sensitive to differences in the density distributions provided by different models. As a consequence, information on the discrepancies between the various models and the true density can be derived from SLR measurements. Moreover, it is found that SLR observations to LEO satellites at very low altitudes are capable to estimate corrections to (scale factors of) the integrated thermospheric density if all other perturbing accelerations are modelled with sufficient accuracy. We derived time series of estimated scale factors of the thermospheric densities provided by the models and obtained the following mean values during the processed period of time at the ANDE-P altitude: 0.65 ± 0.26 for CIRA86, 0.65 ± 0.25 for NRLMSISE00, 0.79 ± 0.24 for DTM2013, and 0.89 ± 0.27 for JB2008. This suggests that all models overestimate the true thermospheric density along the ANDE-P trajectory during the processed period to a certain extent. The thermospheric densities need to be scaled downwards to fit ANDE-P SLR observations with JB2008 requiring the least amount of scaling. Keywords Thermosphere · SLR · ANDE-P · Atmospheric drag · Balistic coefficient

1 Introduction The Earth’s upper atmosphere comprising the thermosphere and the ionosphere exhibits a dynamically coupled non-linear system in terms of chemical and physical processes. The

B 1

Mathis Bloßfeld [email protected] Deutsches Geodätisches Forschungsinstitut of the Technische Universität München (DGFI-TUM), Arcisstr. 21, 80333 Munich, Germany

system also interacts with the magnetosphere as well as the lower atmosphere (Heelis 2011). Several stand-alone or coupled models have been developed to reveal the behaviour of atmospheric target parameters and their interactions (e.g. the neutral and charged particle density of the thermosphere and ionosphere) from different perspectives which are, for instance, based on pure physical or (semi-) empirical models as well as data assimilative approaches combining available models with new sets of observations. For a review of atmospheric models we refer to, e.g., Akmaev (2011), Qian and Solomon (2012), Emmert (2015), and references therein. In

123

F. Panzetta et al.

this contribution, empirical thermospheric models are taken into account and are investigated by assimilating SLR observations through a precise orbit determination (POD). The thermospheric neutral density plays a crucial role within the equation of motion of low (altitude) Earth orbiting (LEO) objects since the drag force is one of the largest and worst modelled non-gravitational perturbations acting on these objects. It is a function of the thermospheric integral density that is still not precisely known. Besides, the density estimation is of critical consideration for re-entry operations, manoeuvre planning, collision avoidance, POD, and satellite lifetime planning (Doornbos 2011; Emmert 2015). In the past, accelerometer instruments on, e.g., the CHAllenging Minisatellite Payload (CHAMP, mission time: 2000–2010; Bruinsma and Biancale 2003; Liu et al. 2005) satellite mission and on the Gravity Recovery And Climate Experiment (GRACE, mission time: 2002–2017) satellite mission have provided thermospheric density data at an unprecedented accuracy and resolution (Sutton et al. 2007; Doornbos et al. 2010; Doornbos 2011). Moreover, the gradiometer measurements of the Gravity field and steady-state Ocean Circulation Explorer (GOCE, mission time: 2009–2013; Doornbos et al. 2013) satellite mission together with the use of precision orbit ephemerides (POE) as observations from the TERRASARX satellite mission (since 2007; McLaughlin et al. 2012) and from the GRACE satellites (Lechtenberg et al. 2013) will further improve the thermospheric density database. The major limiting factor for the accuracy of thermospheric density models derived from POE and accelerometer data is that uncertainties in the complex satellite shape and dimensions and the computation of the satellite-specific drag coefficient directly impact the density estimates. For example, Doornbos (2011) compared different geometrical models for CHAMP and GRACE and found differences between a common simple panel model and the fully 3D CAD (computer-aided design) model ANGARA (Fritsche et al. 1998) of up to 37% for surfaces along certain axes of the satellites. Therefore, the intention of this paper is to make use of LEO satellites with much simpler shape like purely spherical ones. Such satellites have normally no accelerometers on-board. This means other types of observations have to be taken into account. For completeness, an investigation of Bowman (2004) has to be mentioned who used historical radar observations to 13 simply shaped satellites (cubes and spheres) in order to characterize mainly semi-annual thermospheric density variations between 1970 and 2002 at altitudes between 200 and 1100 km. Bowman and Moe (2005) extended the previous investigation determining on a daily basis atmospheric temperature and density corrections to a reference thermospheric model using up to 79 calibration satellites at the altitude range between 150 and 500 km. Satellite laser ranging (SLR) is a geodetic tracking technique which can be used for POD of near-Earth satellites.

123

Fig. 1 The ANDE-2 spherical satellites Castor (left) and Pollux (right), image credit: NRL

SLR provides highly accurate travel time measurements of laser pulses reflected at laser retro-reflector arrays (RRA) mounted on the satellite surface which have been emitted from telescopes on the Earth’s surface. Due to the high precision and accuracy of such measurements, SLR observations are highly sensitive to any perturbing acceleration acting on a satellite. The scope of this investigation is to study the ability of SLR measurements to very low orbiting satellites to derive scale factors of thermospheric densities. One of the first experiments using SLR for thermospheric density determination was performed by Gaposchkin and Coster (1988). They used precise radar tracking data to two spherical satellites and laser ranging data to a third spherical satellite to evaluate the thermospheric models Jacchia77, MSIS86 (MSIS: Mass Spectrometer Incoherent Scatter), CIRA72 (CIRA: COSPAR International Reference Atmosphere), and DTM (Drag Temperature Model). Doornbos et al. (2007) used SLR observations to evaluate the thermospheric densities in the orbit determination of ERS-2 and Envisat at about 800 km altitude, comparing the results for several thermospheric models. Jeon et al. (2011) estimated upper atmosphere (> 800 km) mass density from Starlette using SLR observations. In their investigation, the authors validated the thermospheric models MSIS86 and NRLMSISE00. Also Bruinsma et al. (2012) used SLR observations to Stella and Starlette to evaluate the DTM-09 thermospheric model at altitudes of around 800 km. In this paper, we focus on much lower altitudes using a specially designed satellite mission. The satellite “Atmospheric Neutral Density Experiment-Pollux” (ANDE-P) is one of two spherical satellites of the ANDE-2 mission (see Fig. 1) launched in 2009 as a follow-on mission of ANDERR (launched in 2006; Nicholas et al. 2007). The ANDE-2 low-cost mission of Naval Research Laboratory (NRL) had the aim to monitor the thermospheric density from an initial orbital altitude of about 350 km downwards (see Fig. 2). The other orbital characteristics of the ANDE-2 satellites (“Pollux” and “Castor”) are a mean inclination of 51.6◦ and an eccentricity of 0.0007, provided together with further

19125 19121 5892 5784 1485 1442 957 835 817 798 681 425 400 400 400 350 350 250

Etalon-1 Etalon-2 LAGEOS-1 LAGEOS-2 Ajisai LARES Starlette Westpac BLITS Stella Larets SpinSat GFZ-1 ANDE-RR-A ANDE-RR-P ANDE-Pollux ANDE-Castor LRE

1970

1980

thermospheric drag modelled

satellite mean altitude [km]

Towards thermospheric density estimation from SLR observations of LEO satellites: a case study…

1990

2000

2010

2020

Fig. 2 Mission lifetime for different spherical SLR satellites. The grey box highlights the missions where thermospheric drag has to be considered in the POD. The blue bar indicates the data used in this investigation whereas the yellow bars indicate the individual mission lifetimes

information at https://directory.eoportal.org/web/eoportal/ satellite-missions/a/ande-2. Figure 2 shows a constellation of currently used spherical satellites within the International Laser Ranging Service (ILRS; Pearlman et al. 2002) community. It clearly highlights the potential of SLR observations to supplement existing observations of the integrated neutral thermospheric density. Moreover, it emphasizes that this investigation is only a case study to derive algorithms for the scaling of the density which might be applicable to other satellite missions in future studies. In particular, ANDE-P is a spherical micro-satellite with a diameter of 48.26 cm and a mass of 27.442 kg, equipped with 30 optical retro-reflectors for SLR applications (Nicholas et al. 2009). The high area-to-mass ratio of this satellite (0.006666 m2 /kg) caused its fast orbit decay which results in a short mission lifetime from August 2009 until March 2010. Recently, Lechtenberg (2015) used pre-processed SLR range data of different ANDE satellites for the determination of integral thermospheric densities and compared them with CHAMP and GRACE derived estimates. Another investigation was done by Flanagan (2015) who used the SLR orbit predictions for ANDE-C and ANDE-P in order to obtain integrated thermospheric densities. In this paper, for the first time we present the estimation of scale factors of the integrated thermospheric densities using actual SLR observations (not pre-processed SLR range data or SLR orbit predictions) to the ANDE-P satellite in combination with a full POD, taking advantage of the high precision

and accuracy of the SLR technique in the orbit determination of spherical satellites. A further benefit of this approach is that we can monitor all possible correlations between the estimated parameters in order to have a deep understanding of the system behaviour. The rest of the paper is organized as follows. In Sect. 2, the theoretical background of thermospheric drag modelling is explained in detail. Section 3 describes the solution setup of our experiment. In Sect. 4, the obtained results of the experiment and a refined thermospheric modelling are presented. Section 5 provides our results on stability and reliability of the obtained results. Finally, Sect. 6 concludes the paper and gives a short outlook on future work.

2 Theoretical background The assumptions and considerations described in detail in the following (e.g. those for the drag coefficient CD ) ensure that, in our interpretation of the results, the finally obtained (estimated) scale factor of the drag acceleration aaero can be interpreted as a pure scaling of the integrated neutral thermospheric density. We are aware of the fact that any changes would result in another scale factor. Therefore, we did all considerations described in the following in all conscience. SLR is a space geodetic technique measuring the two-way travel time of a short laser pulse from a ground station to an Earth orbiting satellite equipped with retro-reflectors. SLR is, among current space geodetic techniques, the most accurate technique available to determine the geocentric position of near-Earth satellites since it provides instantaneous measurements of the satellite range of millimetre level precision (distance between the satellite retro-reflector and the observation telescope) from a global network of stations. In general, the one-way range ρ can be expressed in terms of the roundtrip travel time Δτ of light as ρ+ =

1 cΔτ, 2

(1)

where c represents the speed of light in vacuum and  is the measurement error. However, the obtained range ρ is not equal to the purely geometrical distance between the positions of the station and the satellite. Thus, some corrections have to be applied. Within the DGFI Orbit and Geodetic parameter estimation Software (DOGS), developed at Deutsches Geodätisches Forschungsinstitut of the Technische Universität München (DGFI-TUM), the computed one-way range ρ is modelled as (Gerstl 1997; Bloßfeld 2015) ρ +  = r sat (tM + Δt) − r sta (tM + Δt) + Δρ + ctrop (1 + δr ) + crel + csta + cmasc + cmesc , (2)

123

F. Panzetta et al.

where r sat is the position vector of the satellite centre of mass in the Geocentric Celestial Reference System (GCRS), tM is the approximated time epoch of reflection of the laser pulse at the satellite, Δt is the time bias of the measurement, r sta is the position vector of the station converted from the International Terrestrial Reference System (ITRS) to the GCRS using the Earth orientation parameters (EOP), Δρ is the range bias of measurement, ctrop is the tropospheric range correction, δr is the bias of tropospheric refraction, crel is the relativistic range correction, csta is the station-dependent SLR correction, cmasc is the satellite-specific centre of mass correction, and cmesc is the SLR array-dependent correction (e.g. phase centre offset in measurement direction). The position of the satellite r sat at a specified epoch is affected by gravitational and non-gravitational accelerations aG and aNG , respectively. In general, the total acceleration asat acting on a LEO satellite in the GCRS is modelled as a sum of the following contributions (Montenbruck and Gill 2000): aG

   asat = akep + ankep + a3body + atid + antid aNG

   + aaero + asrp + aerp + aother +aemp ,

(3)

where besides the acceleration due to the Keplerian (central) gravitational field akep , the gravitational perturbations include the non-Keplerian (non-central or non-spherical) gravitational acceleration ankep , the acceleration due to third body effects like the Sun, the moon, and the planets a3body , the tidal acceleration atid due to solid Earth and ocean tides, and the acceleration due to non-tidal loading effects a ntid (e.g. atmospheric, hydrological, oceanic). The non-gravitational perturbations comprise the aerodynamic acceleration aaero due to the atmosphere, the acceleration due to solar radiation pressure asrp , the acceleration due to radiation pressure aerp caused by Earth albedo and infrared radiation as well as other perturbing accelerations aother (e.g. satellite thermal emission, Earth’s magnetic fields). In addition, one cycleper-revolution (CPR) empirical accelerations aemp absorb so far non-modelled accelerations. One of the major problems in satellite orbit determination is the assessment of the non-gravitational perturbations since they are essentially generated by complex interactions of the satellite surface with the surrounding space environment and depend on uncertainly determined physical parameters. Among the non-gravitational perturbations, the aerodynamic acceleration is the largest one acting on LEO satellites at altitudes lower than 1000 km (thermosphere). Thus, the aerodynamic perturbation represents the main error source within LEO satellite POD. All spherical satellites like, e.g., ANDE-P are ideally affected only by one component of the total aerodynamic

123

perturbation, namely the atmospheric drag. The drag acceleration is always pointing in opposite direction to the relative velocity vector v rel of the satellite w.r.t. the atmosphere and causes a (continuous) deceleration of the satellite. As a consequence, the spacecraft looses kinetic energy due to the friction which itself results in a decrease in the semi-major axis (orbit decay) and the eccentricity (orbit becomes more circular). Moreover, it has some periodic effects on the other orbital elements (Vallado 2007). For a perfectly spherical satellite, the other two components of the aerodynamic perturbation, namely the lift and side accelerations, disappear because the normal of a sphere (usually modelled as a plane) is always parallel to the drag direction. Therefore, for a spherical satellite like ANDE-P, we can consider these two components to be negligible and assume that the approximation aaero  adrag is valid (Doornbos 2011). In general, the drag acceleration can be modelled as adrag = −

1 Aref 2 CD ρM vrel uˆ D , 2 m

(4)

where uˆ D is the drag unit vector computed as v rel /||v rel ||, Aref is the effective cross-sectional area of the satellite interacting with the atmosphere, m is the satellite mass, CD is the dimensionless aerodynamic drag coefficient describing the interaction of the atmosphere with the satellite surface, and ρM is the thermospheric neutral density (provided by a thermospheric model). In the literature, the term Aref /m · CD is called the satellite-specific ballistic coefficient. To a first approximation, the velocity of the atmosphere can be identified with the Earth co-rotation velocity, given as ωe × r sat , where ωe is the Earth angular velocity vector in the GCRS. The atmospheric co-rotation velocity is proportional to the cosine of the latitude, ranging from a maximum of 483–502 m/s at the equator (at altitudes of 250–500 km) to zero at the poles. This means the deviation from co-rotation becomes significant in the polar regions and, furthermore, at high solar and geomagnetic activity (Doornbos 2011). Therefore, in order to include a more realistic velocity of the atmosphere, it is important to consider also the wind velocity v w w.r.t. a co-rotating atmosphere. Atmospheric winds can reach velocities of hundreds of m/s in the thermosphere with peaks exceeding 600 m/s in the polar regions (Lühr et al. 2007). The wind velocity is usually defined in local coordinates (East, North, Vertical) corresponding to zonal, meridional, and vertical winds, respectively. The transformation of these thermospheric winds from local to inertial coordinates requires the rotation matrix REI from Earth-fixed to inertial coordinates using the EOP. In addition, the rotation matrix RLE from a local to an Earth-fixed coordinate system using satellite’s latitude and longitude (Montenbruck and Gill 2000) is necessary. Finally, the relative velocity of the satellite w.r.t. the atmosphere is computed as the differ-

Towards thermospheric density estimation from SLR observations of LEO satellites: a case study…

ence between the satellite’s orbital velocity v sat and the sum of co-rotation and wind velocities as (Doornbos 2011) v rel = v sat − (ωe × r sat + REI RLE v w ).

(5)

As shown in Eq. (4), the drag acceleration results from the product of four main terms: the area-to-mass ratio Aref /m, the drag coefficient CD , the thermospheric density ρM , and the relative velocity vrel of the satellite w.r.t. the atmosphere. Since none of these quantities is known exactly, an error in the atmospheric drag computation can be due to an error in the determination of any of these four parameters or is very probably due to a combination of their errors [see case (D) in Fig. 10]. A discussion of the inaccuracies caused while modelling satellite drag is given by Marcos (2003); Marcos et al. (2010). For spherical satellites, the area-to-mass ratio is probably the most accurate parameter because the satellite mass is well-known and constant since no mobile parts are mounted on the (passive) satellite. Moreover, the crosssectional area is constant and easy to determine for a sphere. Concerning the contributions to the relative velocity, the orbital and co-rotation velocities are known at an accuracy, much higher than the wind velocity (Doornbos 2011), which has an uncertainty ranging from 10 to 60 m/s as for the Horizontal Wind Model 2014 (HWM14; Drob et al. 2015). On the other hand, the uncertainties in the thermospheric density and the drag coefficient are significant under all conditions of solar and geomagnetic activity and at all latitudes and thermospheric altitudes. Therefore, for satellites with accurately known shape and attitude, the main sources of error in the drag computation are the thermospheric density and the drag coefficient. This is why spherical satellites can be seen as so-called calibration satellites (Marcos 2003) since their satellite-specific parameters (such as area-to-mass ratio) are well known and can be modelled very accurately, as compared to satellites with a complex shape or mass changes due to, e.g. fuel consumption. At the very low altitude of ANDE-P of about 350 km (during September 2009) the thermospheric density is assumed to be the least accurate parameter (Marcos et al. 2010). The density is modelled by a complex function depending on the instantaneous satellite position (altitude, latitude, longitude), the local solar time, solar and geomagnetic activity indices and the harmonics of the day of year, representing seasonal variations. In order to minimize any drag error due to a wrong drag coefficient, it has to be physically modelled as accurately as possible. Historically, a fixed value CD = 2.2 was usually adopted for all satellites with compact shape (Cook 1965).

Nowadays, the drag coefficient is computed analytically, based on physical principles described in a gas–surface interaction (GSI) model which studies how gas particles exchange energy and momentum with the surface of an object. In this sense, the drag coefficient is a macroscopic parameter enclosing the physical properties of the GSI model since it depends on the material and temperature of the satellite surfaces, shape and orientation of the satellite surfaces w.r.t. atmospheric particles and composition of the thermosphere and its temperature. A gas particle, colliding with a spacecraft surface, can be either adsorbed and then potentially re-emitted or directly reflected (specularly or diffusely). These interactions are described by a fundamental parameter in the GSI models, namely the energy accommodation coefficient α, formally defined as the average fraction of energy lost by the molecules impacting the satellite’s surface (Cook 1965) α=

Ei − Er , Ei − Ew

(6)

where E i and E r are, respectively, the average kinetic energies of the incident and the re-emitted molecules, while E w is the average kinetic energy of a molecule re-emitted with a velocity corresponding to the surface (wall) temperature Tw . In the case of a specular reflection, i.e. α = 0 (no accommodation), there is no energy exchange between the molecules and the surface (E i = E r ). For an entirely diffuse reflection, α = 1 (complete accommodation) since the incident molecules and the surface reach thermal equilibrium before the molecules are re-emitted (E r = E w ). The case 0 < α < 1 is referred to as quasi-specular reflection. Another important parameter for the drag coefficient computation is the molecular speed ratio sk which is defined for each thermospheric species k as (Cook 1965) sk =

vrel , vm,k

(7)

being vm,k the most probable molecular velocity at any altitude derived from gas kinetic theory with a Maxwell– Boltzmann distribution of velocity  vm,k =

2RT∞ . mk

(8)

In Eq. (8), R is the universal gas constant, m k is the molecular mass of the gas species and T∞ is the atmospheric translational temperature which can be obtained from thermospheric models. Several simplifying assumptions, briefly reported here, are adopted in order to analytically compute the physical drag coefficient for the ANDE-P satellite.

123

F. Panzetta et al.

C D,k =

4sk4 + 4sk2 − 1 2sk4

2s 2 + 1 2 erf (sk ) + √k 3 e−sk π sk

√  2 π Tk,r + , 3sk T∞

(9)

where erf (sk ) is the error function defined as

Fig. 3 Completely diffuse reflection of the gas particles from the satellite’s surface; fundamental assumption of Sentman’s model (Mehta et al. 2014)

1. The fundamental assumption, common to all GSI models, is the condition of “free molecular flow” of the atmosphere. At altitudes above approximately 140 km, the mean free path of the molecules is much larger than the characteristic satellite dimension. This flow regime is dominated by molecule-surface collisions with almost no intermolecular collisions and is characterized by a Knudsen number K n > 10 (King-Hele 1987). 2. In our case, the adopted GSI model is Sentman’s model (Sentman 1961), valid at altitudes lower than 500 km where the thermosphere composition is rich in atomic oxygen. 3. The model is based on the assumption of thermal flow which describes the incident flow at every point at the satellite’s surface as a superposition of the Maxwell– Boltzmann molecular velocity distribution and the incident velocity vector of the atmosphere relative to the spacecraft. 4. The gas molecules experience a fully diffuse reflection (see Fig. 3) with complete accommodation (α = 1). This effect is caused by a layer of adsorbed atomic oxygen that is supposed to cover the satellite’s surface at altitudes lower than 500 km (Moe and Moe 2005). 5. The re-emitted particles have a Maxwell–Boltzmann velocity distribution. 6. Usually, the satellite’s surface temperature Tw varies in a range between 100 and 500 K. Due to the fact that the sensitivity of the drag coefficient to this parameter is quite small for various satellite geometries (Mehta et al. 2014), we assume an average value of Tw = 300 K for this investigation. It has to be pointed out here once again that these assumptions are crucial for the analytical computation of CD . Any obtained CD relies (to different extents) on these assumptions. The physical drag coefficient CD of a spherical satellite can be computed for each thermospheric species k as (Pilinski et al. 2010)

123

2 erf (x) = √ π



x

e−t dt. 2

(10)

0

The kinetic temperature Tk,r of the reflected particles can be expressed as (Pilinski et al. 2010) Tk,r = Tk,i (1 − α) + αTw =

m k vi2 (1 − α) + αTw , 3kB

(11)

being a function of the energy accommodation coefficient α, the satellite surface temperature Tw and the kinetic temperature Tk,i of the incident particles, where vi is the incident velocity that can be approximated as vi = vrel for all the gas species and kB is the Boltzmann constant. Therefore, Sentman’s analytical solution for computing the physical drag coefficient of a spherical satellite is given by Eq. (9) with Tk,r = Tw since the accommodation coefficient is assumed to be α = 1. Finally, the total drag coefficient CD for the considered mixture of thermospheric species is given by the sum over the thermospheric constituents of the partial drag coefficients CD,k , weighted by the relative mass densities (Walker et al. 2014)  k ρk C D,k , (12) CD =  k ρk where ρk is the mass density of each thermospheric species k. As an example, Fig. 4 shows the total drag coefficient computed using Sentman’s model with the atmospheric translational temperature T∞ and the relative mass concen tration ρk /( k ρk ) of each gas species provided by the model JB2008 (Bowman et al. 2008) over the processed ANDE-P time period. The average value for the modelled CD of ANDE-P is 2.1150 (its standard deviation of 0.0012 is much less than 1% of the value) which is comparable with the value obtained for ANDE-RR satellites (Nicholas et al. 2007).

3 Description of the experiment The ANDE-P satellite is difficult to track because fastrotating SLR telescopes are necessary (low altitude means large angular velocity and short fly-over time). We use SLR

Towards thermospheric density estimation from SLR observations of LEO satellites: a case study…

Fig. 4 Physical drag coefficient computed using Sentman’s model with input from the model JB2008 (Bowman et al. 2008) over the processed ANDE-P time period

data available for ANDE-P from the ILRS. That’s why a, rather limited amount of SLR observations (see Fig. 5) is available for this satellite. The month with the largest number of observations during the ANDE-P mission lifetime is September 2009 which has in total about 1400 normal points (mean values over a certain number of observations). In this investigation, we processed 49 days from 16 August 2009 until 3 October 2009 (blue bar in Fig. 2 and red box in Fig. 5) with a 3.5-day arc resolution. The temporal distribution of ANDE-P SLR observations used per arc is reported in Fig. 6, and the number of SLR tracking stations per arc is shown in Fig. 7, respectively. In fact, SLR observations of ANDEP have globally a very sparse geographic distribution (see Fig. 8) being mostly concentrated in the USA, Australia, Europe, and Asia, with a few observations in South Africa and at Hawaii. The impact of this sparse distribution on the experiment described in this paper is discussed in the next section. The approach selected to estimate the total thermospheric densities from ANDE-P SLR observations consists in performing a fully dynamic POD of the ANDE-P satellite using DOGS. All a priori models are based on the recommendations of the IERS Conventions 2010 (Petit and Luzum 2010). In the analysis, the Earth’s gravitational field model EIGEN-6S (Förste et al. 2011) up to degree/order 120 and the ocean tide model EOT11a (Savcenko and Bosch 2012) up to degree/order 30 are used. A summary of the different dynamical models adopted for ANDE-P orbit determination is given in Table 1. During the processing of ANDE-P satellite data we did not account for the satellite thermal radiation forces (Yarkovsky effect and Yarkovsky–Schach effect), which are caused by the satellite rotational dynamics (Andrés et al. 2004), since no orientation and precise rotation data were available to us. In fact, the satellite was deployed with an intentional initial spin of 2.4 deg/s (0.40 revolutions per minute) (Nicholas et al. 2009), which was necessary for the thermal control of on-board sensors. According to the fully dynamic POD approach, the satellite’s equation of motion can be numerically integrated

Fig. 5 Temporal distribution of SLR observations to ANDE-P during the different months of the mission lifetime

Fig. 6 Temporal distribution of SLR observations to ANDE-P used at the different 3.5-day arcs of the processed data period. The red line indicates the number of unknowns during the POD process

Fig. 7 Number of SLR tracking stations during the different 3.5-day arcs of the processed data period for ANDE-P

Fig. 8 Ground track distribution of SLR observations to ANDE-P during the processed data period

123

F. Panzetta et al. Table 1 Dynamical models adopted for the SLR POD of ANDE-P Dynamical model

Description

Earth gravity field

EIGEN-6S up to degree/order 120 (Förste et al. 2011)

Solid earth tides

IERS-2010 Conventions

Ocean tides

EOT11a up to degree/order 30 (Savcenko and Bosch 2012)

Ocean loading tides

EOT11a (Savcenko and Bosch 2012)

Atmospheric tides

IERS-2010 Conventions

Atmospheric loading

IERS-2010 Conventions

tides

(Ray and Ponte 2003)

Pole tides

IERS-2010 Conventions

Lunar gravity field

(Konopliv et al. 2001) up to degree/order 150

General relativistic

Schwarzschild, deSitter,

correction

Lense-Thirring (IERS-2010 Conventions)

Third-body effect

DE-421: Mercury, Venus, Mars, Jupiter, Saturn, Sun, Moon (Folkner et al. 2008)

Non-tidal gravity perturbation

None

Solar radiation pressure

Cannonball model

Earth radiation pressure

Albedo and infrared (Knocke et al. 1988)

Atmospheric drag

Models alternately used: CIRA86, NRLMSISE00, DTM2013, JB2008 (see Table 2)

Thermal radiation

None

Table 2 Adopted empirical thermospheric models Model

References

Species

JB2008

Bowman et al. (2008)

H, He, O, N2 , O2 , Ar

DTM2013

Bruinsma (2015)

H, He, O, N2 , O2

CIRA86

Hedin et al. (1988)

H, He, O, N2 , N, O2 , Ar

NRLMSISE00

Picone et al. (2002)

H, He, O, N2 , N, O2 , Ar

starting from an initial satellite state vector (position, velocity and other possible model parameters). Gravitational, nongravitational, and empirical force models are used to integrate the satellite orbit. Since the initial satellite state vector is usually not accurate enough and the force models might contain errors, the estimated satellite state vector is adjusted within some iterations in order to minimize the observation residuals in a least-squares sense. The integrator used within DOGS is based on a Gauss-Jackson predictor-corrector algorithm (Jackson 1924) of seventh order with a relative error bound of 10−9 . The step size used for the orbit integration is 30 s. Four empirical thermospheric models, reported in Table 2, are considered in our ANDE-P solution set up, namely

123

Fig. 9 Thermospheric densities of the models NRLMSISE00 (blue), CIRA86 (green), DTM2013(orange) and JB2008 (magenta)

JB2008, DTM2013, CIRA86, and NRLMSISE00. They include different gas species and provide thermospheric temperature and density as functions of the instantaneous satellite position (altitude, latitude, longitude), the local solar time, solar and geomagnetic storm indices and harmonics of the year’s fraction. As an example, the mean thermospheric densities of each species during the processing period provided by the model NRLMSISE00 in kg/m3 are 5.63 × 10−16 for H, 2.81 × 10−14 for He, 3.28 × 10−12 for O, 4.17 × 10−13 for N2 , 3.76 × 10−14 for N, 1.61 × 10−14 for O2 , and 4.34 × 10−17 for Ar. The total thermospheric densities of these models are shown in Fig. 9, where JB2008 shows in general the lowest densities: the mean values are, respectively, 3.78 × 10−12 kg/m3 for NRLMSISE00, 3.80 × 10−12 kg/m3 for CIRA86, 2.87 × 10−12 kg/m3 for DTM2013, and 2.46 × 10−12 kg/m3 for JB2008. Large density differences are visible between JB2008 and both CIRA86 and NRLMSISE00, while smaller differences (and not over all the processed period) exist between JB2008 and DTM2013. In the NRLMSISE00 model, the presence of the so-called anomalous oxygen (consisting of “hot” atomic oxygen and ionospheric atomic oxygen ions O+ ) is taken into account as a significant contributor to the satellite drag during summer at high latitudes and altitudes above 500 km (Picone et al. 2002). Since ANDE-P is orbiting the Earth at an approximative altitude of 350 km, this thermospheric constituent is not included in our thermospheric model used for the ANDE-P analysis. In addition to each thermospheric model the HWM14 is used which provides zonal and meridional winds. To account for the different magnitudes of the thermospheric densities, as shown in Fig. 9, Eq. (4) is extended by a scale factor f s : adrag = −

Aref 1 2 fs CD ρM vrel uˆ D , 2 m

(13)

As pointed out in Sect. 2, the largest error sources of the drag perturbation for spherical LEO satellites are the drag

Towards thermospheric density estimation from SLR observations of LEO satellites: a case study…

Fig. 10 Different possibilities for applying the estimated scale factor in our analysis

coefficient and the thermospheric density. This means the estimated scale factor f s compensates either the error of one of these two quantities or of a combination of both. Therefore, some considerations become necessary in order to decide which of the two parameters might be scaled by f s . Figure 10 shows the possibility to either model the drag coefficient and scale the density [case (A) and (B)], to scale the drag coefficient and model the density [case (C)], or to scale a linear combination of the drag coefficient and the thermospheric density [case (D)]. In this investigation, we scale the thermospheric density by the factor f s since CD is analytically computed using a GSI model (see Sect. 2). Nevertheless, we have to mention that any value for CD relies (to different extents) on the aforementioned assumptions and considerations concerning the atmospheric drag computation made in Sect. 2 and therefore, also the estimated scale factor f s is influenced by them (see also Bruinsma et al. 2012). For a thorough interpretation of CD and/or f s , all assumptions have to be considered. Nevertheless, the physical CD coefficient described in Sect. 2 is the best (most precise) one can currently obtain and therefore, the estimated factors f s are assumed to primarily scale the thermospheric densities. The unknown parameters estimated during the POD of ANDE-P are listed in Table 3 (in total 27 unknowns each 3.5 days). For each parameter, the selected temporal resolution and a priori sigma are specified. In the present study, the scale factor f s has been estimated with a 6-h resolution. Table 3 Estimated parameters, their temporal resolutions and a priori standard deviations (STD)

The satellite orbital period around the Earth is about 90 min, corresponding approximately to 1/4 of the validity interval of a scale factor value. Since the geographical distribution of SLR tracking stations is not uniform, the estimated scale factors are assumed valid also over areas where no observations are registered. In particular, the scale factors are related to densities corresponding to the thermospheric volume swept by the satellite during about four orbits completed in the 6-h validity interval. This means that the scale factors f s can be only seen as a global mean value of four satellite revolutions which are observed from a rather sparse ground network. The impact of this fact is discussed in the framework of Fig. 20.

4 Analysis of the results The results concerning ANDE-P POD are reported in terms of root-mean-square (RMS) values of SLR observation residuals shown in Fig. 11 for each thermospheric model tested in this investigation. Notice that the wind model HWM14 is applied to all solutions and therefore, will not be mentioned any more. The maximum SLR RMS reaches 6 cm for two arcs which correspond to the first half of the GPS weeks 1549 and 1550, while for all other arcs the RMS is lower than 4 cm. If the RMS values within the arcs are compared, no systematics can be identified. This means that the estimated scale factor values compensate most of the density differences of the thermospheric models. This is confirmed also by the binned SLR observation residuals, which are reported for all the models with the mean and standard deviation values shown in Fig. 12. RMS differences within the arcs are caused by the POD in the iterative least-squares adjustment process since the outlier rejection relies, to a certain extent, on the a priori thermospheric model used. The NRLMSISE00 model provides the smallest average values of the standard deviation and absolute mean of the SLR observation residuals among the four thermospheric models tested. The estimated 6-h scale factors for the different models are shown in Fig. 13. The values are always in the positive domain (since no negative densities are possible) and range from about 0.1 to 1.8. Another effect shown in Fig. 13

Estimated parameters

Temporal resolution

A priori STD

Keplerian elements

One set per arc (initial epoch)

0

Solar radiation pressure scale coefficient Cr

One per arc

0.1

Albedo scale coefficient

One per arc

0.1

CPR empirical coefficients

One set per arc (sin/cos; along-/cross-track) 0.1 (m/day2 )

Scale coefficients

Four per day

(for thermospheric density)

0.1

(6-h resolution)

The arc length selected to process ANDE-P SLR data is 3.5 days

123

F. Panzetta et al.

Fig. 11 RMS of ANDE-P SLR observation residuals using different thermospheric models

Fig. 12 Histograms of the SLR observation residuals for four thermospheric models. The circles represent the SLR observation residuals binned at 1 cm for the models NRLMSISE00 (blue), CIRA86 (green), DTM2013 (orange) and JB2008 (magenta). The corresponding solid curves are the splines which interpolate the circles. Mean and standard deviations of the residuals are reported for each model

Fig. 13 Estimated scale factors of the models NRLMSISE00 (blue), CIRA86 (green), DTM2013 (orange) and JB2008 (magenta)

is sequences of scale factors equal to 1.0. These scale factors correspond to intervals without any SLR observation of ANDE-P. Therefore, at these time intervals, no scaling of the a priori thermospheric density is done. It also has to be mentioned here that the estimated scale factors are affected by the subjective outlier detection (analyst influences outlier detection to a certain extent), the arc length, and the selected temporal resolution. In future investigations, the temporal

123

resolution might be varied in order to obtain more stable results (e.g. half-daily, daily or arc-wise resolution). One of the main results of this investigation, namely the mean, the RMS and the standard deviation of the estimated time series of scale factors f s , is summarized in Table 4. The values are computed using only estimated scale factors. The mean values of the scale factors f s of all four models are below 1. Therefore, it can be concluded that the applied thermospheric models overestimate the integral density to some extent (Doornbos 2011). The mean values of the scale factors f s related to the models NRLMSISE00 and CIRA86, namely around 0.65, with a comparable scattering are the lowest and are very similar to each other since both models belong to the same class of thermospheric models called MSIS (Doornbos 2011). The integrated density (total density for all k constituents) of JB2008 is scaled the least, since its mean scale factor (0.89) is the closest to 1. Therefore, the JB2008 model seems to be just slightly overestimating the thermospheric density (Solomon et al. 2015). The standard deviations of the estimated scale factors are similar for all the models, being around 0.25. According to Solomon et al. (2015), the thermospheric density changes by about − 5% per decade due to climate change. At this point, the estimated scale factors shown in Fig. 13 are applied to the corresponding modelled densities illustrated in Fig. 9. A linear interpolation between the 6-h scale factors is adopted when we apply them to the thermospheric densities. The statistics of these scaled densities is reported in Table 5. The mean values of the scaled densities are much closer to each other with respect to the mean values of the modelled densities, with JB2008 having both the lowest mean and RMS of the scaled densities. The scaled thermospheric density distributions given in Fig. 14 show smaller differences between the models than the original ones given in Fig. 9. A zoom over 1 day of the modelled densities, the estimated scale factors and the corresponding scaled densities are shown in Fig. 15. As anticipated by the statistics of the scale factors listed in Table 4, JB2008 is the least scaled model. Indeed, the values of the estimated scale factors strongly depend on the particular satellite orbit (especially altitude) and on the geomagnetic and solar activity at the time of the processed ANDE-P SLR data. As shown in Fig. 16 and reported by Marcos (2003) and Doornbos (2011), ANDE-P data in August, September and October 2009 correspond to a period of low solar and geomagnetic activity. Nevertheless, the obtained scaled integrated thermospheric densities might be used to calibrate empirical density models derived from satellite accelerometry. Together with the thermospheric scale factors, also other parameters are estimated during ANDE-P POD. The frequency of correlations above 0.5 between the estimated parameters is provided in Fig. 17 for the model JB2008-

Towards thermospheric density estimation from SLR observations of LEO satellites: a case study… Table 4 Mean, RMS and standard deviations (STD) of the 6-h estimated scale factors f s for each thermospheric model (the wind model HWM14 is always included)

Table 5 Mean, RMS and standard deviations (STD) of the scaled thermospheric densities for all the models

Mean ( f s )

RMS ( f s )

STD ( f s )

NRLMSISE00

0.651

0.697

0.250

CIRA86

0.646

0.697

0.263

DTM2013

0.790

0.825

0.236

JB2008

0.888

0.928

0.270

Thermospheric model

Mean (kg/m3 )

RMS (kg/m3 )

STD (kg/m3 )

NRLMSISE00

2.44 × 10−12

2.80 × 10−12

1.36 × 10−12

CIRA86

2.19 × 10−12

2.45 × 10−12

1.10 × 10−12

× 10−12

2.22

× 10−12

0.92 × 10−12

2.22

× 10−12

0.99 × 10−12

Thermospheric model

DTM2013

2.02

JB2008

1.99 × 10−12

Fig. 14 Thermospheric total densities scaled using the 6-h estimated scale factors of the models NRLMSISE00 (blue), CIRA86 (green), DTM2013 (orange) and JB2008 (magenta)

Fig. 17 Frequency of correlations above 0.5 for pairs of estimated parameters during the processed period, while using the model JB2008HWM14

Fig. 15 Zoom of the modelled thermospheric densities (top), the thermospheric scale factors (middle) and the corresponding scaled densities (bottom) over 1 day for the models NRLMSISE00 (blue), CIRA86 (green), DTM2013 (orange) and JB2008 (magenta)

Fig. 16 Solar radio flux index F10.7 and daily geomagnetic index A p covering the period 2000–2010 (Doornbos 2011)

HWM14 and is very similar for the other models NRLMSISE00, CIRA86, DTM2013. The highest frequency that can be observed for each pair of parameters is 14, because 14 arcs of 3.5 days are processed in total. The estimated parameters are indicated with numbers from 1 to 27 according to the following sequence: – parameters in the range 1–6 are the orbital elements (semi-major axis, eccentricity, inclination, argument of perigee, longitude of ascending node and mean anomaly); – parameter 7 is the solar radiation pressure scale coefficient; – parameter 8 is the Earth albedo scale coefficient; – parameters in the range 9–23 are the thermospheric scale factors;

123

F. Panzetta et al.

– parameters in the range 24–27 are the empirical acceleration coefficients (CPRs). It can be observed that the estimated Earth albedo scale coefficient does not have correlations above 0.5 with the estimated thermospheric scale factors. This result is expected because the albedo radiation pressure forcing is directed radially away from the Earth, so it is approximately perpendicular to the satellite velocity and to the atmospheric drag forcing for the near-circular ANDE-P orbit. On the other hand, the estimated solar radiation pressure scale coefficient shows sparse correlations above 0.5 with the estimated thermospheric scale factors. Considering that these scale factors are 14 per arc and that we processed 14 arcs, there is a maximum of 30 scale factors (over a total of 196) which are correlated to the solar radiation pressure scale coefficient when using the model JB2008-HWM14. There is a higher number of correlations above 0.5 between the orbital elements and the thermospheric scale factors, especially at the beginning of an orbital arc, and between the thermospheric scale factors themselves, especially between factors close in time, i.e. close in sequence. Finally, there are no correlations above 0.5 of any parameter with the CPR empirical acceleration coefficients. Other parameters are also investigated, besides the thermospheric scale factors estimated at each arc using four thermospheric density models. In Figs. 18 and 19, the estimated solar radiation pressure and the Earth albedo scale coefficients are shown for the different thermospheric models, respectively. It can be observed that there are no relevant discrepancies among the albedo scale coefficients estimated using different models for the same arc. In case of the solar radiation pressure scale coefficients Cr , however, significant differences can be detected. This can be explained by correlations of the Keplerian elements (initial state vectors) with the Cr coefficient (see Fig. 17). Besides these correlations, the Cr coefficient is also correlated with the thermospheric scale coefficients. Their reliability and stability are discussed in more detail in Sect. 5.

Fig. 18 Solar radiation pressure scale coefficient estimated once per arc using different thermospheric models

123

Fig. 19 Earth albedo scale coefficient estimated once per arc using different thermospheric models

Fig. 20 Differences between the semi-major axes estimated at the initial epoch and the minimum among them within each arc, using different thermospheric models

In Fig. 20, the variations of the estimated semi-major axis are shown for each thermospheric model. For each arc, the minimum semi-major axis is selected and subtracted from all other semi-major axes of that arc. That means in each arc the bar corresponding to the thermospheric model having the minimum semi-major axis is zero. In general, the estimated semi-major axes of CIRA86 and NRLMSISE00 often show significant discrepancies, even of several tens of metres, w.r.t. those of DTM2013 and JB2008. The differences can be explained by the fact that the thermospheric densities are only partly scaled at the different arcs processed in this experiment (see Fig. 13). Since different thermospheric densities are used a priori, those 6-h intervals where ρM is not scaled at all √ ( f s = 1.0) result in a wrong mean motion of the satellite ( (G M/a 3 )). If the satellite is not observed in the perigee or apogee, this wrong mean motion (caused by a wrong ρM ) is compensated by an offset in the semi-major axis a and eccentricity e. During these time periods, the models CIRA86 and NRLMSISE00 are overestimating the thermospheric density compared to JB2008 and DTM2013 (see Fig. 14). Due to this density offset, the mean velocity of ANDE-P during an arc is reduced. Since SLR observations to ANDE-P are rather sparse, the change in the mean velocity can be compensated

Towards thermospheric density estimation from SLR observations of LEO satellites: a case study…

Fig. 21 Transverse (along-track) cosine and sine terms (CPR-T) of the empirical acceleration estimated once per revolution using different thermospheric models

Fig. 22 Normal (cross-track) cosine and sine terms (CPR-N) of the empirical acceleration estimated once per revolution using different thermospheric models

by a change in the semi-major axis a and the eccentricity e of the ANDE-P orbits. The other estimated orbital elements, namely inclination, argument of perigee, longitude of the ascending node and mean anomaly, are not shown here explicitly since no relevant or systematic differences caused by the thermospheric models can be seen. The cosine and sine terms of the empirical acceleration in the transverse (along-track, CPR-T) and normal (crosstrack, CPR-N) directions are plotted in Fig. 21 and in Fig. 22, respectively. These values remain small for the most of the arcs. Moreover, they do not show a notable trend caused by the thermospheric models.

the first one, we fix the solar radiation pressure scale coefficient Cr to 1.0 for all processed orbital arcs with the four thermospheric models. This causes rather small differences (0.001–0.012, i.e. 0.1–1.8%) of the mean values of the estimated thermospheric scale factors (see Table 6), as compared to the nominal case described in Sect. 4. In Table 6, we provide just the results for the NRLMSISE00 and DTM2013 models, which experience, respectively, the major and minor impact on the f s among the four models used. Fixing Cr increases the standard deviation of the estimated thermospheric scale factors for three models except for NRLMSISE00. Therefore, not estimating the solar radiation pressure scale coefficient has generally a minor impact on the mean, RMS, and standard deviation of the estimated scale factors and does not change the conclusions of our investigation. On the other hand, fixing Cr increases the RMS of SLR observation residuals, of cosine and sine terms of the normal and transverse empirical accelerations, the estimated value of the albedo scale coefficient and its standard deviation, but reduces the mean value of SLR observation residuals, except those of the NRLMSISE00 model. In the second experiment, we do not estimate the CPR empirical accelerations and fix them to zero. This has an even smaller effect on the mean, RMS and standard deviation values of the estimated scale factors f s (see Table 6). As an example, Fig. 23 provides time series of the estimated scale factor f s for the NRLMSISE00 model for the three test cases. In view of these results, we conclude that our solution in terms of thermospheric scale factors is not systematically or significantly affected by the estimation of the solar radiation pressure scale coefficient and CPR empirical acceleration parameters, leading to reliable and realistic estimated thermospheric scale factors. In the third experiment, we fixed the scale factors for each arc to the computed mean value documented in Table 4. This resulted in a large increase in the empirical CPR accelerations for the orbits where the 6-hourly estimated scale factors shown in Fig. 13 differ significantly from the mean value. This indicates that the mean values of the estimated scale factors should not be used for POD of ANDE-P and should not be over-interpreted. One should estimate a time series of scale factors during a POD to get a reliable orbit and a realistic estimation of the thermospheric density.

6 Conclusions and outlook 5 Reliability of the estimated thermospheric density scale factors In order to investigate the stability and reliability of the estimated scale factors f s which might be affected by correlations with other estimated parameters, we perform three experiments for each considered thermospheric model. In

In this investigation, we tested the sensitivity of SLR observations to different thermospheric models in the POD of the spherical ANDE-P satellite. Since the thermospheric densities computed using these models show different values, we refined the atmospheric drag perturbation modelling and estimated scale factors in order to compensate possi-

123

F. Panzetta et al. Table 6 Solution characteristics for two selected models and three test cases: no solar radiation pressure scale coefficient Cr is estimated, no CPR empirical accelerations are estimated, and the nominal case, when all these parameters are estimated Thermospheric model

Case

f s Mean (–)

f s RMS (–)

f s STD (–)

RMS SLR Mean SLR Albedo residuals residuals coefficient (cm) (cm) (–)

NRLMSISE00

No Cr estim.

0.639

0.685

0.247

3.428

0.416

2.296 ± 3.154 1.4448

NRLMSISE00

No CPR estim. 0.649

0.695

0.249

1.833

0.052

1.007 ± 0.056 0.0000

NRLMSISE00

Nominal

0.651

0.697

0.250

1.832

0.051

1.007 ± 0.056 0.0313

DTM2013

No Cr estim.

0.791

0.828

0.243

2.230

0.210

1.171 ± 0.342 0.0400

DTM2013

No CPR estim. 0.790

0.825

0.236

2.126

0.237

0.957 ± 0.117 0.0000

DTM2013

Nominal

0.825

0.236

2.126

0.237

0.957 ± 0.116 0.0295

0.790

Fig. 23 Estimated scale factor of the NRLMSISE00 model for three cases: nominal, no CPR empirical accelerations are estimated and no solar radiation pressure scale coefficient is estimated. The blue curve is behind the green curve

ble integral density differences. ANDE-P was chosen for this experiment since its spherical shape makes this satellite motion easier to model than of satellites with complex shape and since it was flying at an altitude of about 350 km where the thermospheric drag is the largest non-gravitational perturbation acting on LEO satellites. Due to this high sensitivity, reliable scale factors for the thermospheric density could be estimated. The investigation was done for four different thermospheric models, namely JB2008, DTM2013, NRLMSISE00 and CIRA86, in order to be able to intercompare the obtained results. We derived time series of scale factors of the thermospheric density at the time interval studied (16 August 2009 until 3 October 2009) and obtained their following mean values: 0.65 ± 0.25 for NRLMSISE00, 0.65 ± 0.26 for CIRA86, 0.79 ± 0.24 for DTM2013 and 0.89 ± 0.27 for JB2008 model. The obtained results establish that JB2008 is the least scaled model, followed by DTM2013, while CIRA86 and NRLMSISE00 have signifi-

123

RMS of CPRT cosine term (m/day2 )

cant scale factors. This suggests that all the models, to a larger or lower extent, are in general overestimating the thermospheric densities along ANDE-P orbit during the analyzed period. We studied the correlations between the estimated parameters and showed that the estimated scale factors of the thermospheric density are not systematically or significantly affected by our setup of the estimated parameters, namely, by the estimation of the solar radiation pressure coefficient and once-per-revolution empirical accelerations. All thermospheric densities, when scaled using derived scale factors, show similar mean densities for four models (see Table 5). This is also a clear indication that the scaled densities are reliable. Besides the satellite altitude, the results depend on the geomagnetic and solar activity which is at a very low level during the processed time interval. The approach presented in this paper provides the possibility to calibrate empirical models like those obtained from satellite accelerometry measurements. Nevertheless, one has to keep in mind that the estimated densities are not purely scaled copies of the background models used since the model of the drag coefficient also holds small latitudinal variations due to its dependence on the atmosphere composition and temperature (see also Bruinsma et al. 2012). One should derive a time series of scale factors to obtain realistic scaled thermospheric density. In addition, one always has to keep the key parameters and assumptions in mind which are crucial in order to reliably interpret the obtained scaled densities. In addition to the empirical models investigated in the previous sections, thermospheric densities can also be obtained from physics-based models which have been developed to simulate the thermosphere behaviour, for instance, the Thermosphere Ionosphere Electrodynamic General Circulation Model (TIEGCM) (Roble et al. 1988; Richmond et al. 1992; Qian et al. 2014). Their performance against the empirical ones can be revealed in a similar way by estimating thermospheric scale factors in the SLR-based orbit determination. Beyond that, physics-based models can be extended by incorporation which paves the way of establishing a physics-based

Towards thermospheric density estimation from SLR observations of LEO satellites: a case study…

data assimilation approach. One should consider that the thermospheric densities can be improved by injecting ionospheric observations into the data assimilation system (see, e.g. Lee et al. 2012). The obtained results can be validated again using SLR observations. Validation and comparison of physics-based approaches with empirical ones including the present results will be considered in future work. As a future work, further tests can be performed using different temporal resolutions and different a priori values (not equal to 1) of the estimated scale factors (e.g. mean values from Table 4). Moreover, it would be interesting to process the complete time series of ANDE-P and ANDE-C together with the other two ANDE-RR satellites and, e.g. Spinsat. This investigation would provide the possibility to compare the estimated thermospheric density distributions for satellites flying at a similar altitude. Acknowledgements The work described in this paper was carried out within the project “Interactions of Low-orbiting Satellites with the Surrounding Ionosphere and Thermosphere (INSIGHT)” funded by the German Research Foundation (DFG) through the Special Priority Program 1788 “Dynamic Earth”. The authors also want to thank the editor-in-chief, Prof. Jürgen Kusche, and the three anonymous reviewers for their comments and suggestions which clearly improved the quality and readability of the paper.

References Akmaev RA (2011) Whole atmosphere modeling: connecting terrestrial and space weather. Rev Geophys 49(4):4004. https://doi.org/10. 1029/2011RG000364 Andrés JI, Noomen R, Bianco G, Currie DG, Otsubo T (2004) Spin axis behavior of the LAGEOS satellites. J Geophys Res 109:B06403. https://doi.org/10.1029/2003JB002692 Bloßfeld M (2015) The key role of Satellite Laser Ranging towards the integrated estimation of geometry, rotation and gravitational field of the Earth. PhD Dissertation, Technische Universität München (TUM), Munich Bowman BR (2004) The semiannual thermospheric density variation from 1970 to 2002 between 200–1000 km. Adv Astronaut Sci 119:04–174 Bowman BR, Moe K (2005) Drag coefficient variability at 175-500 km from the orbit decay analyses of spheres. In: AAS/AIAA astrodynamics specialist conference, Lake Tahoe, CA, AAS 05-257 Bowman BR, Tobiska WK, Marcos FA, Huang CY, Lin CS, Burke WJ (2008) A new empirical thermospheric density model JB2008 using new solar and geomagnetic indices. AIAA/AAS astrodynamics specialist conference and exhibit, AIAA 2008-6438 Bruinsma SL, Biancale R (2003) Total densities derived from accelerometer data. J Spacecr Rockets 40(2):230–236. https://doi. org/10.2514/2.3937 Bruinsma SL, Sánchez-Ortiz N, Olmedo E, Guijarro N (2012) Evaluation of the DTM-2009 thermosphere model for benchmarking purposes. J Space Weather Space Clim 2:A04. https://doi.org/10. 1051/swsc/2012005 Bruinsma SL (2015) The DTM-2013 thermosphere model. J Space Weather Space Climate 5:A1. https://doi.org/10.1051/swsc/ 2015001 Cook GE (1965) Satellite drag coefficients. Planet Space Sci 13:929– 946. https://doi.org/10.1016/0032-0633(65)90150-9

Doornbos E, Klinkrad H, Scharroo R, Visser PNAM (2007) Thermosphere density calibration in the orbit determination of ERS-2 and Envisat. In: Proceedings of envisat symposium, Montreux, ESA SP-636 Doornbos E, Van Den Ijssel J, Lühr H, Förster M, Koppenwallner G (2010) Neutral density and crosswind determination from arbitrarily oriented multiaxis accelerometers on satellites. J Spacecr Rockets 47(4):580–589. https://doi.org/10.2514/1.48114 Doornbos E (2011) Thermospheric density and wind determination from satellite dynamics. PhD Dissertation, Technische Universiteit Delft, The Netherlands, https://doi.org/10.1007/978-3-64225129-0 Doornbos E, Bruinsma S, Fritsche B, Visser PNAM, Van Den IJssel J, Teixeira da Encarnação J, Kern M (2013) Air density and wind retrieval using GOCE data. In: Proceedings of the ESA living planet symposium, Edinburgh (United Kingdom), ESA SP-722 Drob DP, Emmert JT, Meriwether JW, Makela JJ, Doornbos E, Conde M, Hernandez G, Noto J, Zawdie KA, McDonald SE, Huba JD, Klenzing JH (2015) An update to the Horizontal Wind Model (HWM): the quiet time thermosphere. Earth and Space Science 2:301–319. https://doi.org/10.1002/2014EA000089 Emmert JT (2015) Thermospheric mass density: a review. Adv Space Res 56(5):773–824 Flanagan HP (2015) Improved atmospheric density estimation for ANDE-2 satellites using drag coefficients obtained from gassurface interaction equations. Master of Science Thesis, Aerospace Engineering, University of Kansas, USA, https://kuscholarworks. ku.edu/bitstream/handle/1808/19392/Flanagan_ku_0099M_ 14309_DATA_1.pdf Folkner WF, Williams JG, Boggs DH (2008) Planetary ephemeris DE421 for Phoenix navigation, IOM 343R-08-002 Förste C, Bruinsma S, Shako R, Marty JC, Flechtner F, Abrikosov O, Dahle C, Lemoine JM, Neumayer KH, Biancale R, Barthelmes F, König R, Balmino G (2011) EIGEN-6—a new combined global gravity field model including GOCE data from the collaboration of GFZ-Potsdam and GRGS-Toulouse. Geophysical Research Abstracts, 13, EGU2011-3242-2, EGU General Assembly, 2011 Fritsche B, Ivanov M, Kashkovsky A, Koppenwallner G, Kudryavtsev A, Voskoboinikov U, Zhukova G (1998) Radiation pressure forces on complex spacecraft, final report, ESOC contract 11908/96/D/IM, HTG, Germany and ITAM, Russia Gaposchkin EM, Coster AJ (1988) Analysis of satellite drag. Linc Lab J 1(2), 203–224 1(2):203–224 Gerstl M (1997) Parameterschätzung in DOGS-OC, In: DGFI Interner Bericht, MG/01/1996/DGFI, 2nd edn (in German) Hedin AE, Spencer NW, Killeen TL (1988) Empirical global model of upper thermosphere winds based on Atmosphere and Dynamics Explorer satellite data. J Geophys Res 93:9959–9978. https://doi. org/10.1029/JA093iA09p09959 Heelis RA (2011) Aspects of Coupling Processes in the Ionosphere and Thermosphere. In: Huba J, Schunk R, Khazanov G (eds) Modeling the Ionosphere-Thermosphere System. Wiley, Chichester Jackson J (1924) Note on the numerical integration of d2 x/dt2 =f(x, t). Mon Not R Astron Soc 84:602–606 Jeon HS, Cho S, Kwak YS, Chung JK, Park JU, Lee DK, KuzmiczCieslak M (2011) Mass density of the upper atmosphere derived from Starlette’s Precise Orbit Determination with Satellite Laser Ranging. Astrophys Space Sci 332(2):341–351. https://doi.org/10. 1007/s10509-010-0528-2 King-Hele D (1987) Satellite orbits in an atmosphere: theory and applications. Blackie and Son Ltd., Glasgow Knocke PC, Ries JC, Tapley BD (1988) Earth radiation pressure effects on satellites. In: AIAA/AAS astrodynamics conference, AIAA 88– 4292. Minneapolis, MN. https://doi.org/10.2514/6.1988-4292

123

F. Panzetta et al. Konopliv AS, Asmar SW, Carranza E, Sjogren WL, Yuan D-N (2001) Recent gravity models as a result of the Lunar Prospector Mission. Icarus 150(1):1–18 Lechtenberg TF, McLaughlin CA, Locke T, Krishna DM (2013) Thermospheric density variations: observability using precision satellite orbits and effects on orbit propagation. Space Weather 11:34–45. https://doi.org/10.1029/2012SW000848 Lechtenberg TF (2015) Density model corrections derived from orbit data to characterize upper atmospheric density variations. PhD Dissertation, Aerospace Engineering, University of Kansas, USA, https://kuscholarworks.ku.edu/bitstream/handle/ 1808/18430/Lechtenberg_ku_0099D_13830_DATA_1.pdf Lee IT, Matsuo T, Richmond AD, Liu JY, Wang W, Lin CH, Anderson JL, Chen MQ (2012) Assimilation of FORMOSAT-3/COSMIC electron density profiles into a coupled thermosphere/ionosphere model using ensemble Kalman filtering. J Geophys Res 117(A10):A10318. https://doi.org/10.1029/2012JA017700 Liu H, Lühr H, Henize V, Köhler W (2005) Global distribution of the thermospheric total mass density derived from CHAMP. J Geophys Res 110:A04301. https://doi.org/10.1029/2004JA010741 Lühr H, Rentz S, Ritter P, Liu H, Häusler K (2007) Average thermospheric wind patterns over the polar regions, as observed by CHAMP. Ann Geophys 25:1093–1101. https://doi.org/10.5194/ angeo-25-1093-2007 Marcos FA (2003) New satellite drag modeling capabilities. In: 44th AIAA aerospace sciences meeting and exhibit. https://doi.org/10. 2514/6.2006-470 Marcos FA, Lai ST, Huang CY, Lin CS, Retterer J, Delay S, Sutton E (2010) Towards next level satellite drag modeling. In: AIAA atmospheric and space environments conference, Toronto, Ontario, Canada. https://doi.org/10.2514/6.2010-7840 Mehta PM, Walker A, McLaughlin CA, Koller J (2014) Comparing physical drag coefficients computed using different gas-surface interaction models. J Spacecr Rockets 51(3):873–883. https://doi. org/10.2514/1.A32566 McLaughlin CA, Lechtenberg T, Fattig E, Krishna DM (2012) Estimating density using precision satellite orbits from multiple satellites. J Astronaut Sci 59(1):84–100. https://doi.org/10.1007/s40295-0130007-4 Moe K, Moe MM (2005) Gas-surface interactions and satellite drag coefficients. Planet Space Sci 53:793–801. https://doi.org/10. 1016/j.pss.2005.03.005 Montenbruck O, Gill E, Montenbruck O, Gill E (2000) Satellite Orbits. Springer, Berlin. https://doi.org/10.1007/978-3-642-58351-3 Nicholas AC, Picone JM, Emmert J, DeYoung J, Healy L, Wasiczko L, Davis MA, Cox C (2007) Preliminary results from the atmospheric neutral density experiment risk reduction mission. In: Proceedings of the AAS/AIAA astrodynamics specialist conference, AAS 07265. http://ilrs.gsfc.nasa.gov/docs/ANDE_AAS_07-265.pdf

123

Nicholas AC, Finne T, Davis MA, Kessel R (2009) Atmospheric neutral density experiment (ANDE-2) flight hardware details, May 26, 2009. http://ilrs.gsfc.nasa.gov/docs/andehw.pdf Pearlman MR, Degnan JJ, Bosworth JM (2002) The international laser ranging service. Adv Space Res 30(2):135–143. https://doi.org/ 10.1016/S0273-1177(02)00277-6 Petit G, Luzum B (2010) IERS conventions (2010), technical note 36. Verlag des Bundesamts für Kartographie und Geodäsie, Frankfurt, Germany. ISBN 3-89888-884-3 Picone JM, Hedin AE, Drob DP, Aikin AC (2002) NRLMSISE-00 empirical model of the atmosphere: statistical comparisons and scientific issues. J Geophys Res 107(A12):1468. https://doi.org/ 10.1029/2002JA009430 Pilinski MD, Argrow BM, Palo SE (2010) Semi-empirical model for satellite energy-accommodation coefficients. J Spacecr Rockets 47(6):951–956. https://doi.org/10.2514/1.49330 Qian L, Burns AG, Emery BA, Foster B, Lu G, Maute A, Richmond AD, Roble RG, Solomon SC, Wang W (2014) The NCAR TIE-GCM: a community model of the coupled thermosphere/ionosphere system. In: Huba J, Schunk R, Khazanov G (eds) Modelling the ionosphere-thermosphere system. Wiley, Chichester Qian L, Solomon SE (2012) Thermospheric density: an overview of temporal and spatial variations. Space Sci Rev 168(1):147–173. https://doi.org/10.1007/s11214-011-9810-z Ray RD, Ponte RM (2003) Barometric tides from ECMWF operational analyses. Ann Geophys 21:1897–1910 Richmond AD, Ridley EC, Roble RG (1992) A thermosphere/ionosphere general circulation model with coupled electrodynamics. Geophys Res Lett 19(6):601604. https://doi. org/10.1029/92GL00401 Roble RG, Ridley EC, Richmond AD, Dickinson RE (1988) A coupled thermosphere/ionosphere general circulation model. Geophys Res Lett 15:1325–1328. https://doi.org/10.1029/GL015i012p01325 Savcenko R, Bosch W (2012) EOT11A—Empirical Ocean Tide model from multi-mission satellite altimetry, Report No. 89, Deutsches Geodätisches Forschungsinstitut (DGFI), München, Germany Sentman LH (1961) Free molecule flow theory and its application to the determination of aerodynamic forces. Lockheed Missile and Space Co., LMSC-448514, AD265-409 Solomon SC, Qian L, Roble RG (2015) New 3-D simulations of climate change in the thermosphere. J Geophys Res Space Phys 120(3):2183–2193. https://doi.org/10.1002/2014JA020886 Sutton EK, Nerem RS, Forbes JM (2007) Density and winds in the thermosphere deduced from accelerometer data. J Spacecr Rockets 44(6):1210–1219. https://doi.org/10.2514/1.28641 Vallado DA (2007) Fundamentals of Astrodynamics and Applications, Third edn. Space Technology Library Walker A, Mehta P, Koller J (2014) Drag coefficient model using the Cercignani–Lampis–Lord gas-surface interaction model. J Spacecr Rockets 51(5):1544–1563. https://doi.org/10.2514/1. A32677