New modes and mechanisms of thermospheric mass density variations from GRACE accelerometers

New modes and mechanisms of thermospheric mass density variations from GRACE accelerometers Andres Calabia 1,2 and Shuanggen Jin 1 1 Key Laboratory of Planetary Sciences, Shanghai Astronomical Observatory, Chinese Academy of Sciences, Shanghai, China, 2 University of Chinese Academy of Sciences, Beijing, China Abstract Monitoring and understanding the upper atmosphere processes is important for orbital decay and space physics. Nowadays, Low Earth Orbit (LEO) accelerometers provide a unique opportunity to study thermospheric density variations with unprecedented details. In this paper, thermospheric mass densities variations from Gravity Recovery and Climate Experiment (GRACE) accelerometers are investigated for the period 2003–2016 using the principal component analysis (PCA). The resulting modes are analyzed and parameterized in terms of solar and magnetospheric forcing, local solar time (LST), and annual variations. A better understanding of global thermospheric air density variations is presented, which validates the suitability of our technique and model. The parameterization of the subsolar-point annual variation shows two maxima around June and only one in December. The LST parameterization shows a new fluctuation controlling a middle latitude four-wave pattern, with two maxima at 12 h and 21 h LST and two minima at 1 h and 17 h LST. Our parameterizations are suitable to represent small-scale variations including, e.g., the equatorial mass density anomaly (EMA) and the midnight density maximum (MDM). Finally, the residuals are analyzed in the spectral domain, and additional contributions are found at the frequencies of the radiational tides and at the periods of 83, 93, 152, and 431 days. 1. Introduction The thermosphere is highly variable in space and time, and its physical processes are still not well understood [e.g., Jin et al., 2008, 2013]. On the one hand, the extreme ultraviolet irradiance is absorbed at the thermosphere and ionizes, dissociates, and excites the neutral species. On the other hand, the solar wind plasma, combined with a favorable alignment of the interplanetary magnetic field, produces thermospheric Joule heating and particle precipitation along the Earth’s magnetic field lines. Subsequently, since the density at this altitude is too low for molecular interaction, the dynamics act through the collisions between the ionospheric plasma and the neutrals, creating variations in the thermospheric mass density distribution. Accurate knowledge of thermospheric density variations is essential for Precise Orbit Determination (POD) in a Low Earth Orbit (LEO). Currently, half of the world’s active satellites operate in LEO (http://www.ucsusa.org/ nuclear-weapons/space-weapons/satellite-database#.V3jtOPl95QI), and more than 20,000 inactive objects are tracked for collision avoidance (http://www.orbitaldebris.jsc.nasa.gov/index.html). In a LEO environment, atmospheric drag is a major cause of orbital decay and perturbations, limiting the mission life span with a destructive reentry. Therefore, it is clear that an accurate air density model is essential for the upper atmosphere research and applications. However, in view of the high expected accuracy from present and future LEO-POD products, the current thermospheric air density models are unable to predict the variability as accurately and efficiently required. During the last decade, the technique of equaling the drag-force formula to aerodynamic accelerations has been providing an unprecedented detail when satellite-accelerometer measurements are employed [e.g., Bruinsma et al., 2004; Doornbos, 2011]. Consequently, the technical advance in this area has enlarged the research efforts on thermospheric mass density distribution and variations. For instance, Liu et al. [2005] derived a global distribution of the thermospheric total mass density from 1 year (2002) of Challenging Minisatellite Payload (CHAMP) accelerometer measurements. Several features were shown to be related to the equatorial ionization anomaly and to the midnight temperature maximum. However, only a general characterization in local solar time (LST) coordinates was provided, and no temporal variations were analyzed due to limited data. Then the concept of equatorial thermospheric mass density anomaly (EMA) was introduced in Liu et al. [2007], where the variations with season, geomagnetic activity, and solar flux levels were CALABIA AND JIN THERMOSPHERIC MASS DENSITY VARIATIONS 11,191 PUBLICATION S Journal of Geophysical Research: Space Physics RESEARCH ARTICLE 10.1002/2016JA022594 Key Points: •Thermospheric density variations are investigated from GRACE accelerometers •New periodic contributions are found at the frequencies of the radiational tides •Thermospheric mass density variations are mainly driven by radiational waves Supporting Information: •Supporting Information S1 •Data Set S1 Correspondence to: S. Jin, [email protected]; [email protected] Citation: Calabia, A., and S. Jin (2016), New modes and mechanisms of thermospheric mass density variations from GRACE accelerometers, J. Geophys. Res. Space Physics,121, 11,191–11,212, doi:10.1002/2016JA022594. Received 23 FEB 2016 Accepted 17 SEP 2016 Accepted article online 21 SEP 2016 Published online 5 NOVEMBER 2016 Corrected 24 NOV 2016 This article was corrected on 24 NOV 2016. See the end of the full text for details. ©2016. American Geophysical Union. All Rights Reserved. investigated using 4 years of CHAMP measurements (2002–2005). The EMA was defined as a minimum on the dayside, clearly aligned to the geomagnetic equator, and with two maxima at ±20° geomagnetic latitude. However, only the geographical distribution of the dayside mass density was depicted, and no variations at high latitude (above 60°) were studied. Later, Lei et al. [2012b] showed that the equatorial thermosphere anomaly (ETA) can be attributed to energy transfer from thermal plasma to the neutrals through collisions, due to their temperature differences. Recently, an anomalous behavior of the equatorial anomaly has been investigated in the nightside mass density distribution [e.g., Ma et al., 2010; Akmaev et al., 2010; Ruan et al., 2014, 2015], i.e., the midnight density maximum (MDM). It seems that the MDM has a maximum, instead of minimum around the equator, and lower values instead of higher values, at middle latitudes. More complete studies on changes with geophysical conditions have been presented, e.g., in Müller et al. [2009], Guo et al. [2007], and Lathuillère et al. [2008]. For instance, Müller et al. [2009] studied the solar and magnetospheric forcing of the low-latitude thermospheric mass density as observed by CHAMP accelerometers, and provided a detailed study on the day-to-night mass density ratio and the seasonal variations. Unfortunately, their day-to-night study was insufficient to properly describe local time variations, but the seasonal dependence was better characterized than that of previous studies [e.g., Bowman et al., 2008a; Guo et al., 2008]. The improvement was achieved by removing the local time, solar flux, and magnetic activity effects before the analysis. In recent years, several other studies have investigated the seasonal variation of thermospheric neutrals using different and innovative techniques. For instance, Guo et al. [2008] deduced the seasonal variations from CHAMP accelerometer measurements (2002–2005), using daily averages compared to the JB2006 model [Bowman et al., 2008b], in which the seasonal dependence was deactivated. Unfortunately, their results showed strong dependence on the performance of the empirical model JB2006. Several studies have detected that the amplitude of seasonal variations increases with solar activity [e.g., Emmert and Picone, 2010; Matsuo and Forbes, 2010], and a further study of on the latitudinal dependence of seasonal variations was performed in Lei et al. [2012a]. Concerning longitudinal variations, Xu et al. [2013] revealed from CHAMP and GRACE accelerometers a hemispherically asymmetric distribution, positive always near the magnetic poles and with strong seasonal variations. From the above studies, it is clear that more comprehensive studies are required to understand all involved parameters and the interconnections between the functional dependences. In addition, the mechanisms related to equatorial anomalies are still not well defined [e.g., Ruan et al., 2014, 2015; Akmaev et al., 2010; Ma et al., 2010; Liu et al., 2007]. It is clear that the above research has been addressed to different aspects of thermospheric density variations, but the scientific community still needs more complete studies to characterize and predict the thermospheric density variability in an accurate and efficient manner. In the spatiotemporal data analysis of satellite measurements, the orbital precession represents a major limitation. For that reason, the studies of Matsuo and Forbes [2010] and Lei et al. [2012a] have provided the first attempts to characterize the dominant modes of global thermosphere density variability, with the use of the principal component analysis (PCA) technique. In their work, the principal components were estimated with a sequential nonlinear regression analysis of one-dimensional measurements along satellite tracks, using a time-invariant orthogonal basis at every orbit time. Providing similar results, Calabia and Jin [2016] have recently introduced a more feasible technique for the PCA of sparse measurements along satellite tracks. The method is based on the eigenvalue solution of the covariance matrix of static grids, which are obtained by interpolation of consecutive arcs of orbits. In this scheme, the spatiotemporal behavior of thermospheric air density variations with solar and magnetospheric forcing, as well as LST and annual fluctuations, can be investigated together with all driver-response relationships. In this study, thermospheric neutral density variations from Gravity Recovery and Climate Experiment (GRACE) accelerometry are investigated for the period 2003–2016 using the PCA. The structure of this paper is presented as follows. In section 2, a brief overview of neutral density retrieval from GRACE measurements is given; the detailed technique for the analysis is presented in section 3; the results are analyzed and discussed in section 4; and finally, conclusions are summarized in section 5. 2. Observations and Data Processing 2.1. Observations and Models GRACE (Gravity Recovery and Climate Experiment) is a joint mission between the National Aeronautics and Space Administration (NASA) and the Deutsches Zentrum Für Luft und Raumfahrt (DLR). GRACE’s Journal of Geophysical Research: Space Physics 10.1002/2016JA022594 CALABIA AND JIN THERMOSPHERIC MASS DENSITY VARIATIONS 11,192 21699402, 2016, 11, Downloaded from https://agupubs.onlinelibrary.wiley.com/doi/10.1002/2016JA022594 by Spanish Cochrane National Provision (Ministerio de Sanidad), Wiley Online Library on [28/11/2025]. See the Terms and Conditions (https://onlinelibrary.wiley.com/terms-and-conditions) on Wiley Online Library for rules of use; OA articles are governed by the applicable Creative Commons License Level 1B format record files can be downloaded from the Information System and Data Center (ISDC) GeoForschungsZentrum (GFZ) website (http://isdc.gfz-potsdam.de/) in binary big-endian byte-ordering format. The twin satellites of the GRACE mission are equipped with three-axis capacitive Super-STAR accelerometers to measure the nongravitational acceleration acting on the satellites. Plus and minus drive voltages are applied to electrodes with respect to opposite sides of a proof mass, whose electrical potential is maintained at a dc biasing voltage. Unfortunately, this dc level is the source of bias and bias fluctuations which must be corrected, e.g., using the parameterizations of Calabia et al. [2015]. The star cameras mounted on GRACE satellites provide the precise attitude references to determine the satellite’s absolute orientation with respect to the International Celestial Reference System. Precise orbit ephemeris (position and velocity) are generated by processing 5 min GPS tracking data, through reduced-dynamic POD with stochastic acceleration estimated every 5 min [Bertiger et al., 2002]. Precise orbit ephemeris, thruster activation times, and satellite mass records are also included within the record files. NRLMSISE-00 is an empirical, global model of the Earth’s atmosphere, from ground to space, and is the standard for international space research. A primary use is to aid predictions of satellite orbital decay due to atmospheric drag. The earlier models MSIS-86 and MSISE-90 are based on Mass Spectrometer and Incoherent Scatter Radar measurements, and the current model has been updated with satellite drag data [Picone et al., 2002]. 2.2. Thermospheric Air Mass Density Estimation In this study, thermospheric mass densities are derived from accelerometer measurements of GRACE (GFZ database) using the drag-force formula: FD¼1 2CAρv2 r(1) where Cis the drag coefficient vector, Ais the cross-sectional area perpendicular to v r ,ρis the atmospheric mass density, and v r is the relative velocity of the atmosphere with respect to the spacecraft, including corotating atmosphere and horizontal winds [Doornbos, 2011]. Horizontal wind velocities have been computed from the HWM07 model [Drob et al., 2008] in the local ENU (East, North, Up) reference frame, and transformed to the terrestrial reference frame (International Terrestrial Reference System) to be combined with satellite’s velocities. The velocity of the corotating atmosphere has been computed as the vector product between the Earth’s angular rotation and the satellite’s position vector. The forces measured by satellite accelerometers mainly include atmospheric drag (F D ), acceleration due to radiation pressure, and thruster firings. While thruster firings have been easily removed by computing the activation times, irradiative accelerations comprise a more elaborated modeling. Basically, the three main sources of irradiative accelerations are direct solar radiation pressure, reflected solar radiation pressure, and terrestrial infrared radiation. While terrestrial infrared radiation is almost independent from illumination conditions, direct and reflected solar radiations have to be multiplied by a ratio to account for planetary eclipses [Montenbruck and Gill, 2013]. On the plates of the user’s satellite, one part of the incoming radiation is absorbed and the other is reflected diffusely and specularly. The equation to determine the entire resultant force on the satellite due to solar radiation [Luthcke et al., 1997] accounts for the plate areas and their orientation, their coefficients of diffusive and specular reflectivity, and the mass of the satellite. GRACE’s panel model properties [Bettadpur, 2007] provide the coefficients of reflectivity for both the visible (VIS) and the infrared (IR) part of the electromagnetic spectrum. We have combined these values by weighting them in function of the amount of solar flux given in each spectral window (43% for VIS and 53% for IR). In a similar way, Earth albedo has been computed accounting for the reflected solar radiation and the terrestrial infrared radiation. For the terrestrial infrared radiation, Knocke and Ries [1987] modeled seasonal and latitudinal variations of a blackbody with a surface temperature of 288°K. The reflected radiant flux is a fraction of the incoming flux, and it has been computed from the reflectivity index available from the Total Ozone Mapping Spectrometer. GRACE’s drag coefficient vector Chas been computed from the energy accommodation coefficient model generated by Mehta et al. [2013]. For the molecular mass ratio between the mean atmosphere and each plate, note that the surface of GRACE is covered with Kapton (382 g/mol), solar array glass (144 g/mol), and Teflon (100 g/mol). For the mean molecular density (g/mol) of the atmosphere, we employ the sum of partial number densities (m 3 ) from the MSIS model, multiplied by each corresponding Journal of Geophysical Research: Space Physics 10.1002/2016JA022594 CALABIA AND JIN THERMOSPHERIC MASS DENSITY VARIATIONS 11,193 21699402, 2016, 11, Downloaded from https://agupubs.onlinelibrary.wiley.com/doi/10.1002/2016JA022594 by Spanish Cochrane National Provision (Ministerio de Sanidad), Wiley Online Library on [28/11/2025]. See the Terms and Conditions (https://onlinelibrary.wiley.com/terms-and-conditions) on Wiley Online Library for rules of use; OA articles are governed by the applicable Creative Commons License atomic mass (g/mol), and finally divided by the total volume. Detailed study on using different energy accommodation coefficients and the repercussions on GRACE’s density estimation can be found in Mehta et al. [2013]. After solving the drag-force equation for neutral densities at each satellite position, the derived values have been normalized to a common height of 475 km with the MSIS model, following the indications of Rentz and Lühr [2008]. Since the GRACE orbit altitude varies from 500 km in 2003 to 400 km in 2015 (scale height of 100 km), the errors caused by the normalization are expected to be small enough (within 5%), as discussed in Bruinsma et al. [2006]. Altitude information for the GRACE mission can be found at the website of the University of Texas, Center for Space Research (http://www.csr.utexas.edu/ grace/operations/configuration.html). 3. Analysis Methods 3.1. Spatiotemporal Analysis In order to derive a complete grid of densities, the orbital trajectory has been divided in ascending and descending legs. For the development in longitude, the values have been biased 360° after a grid is completed. Then the interpolation and clipping of each grid has been conducted. Taking into account time processing and required physical memory for the PCA analysis, we achieved an optimal resolution of 3° (120 × 60). Detailed procedure for the computation can be found in Calabia and Jin [2016]. Following to this approach, the variability of each grid is synthesized in two temporal dimensions: the latitudinal variation, defined by the half orbital period (3.99 min per westward degree), and the longitudinal variation, defined by the equatorial orbit shift (31.37 s per latitudinal degree). The spatial patterns of variability, their time variation, and the measure of their importance are conducted via the PCA. The aim of this technique is to determine a new set of variables that capture most of the observed variance through a linear combination of the original variables, based on eigen decomposition. The analysis is presented as a series of modes which can be reconstructed multiplying the spatial patterns of variability with their corresponding temporal components. Detailed analyses and the selection of retained modes for static grids can be found in Preisendorfer [1988] and Wilks [1995], and a readily computable algorithm in Bjornsson and Venegas [1997]. The method requires the use of a two-dimensional matrix where each row represents a given instant of time, and each column a time series for a given location. In this study, each grid has been expanded into one-dimensional array (120 × 60 = 7200), and a matrix of 7200 × 4762 values has been created from the 4762 grids to analyze (grids from 2003 to 2016). Note that each grid corresponds to 23 h 56 min, so we obtain 2 × 17 grids extra in the period 2003–2016. In addition, ascending and descending orbits have been concatenated to obtain a singular solution, so the final matrix becomes 7200 × 9524. Then the covariance matrix has been computed for the eigenvalue problem [Bjornsson and Venegas, 1997]. The solution is a diagonal matrix containing the eigenvalues of the covariance matrix and a matrix of the corresponding eigenvectors. Each eigenvectors is regarded as a map, and its contribution to the total variability is given by the corresponding eigenvalue. For the purpose of this research, subdaily variations are not investigated, and since all-derived grids are equally affected by longitudinal and latitudinal variations, the resulting PCA components will follow the same pattern. Additional details about the features of this method can be found in Calabia and Jin [2016]. 3.2. Radiational Constituents In Munk and Cartwright [1966, Figure 8], the authors introduced the concept of “radiational tides”to associate their residuals with the tidal phenomena derived from solar heating. Therefore, since the driving force of radiational variations in the thermosphere should contain the same spectral structure as the main radiation tides, it seems to be meaningful to look into the thermospheric air mass density time series for peaks at frequencies, which are marked in the tidal spectrum by the most energetic radiational constituents (Table 1). From the perspective of the GRACE satellites (subscript sat), the wave frequencies associated to each radiational constituent vary depending on the satellite’s equatorial orbit shift, i.e., GRACE’s precession rate σSun sat ¼1=322 cycle=dayðÞ. From the perspective of a ground-based station (subscript Earth), the Earth’s self-rotation frequency with respect to the Sun is σSun Earth ¼1 cycle=dayðÞand with respect to, e.g., the P1 radiational constituent is σP1 Earth ¼σSun Earth 1=365:2434, so 1=σP1 Earth ¼1:002745 or 24.06589 h, which is confirmed as the solar diurnal P1 period [Petit and Luzum, 2010]. Since the phase speed of the radiational Journal of Geophysical Research: Space Physics 10.1002/2016JA022594 CALABIA AND JIN THERMOSPHERIC MASS DENSITY VARIATIONS 11,194 21699402, 2016, 11, Downloaded from https://agupubs.onlinelibrary.wiley.com/doi/10.1002/2016JA022594 by Spanish Cochrane National Provision (Ministerio de Sanidad), Wiley Online Library on [28/11/2025]. See the Terms and Conditions (https://onlinelibrary.wiley.com/terms-and-conditions) on Wiley Online Library for rules of use; OA articles are governed by the applicable Creative Commons License constituent (e.g., in the P1 case C P1 = 1/365.2434) is constant for each sinusoidal mode [Forbes et al., 2013], by subtracting the equations derived from both perspectives ref ={sat, Earth}, i.e., σP1 ref ¼σSun ref CP1,we can isolate the satellite-associated wave frequencies, e.g., equation (2). Note that equation (2) is given as example to calculate the wave frequency associated to the P1 radiational constituent. In the case of GRACE, the precession rate has a negative value, reflecting the opposite precession of the satellite than that given by the Earth, both with respect to the Sun. σP1 sat ¼σSun sat þσP1 Earth σSun Earth ¼0:005843 ¼ 1 171:1306 (2) The corresponding periods for the radiational constituents are extracted from Petit and Luzum [2010] and listed in Table 1, where we also include the corresponding alias period for GRACE. 4. Results and Discussion 4.1. Thermospheric Air Mass Density Estimation Since we have generated a new density data set, we provide a brief statistical comparison with respect to a previous density derivation [Sutton, 2011], as well as the differences to NRLMSISE-00. We employ density estimations derived from Sutton [2011], which have been downloaded from the University of Colorado at Boulder website (http://sisko.colorado.edu/sutton/data/ver2.2/) and have been normalized to 475 km with the MSIS model, following the indications of Rentz and Lühr [2008]. The estimates of Sutton [2011] are the solution of Sentman’s [1961] drag coefficient formulas, which include the thermal drag effect (phenomenon that increases the coefficient of drag for long satellites). The resulting statistical comparisons are presented in Figure 1, with correlation coefficients, relative error (with respect to the background density), and relative standard deviation (SD) every 15 day from 2003 to 2009. Note these relative statistics are strongly dependent on the satellite’s LST, as well as solar flux conditions (i.e., satellite positions at around 6 or 18 h LST or low values in the background density will influence negatively in the relative statistics). With a mean correlation of 94%, the densities of Sutton [2011] have 10% better correlation to our estimates than the MSIS model. In Figure 1c, while the mean difference with respect to Sutton [2011] is centered at zero, the MSIS estimates are positively deviated in 25% (with respect to the background density). Concerning the relative SD (Figure 1d), with a mean value of 16% with respect to the background density, the differences with respect to Sutton [2011] are 11% smaller than the differences to MSIS. Plotted in geographical coordinates, Figure 2 shows the relative differences between the GRACE measurements and the MSIS model during the period 2003–2016. In Figure 2a, the median average of relative differences, i.e., (MSIS-GRACE/GRACE), shows a mean global overestimation of about 0.15, being less pronounced at the northern polar cap and at the Cape Horn region, but more prominent at the equator (~0.3). The SD of relative differences (Figure 2b) shows a mean value of 0.35, showing marked values at the southern polar Table 1. GRACE Alias Periods for the Main Radiation Waves Darwin Symbol Period (day) Alias Period (day) P1 1.0027454 171.13 S1 (LST) 1.000000000 322.00 K1 0.9972696 2719.68 K+ 0.9971964 3399.74 T2 0.5006854 111.74 R2 0.4993165 287.89 Figure 1. Statistical comparison with respect the estimates of Sutton [2011] in black and with respect the NRLMSISE00 empirical model in green. (a) Background density, (b) correlation coefficients, (c) relative error, and (d) relative SD. Journal of Geophysical Research: Space Physics 10.1002/2016JA022594 CALABIA AND JIN THERMOSPHERIC MASS DENSITY VARIATIONS 11,195 21699402, 2016, 11, Downloaded from https://agupubs.onlinelibrary.wiley.com/doi/10.1002/2016JA022594 by Spanish Cochrane National Provision (Ministerio de Sanidad), Wiley Online Library on [28/11/2025]. See the Terms and Conditions (https://onlinelibrary.wiley.com/terms-and-conditions) on Wiley Online Library for rules of use; OA articles are governed by the applicable Creative Commons License region and along the equator (~0.4), and lower values at the northern polar region (~0.3). 4.2. Spatiotemporal Analysis of Thermospheric Densities From PCA The spatial patterns of density variability, their time variation, and the measure of their importance have been conducted via the PCA. The eight main eigenvectors together account for 98.6% of the total variance and, individually, explain the 90.3%, 3.5%, 2.9%, 1.2%, 0.3%, and 0.1% for each of the last three ones. The fourth, fifth and sixth modes have been recognized as the gradient of density variations driven by energetic geomagnetic contributions and therefore excluded from our analysis (no clear periodicities to investigate were detected). We decided to study these nonperiodic components in future research. Figure 3 show the five selected modes of thermospheric density variability in geographic coordinates (i.e., to represent a specific contribution, e.g., local time or day of the year, each map must be multiplied by the queried value, which is given by the corresponding time expansion coefficient). The first PCA component (Figure 3a) shows higher values in the Southern Hemisphere, a trough that follows the geomagnetic equator, and two peaks located at the cusps. Maxima values are located at the southern cusp and minima at the northern cusp. We might induce that the distribution is coupled with the orientation of the dip angle, producing positive density contributions in the southern cusp and negative in the northern cusp, both with respect to the variable global mean value. Similar results have been pictured in Xu et al. [2013], with a hemispherically asymmetric distribution, showing higher values near the magnetic poles. Previous studies using PCA techniques [Matsuo and Forbes, 2010; Lei et al., 2012a] provided no clear descriptions for these features. In comparison with Liu et al. [2005, 2007, 2009], our PCA1 shows an improved characterization of the EMA, which shows the clear relation with the ETA [Lei et al., 2012b]. Since the PCA technique separates in a different component the variations due to the movement of the subsolar point, the EMA structure, and its fluctuations in time can be better analyzed. For instance, Figure 3a shows more pronounced and symmetric crests where the geomagnetic equator is close to the geographical equator (i.e., 90° to 180°W); while in the other side (i.e., 0° to 180°E) the northern crest becomes less prominent. The second and third components (Figures 3b and 3c) show a clear relation to latitudinal variations following the subsolar point, plus changes in latitudinal curvature due to LST variations. Figures 3d and 3e show more detailed variations in latitude, and a clear four-peak longitudinal pattern, which has previously reported in Liu et al. [2009]. A detailed analysis of the frequency spectrum of time domain signals (right column in Figure 3) is given in the followings sections. 4.3. Parameterization of PCA Time Expansion Coefficients Since the strongest contribution on thermospheric density variations is the solar and magnetospheric forcing, the main time expansion coefficient (PCA1, Figure 3a) has been simultaneously modeled in a Figure 2. (a) Median average and (b) SD of relative differences between the NRLMSISE00 and the GRACE density estimates, i.e., (NRLMSISE00 GRACE/ GRACE). Dip isoclinic lines are plotted in dash-dotted gray lines to show the alignments. Values are dimensionless. Journal of Geophysical Research: Space Physics 10.1002/2016JA022594 CALABIA AND JIN THERMOSPHERIC MASS DENSITY VARIATIONS 11,196 21699402, 2016, 11, Downloaded from https://agupubs.onlinelibrary.wiley.com/doi/10.1002/2016JA022594 by Spanish Cochrane National Provision (Ministerio de Sanidad), Wiley Online Library on [28/11/2025]. See the Terms and Conditions (https://onlinelibrary.wiley.com/terms-and-conditions) on Wiley Online Library for rules of use; OA articles are governed by the applicable Creative Commons License Figure 3. Main PCA components of thermospheric air mass density variability from GRACE measurements (475 km altitude) (a–e) Each component individually explains the 90.3%, 3.5%, 2.9%, 0.1, and 0.1% of the total variance. Dip isoclinic lines are plotted in dash-dotted gray lines to show the alignments. Values are dimensionless. Journal of Geophysical Research: Space Physics 10.1002/2016JA022594 CALABIA AND JIN THERMOSPHERIC MASS DENSITY VARIATIONS 11,197 21699402, 2016, 11, Downloaded from https://agupubs.onlinelibrary.wiley.com/doi/10.1002/2016JA022594 by Spanish Cochrane National Provision (Ministerio de Sanidad), Wiley Online Library on [28/11/2025]. See the Terms and Conditions (https://onlinelibrary.wiley.com/terms-and-conditions) on Wiley Online Library for rules of use; OA articles are governed by the applicable Creative Commons License two-dimensional degree-2 polynomial fitting (Figure 4). Instead the lineal response given by previous studies [e.g., Müller et al., 2009, Figure 3; Ruan et al., 2014 Figure 2], Figure 4 depicts a clear quadratic response. The fitting is based on applying a robust least squares regression analysis, where the robust linear fitting M-estimator method Tukey’s biweight (also known as bisquare) is employed to avoid possible outliers. The polynomial fitting uses the coefficients given in Table 2 in the following form: FAm;P10:7 ðÞ¼p00 þp10 Am þp01 P10:7þp20 Am2þp11 Am P10:7þp02 P10:72(3) where Am is the k-derived planetary geomagnetic index (http://isgi.unistra.fr/data_download.php), and P 10.7 is the composite index P 10.7 = 0.5 (F 10.7 +(F 10.781 days )) from the previous day, which is derived from the F 10.7 solar radio flux at 10.7 cm (http://omniweb.gsfc.nasa.gov/form/dx1.html). Previous authors have shown that the P 10.7 (previous day) and the Am indices better correlate with the thermospheric density variations [e.g., Müller et al., 2009; Guo et al., 2007]. After removing solar and magnetospheric contributions from the first time expansion PCA coefficient, the residuals have shown variations at LST and at annual periods, which are modulated in amplitude by a factor strongly correlated to the solar flux. These variations have been extracted from the signals, and parameterized in terms of LST and DOY (day of year), in sinusoidal functions modulated by the P 10.7 index. Note throughout this manuscript that each function has an additive contribution to the total variance, and the linear combination of all the functions associated to a specific contribution, e.g., LST, represents the complete spectrum of that contribution. In order to reduce the number of variables in the fitting process, an iterative two-step process has been conducted as follows: First, data normalization to common flux has been implemented following the indications in Müller et al. [2009], then a sinusoidal fitting based on Fourier functions has been conducted as done in Jin et al. [2007], and finally, a polynomial fitting modulates the flux-dependent amplitude of the functions computed in the previous step. Note that all cases employ the bisquare least squares fitting (same as for equation (3)). The resulting functions for the main secular variations are shown in Figures 5, 6, and 8 to 11. For these variations, the coefficients given in Table 3 are to be used in the following form: Gχ;P10:7 ðÞ¼1015 10aP10:7ba0þX n i¼1 an cos nχðÞþbn sin nχðÞ½ ! (4) where a,b,a0, an,andbn are the constant and amplitudes, and χtakes different values for each periodic case: For the annual variations, χis DOY · 2π/365; for the radiational variations (e.g., LST), χis the argument defined by the vector product between the six vectors of Doodson’s fundamental arguments and their corresponding six vectors of multipliers [Petit and Luzum, 2010]; for the other periodic cases, χis the product Figure 4. Fitting of the first PCA time expansion coefficient (Figure 3a) in terms of P 10.7 and Am indices. Parameterization is given in Table 2. Corresponding spatial pattern is shown in Figure 3a and given in Table 4. Yaxis values are dimensionless. Table 2. Coefficients for Solar and Magnetospheric Variations Used in Equation (3) PCA 1 P 10.7 and Am p00 1.98E 14 p10 5.55E 16 p01 1.62E 17 p20 6.09E 18 p11 1.45E 17 p02 7.10E 18 Journal of Geophysical Research: Space Physics 10.1002/2016JA022594 CALABIA AND JIN THERMOSPHERIC MASS DENSITY VARIATIONS 11,198 21699402, 2016, 11, Downloaded from https://agupubs.onlinelibrary.wiley.com/doi/10.1002/2016JA022594 by Spanish Cochrane National Provision (Ministerio de Sanidad), Wiley Online Library on [28/11/2025]. See the Terms and Conditions (https://onlinelibrary.wiley.com/terms-and-conditions) on Wiley Online Library for rules of use; OA articles are governed by the applicable Creative Commons License between the corresponding angular frequency (= 2π/period) and the days since J2000.0 (days since 51544.5 in modified Julian date format). Note that Doodson’s fundamental arguments are computed at the times when the satellite crosses the Greenwich meridian. Fitted parameters and correlation coefficients between the time expansion PCA series and the parameterized ones are given in Table 3. The parameterizations achieve 96%, 92%, 91%, 75%, and 75% correlation with respect to the 1 to 5 PCA components, respectively. Figure 5. Fitting of PCA time expansion coefficients, as seen from periodic variations of the S1 constituent (represented in LST). Data normalized to P 10.7 = 110 common flux. Parameterizations are given in Table 3, including the modulation in amplitude. Corresponding spatial patterns are shown in Figure 3 and given in Table 4. Yaxis values are dimensionless. Journal of Geophysical Research: Space Physics 10.1002/2016JA022594 CALABIA AND JIN THERMOSPHERIC MASS DENSITY VARIATIONS 11,199 21699402, 2016, 11, Downloaded from https://agupubs.onlinelibrary.wiley.com/doi/10.1002/2016JA022594 by Spanish Cochrane National Provision (Ministerio de Sanidad), Wiley Online Library on [28/11/2025]. See the Terms and Conditions (https://onlinelibrary.wiley.com/terms-and-conditions) on Wiley Online Library for rules of use; OA articles are governed by the applicable Creative Commons License of 152 day period can be seen in PCA1 and PCA4; the presence of 171 day period (P) is clear in all the PCA components; the period of 290 days (R) is clear in the PCA4; a 431 day periodicity can be seen in PCA1; and additional periodicities are found at 1700 (K+) and 2719 days (K). In order to explain the peaks of periodicities in a scientific manner, we studied if a possible geophysical explanation could be given with the constituents of the theory of tides. The suppositions are based on driving forces of radiational origin, which should contain the same spectral signature as the main radiation tides. In Munk and Cartwright [1966, Figure 8], the most energetic radiation was depicted by the T, P, R, S, and K constituents. From the perspective of GRACE, the wave periods associated to these radiational constituents are 112, 171, 288, 322, and 2720 days, respectively (shown in Table 1). Validating this hypothesis with our Figure 9. Fitting of PCA time expansion coefficients, as seen from periodic 2731 day variations (Kwave). Data are normalized to P 10.7 = 110 common flux. Parameterizations are given in Table 3, including the modulation in amplitude. Corresponding spatial patterns are shown in Figure 3 and given in Table 4. Yaxis values are dimensionless. Figure 10. Fitting of PCA time expansion coefficients, as seen from periodic 283 day variations (R2 wave). Data are normalized to P 10.7 = 110 common flux. Parameterizations are given in Table 3, including the modulation in amplitude. Corresponding spatial patterns are shown in Figure 3 and given in Table 4. Yaxis values are dimensionless. Journal of Geophysical Research: Space Physics 10.1002/2016JA022594 CALABIA AND JIN THERMOSPHERIC MASS DENSITY VARIATIONS 11,206 21699402, 2016, 11, Downloaded from https://agupubs.onlinelibrary.wiley.com/doi/10.1002/2016JA022594 by Spanish Cochrane National Provision (Ministerio de Sanidad), Wiley Online Library on [28/11/2025]. See the Terms and Conditions (https://onlinelibrary.wiley.com/terms-and-conditions) on Wiley Online Library for rules of use; OA articles are governed by the applicable Creative Commons License results, the 322 day period (S) has been already removed by LST variations in previous section, and the others are clearly present in the residuals (Figure 7). The resulting sinusoidal functions for the P1, K, R2, and T2 variations are plotted in Figures 8 to 11, respectively. Note that the number next to each period (e.g., R2) indicates diurnal (1) or semidiurnal (2) species. Concerning the period of 93 days, this contribution could be originated by the GRACE’s drift of perigee (having the same period). We note that unlike for the radiational periods, NRLMSISE-00 did not include any periodicity at 93 days. We might induce that scaling the densities at 475 km has not completely removed the variations produced by orbiting at different altitudes (see last paragraph in section 2.2). The period of 431 days could be related with the free-core nutation theory. We suspect that the core-mantle electromagnetic dynamo (interaction of the mantle and the fluid, ellipsoidal core as it rotates) produces changes in the magnetic field and therefore induces variations through the magnetosphere-ionosphere-thermosphere (MIT) coupling (the collisions between the ionospheric plasma and the neutrals fluctuate in accordance to the 431 day induced magnetic field). Concerning the 83 and 152 day periods, several studies have reported similar periodicities in solar activity [e.g., Cane et al., 1998; Joshi and Joshi, 2005], which could be the precursors of thermospheric density variations at these periods. 4.6. Parameterized Results Parameterizations for the main variations are given in Tables 2 and 3, to be used with equations (3) and (4), as time-varying coefficients for the main spatial patterns given in Table 4. A feasible MATLAB computation to obtain the grids of these spatial patterns is given in Appendix A. Figure 14 shows the global averaged relative residuals (i.e., residual/density) for both the MSIS model and our PCA parameterizations. The relative residuals for the PCA parameterizations are ~60% smaller than those of NRLMSISE-00 and show values below the ~20%. 5. Summary In this paper, a 13 year time series of thermospheric mass densities inferred from the GRACE accelerometry is investigated via PCA. Our results are in agreement with that given by previous studies and show a better characterization of the global thermospheric air mass density distribution and variations than that represented by the NRLMSISE00 empirical model. Overall, the differences between the analyzed densities and the MSIS model are significant in both structure and magnitude. All selected PCAs have shown to be strongly controlled by solar flux and geomagnetic activity, and modulated by radiation waves, and annual fluctuations. At first, the selected PCA time expansion coefficients have been parameterized in terms of LST, DOY, and solar P 10.7 flux and geomagnetic Am indices. The modeling of all periodic variations is the modulation in amplitude of sinusoidal functions, which have been previously fitted into a common flux normalized data. The parameterization of the subsolar-point annual variation shows two maxima around June (one at early May and a reduced at the end of July) and only one in December. We found a middle latitude LST fluctuation Figure 11. Fitting of PCA time expansion coefficients, as seen from periodic 112 day variations (T2 wave). Data are normalized to P 10.7 = 110 common flux. Parameterizations are given in Table 3, including the modulation in amplitude. Corresponding spatial patterns are shown in Figure 3 and given in Table 4. Yaxis values are dimensionless. Journal of Geophysical Research: Space Physics 10.1002/2016JA022594 CALABIA AND JIN THERMOSPHERIC MASS DENSITY VARIATIONS 11,207 21699402, 2016, 11, Downloaded from https://agupubs.onlinelibrary.wiley.com/doi/10.1002/2016JA022594 by Spanish Cochrane National Provision (Ministerio de Sanidad), Wiley Online Library on [28/11/2025]. See the Terms and Conditions (https://onlinelibrary.wiley.com/terms-and-conditions) on Wiley Online Library for rules of use; OA articles are governed by the applicable Creative Commons License Figure 12. Thermospheric density variations at 475 km altitude, in function of LST and latitude (fixed at Greenwich meridian), and for different season and solar flux conditions. (a) P 10.7 = 80 sfu and (b) P 10.7 = 120 sfu. From top to bottom, it plots at ME, JS, SE, and DS. In each scenario (eight cases), from left to right, averaged accelerometer and MSIS-based densities along GRACE orbits (raster plots) and direct calculations (contour plots) from our parameterization and from MSIS. Each bar graph corresponds to the occurrences between ±15 sfu and ±3 h LST with respect to each scenario (for raster plots only). The direct calculations using the models have been computed with the representative parameters of each scenario. Journal of Geophysical Research: Space Physics 10.1002/2016JA022594 CALABIA AND JIN THERMOSPHERIC MASS DENSITY VARIATIONS 11,208 21699402, 2016, 11, Downloaded from https://agupubs.onlinelibrary.wiley.com/doi/10.1002/2016JA022594 by Spanish Cochrane National Provision (Ministerio de Sanidad), Wiley Online Library on [28/11/2025]. See the Terms and Conditions (https://onlinelibrary.wiley.com/terms-and-conditions) on Wiley Online Library for rules of use; OA articles are governed by the applicable Creative Commons License Figure 13. Thermospheric density variations at 475 km altitude, in function of DOY and latitude (fixed at Greenwich meridian), and for different LST and solar flux conditions. (a) P 10.7 = 80 sfu and (b) P 10.7 = 120 sfu. From top to bottom, plots at 5, 11, 17, and 23 h LST, In each scenario (eight cases), from left to right, averaged accelerometer and MSIS-based densities along GRACE orbits (raster plots) and direct calculations (contour plots) from our parameterization and from MSIS. Each bar graph corresponds to the occurrences between ±15 sfu and ±3 h LST with respect to each scenario (for raster plots only). The direct calculations using the models have been computed with the representative parameters of each scenario. Journal of Geophysical Research: Space Physics 10.1002/2016JA022594 CALABIA AND JIN THERMOSPHERIC MASS DENSITY VARIATIONS 11,209 21699402, 2016, 11, Downloaded from https://agupubs.onlinelibrary.wiley.com/doi/10.1002/2016JA022594 by Spanish Cochrane National Provision (Ministerio de Sanidad), Wiley Online Library on [28/11/2025]. See the Terms and Conditions (https://onlinelibrary.wiley.com/terms-and-conditions) on Wiley Online Library for rules of use; OA articles are governed by the applicable Creative Commons License with two maxima at 12 h and 21 h LST, and two minima at 1 h and 17 h LST. Our model is suitable to represent small-scale variations including, e.g., EMA and MDM, and can be used to improve the current MIT modeling. Finally, the residuals have been analyzed in the spectral domain, and additional periodic contributions have been found at the frequencies of the radiational tides. In addition, periodic contributions are found at the periods of 83, 93, 152, and 431 days. The 93 day period could be caused by the satellite’s drift of perigee, and the 83 and 152 day periods might be attributed to solar activity. Variations at the free-core nutation frequency (431 days) suggest a possible core-magnetosphere-ionosphere-thermosphere (CMIT) coupling. Appendix A: MATLAB Computation to Obtain Grids From Coefficients Given in Table 4 function [DATA] = Stokes2grid(Alm,Blm) degree = 8; phi = 90:-1:-90; phi = phi’*pi./180; [plm ~] = legnorm(phi,degree); rho = 180/pi; for lon = 0:359 DATA(:,lon + 1) = zeros(181,1); for n = 0:degree for m = 0:n Pl(:,:) = plm(n + 1,m + 1,:); aux = Pl.*cos(m*lon/rho)*Alm(n + 1,m + 1) + Pl.*sin(m*lon/rho)*Blm(n + 1,m + 1); DATA(:,lon + 1) = DATA(:,lon + 1) + aux; end end end References Akmaev, R. A., F. Wu, T. J. Fuller-Rowell, H. Wang, and M. D. Iredell (2010), Midnight density and temperature maxima, and thermospheric dynamics in Whole Atmosphere Model simulations, J. Geophys. Res.,115, A08326, doi:10.1029/2010JA015651. Bertiger, W., et al. (2002), GRACE: Millimeters and microns in orbit, In: Proceedings of ION GPS 2002, Portland, OR, September 24–27, 2002, pp. 2022-2029. Bettadpur, S. (2007), Product specification document, Technical Report CSR-GR-03-02, Center for Space Research, The Univ. of Texas at Austin. Bjornsson, H., and S. A. Venegas (1997), A manual for EOF and SVD analyses of climatic data, MCGill Univ., CCGCR Report No. 97-1 Montréal, Québec, 52 pp. Bowman, B. R., W. K. Tobiska, and M. J. Kendra (2008a), The thermospheric semiannual density response to solar EUV heating, J. Atmos. Sol. Terr. Phys.,70(11–12), 1482–1496, doi:10.1016/j.jastp.2008.04.020. Bowman, B. R., W. K. Tobiska, F. A. Marcos, and C. Valladares (2008b), The JB2006 empirical thermospheric density model, J. Atmos. Sol. Terr. Phys.,70(5), 774–793, doi:10.1016/j.jastp.2007.10.002. Bruinsma, S., D. Tamagnan, and R. Biancale (2004), Atmospheric densities from CHAMP/STAR accelerometer observations, Planet. Space Sci., 52, 297–312, doi:10.1016/j.pss.2003.11.004. Bruinsma, S., J. M. Forbes, R. S. Nerem, and X. Zhang (2006), Thermosphere density response to the 20–21 November 2003 solar and geomagnetic storm from CHAMP and GRACE accelerometer data, J. Geophys. Res.,111, A06303, doi:10.1029/2005JA011284. Figure 14. Global averaged relative residuals (i.e., residual/density) for (a) ascending orbits and (b) descending orbits, for the parameterized PCA and for the MSIS empirical model. Journal of Geophysical Research: Space Physics 10.1002/2016JA022594 CALABIA AND JIN THERMOSPHERIC MASS DENSITY VARIATIONS 11,210 Acknowledgments This work was supported by the National Keystone Basic Research Program (MOST 973) (grant 2012CB72000), National Natural Science Foundation of China (NSFC) Project (grant 11573052), and Shanghai Science and Technology Commission Project (grant 12DZ2273300). Great appreciation is extended to ISDC for providing the data access. The authors declare that there is no conflict of interest regarding the publication of this paper. The GRACE data were obtained from the Information System and Data Center (ISDC) GeoForschungsZentrum (GFZ) website (http://isdc.gfz-potsdam.de/). Supporting data are included in a SI file; any additional data may be obtained from A. Calabia (email: andres@calabia. com). Thanks to the two reviewers for their constructive comments and helpful suggestions on an earlier version of the manuscript. 21699402, 2016, 11, Downloaded from https://agupubs.onlinelibrary.wiley.com/doi/10.1002/2016JA022594 by Spanish Cochrane National Provision (Ministerio de Sanidad), Wiley Online Library on [28/11/2025]. See the Terms and Conditions (https://onlinelibrary.wiley.com/terms-and-conditions) on Wiley Online Library for rules of use; OA articles are governed by the applicable Creative Commons License Calabia, A., and S. G. Jin (2016), Assessment of conservative force models from GRACE accelerometers and precise orbit determination, Aerosp. Sci. Technol.,49,80–87, doi:10.1016/j.ast.2015.11.034. Calabia, A., S. G. Jin, and R. Tenzer (2015), A new GPS-based calibration of GRACE accelerometers using the arc-to-chord threshold uncovered sinusoidal disturbing signal, Aerosp. Sci. Technol.,45, 265–271, doi:10.1016/j.ast.2015.05.013. Cane, H. V., I. G. Richardson, and T. T. von Rosenvinge (1998), Interplanetary magnetic field periodicity of ∼153 days, Geophys. Res. Lett., 25(24), 4437–4440, doi:10.1029/1998GL900208. Doornbos, E. (2011), Thermospheric density and wind determination from satellite dynamics, PhD thesis. TU Delft. Drob, D. P., et al. (2008), An empirical model of the Earth’s horizontal wind fields: HWM07, J. Geophys. Res.,113, A12304, doi:10.1029/ 2008JA013668. Emmert, J. T., and J. M. Picone (2010), Climatology of globally averaged thermospheric mass density, J. Geophys. Res.,115, A09326, doi:10.1029/2010JA015298. Forbes, J. M., X. Zhang, S. Bruinsma, and J. Oberheide (2013), Lunar semidiurnal tide in the thermosphere under solar minimum conditions, J. Geophys. Res. Space Physics,118, 1788–1801, doi:10.1029/2012JA017962. Guo, J., W. Wan, J. M. Forbes, E. Sutton, R. S. Nerem, T. N. Woods, S. Bruinsma, and L. Liu (2007), Effects of solar variability on thermosphere density from CHAMP accelerometer data, J. Geophys. Res.,112, A10308, doi:10.1029/2007JA012409. Guo, J., W. Wan, J. M. Forbes, E. Sutton, R. S. Nerem, and S. Bruinsma (2008), Interannual and latitudinal variability of the thermosphere density annual harmonics, J. Geophys. Res.,113, A08301, doi:10.1029/2008JA013056. Jin, S., J.-U. Park, J.-H. Cho, and P.-H. Park (2007), Seasonal variability of GPS-derived zenith tropospheric delay (1994–2006) and climate implications, J. Geophys. Res.,112, D09110, doi:10.1029/2006JD007772. Jin, S. G., O. F. Luo, and P. H. Park (2008), GPS observations of the ionospheric F 2 -layer behavior during the 20th November 2003 geomagnetic storm over South Korea, J. Geod.,82(12), 883–892, doi:10.1007/s00190-008-0217-x. Jin, S. G., T. van Dam, and S. Wdowinski (2013), Observing and understanding the Earth system variations from space geodesy, J. Geodyn.,72, 1–10, doi:10.1016/j.jog.2013.08.001. Joshi, B., and A. Joshi (2005), Intermediate-term periodicities in soft X-ray flare index during solar cycles 21, 22 and 23, Sol. Phys.,226, 153, doi:10.1007/s11207-005-5716-9. Knocke, P., and J. Ries (1987), Earth radiation pressure effects on satellites, Center for Space Research, The Univ. of Texas at Austin, Technical Memorandum, CSR-TM-87-01. Lathuillère, C., M. Menvielle, A. Marchaudon, and S. Bruinsma (2008), A statistical study of the observed and modeled global thermosphere response to magnetic activity at middle and low latitudes, J. Geophys. Res.,113, A07311, doi:10.1029/2007JA012991. Lei, J., T. Matsuo, X. Dou, E. Sutton, and X. Luan (2012a), Annual and semiannual variations of thermospheric density: EOF analysis of CHAMP and GRACE data, J. Geophys. Res.,117, A01310, doi:10.1029/2011JA017324. Lei, J., J. P. Thayer, W. Wang, X. Luan, X. Dou, and R. Roble (2012b), Simulations of the equatorial thermosphere anomaly: Physical mechanisms for crest formation, J. Geophys. Res.,117, A06318, doi:10.1029/2012JA017613. Liu, H., H. Lühr, V. Henize, and W. Köhler (2005), Global distribution of the thermospheric total mass density derived from CHAMP, J. Geophys. Res.,110, A04301, doi:10.1029/2004JA010741. Liu, H., H. Lühr, and S. Watanabe (2007), Climatology of the equatorial thermospheric mass density anomaly, J. Geophys. Res.,112, A05305, doi:10.1029/2006JA012199. Liu, H., M. Yamamoto, and H. Lühr (2009), Wave-4 pattern of the equatorial mass density anomaly: A thermospheric signature of tropical deep convection, Geophys. Res. Lett.,36, L18104, doi:10.1029/2009GL039865. Luthcke, S. B., J. A. Marshall, S. C. Rowton, K. E. Rachlin, C. M. Cox, and R. G. Williamson (1997), Enhanced radiative force modeling of the tracking and data relay satellites, J. Astronaut. Sci.,45(3), 349–370. Ma, R., J. Xu, W. Wang, J. Lei, H.-L. Liu, A. Maute, and M. E. Hagan (2010), Variations of the nighttime thermospheric mass density at low and middle latitudes, J. Geophys. Res.,115, A12301, doi:10.1029/2010JA015784. Matsuo, T., and J. M. Forbes (2010), Principal modes of thermospheric density variability: Empirical orthogonal function analysis of CHAMP 2001–2008 data, J. Geophys. Res.,115, A07309, doi:10.1029/2009JA015109. Mehta, P. M., C. A. McLaughlin, and E. K. Sutton (2013), Drag coefficient modeling for grace using direct simulation Monte Carlo, Adv. Space Res.,52, 2035–2051, doi:10.1016/j.asr.2013.08.033. Montenbruck, O., and E. Gill (2013), Satellite Orbits: Models, Methods and Applications, Springer, Berlin. Müller, S., H. Lühr, and S. Rentz (2009), Solar and geomagnetic forcing of the low latitude thermospheric mass density as observed by CHAMP, Ann. Geophys.,27, 2087–2099, doi:10.5194/angeo-27-2087-2009. Munk, W. H., and D. E. Cartwright (1966), Tidal spectroscopy and prediction, Phil. Trans. R. Soc. Lond.,A259, 533–581. Petit, G., and B. Luzum (2010), IERS Conventions (2010), IERS Tech. Note, vol 36. Verlag des Bundesamts für Kartogr. und Geod., Frankfurt am. Picone, J. M., A. E. Hedin, D. P. Drob, and A. C. Aikin, (2002), NRLMSISE-00 empirical model of the atmosphere: Statistical comparisons and scientific issues, J. Geophys. Res.,107(A12), 1468, doi:10.1029/2002JA009430. Preisendorfer, R. (1988), Principal Component Analysis in Meteorology and Oceanography, Elsevier, Amsterdam. Rentz, S., and H. Lühr (2008), Climatology of the cusp-related thermospheric mass density anomaly, as derived from CHAMP observations, Ann. Geophys.,26, 2807–2823, doi:10.5194/angeo-26-2807-2008. Ruan, H., J. Lei, X. Dou, W. Wan, and Y. C.-M. Liu (2014), Midnight density maximum in the thermosphere from the CHAMP observations, J. Geophys. Res. Space Physics,119, 3741–3746, doi:10.1002/2013JA019566. Ruan, H., J. Du, M. Cook, W. Wang, J. Yue, Q. Gan, X. Dou, and J. Lei (2015), A numerical study of the effects of migrating tides on thermosphere midnight density maximum, J. Geophys. Res. Space Physics,120, 6766–6778, doi:10.1002/2015JA021190. Sentman, L. H. (1961), Comparison of the exact and approximate methods for predicting free-molecular aerodynamic coefficients, Am. Rocket Soc. J.,31, 1576–1579. Sutton, E. K. (2011), Accelerometer-derived atmospheric densities from the CHAMP and GRACE accelerometers: Version 2.3, Tech.Memo, Air Force Res. Lab., Kirtland Air Force Base. Wilks, D. S. (1995), Statistical Methods in the Atmospheric Sciences, Academic Press, San Diego, Calif. Xu, J., W. Wang, and H. Gao (2013), The longitudinal variation of the daily mean thermospheric mass density, J. Geophys. Res. Space Physics, 118, 515–523, doi:10.1029/2012JA017918. Journal of Geophysical Research: Space Physics 10.1002/2016JA022594 CALABIA AND JIN THERMOSPHERIC MASS DENSITY VARIATIONS 11,211 21699402, 2016, 11, Downloaded from https://agupubs.onlinelibrary.wiley.com/doi/10.1002/2016JA022594 by Spanish Cochrane National Provision (Ministerio de Sanidad), Wiley Online Library on [28/11/2025]. See the Terms and Conditions (https://onlinelibrary.wiley.com/terms-and-conditions) on Wiley Online Library for rules of use; OA articles are governed by the applicable Creative Commons License Erratum In the originally published version of this article, the legends of Figure 7 through 13 were mistakenly transposed. They are now represented accurately. In addition, a small correction to an entry in Table 4 was made and a “#”was removed from the axis in Figure 13. The present version may be considered the authoritative version of record. Journal of Geophysical Research: Space Physics 10.1002/2016JA022594 CALABIA AND JIN THERMOSPHERIC MASS DENSITY VARIATIONS 11,212 21699402, 2016, 11, Downloaded from https://agupubs.onlinelibrary.wiley.com/doi/10.1002/2016JA022594 by Spanish Cochrane National Provision (Ministerio de Sanidad), Wiley Online Library on [28/11/2025]. See the Terms and Conditions (https://onlinelibrary.wiley.com/terms-and-conditions) on Wiley Online Library for rules of use; OA articles are governed by the applicable Creative Commons License