Research Article
Semiempirical thermosphere model evaluation at low altitude with GOCE densities
^{1}
OMPGET/CNES, Department of Space Geodesy, 18 Avenue E. Belin, 31401
Toulouse Cedex 4, France
^{2}
Astronomical Institute of the University of Bern, Hochschulstrasse 6, 3012
Bern, Switzerland
^{3}
Elecnor Deimos, Ronda de Poniente 19, 28760
Tres Cantos, Spain
^{*} Corresponding author: sean.bruinsma@cnes.fr
Received:
27
July
2016
Accepted:
21
December
2016
Aims: The quality of the Committee on Space Research (COSPAR) International Reference Atmosphere models NRLMSISE00, JB2008, and DTM2013 in the 150–300 km altitude range has never been thoroughly evaluated due to a lack of good density data. This study aims at providing the model accuracies thanks to the recent highresolution highaccuracy Gravity field and steadystate Ocean Circulation Explorer (GOCE) density dataset. The evaluation was performed on yearly, monthly, and daily time scales, which are important for different applications such as mission design, mission operation, or reentry predictions.
Methods: The accuracy of the models was evaluated by comparing to the GOCE density observations of the Science Mission (1 November 2009–20 October 2013) and new density data at the lowest altitudes derived for the last weeks before the reentry (22 October–8 November 2013) according to a metric, which consists of computing mean, standard deviation and root mean square (RMS) of the observedtomodel ratios, and correlation. Mean statistics are then calculated over the three time scales.
Results: The range of model biases, standard deviations, and correlations becomes larger when the time interval decreases, and this study provides COSPAR International Reference Atmosphere (CIRA) model statistics in the altitude range of 275–170 km. DTM2013 is the least biased and most accurate model on all time scales, essentially thanks to the database, notably containing two years of GOCE densities, to which it was fitted. NRLMSISE00 performs worst, with considerable bias of about 20% in 2009 and 2013, and systematically higher standard deviations (lower correlations) than JB2008 and DTM2013. The performance of JB2008 is presently only slightly behind DTM2013, thanks to the new release 4_2g solar activity proxies. However, it still presents some weakness under the lowest solar activity conditions in 2009 and 2010. Comparison to Challenging MiniSatellite Payload (CHAMP) density data showed that the results based on GOCE densities, despite limited local solar time coverage of 6–8 am & pm, are representative of model performance.
Key words: Thermosphere model / CIRA models / GOCE density
© S. Bruinsma et al., Published by EDP Sciences 2017
This is an Open Access article distributed under the terms of the Creative Commons Attribution License (http://creativecommons.org/licenses/by/4.0), which permits unrestricted use, distribution, and reproduction in any medium, provided the original work is properly cited.
1. Introduction
Semiempirical thermosphere models are used in the computation of the atmospheric drag force in satellite orbit determination and prediction, as well as in atmospheric studies. They predict pointwise temperature and (partial) density as a function of location (altitude, latitude, longitude, local solar time), solar and geomagnetic activities, and season. The Committee on Space Research (COSPAR) selected three thermosphere models in 2012, and they are described in the COSPAR International Reference Atmosphere (CIRA2012) report. These models, NRLMSISE00 (Picone et al. 2002), JB2008 (Bowman et al. 2007), and DTM2009 (Bruinsma et al. 2012), are also the ISO (International Organization for Standardization) models for the neutral upper atmosphere. The Drag Temperature Model (DTM) was updated recently, and the current version is DTM2013 (Bruinsma 2015). It was developed in the framework of the Advanced Thermosphere Modelling and Orbit Prediction project (ATMOP; http://www.atmop.eu), which was a European Union 7th Framework project. DTM2013 and the preATMOP benchmark DTM2009 were evaluated by comparing to total density data in order to demonstrate and quantify the significant improvements made, especially in the 250–500 km altitude range (Bruinsma et al. 2012; Bruinsma 2015).
The satellite GOCE (Gravity field and steadystate Ocean Circulation Explorer; ESA 1999), was the first Earth Explorer mission of the European Space Agency (ESA). It was launched in March 2009 in a 96.5° inclination, quasi duskdawn orbit. The Science Mission phase, maintained at constant mean altitudes between 275 and 230 km thanks to drag compensation using ion propulsion, lasted from 1 November 2009 through 20 October 2013. GOCE neutral densities, with a resolution of 80 km along the orbit and a precision of a few percent (Bruinsma et al. 2014), are made available by ESA (https://earth.esa.int/web/guest/missions/goce/gocethermosphericdata). Whereas the Bruinsma et al. (2014) paper compared GOCE densities (data version 1_0) to models mainly for validation purposes of the density data up to May 2012, in this study the validated and complete GOCE density dataset (data version 1_4) is exploited to evaluate the models thoroughly.
The GOCE dataset is the only highresolution, highaccuracy dataset at low altitude. As such, it is the best dataset for thermosphere model evaluation at low altitude. The Xenon for the ion propulsion was exhausted on 21 October 2013, and GOCE reentered the atmosphere three weeks later on 11 November 2013. The Global Positioning System (GPS) receiver operated almost to the end, and the accelerometers measured the drag acceleration until they saturated on 8 November 2013. Density was derived from these measurements, which provides additional information for model testing at altitudes between 230 and 170 km. This density dataset was compared to the models separately from the ESA density dataset (described in Sect. 2.1.1), because it was computed with different software as part of this study (see Sect. 2.1.2), and because the GOCE satellite was no longer in Science mode.
The objective of this study is to quantify the performance of the CIRA models at low altitudes of 270 km and less, i.e. when reentry of objects in nearcircular orbits becomes imminent, using GOCE density data. Performance is evaluated on long time scales of the order of years, which is important in case of mission design and lifetime (i.e. required thrust and fuel), to typical operational time scales of a month, as well as time scales of the order of days, which is important in case of mission operations and reentry predictions. Model evaluations have been published in the past, but they may no longer be pertinent because the CIRA models have since been updated and because the quality of the density data is not comparable to GOCE (e.g. Marcos et al. 1983; Hedin 1987; Marcos 1990). Model evaluations are also given in the NRLMSISE00, DTM2009, and DTM2013 papers, but not in such detail at low altitude. Other studies focused on correcting model error (model calibration), which required quantifying the error, i.e. model evaluation (e.g. Doornbos et al. 2008).
The next section briefly reviews the models and the density data used in the study. The evaluation results are presented in Section 3. A summary of findings and the conclusions of the analyses are given in Section 4.
2. Description of models and data
2.1. Relevant information on NRLMSISE00, JB2008, and DTM2013
The thermosphere varies on global scales with annual, semiannual, solar rotation (about 27 days), solar active region (a few months), and solar cycle (11 years on average) periodicities. These variations are represented in current models of the thermosphere. The CIRA models were constructed by fitting to their respective underlying density databases as well as possible in the leastsquares sense. They are climatology (or “specification”) models of the thermosphere with a low spatial and temporal resolution of the order of thousands of kilometers and hours, respectively. Density variations with smaller spatial scales or shorter time scales are sources of geophysical noise (e.g. gravity waves, upward propagating tides; Bruinsma & Forbes 2008), which, in the absence of all other sources of error, represent the ultimately achievable model accuracy. The solar and geomagnetic indices limit the temporal resolution of these models to 24 and 1–3 h, respectively. The main data sources and solar and geomagnetic indices of each model are listed in Table 1; more detailed information for each model is given in Picone et al. (2002), Bowman et al. (2007), and Bruinsma (2015). The dates and cadences for some of the NRLMSISE00 datasets could not be retrieved, but the purpose of this table is to show the differences in the databases. DTM2013 is the only model that has partly assimilated the GOCE density dataset, and consequently it is expected to compare best. NRLMSISE00 has not assimilated any recent total density data. Figure 1 shows predictions of the CIRA models at 250 km altitude for low (mean F10.7 = 75 sfu; 14 December 2009) and medium (mean F10.7 = 144 sfu; 14 December 2011) solar activity in the Science Mission phase of GOCE. The geomagnetic activity kp was less than 2 on both dates. It is a qualitative example of model predictions from which typical differences can be gleaned that are due to the underlying database of each model. DTM2013 and JB2008 predictions are rather similar in shape, with differences in amplitude of about 10–15%. The amplitudes predicted with NRLMSISE00 lead to about the same differences, but its structure is rather different. This is due to the very large amount of recent satellitedrag inferred densities assimilated in DTM2013 and JB2008 as opposed to incoherent scatter radar data (i.e. a temperature measurement) combined with a small amount of satellitedrag inferred densities assimilated in NRLMSISE00, and a more (NRLMSISE00) or less (JB2008) sophisticated modeling of the diurnal variation, e.g. with or without cross coupling with solar and geomagnetic activities.
Fig. 1.
Model predictions, in g/cm^{3}, at 250 km altitude close to winter solstice for low (left frames) and medium solar activity. 
The main data used in the construction of the CIRA models.
DTM2013 and NRLMSISE00 can be used for the entire space age (ap, F10.7 and F30 measurements start before 1957), and they predict temperature, density, and composition. JB2008 is constructed using a combination of solar and geomagnetic proxies and indices that have been available since 1998 (i.e. it cannot be used before that year), and it predicts temperature and density. An additional weakness of JB2008 is that atmospheric composition, which is required together with temperature to calculate the aerodynamic coefficient of a satellite, is not standard output. Figure 2 shows the 81day mean F10.7 (NRLMSISE00 and JB2008), F30 (DTM2013), and S81c (JB2008; integrated 26–34 nm) proxies.
Fig. 2.
81day mean of the solar activity proxies F10.7, F30, and S81c (left axis), and planetary geomagnetic index ap (right axis). 
The way the solar and geomagnetic indices should be used (NRLMSISE00 and DTM2013) or are used (JB2008) in the models is not entirely clear. Specifically, the models can be run with indices that are constant or interpolated over the day in case of solar proxies, or constant or interpolated over the 3hour intervals in case of planetary geomagnetic indices ap. There is no consensus or standard method, and this may lead to differences in model predictions of a few percent. In this study, mean (daily) solar proxies were interpolated linearly over the current day (the current day delayed by 24 h) for DTM2013 and NRLMSISE00. For DTM2013, the ap, delayed by 3 h and the mean of the last 24 h, were interpolated linearly in the 3hour intervals and then converted to kp. NRLMSISE00 was run in precise mode, with seven (interpolated) values of ap as specified in the subroutine header. JB2008 code was run with the proxy files (SOLFSMY and DTCFILE) made available by Space Environment Technologies. The user only has to specify the date, and the JB2008 subroutine then returns a density.
2.2. The total density data
The neutral density data that were used in this study are inferred from GOCE data, with additional highresolution density data from Challenging MiniSatellite Payload (CHAMP) for specifically testing the performance as a function of local solar time.
2.2.1. Total density dataset of GOCE in Science Mode
Densities were inferred from thruster data for the Science Mission (Doornbos et al. 2014), which lasted from 1 November 2009 through 20 October 2013. GOCE neutral densities have a 0.1 Hz sampling, which leads to approximately 80 km resolution along track, with a precision of a few percent at worst and better than 1% after 2010 (Bruinsma et al. 2014). Due to the quasi duskdawn orbit, the local solar time coverage of GOCE densities is limited to 6–8 am & pm (the orbit precessed from 6 to 8 at the end of the mission). The density data is freely available on the dedicated ESA GOCE website (https://earth.esa.int/web/guest/missions/goce/gocethermosphericdata). In this study, the GOCE data were scaled by a factor of 1.25, as was done in the ESA Final Report (Doornbos et al. 2014), effectively assuming that the High Accuracy Satellite Drag Model (HASDM) of the US Air Force Space Command (Storz et al. 2005) is unbiased. This is not necessarily true because the scale of a thermosphere model is closely related to the aerodynamic coefficients that were used in the derivation of atmospheric density. As a consequence, DTM2013 and JB2008 (which can be considered the climatology version of the HASDM model) are expected to be unbiased with respect to the scaled GOCE densities. Applying a different scale factor than 1.25 will lead to different results of the evaluation, because values of observed and modeled density are analyzed in the form of mean of density ratios (i.e. “model bias”), and not only as standard deviation and correlation. The dailymean densities and the spherical (above 6378 km) altitude are shown in Figure 3 (the slight scatter in the altitude is due to averaging over incomplete days). The increasing density in 2012–2013 is due to the staggered lowering of the orbit altitude as well as solar activity.
Fig. 3.
Dailymean GOCE densities and mean spherical altitude above 6378 km (red; right axis). 
2.2.2. Total density dataset of GOCE in the reentry phase
The last three weeks of the GOCE mission are not part of the Science Mission, and the GOCE accelerometer measurements in a lower resolution mode were provided by ESA as a special dataset. Densities were derived from the special dataset using the methodology described in Bruinsma et al. (2004) from 22 October to 8 November 2013, and they were also scaled to HASDM for consistency. Their relative precision was estimated at a few percent, i.e. comparable to the thrusterinferred densities. This high precision was possible thanks to the everincreasing drag acceleration experienced as GOCE decayed, which offset the lower accelerometer resolution in the reentry phase. Precise kinematic orbit positions were used in a leastsquares dynamic orbit adjustment procedure in which the accelerometer bias was estimated every 12 h. The scale factors have been fixed to the values provided in Visser & van den IJssel (2015). Actual scale factors in the reentry phase are likely different by a few percent, but scaling the densities to HASDM minimizes the already small effect of the error. The satellite macro model for GOCE was simply a rectangular box with a frontal area of 0.7 m^{2} and sides of 10.77 and 5.90 m^{2} in the cross track and radial directions, respectively, and a mass of 1000 kg. The drag coefficient was computed using the Sentman model (Sentman 1961) for flat plates. The macro model was correctly oriented in inertial space thanks to the level1b attitude quaternions (GO_CONS_EGG_IAQ_2C; Gruber et al. 2012).
The precise kinematic GOCE orbits were obtained by processing the dualfrequency GPS observations collected by the onboard 12channel Lagrange GPS receiver (Intelisano et al. 2008). The undifferenced carrier phase data were used in the ionospherefree linear combination for a leastsquares adjustment of the epochwise kinematic positions and receiver clock corrections, as well as the carrier phase ambiguities (Svehla & Rothacher 2005). The pseudorange observations were used only for the initial receiver clock synchronization and to obtain a very first a priori orbit from a kinematic code point positioning. Using the latest version of the Bernese GNSS Software v5.3 (Dach et al. 2015), the GPS data were processed at the full sampling of 1 Hz, together with the attitude data delivered by the onboard star cameras according to Bock et al. (2014). The GPS orbits, the 5 s GPS clock corrections, and the Earth rotation parameters of the Center for Orbit Determination in Europe (CODE; Bock et al. 2009; Dach et al. 2009) were used, to guarantee the highest consistency between the GOCE orbit determination and the GPS orbit determination, which is also performed with the Bernese GNSS Software.
For the performance of the GOCE POD during the Science Mission phase, the reader is referred to Bock et al. (2014). In particular, the independent orbit validation by means of Satellite Laser Ranging (SLR) yielded a mean value of 1 mm and a root mean square (RMS) of 24.2 mm for the SLR residuals of the kinematic orbits over the entire Science Mission phase. During the last three weeks, GOCE was tracked by SLR on three days only, due to the challenge imposed by the low orbital altitude to provide reliable orbit predictions to the tracking network (Jäggi et al. 2011). Table 2 shows the daily statistics of the SLR validation for the kinematic GOCE orbits of these three days. It has to be noted, however, that the number of SLR observations is very limited and the statistics, therefore, not very robust. Figure 4 shows the daily RMS values of the ionospherefree carrier phase residuals for the last three weeks of GOCE, which is comparable to results of the Science Mission. The quality of the GOCE orbits in the last three weeks is only very slightly lower than during the Science Mission; it allows accurate accelerometer bias estimation to the end of the mission. As a result, the densities are only minimally affected by calibration error, and negligibly by orbit position error.
Fig. 4.
RMS of the GOCE GPS ionospherefree phase residuals. 
Number of SLR observations, dailymean values, standard deviations, and RMS of SLR residuals of the kinematic orbits of the last three weeks of GOCE. SLR stations: 7810 (Zimmerwald, Switzerland), 7110 (Monument Peak, CA, USA), and 7090 (Yarragadee, Australia).
2.2.3. Total density inferred from accelerometer measurements on CHAMP
Neutral densities were derived from accelerometer measurements on the Challenging MiniSatellite Payload (CHAMP; Reigber et al. 1996) in the altitude range 450–250 km using the methodology described in Bruinsma et al. (2004). The CHAMP dataset used in this study covers the period 20 May 2001 through 2 September 2010, i.e. up to 17 days before its reentry in the atmosphere. Due to its precessing orbit plane, 24hour local solar time sampling was achieved approximately every four months. The accelerometer provided highresolution measurements from which densities were inferred with about 80 km along track resolution and a precision of a few percent (Bruinsma et al. 2004) nearly from poletopole (87° inclination).
3. Results
The models are evaluated by computing the mean of the density ratios, defined as observedtocalculated ratios (“O/C”; unity for an unbiased model), the standard deviation and RMS of the density ratios minus one (simply “standard deviation of the density ratios” and “RMS of the density ratios” in the following). The numbers reflect the relative precision of the models; the absolute precision, i.e. in kg/m^{3}, is rather difficult to interpret due to its (orders of magnitude) difference as a function of altitude or solar cycle. Model bias, i.e. the mean of the density ratios minus one, is most damaging in orbit extrapolation because it causes position errors that increase in time. The standard deviation represents a combination of the ability of the model to reproduce the observed variations and the geophysical and instrumental noise in the observations, inclusive of bias in case of RMS. The correlation coefficients R are also computed. Contrary to RMS, they are insensitive to model bias and R ^{2} represents the fraction of observed variance captured by the model. In the next sections, the relative precisions of the models are evaluated on three typical time scales, ranging from years (i.e. solar cycle or part of), to months (i.e. solar rotation), to days.
3.1. Performance on long time scales
The model performance on long time scales is evaluated by computing the density ratios averaged over the entire GOCE Science Mission of nearly 4 years, as well as per calendar year. The results are listed in Table 3. The bias of the models over the entire Science Mission is small, less than 5% for all three models. This remains true only for DTM2013 when the averaging is done over one year intervals (its bias is smallest, 2% or less), whereas the NRLMSISE00 bias is 11% for 2013. The number of observations varies per year due to data gaps, which is why the results for the entire Science Mission are not the average of the annual results. The standard deviations of the density ratios are smallest for DTM2013, notably in 2010 during the low solar activity conditions and the onset of solar cycle 24; it is less than 10% in 2012 and 2013. JB2008 performs worst when solar activity was low, with largest bias and standard deviation in 2010, and smallest standard deviation in 2013. NRLMSISE00 on the other hand performs worst when solar activity was highest, with largest bias and standard deviation in 2013. The standard deviation of the density ratios is expected to decrease as solar activity increases (Bruinsma & Forbes 2008). The correlation coefficient of the models is 0.97 or larger for the entire Science Mission. DTM2013 has rather similar correlations per year, whereas NRLMSISE00 presents the largest variations in annual correlations. JB2008 has the smallest correlation in 2010, confirming its worst performance when solar activity was low. The three models have the highest correlation in 2011 for unexplained reasons.
Model performance on long time scales. Results for all GOCE Science Mission data (top) and per year (bottom: 2010–2013). NB: JB2008 proxies release 4_2g were used (June 2016).
3.2. Performance on monthly time scales
The monthly evaluation is done because it is an operational time scale (mission planning) as well as a physical one, close to the solar rotation period of approximately 27 days. The performance of the models is expected to decrease when averaging is done over shorter time spans. This is confirmed in Figure 5, which displays the monthly mean density ratios (Fig. 5a), standard deviations and correlations (bottom). The DTM2013 (JB2008; NRLMSISE00) bias shown in Figure 5a, which was 2% (6.3%; 11.5%) or less for the yearly averages, increases to 8.4% (12.3%; 20.1%) or less. It does not present a clear variation that hints at model error. The same holds for JB2008 except for 2010, which displays a semiannual variation in the density ratios. NRLMSISE00 displays an increasing bias in 2013, which is probably due to the significantly weaker F10.7 proxy for that year compared to F30 and S81c (Fig. 2). The arithmetic mean and standard deviation of the time series of the monthly values is given for each model, both of which are best for DTM2013.
Fig. 5.
Monthly mean density ratios (a), standard deviation (std) of the density ratios (b), and correlations (c). The last (solid) symbol is with the reentry dataset. JB2008 computed with proxies from June 2016. 
The monthly mean standard deviations of the density ratios (Fig. 5b) all present a negative slope in 2010–2011, after which they stabilize for DTM2013 and JB2008 but increase again for NRLMSISE00 in 2013. DTM2013 and JB2008 achieve very low standard deviations of less than 8% for several months. For the low solar activity conditions in the first few months of 2010, DTM2013 and NRLMSISE00 have much smaller standard deviations than JB2008, and this may be due to the different variation in activity between S81c, F10.7, and F30 that is visible in Figure 2.
The monthly correlation coefficients (Fig. 5c) are highest with DTM2013 for most of the months. NRLMSISE00 has the least high correlations, and an annual minimum around January is visible. JB2008 correlations become higher when solar activity increased and orbit altitude decreased. It has the smallest correlation of the three models in the Science Mission for the month January 2010, when it drops to just below 0.86. Periods of diminished correlation are visible for all three models. This will be shown in more detail in the next section.
The last month (#47), plotted in solid symbols, is actually the 3week reentry period dataset described in Section 2.1.2. The models achieve the highest correlation (0.98) in that period. JB2008 has the smallest bias and standard deviation with the reentry dataset, followed by DTM2013 and NRLMSISE00. The models’ performances compared with the Science Mission and reentry datasets are the same, i.e. they are the same in the altitude range 275–170 km; the best results are obtained with DTM2013 and JB2008.
3.3. Performance on daily time scales
The model performance on daily time scales is most important for (imminent) reentry predictions, as well as spacecraft conjunction predictions. The performance of the models decreases again compared to the monthly averages. This is shown in Figure 6, which displays the dailymean density ratios with DTM2013 (Fig. 6a), NRLMSISE00 and JB2008 (Figs. 6b and 6c). The reentry data is plotted in black for the three models. The mean and standard deviation of the time series of dailymean density ratios are also given for each model. The maximum model biases, and the standard deviations of the time series are much larger than the corresponding numbers of the monthly averages (Fig. 5a). This is mainly due to weaker model performance during periods of enhanced geomagnetic activity or storms, which typically have a duration of hours to a few days; these events are the main cause of the spikes in Figure 6. Inspection of time series of density ratios reveals that the largest biases (spikes) for DTM2013 (JB2008; NRLMSISE00) on daily time scales are 16–24% (20–28%; 30–38%), and the standard deviation of its time series increases by 59% (57%; 46%) compared to monthly averages. However, the standard deviations are still rather small for DTM2013 and JB2008 (5.9% and 6.4%, respectively), and larger for NRLMSISE00 (9.4%) due to the drift in the density ratios. The density ratios of the three models computed with the reentry dataset present large variations. It is largest with DTM2013, which has a clear step about in the middle of the period (31 October), and smallest with NRLMSISE00 (but biased by 15% on average). The performance of the three models is (visually as well as in terms of mean and standard deviation of the time series, which are given for each model in Fig. 6) in agreement with results for the entire Science Mission.
Fig. 6.
Dailymean density ratios with DTM2013 (a), NRLMSISE00 (b), and JB2008 (c) computed with proxies from June 2016. The reentry dataset is plotted in black for all three models. 
Due to the small number of geomagnetic storms over the GOCE mission, we cannot glean robust and accurate statistical results on model performance with respect to geomagnetic activity level. DTM2013 has a tendency to overestimate GOCE densities for enhanced geomagnetic and storm conditions by 10–20% (not shown), whereas NRLMSISE00 by about the same amount underestimates them; JB2008 does not appear to have a bias that is proportional to geomagnetic activity. The reduced model performance with geomagnetic activity can however be demonstrated by sorting the standard deviations of the density ratios and correlations over kp (the daily mean). This is displayed in Figure 7 (a zoom on days for which kp > 2), which shows that for enhanced and storm conditions, the standard deviations of the density ratios increase while correlations decrease. Averages for low activity and enhanced to storm activity were calculated and provided in the figure frames.
Fig. 7.
Dailymean standard deviations (std) of the density ratios (a) and correlations (b) sorted over kp (black; right axis). JB2008 computed with proxies from June 2016. 
The dailymean RMS of the density ratios of the three models and the 81day mean F30 flux are displayed in Figure 8. The time series are shown for NRLMSISE00 (Fig. 8a), JB2008 (Fig. 8b), and DTM2013 (Fig. 8c); the average daily RMS of the density ratios of DTM2013 (JB2008 and NRLMSISE00) is 0.10 (0.12 and 0.14). The relative variability of the thermosphere is highest under low solar activity conditions (Bruinsma & Forbes 2008), which is confirmed by the high RMS in 2009–2010 (for which model biases are small, i.e. RMS is comparable to standard deviation). In 2011, when the solar activity increased sharply, DTM2013 and JB2008 become significantly more accurate, their RMS dropping below 0.10. It is less pronounced for NRLMSISE00 due to its bias (Fig. 6b); in 2013, the bias becomes larger than 20% and the RMS is largest here too. Solar cycle 24 is weak with a small maximum reached in April 2014. Model performance (relative) is expected to be even better for average and strong cycles, for which the standard deviation (and RMS in case bias is of the order of a few percent) of the density ratios is expected to drop again by 0.01–0.03 at GOCE altitude.
Fig. 8.
Dailymean RMS of the density ratios of NRLMSISE00 (a), JB2008 (b), and DTM2013 (c), and the 81day mean F30 (black; right axis). JB2008 computed with proxies from June 2016. 
The daily correlations, shown for DTM2013 in Figure 9, reveal several 2–4 month intervals that stand out. The most important ones are centered on December–January (in the black ellipses) and present a clear minimum. The same intervals are detected in the NRLMSISE00 and JB2008 daily correlations (not shown). The 81day mean F30 and the daily minus the mean F30 (i.e. how the proxies are actually used in the models) are also plotted in Figure 9 in order to visually verify a possible relation with solar activity. None is evident except for the peak in activity at the end of 2011. A similar inspection with regard to geomagnetic activity (kp; not shown) also does not indicate a relation with the intervals. The minimum reached at the end of 2011 is less profound with JB2008 than with DTM2013 and NRLMSISE00, which means that at least part of the effect is due to the solar activity proxies used. A spacecraftrelated issue (thruster, attitude, and/or accelerometers) affecting the densities can be ruled out because it would have had direct consequences on the orbit altitude for example. That leaves another geophysical source of perturbations, which is waves generated in the lower atmosphere that propagate to the thermosphere. Verifying this hypothesis requires a Whole Atmosphere Model, which is beyond the scope of this study. The effect of waves and tides in the GOCE data has been analyzed (e.g. Gasperini et al. 2015; Häusler et al. 2015), although not in the context of a modification of the mean state of the thermosphere.
Fig. 9.
Dailymean correlations with DTM2013, and the 81day mean F30 and daily minus the mean F30 (purple; right axis). 
3.4. Evaluation over 24hour local solar time with CHAMP
GOCE was in a quasi Sun synchronous duskdawn orbit and, as a consequence, the densities covered only 4 h in local solar time from 6 to 8 am & pm. Therefore, the model evaluations presented thus far may not be representative of performance for all local solar times and consequently be overly optimistic, or pessimistic. As stated earlier, GOCE delivered the only highresolution and highprecision dataset at low altitude and a dataset with full local time coverage and comparable quality is only available at higher altitudes. The CHAMP densities, in the altitude range 450–250 km (but only the last month below 275 km), were used to verify the model error as a function of local solar time. Figure 10 presents the mean density ratios and the correlations using the nineyear CHAMP dataset binned in local solar time. Except for considerable but quite stable biases of JB2008 and NRLMSISE00, a typical modeling error (e.g. a diurnal variation) is not evident in the density ratios. The correlations reveal minima for the GOCE orbit configuration and are smallest in the dawn sector. Therefore, we conclude that model evaluation with GOCE densities provides performance results that are representative or even slightly pessimistic, when applied to all local solar times.
Fig. 10.
CHAMP density ratios (a) and correlations (b) binned in local solar time. JB2008 computed with proxies from June 2016. 
4. Summary and conclusions
GOCE thruster (Science Mission phase) and accelerometerinferred (3week reentry phase) densities were used to evaluate the performance of the CIRA models DTM2013, JB2008 (using the June 2016 version of solar activity proxies) and NRLMSISE00 in the 275–170 km altitude range. Three typical time scales – annual and longer, monthly and daily – were used in the evaluation in order to address different applications of the models. The range of possible values for bias, standard deviation, RMS, and correlation becomes larger when the evaluation interval decreases; this study provides these statistics based on data in the altitude range of 275–170 km (and probably applicable to a slightly wider range) that heretofore were not documented. DTM2013 performed best of the three models on all time scales. This was expected because it is the only model that has assimilated (part of) the GOCE density dataset as well as the entire CHAMP dataset. NRLMSISE00 performs worst, with considerable bias in 2009 and 2013, and systematically higher standard deviations (lower correlations) than the other two models. Nonetheless, its performance is exceptional taking the absence of recent and high accuracy density data in its construction into account. JB2008 performance is presently, thanks to the release 4_2g solar activity proxies, very close to that of DTM2013. However, it still presents some weakness under the lowest solar activity conditions in 2009 and 2010.
On the longest (4 years) time scale, all three models have small biases of less than 5% and correlations of 0.97. When the evaluation is performed on annual time scales, the range of calculated biases increases by a few percent, and the correlations remain high but decrease to 0.93–0.96. The standard deviation of the density ratios becomes smaller than 10% with DTM2013 for the years 2012 and 2013.
On monthly and daily time scales, the model biases cover an everwider range, roughly doubling from about 10% to 20%, respectively; standard deviation of the density ratios reaches about 8% (6%) for the best months (days) and about 15% (30%) for the worst with DTM2013 and JB2008; the minimum correlations drop below 0.90 and 0.75, respectively. The sometimes significantly diminished model performance on daily time scales is due to geomagnetic enhanced and storm conditions, which have typical durations of hours to a few days.
Model comparisons using the entire Science Mission dataset or the reentry dataset, with densities from 230 to 170 km, showed statistically the same performance according to our metric. While overall most accurate, the DTM2013 performance is trailing JB2008 in the reentry phase.
The CHAMP densities were used here to establish the model errors as a function of local solar time, which is not possible with only GOCE due to its quasi Sunsynchronous duskdawn orbit with 4 h of local time coverage (6–8 am & pm). The correlations were smallest for the GOCE orbit configuration, i.e. the model errors are not underestimated when comparing to the GOCE densities.
The CIRA model accuracy at low altitude from 275 to 170 km is presently quantified.
Acknowledgments
This study received funding from ESA under Contract 4000115173/15/F/MOS. The authors wish to thank the reviewers for their suggestions to improve the paper. The editor thanks Eelco Doornbos and an anonymous referee for their assistance in evaluating this paper.
References
 Bock, H., R. Dach, A. Jäggi, and G. Beutler. Highrate GPS clock corrections from CODE: support of 1 Hz applications. J. Geod., 83 (11), 1083–1094, 2009, DOI: 10.1007/s0019000903261. [CrossRef] (In the text)
 Bock, H., A. Jäggi, G. Beutler, and U. Meyer. GOCE: precise orbit determination for the entire mission. J. Geod., 88 (11), 1047–1060, 2014, DOI: 10.1007/s0019001407428. [CrossRef] (In the text)
 Bowman, B.R., W.K. Tobiska, F.A. Marcos, and C. Valladares. The JB2006 empirical thermospheric density model. J. Atmos. Sol. Terr. Phys., 70, 774–793, 2007, DOI: 10.1016/j.jastp.2007.10.002. [CrossRef] (In the text)
 Bruinsma, S., D. Tamagnan, and R. Biancale. Atmospheric densities derived from CHAMP/STAR accelerometer observations. Planet. Space Sci., 52, 297–312, 2004. [CrossRef] (In the text)
 Bruinsma, S.L. The DTM2013 thermosphere model. J. Space Weather Space Clim., 5, A1, 2015, DOI: 10.1051/swsc/2015001. [CrossRef] [EDP Sciences] (In the text)
 Bruinsma, S.L., and J.M. Forbes. Medium to largescale density variability as observed by CHAMP. Space Weather, 6, S08002, 2008, DOI: 10.1029/2008SW000411. [CrossRef] (In the text)
 Bruinsma, S.L., N. SánchezOrtiz, E. Olmedo, and N. Guijarro. Evaluation of the DTM2009 thermosphere model for benchmarking purposes. J. Space Weather Space Clim., 2, A04, 2012, DOI: 10.1051/swsc/2012005. [CrossRef] [EDP Sciences] (In the text)
 Bruinsma, S.L., E. Doornbos, and B.R. Bowman. Validation of GOCE densities and thermosphere model evaluation. Adv. Space Res., 54, 576–585, 2014, DOI: 10.1016/j.asr.2014.04.008. [CrossRef] (In the text)
 Dach, R., E. Brockmann, S. Schaer, G. Beutler, M. Meindl, L. Prange, H. Bock, A. Jaeggi, and L. Ostini, GNSS processing at CODE: status report. J. Geod., 83, 3–4, 353–365, 2009, DOI: 10.1007/s0019000802812 [CrossRef] (In the text)
 Dach, R., S. Lutz, P. Walser, and P. Fridez, Editors. Bernese GNSS software version 5.2, Astronomical Institute, University of Bern, Switzerland, 2015, DOI: 10.7892/boris.72297. (In the text)
 Doornbos, E., S. Bruinsma, B. Fritsche, G. Koppenwallner, P. Visser, J. van den IJssel, and J. de Teixeira de Encarnação. ESA contract 4000102847/NL/EL, GOCE+ Theme 3: air density and wind retrieval using GOCE data – final report, TU Delft, Delft, The Netherlands, 2014. (In the text)
 Doornbos, E., H. Klinkrad, and P. Visser. Use of twoline element data for thermosphere neutral density model calibration. Adv. Space Res., 41, 1115–1122, 2008. [CrossRef] (In the text)
 European Space Agency. Gravity Field and SteadyState Ocean Circulation Mission – The Four Candidate Earth Explorer Core Missions. Technical Report, ESA SP1233(1), European Space Agency Publications Division, Noordwijk, The Netherlands, 1999. (In the text)
 Gasperini, F., J.M. Forbes, E.N. Doornbos, and S.L. Bruinsma. Wave coupling between the lower and middle thermosphere as viewed from TIMED and GOCE. J. Geophys. Res., 120, 5788–5804, 2015, DOI: 10.1002/2015JA021300. [CrossRef] (In the text)
 T., Gruber, R. Rummel, O. Abrikosov, R. van Hees, Editors. GOCE Level 2 Product Data Handbook. GOMAHPFGS0110, Issue 4, Revision 3, European Space Agency, Noordwijk, The Netherlands, 2012. (In the text)
 Häusler, K., M.E. Hagan, J.M. Forbes, X. Zhang, E. Doornbos, S. Bruinsma, and G. Lu. Intraannual variability of tides in the thermosphere from model simulations and in situ satellite observations. J. Geophys. Res., 120, 751–765, 2015, DOI: 10.1002/2014JA020579. [CrossRef] (In the text)
 Hedin, A.E.. MSIS86 thermospheric model. J. Geophys. Res., 92, 4649–4662, 1987. [CrossRef] (In the text)
 Intelisano, A., L. Mazzini, A. Notarantonio, S. Landenna, A. Zin, L. Scaciga, and L. Marradi. Recent flight experiences of TASI onboard navigation equipments. In: Proceedings of the 4th ESA workshop on satellite navigation user equipment technologies, NAVITEC2008, 10–12 Dec 2008, Noordwijk, The Netherlands, 2008. (In the text)
 Jäggi, A., H. Bock, and R. Floberghagen. GOCE orbit predictions for SLR tracking. GPS Solut., 15 (2), 129–137, 2011, DOI: 10.1007/s1029101001766. [CrossRef] (In the text)
 Marcos, F.A. Accuracy of atmospheric drag models at low satellite altitudes. Adv. Space Res., 10 (3), 417–422, 1990. [CrossRef] (In the text)
 Marcos, F.A., D.F. Gillette, and E.C. Robinson. Evaluation of selected global thermospheric density models during low solar flux conditions. Adv. Space Res., 3, 85–89, 1983. [CrossRef] (In the text)
 Picone, J.M., A.E. Hedin, D.P. Drob, and A.C. Aikin. NRLMSISE00 empirical model of the atmosphere: statistical comparisons and scientific issues. J. Geophys. Res., 107 (A12), 1468, 2002, DOI: 10.1029/2002JA009430. [CrossRef] (In the text)
 Reigber, Ch., R. Bock, Ch. Foerste, L. Grunwaldt, N. Jakowski, H. Luehr, P. Schwintzer, and C. Tilgner. CHAMP phase B executive summary, Scientific Technical Report STR96/13, GeoForschungsZentrum Potsdam, Germany, 1996. (In the text)
 Sentman, L.H. Comparison of the exact and approximate methods for predicting free molecule aerodynamic coefficients. ARS J., 31, 1576–1579, 1961. (In the text)
 Storz, M.F., B.R. Bowman, M.J.I. Branson, S.J. Casali, and W.K. Tobiska. High accuracy satellite drag model (HASDM). Adv. Space Res., 36, 2497–2505, 2005, DOI: 10.1016/j.asr.2004.02.020. [CrossRef] (In the text)
 Svehla, D., and M. Rothacher. Kinematic precise orbit determination for gravity field determination. In: F. Sansò, Editor. A Window on the Future of Geodesy, vol. 128, Springer, Berlin, 181–188, 2005, DOI: 10.1007/3540274324_32. [CrossRef] (In the text)
 Visser, P.N.A.M., and J.A.A. van den IJssel. Calibration and validation of individual GOCE accelerometers by precise orbit determination. J. Geod., 90, 1–13, 2015, DOI: 10.1007/s0019001508500. [CrossRef] (In the text)
Cite this article as: Bruinsma S, Arnold D, Jäggi A & SánchezOrtiz N. Semiempirical thermosphere model evaluation at low altitude with GOCE densities. J. Space Weather Space Clim., 7, A4, 2017, DOI: 10.1051/swsc/2017003.
All Tables
Number of SLR observations, dailymean values, standard deviations, and RMS of SLR residuals of the kinematic orbits of the last three weeks of GOCE. SLR stations: 7810 (Zimmerwald, Switzerland), 7110 (Monument Peak, CA, USA), and 7090 (Yarragadee, Australia).
Model performance on long time scales. Results for all GOCE Science Mission data (top) and per year (bottom: 2010–2013). NB: JB2008 proxies release 4_2g were used (June 2016).
All Figures
Fig. 1.
Model predictions, in g/cm^{3}, at 250 km altitude close to winter solstice for low (left frames) and medium solar activity. 

In the text 
Fig. 2.
81day mean of the solar activity proxies F10.7, F30, and S81c (left axis), and planetary geomagnetic index ap (right axis). 

In the text 
Fig. 3.
Dailymean GOCE densities and mean spherical altitude above 6378 km (red; right axis). 

In the text 
Fig. 4.
RMS of the GOCE GPS ionospherefree phase residuals. 

In the text 
Fig. 5.
Monthly mean density ratios (a), standard deviation (std) of the density ratios (b), and correlations (c). The last (solid) symbol is with the reentry dataset. JB2008 computed with proxies from June 2016. 

In the text 
Fig. 6.
Dailymean density ratios with DTM2013 (a), NRLMSISE00 (b), and JB2008 (c) computed with proxies from June 2016. The reentry dataset is plotted in black for all three models. 

In the text 
Fig. 7.
Dailymean standard deviations (std) of the density ratios (a) and correlations (b) sorted over kp (black; right axis). JB2008 computed with proxies from June 2016. 

In the text 
Fig. 8.
Dailymean RMS of the density ratios of NRLMSISE00 (a), JB2008 (b), and DTM2013 (c), and the 81day mean F30 (black; right axis). JB2008 computed with proxies from June 2016. 

In the text 
Fig. 9.
Dailymean correlations with DTM2013, and the 81day mean F30 and daily minus the mean F30 (purple; right axis). 

In the text 
Fig. 10.
CHAMP density ratios (a) and correlations (b) binned in local solar time. JB2008 computed with proxies from June 2016. 

In the text 