Eigen-Swarm: Swarm's Thermospheric Mass Density Modeling via Eigen-Decomposition

Eigen‐Swarm: Swarm's Thermospheric Mass Density Modeling via Eigen‐Decomposition Charles Owolabi 1 , Hyunju Connor 2 , Don Hampton 1 , Denny M. Oliveira 2,3 , Andres Calabia 4 , V. Sai Gowtam 2,5 , and Eftyhia Zesta 2 1 Geophysical Institute, University of Alaska Fairbanks, Fairbanks, AK, USA, 2 NASA Goddard Space Flight Center, Greenbelt, MD, USA, 3 Goddard Planetary Heliophysics Institute, University of Maryland, Baltimore, MD, USA, 4 Department of Physics and Mathematics, University of Alcala, Madrid, Spain, 5 Department of Physics, The Catholic University of America, Washington, D.C., USA Abstract Precise thermospheric mass density (TMD) prediction is essential for satellite orbital tracking, reentry calculations, and upper atmospheric processes under varying solar and magnetospheric conditions. In this paper, we construct an empirical model of TMD at 450 km altitude with accelerometer‐derived (ACC) TMD observations from the Swarm‐C satellite during 2014–2020. We employ the Eigen‐Decomposition technique to extract dominant spatio‐temporal modes, with the first three capturing 99.12% of the variance, forming the basis of the Swarm‐based Eigen‐Decomposition model. We study the factors controlling the observed TMD variability and investigate their relation to longitude, latitude, local solar time, seasonal effects, solar and geomagnetic indices. The Eigen‐Decomposition model performance is validated by comparison with the Jacchia‐Bowman 2008 (JB2008), Naval Research Laboratory Mass Spectrometer Incoherent Scatter 2.0 (NRLMSIS2.0), and Calabia and Jin (CAJIN) models, as well as TMD data from the Gravity Recovery and Climate Experiment Follow‐On mission during 2018–2020, using root mean square error (RMSE) as the evaluation metric. The Eigen‐Decomposition model achieves an RMSE of 19.45%, outperforming JB2008 (29.83%), NRLMSIS2.0 (65.16%), and CAJIN (45.25%). Additional metrics, including correlation coefficient (R), mean (μ), and variance (σ 2 ), further confirm the improved accuracy and fidelity of our approach across different solar activity conditions. This work demonstrates the effectiveness of data‐driven techniques in capturing TMD dynamics and deepening our understanding of the thermospheric response to space weather conditions. Plain Language Summary Accurate prediction of thermospheric mass density (TMD) is vital for safeguarding satellite operations from space weather impacts and also provides critical information about satellite drags. Recent events, such as the loss of 38 satellites by SpaceX due to a geomagnetic storm, have highlighted the importance of having more accurate estimation and prediction of TMD variability. While several empirical models are used to predict TMD changes, they often struggle to accurately capture the TMD variability, especially during the geomagnetic storms, highlighting the need for continuous validation and improvement. Using the Eigen‐Decomposition and regression techniques, the dominant modes of TMD variability from Swarm‐C observations (2014–2020) were extracted and parameterized with respect to external drivers to develop the Eigen‐Decomposition model. The newly developed model effectively reproduces TMD variability when validated against traditional empirical models and independent TMD observations. Our findings advance our understanding of thermospheric dynamics and enhance predictions crucial for satellite operations and space weather monitoring. 1. Introduction As the number of anthropogenic space objects in the low‐Earth orbit (LEO) grows, accurate modeling of the global thermospheric variability for orbit prediction becomes more important, especially for space weather forecasting and monitoring. Understanding thermospheric mass density (TMD) variability is crucial for comprehending the Sun‐Earth interaction and its detrimental impacts on LEO satellites. The primary source of uncertainty in predicting satellite orbits in the Earth's upper atmosphere is TMD variability, which is characterized by various disturbances that originate from multiple forcing mechanisms (e.g., Calabia et al., 2023; Krauss et al., 2015; Oliveira & Zesta, 2019; Yu et al., 2020). Such mechanisms include the diurnal and seasonal cycles from the variable solar irradiance, the electromagnetic energy deposited in the high‐latitude atmosphere due to RESEARCH ARTICLE 10.1029/2025SW004351 Special Collection: Recent Advances on Modelling and Observations in Space and Earth Sciences (0) Key Points: •An Eigen‐Decomposition model of the thermospheric mass density is constructed based on the Swarm‐C satellite observations •The model characterizes density variability as a function of longitude, latitude, local time, season, and solar and geomagnetic activity •Eigen‐Decomposition model performs better than the JB2008 model and its capture efficiency is higher than that of the NRLMSIS2.0 model Correspondence to: C. Owolabi, [email protected] Citation: Owolabi, C., Connor, H., Hampton, D., Oliveira, D. M., Calabia, A., Gowtam, V. S., & Zesta, E. (2025). Eigen‐Swarm: Swarm's thermospheric mass density modeling via Eigen‐Decomposition. Space Weather,23, e2025SW004351. https://doi. org/10.1029/2025SW004351 Received 14 JAN 2025 Accepted 18 JUL 2025 Author Contributions: Conceptualization: Charles Owolabi Formal analysis: Charles Owolabi Funding acquisition: Hyunju Connor, Don Hampton Investigation: Charles Owolabi, Denny M. Oliveira, Andres Calabia Methodology: Charles Owolabi Resources: Don Hampton Supervision: Hyunju Connor Validation: Hyunju Connor, Denny M. Oliveira, Andres Calabia, V. Sai Gowtam, Eftyhia Zesta Visualization: Hyunju Connor, V. Sai Gowtam © 2025. The Author(s). This is an open access article under the terms of the Creative Commons Attribution License, which permits use, distribution and reproduction in any medium, provided the original work is properly cited. OWOLABI ET AL. 1 of 14 solar wind‐magnetosphere coupling, and effects from the lower and middle atmosphere (e.g., Bowman et al., 2008; Burns et al., 2021; Calabia & Jin, 2016b; Oliveira et al., 2017). At high latitudes, the TMD variability is primarily influenced by Joule heating and energetic electron precipitation, while also being affected by other factors such as winds, plasma turbulence, and variable electric fields and conductivity (e.g., Calabia & Jin, 2019; Oliveira & Zesta, 2024; Sutton et al., 2012; Weimer, 2005). On the other hand, in the lower latitude regions, the TMD variability is more actively driven by the solar extreme ultraviolet (EUV) irradiance and the soft X‐ray flux (e.g., Sutton et al., 2006; Vourlidas & Bruinsma, 2018). These features make the characterization and modeling of TMD dynamics, a key parameter for orbit prediction, difficult, while our modern society strongly depends on space assets (e.g., Acciarini et al., 2024; Oliveira et al., 2021). A new motivator for accurate TMD forecasts is satellite collision avoidance. SpaceX alone performed >100,000 collision avoidance maneuvers in 2024, and all conjunction assessment and collision avoidance capability in LEO is severely degraded during geomagnetic storms because we couldn't accurately predict satellite trajectories due to variable drag (especially during the Gannon storm). The important role of TMD in precise orbit determination and atmospheric research has forced researchers from all over the world to study the spatio‐temporal characteristics of TMD with empirical and physics‐based modeling (e.g., Acciarini et al., 2024; Calabia et al., 2020; Calabia & Jin, 2016b; Gondelach & Linares, 2021; Lei et al., 2012; Liu et al., 2013; Matsuo & Forbes, 2010; Oliveira et al., 2021; Yu et al., 2020). Accurate monitoring and prediction of TMD variability is a challenging task, involving many variables and a multitude of available techniques, and several atmospheric models are currently in use for the dynamics of the thermosphere‐ionosphere system. Some well‐known empirical models developed over the years are the Naval Research Laboratory Mass Spectrometer Incoherent Scatter 2.0 (NRLMSIS2.0) (e.g., Picone et al., 2002), Jacchia‐ Bowman 2008 (JB2008) (e.g., Bowman et al., 2007), High Accuracy Satellite Drag Model (HASDM) (e.g., Tobiska et al., 2021), and Drag‐Temperature Model DTM‐2013 (e.g., Bruinsma, 2015), which are currently in use for practical applications such as the implementation of drag in precise orbit determination. For several years, these models have been improved with respect to earlier versions (e.g., Emmert et al., 2020), since large uncertainties usually occur, especially during geomagnetic storms (e.g., Calabia & Jin, 2016a; Fang et al., 2022; He et al., 2023; Marcos, 1990; Oliveira et al., 2021). Due to the importance of TMD in spacecraft operation, especially for LEO satellites, researchers have explored its spatial and temporal properties of TMD and developed empirical models using methods such as machine learning (e.g., Acciarini et al., 2024; Weng et al., 2020; Zhang et al., 2021), dynamic reduced‐order modeling (e.g., Mehta et al., 2018), data assimilation (e.g., Matsuo et al., 2012) and the Thermosphere–Ionosphere Electrodynamic General Circulation Model, yielding significant results (e.g., Elvidge et al., 2016). Over the past years, substantial research has improved TMD variability estimation during either quiet or stormy time conditions (e.g., Calabia & Jin, 2016b; Lei et al., 2012; Liu et al., 2013), since even moderate geomagnetic storms can significantly impact satellite orbital drag (e.g., Baruah et al., 2024; Lin et al., 2022). The ACC‐derived TMDs from CHAMP (CHAllenging Minisatellite Payload), Swarm and twin‐GRACE (Gravity Recovery and Climate Experiment) showed that the accuracy of existing empirical models is limited, especially during space weather events (e.g., Liu et al., 2005; Oliveira et al., 2021; Sutton et al., 2007). In an effort to account for new TMD data in the models, different modeling approaches have been used, including, for example, the use of the Eigen‐Decomposition analysis, usually known as Principal Component Analysis in geophysics, or Empirical Orthogonal Function analysis in climate studies (e.g., Calabia & Jin, 2016b; Gondelach & Linares, 2021; Lei et al., 2012; Matsuo & Forbes, 2010; Owolabi et al., 2021,2022). Here we will name the method by its original form, namely the Eigen‐Decomposition method. For instance, Lei et al. (2012) investigated annual and semiannual variations using a sequential nonlinear regression analysis of TMD data from CHAMP and GRACE at 400 km altitude during 2002–2010. The authors highlighted the capabilities of their method in capturing significant temporal patterns and spatial structures in TMD data, offering a refined understanding of its behavior under different conditions. The authors discovered a correlation between TMD and NRLMSIS‐00 changes with geographic location, season, and solar‐geomagnetic activity. In Sutton et al. (2012), the authors developed a series of dominant functions to better capture the spatial variability and response of the thermosphere to improve the HASDM model accuracy in TMD prediction. They improved HASDM calibration with sparse data, reducing the root mean square error (RMSE) by 32.9% relative to synthetic data. A study by Liu et al. (2013) employed in situ measurements of TMD from CHAMP satellite and constructed a model that empirically investigates the solar impact of TMD. Their findings indicate a significant correlation between the equinoctial asymmetry of TMD and the solar cycle, wherein the magnitude and phase of the asymmetry are heavily influenced by different levels of solar and geomagnetic activity. Writing – review & editing: Hyunju Connor, Denny M. Oliveira, Andres Calabia, Eftyhia Zesta Space Weather 10.1029/2025SW004351 OWOLABI ET AL. 2 of 14 15427390, 2025, 7, Downloaded from https://agupubs.onlinelibrary.wiley.com/doi/10.1029/2025SW004351 by Spanish Cochrane National Provision (Ministerio de Sanidad), Wiley Online Library on [06/12/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 Using a grid‐based technique, Calabia and Jin (2016b) applied Eigen‐Decomposition to dimensionally reduce 13 years of ACC‐derived TMD from the GRACE mission. They analyzed the resulting principal modes and parameterized them according to solar and magnetospheric forcing, local solar time (LST), and annual variations to create an empirical model available in MATLAB language. This model was later upgraded with a parameterized magnetospheric forcing in Calabia and Jin (2019), hereafter named the CAJIN model. Their parameterizations were suitable to represent small‐scale variations, including, for example, the equatorial mass density anomaly (EMA) and the midnight density maximum, and the residuals showed additional contributions at the frequencies of the radiational tides and at the periods of 83, 93, 152, and 431 days. Then, in Calabia et al. (2020), the authors used Swam‐C data to make several orbit predictions and analyses. They found that the NRLMSISE‐00 model overestimates the in situ TMD by about 20%, while the CAJIN model underestimates it by the same amount. Yang et al. (2022) developed an empirical model to predict TMD variability from 2002 to 2010 using CHAMP and GRACE observations, which showed better agreement with independent GRACE‐A observations than NRLMSIS‐00, accurately capturing EMA and seasonal variations. Similarly, Xiong et al. (2018) introduced CH‐Therm‐2018, based on nine years of CHAMP data (2000–2009), which outperformed NRLMSIS‐00 in capturing TMD variability, particularly during the 2008–2009 solar minimum. Despite considerable advances in TMD model formulation over the years, there is a gap for model improvement in terms of performance and accuracy. This improvement is necessary to surpass the level of agreement achieved by existing models in practical applications (e.g., Calabia et al., 2020; Marcos, 1990). The TMD discrepancy between models and data suggests the need for a better understanding of the Sun–Earth system. Thus, researchers need to keep track of geomagnetic and solar components that affect TMD variability to improve present models and benchmarks. The previous solar cycle (2009–2020) is uniquely quiet compared to other solar cycles and thus the Eigen‐Decomposition model may give insight into thermosphere responses during our recent quiet solar cycle. Here, in this paper, we present a new TMD model based on the Eigen‐Decomposition of ACC‐derived data from 2014 to 2020. The Eigen‐Decomposition model differs from previous work by (a) focusing on a quiet solar cycle, which provides an opportunity to isolate and characterize thermospheric variability under weak solar and geomagnetic forcing conditions; (b) parameterizing TMD by latitude, longitude, LST, and day of year (DOY) to resolve fine‐scale structures; and (c) benchmarking against JB2008, NRLMSIS2.0, and CAJIN under various solar and geomagnetic conditions. Unlike the JB2008 and NRLMSIS2.0 models built on fixed empirical relationships derived from historical data, our adaptive statistical approach can evolve with new data and uncover latent structures. In principle, this comparative approach highlights the strengths and limitations of existing models under quiet conditions. The next section details the data and methods used in this study, followed by the Eigen‐ Decomposition model results and discussion in Section 4. The paper concludes with a summary in the final section. 2. Data and Method 2.1. Data and Empirical Atmospheric Models In this study, we utilized the ACC‐derived TMD from the Swarm‐C satellite to construct an Eigen‐Decomposition model. The ACC‐derived TMD offers higher temporal resolution and captures finer‐scale variability. The Swarm mission comprises constellation of three identical satellites named Alpha (A), Bravo (B), and Charlie (C), launched into near‐polar LEO on 22 November 2013 (e.g., Siemes et al., 2016). These satellites primarily study geomagnetic fields and are equipped with accelerometers positioned at their centers of mass to collect TMD data crucial for understanding Earth's thermosphere dynamics (e.g., Calabia & Jin, 2016b; Kodikara et al., 2018). Swarm‐A and ‐C orbit at an altitude of ∼462 km in near‐polar orbits (inclination 87.35°). The two satellites are separated by approximately 1.4° in longitude at the equator due to their slightly offset orbital planes. Swarm‐B orbits at a higher altitude of ∼510–530 km with an inclination of 87.75°, resulting in a longer orbital period of ∼94.5 min compared to ∼93 min for Swarm‐A and ‐C satellites. Different local time precession rates apply to Swarm‐A, ‐B, and ‐C satellites as follows: The 2 lower satellites precess through 12 hr of local time in 133 days, while the upper spacecraft precesses through 12 hr in 144 days, resulting in a local time separation of 10 hr after 5 years of mission duration. In October 2021, the orbital planes realigned, with Swarm B orbiting in the opposite direction of Swarms A and C satellites. We also used TMD data from the Gravity Recovery and Climate Experiment Follow‐On (GRACE‐FO) mission, a joint NASA and the German Research Centre for Geosciences (GFZ) initiative launched on 22 May 2018, as GRACE's successor. It consists of two identical satellites flying ∼220 km apart in a polar orbit at ∼500 km altitude (e.g., Bandikovaa et al., 2019; Ciracì et al., 2020; Landerer Space Weather 10.1029/2025SW004351 OWOLABI ET AL. 3 of 14 15427390, 2025, 7, Downloaded from https://agupubs.onlinelibrary.wiley.com/doi/10.1029/2025SW004351 by Spanish Cochrane National Provision (Ministerio de Sanidad), Wiley Online Library on [06/12/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 et al., 2020). Global TMD data, available at higher spatial and temporal resolutions, are freely accessible and widely used for TMD estimation (https://swarm‐diss.eo.esa.int). To validate the Eigen‐Decomposition model outputs, we estimated TMD values from JB2008, NRLMSIS2.0 and CAJIN models along the Swarm‐C satellite's orbit (2014–2020). The NRLMSIS2.0 and JB2008 models are empirical and semi‐empirical frameworks that rely on long‐term historical data sets and established physical relationships. NRLMSIS2.0 for instance, is an updated thermospheric model (0–1,000 km) based on satellite and ground‐based radar observations. It characterizes the upper atmosphere, providing altitude‐dependent exospheric and neutral temperatures and TMD. It takes inputs such as altitude, latitude, longitude, DOY, LST, daily and 81‐ day averaged F10.7 flux, and Ap index. On the other hand, the JB2008 is tailored for satellite drag prediction and incorporates solar and geomagnetic indices to model TMD. It includes empirical corrections for storm‐time effects and diurnal variations, making it robust for operational use in space weather forecasting. Compared to NRLMSIS2.0, it incorporates additional solar proxies S10.7 (26–34 nm bandpass EUV chromospheric irradiance index), M10.7 (160 nm FUV Schumann‐Runge photospheric irradiance proxy), Y10.7 (121.6 nm Lyman‐alpha chromospheric region irradiance convolved with the 0.1–0.8 nm X‐ray coronal irradiance), their 81‐day averages, and geomagnetic indices (Ap, Dst) (e.g., Tobiska et al., 2008). Using the NRLMSIS2.0 model, Swarm‐C observations were normalized to a height of 450 km in order to exclude TMD variability caused by changes in satellite orbital altitude, as follows: ρ(450km) = ρ(z)∗ρm(450km) ρm(z). In this equation, ρ m is the NRLMSIS2.0 estimate, and zis Swarm‐C satellite height. The normalization error should be low. Liu et al. (2007) found 5% uncertainty acceptable in multi‐year comparisons, especially when the spacecraft's altitude shift is negligible (e.g., Bruinsma et al., 2023). 2.2. Space Weather Indices and Formulation of the TMD Model To develop a reliable empirical model of TMD, it is essential to incorporate key space weather drivers that influence upper atmospheric variability. In this study, we utilize the Mg II core‐to‐wing index from the Solar Radiation and Climate Experiment (SORCE) as a sensitive proxy for solar activity, along with the planetary geomagnetic Ap index to account for geomagnetic disturbances. The MgII core‐to‐wing index, scaled as M10.7 =8,083 ×MgII–1,154 provides a more reliable proxy for solar activity than the F10.7 flux (e.g., Ruan et al., 2018; Viereck et al., 2001) and has been shown to effectively represent solar activity influences on TMD variability (e.g., Solomon et al., 2011). The Eigen‐Decomposition model is constructed using M10.7p, M10.7A, Ap and ApD indices. Where M10.7prefers to the daily M10.7 index for the day before, M10.7A represents the 81‐ day average value centered on the current day, and ApD stands for the daily averaged 3‐hr Ap index. Figure 1c shows the daily and 81‐day averaged solar M10.7 index, which ranges from 243 to 60 sfu for different solar activity levels. The 3‐hr Ap index obtained from the OMNI webpage, shown in Figure 1d, also varies during this investigation and is higher during high solar activity. We apply the Eigen‐Decomposition approach to decompose the time‐varying TMD(lon,lat,t)into a set of mutually orthogonal spatial base functions ϕn(lon,lat)and time series αn(t)as described in numerous studies (e.g., Calabia & Jin, 2016b; Lei et al., 2012; Ruan et al., 2018; Yu et al., 2020). In this study, we extracted the first three dominant spatio‐temporal modes, which collectively explained 99.12% of the total variance. Mathematically, this can be expressed as: TMD(lon,lat,t)=∑ n=M n=1 αn(t)∗ϕn(lon,lat) (n=1,2,…M)(1) where TMD(lon,lat,t)is the observed TMD for certain longitude, latitude, and time, αn(t)is the component's amplitude that varies with respect to time, t, and ϕn(lon,lat)is the nth Eigen mode that varies as a function of geographic longitude and latitude. Our approach involves grid‐based Eigen‐Decomposition followed by the computation of the covariance matrix across the grid, akin to the method used by Calabia and Jin (2016b) and Ruan et al. (2018). The spatial vector TMD(lon, lat)is represented by M ×N elements on a regular grid of m=1, 2, … M, and n=1, 2, …N is the modal index, which is less than or equal to M. Our method involves using a two‐ dimensional matrix where each row represents LST and each column represents a time series for a specific latitude. Each daily global grid, spanning from 2014 to 2020, has been transformed into a one‐dimensional array (72 ×36 =2,592). These arrays were concatenated from ascending and descending satellite orbits to form a consolidated data set, resulting in a final matrix size of 2,592 ×9,524. By analyzing the covariance structure of Space Weather 10.1029/2025SW004351 OWOLABI ET AL. 4 of 14 15427390, 2025, 7, Downloaded from https://agupubs.onlinelibrary.wiley.com/doi/10.1029/2025SW004351 by Spanish Cochrane National Provision (Ministerio de Sanidad), Wiley Online Library on [06/12/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 TMD, our method identifies dominant modes of variability, essentially capturing how different regions or parameters co‐vary over time. This technique allows for dimensionality reduction and pattern recognition without relying on predefined physical laws. It is particularly powerful for uncovering regional or temporal anomalies and can adapt flexibly to new data sets, making it well‐suited for localized or evolving atmospheric phenomena. We parameterize αn(t)for each mode spanning from 2014 to 2020, incorporating periodic oscillations such as annual, semi‐annual, and terannual variations using harmonic functions. The LST variation includes 24, 12, and 8h components. The change in the harmonic amplitudes is related to the solar flux and geomagnetic indices, as indicated in Equation 2, accounting for solar cycle dependence and geomagnetic activity. The unknown coefficients a,b,c,d,p, and qare derived through linear regression by minimizing the sum of squared differences between the observed values of α and the values predicted by the regression equation. The calculated coefficients are then integrated into the regression equation, together with the input variables (LON, LAT, DOY, LST, M10.7A, M10.7p, Ap, ApD, and first three modes) to reconstruct and provide accurate predictions of TMD values. αn=Xs∗Ys Xs=(ao+a1∗M10.7p+a2∗M10.7A+a3∗M10.7p∗M10.7A) (bo+b1∗Ap +b2∗ApD) Ys=c0+∑ 3 k=1(pkcos (2π∗DOY ∗k 365.25 )+qksin (2π∗DOY ∗k 365.25 )) +∑ 3 k=1(ckcos (2π∗LST ∗k 24 )+dksin (2π∗LST ∗k 24 )) ⎫⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎬ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎭ (2) Figure 1. (a–c) Leading orthogonal components of the thermospheric mass density (PC1–PC3) at 450 km altitude derived from Swarm‐C observations (2014–2020) in unit of 10 −2 as a function of longitude and latitudes and (d) their relative contributions to the total variance. Space Weather 10.1029/2025SW004351 OWOLABI ET AL. 5 of 14 15427390, 2025, 7, Downloaded from https://agupubs.onlinelibrary.wiley.com/doi/10.1029/2025SW004351 by Spanish Cochrane National Provision (Ministerio de Sanidad), Wiley Online Library on [06/12/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 where Xsdepends on solar and geomagnetic activity effects, Ysis referred to as UT/longitude and the seasonality parameter which signifies the inter‐annual variations, c0is the mean value of αn. Here pkand qk(k=1, 2, 3) represent the annual, semiannual and terannual variations of the TMD variations, ckand dk(k=1, 2, 3) indicate the diurnal, semidiurnal and terdiurnal variations, LST stands for current LST in hours, including the diurnal, semidiurnal, and terdiurnal variations; and DOY represents the number of days from January 2014 to December 2020, encompassing the annual, semiannual, and terannual variations. Leap years were accounted for by adding a 0.25‐factor to the 365‐day calculation. To assess performance, we compared TMD predictions from the Eigen‐ Decomposition model against JB2008, MSIS2.0, and CAJIN models using key metrics that quantify how well each model captures TMD variability. These include RMSE, correlation coefficient (R), mean (μ), and variance (σ2). Here are the equations for these metrics: RMSE = 1 h∑ h i=1(TMDobs i−TMDmod i)2 √ √ √ R=Cov(TMDmod i,TMDobs i)  Var(TMDmod i)∗Var(TMDobs i) √ μ=1 h∑ n i TMDobs i σ2=1 h∑ n i(TMDobs i−μ)2 ⎫⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎬ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎭ (3) where his the number of observations, TMDmod iand TMDobs iare the ith actual and predicted TMD values, μ is the mean of TMD values, Cov and Var denote covariance and variance, respectively. We used TMD data from the GRACE‐FO mission spanning 1 May 2018 to 31 December 2020 for model validation. 3. Model Results and Discussion Figures 1a–1c shows the global distribution of the first three dominant modes (ϕ) that determine spatial TMD variability in geographic coordinates, while Figure 1d shows the relative contribution of the top 5 spatial modes. The first orthogonal mode (ϕ1) retains 96.82%, the second orthogonal mode (ϕ2) retains 1.29%, and the third orthogonal mode (ϕ3) retains 1.01% of the total variability. In total, the first 3 modes account for 99.12% of the total variance. It is evident that each mode represents the main spatial features of TMD variability. For example, the equatorial anomaly feature is pictured in the first mode (Figure 1a), with the TMD peak in the southern polar region and other peaks appearing in the equatorial region around 0° longitude. The maximum values are situated at the southern cusp, while the minimum values are found at the northern cusp. It can be inferred that the distribution is associated with the orientation of the dip angle, resulting in positive TMD contributions in the southern cusp and negative contributions in the northern cusp (e.g., Lei et al., 2012). Since the first component contributes 96.82% of overall variance, its global average mode is significant. The second and third modes show north/south dominance of TMD variability. They exhibit different patterns in the northern and southern hemispheres, suggesting that it could be a seasonal variation in the thermosphere. As previously noted in Liu et al. (2009), the second and third modes clearly relate to changes in latitudinal variations caused by LST as well as latitudinal variations following the subsolar point. Figures 2a, 2c, and 2e display the temporal variations of the first 3 modes (α1−3) for the ascending orbit, whereas Figures 2b, 2d, and 2f mirror the same variations for the descending orbit. It is evident that the (α1−3) contains clear annual and semiannual variations that are much more obvious in the higher‐order modes. Notably, their (α1−3) amplitudes gradually increase from 2014 to 2020 and are solar cycle dependent, which correlates with the first three temporal coefficients parameterized in terms of longitude, latitude, LST, DOY, M10.7p, M10.7A, Ap, and ApD with correlation values of 99%, 96%, and 93%, respectively. As seen by the red solid lines and the correlation values, the regression analysis results are well parameterized and aid the TMD reconstruction. Figures 2g and 2h are an identical plot of the solar irradiance proxy M10.7 index and M10.7A. Figures 2i and 2j are also identical plots and Space Weather 10.1029/2025SW004351 OWOLABI ET AL. 6 of 14 15427390, 2025, 7, Downloaded from https://agupubs.onlinelibrary.wiley.com/doi/10.1029/2025SW004351 by Spanish Cochrane National Provision (Ministerio de Sanidad), Wiley Online Library on [06/12/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 show the daily average of the 3‐hourly Ap index during this period. As illustrated in this Figure, the M10.7 flux ranged from 80 solar flux units (sfu) in 2014 to about 140 sfu in 2015. The Ap index varied significantly during 2014–2020, and its variability is more profound during the period of high solar activity in 2014–2016 and becomes lower at low solar activity in 2018–2020. Figure 3shows the TMD obtained from Swarm‐C satellite (magenta dot), alongside statistical results from the JB2008 model (blue line), NRLMSIS2.0 model (gray line), and CAJIN model (black line). As shown in Figure 3a, the Eigen‐Decomposition model (red dot) captures overall TMD variability seen in the Swarm‐C observations by using the following input parameters: M10.7p, M10.7A, Ap, ApD, longitude, latitude, LST, and DOY. Furthermore, the TMDs estimated from the Eigen‐Decomposition model exhibit a temporal variation i.e. closer to the Swarm‐C values than those estimated by the MSIS2.0 and CAJIN models. As shown in Figure 3a, the absolute TMD estimates increase with solar activity because the thermosphere is heated by solar radiation as well as geomagnetic activity and expands to higher altitudes. Note that the Swarm‐C observations were not used when the JB2008 and MSIS2.0 models were developed (e.g., Bowman et al., 2007; Emmert et al., 2020). Figure 3b shows the model‐to‐data ratio of TMD obtained by the Eigen‐Decomposition model (red line), JB2008 (blue line), MSIS2.0 (gray line) and CAJIN (black line). This ratio is calculated along with standard deviation between the TMD from the models and Swarm‐C observations [(TMD data–Model value)/Model value)]. The Eigen‐Decomposition model shows a good model‐to‐TMD data ratio of 0.8–1.1, while the JB2008 and CAJIN models show the ratio of 1.2–1.5, and MSIS2.0 shows the ratio of 1.3–2.5, respectively. Comparatively, the seasonal‐averaged ratios between the Eigen‐Decomposition model and Figure 2. Temporal functions representing the PC 1–3 amplitudes derived from Swarm‐C observations fitting (blue dots) and the parameterization results (red lines) for ascending (a, c, e) and descending (b, d, f) orbits during 2014–2020. (g and h) Time series of the daily M10.7 (or M10.7eff =(M10.7 +M10.7A)/2) index in gray and its 81‐running average (M10.7A) in red. (i and j) The daily average of geomagnetic activity index Ap (i.e., ApD) between 2014 and 2020. Space Weather 10.1029/2025SW004351 OWOLABI ET AL. 7 of 14 15427390, 2025, 7, Downloaded from https://agupubs.onlinelibrary.wiley.com/doi/10.1029/2025SW004351 by Spanish Cochrane National Provision (Ministerio de Sanidad), Wiley Online Library on [06/12/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 Swarm‐C observations are within ±20%, while MSIS2.0 overestimates from 5% for high (2014:2015) and moderate (2015:2016) solar activity to 40% during low solar (2017:2020) activity periods. Another thing worth noting is that the Eigen‐Decomposition model, JB2008 and MSIS2.0 models have a higher precision and applicability in the years 2014–2015 (high solar activity period) than in the year 2020 (low solar activity period), which are consistent with previous studies (e.g., Ruan et al., 2018). We observed that models approaching a quotient near one exhibit lower errors compared to observations. The MSIS2.0 shows poor performance for low‐ density values during 2018–2021, while JB2008 and CAJIN remain within a 0.5 deviation. Notably, CAJIN shows better agreement with Swarm‐C observations under certain conditions, especially during periods of moderate geomagnetic activity. This advantage can be attributed to CAJIN being a recent data‐driven model trained on a more extensive data set, including post‐2014 satellite observations such as Swarm and CHAMP observations. Figure 3c shows the statistical distribution of model‐Swarm‐C daily mean TMD relative errors. The RMSE of the TMD from the Eigen‐Decomposition model, JB2008, MSIS2.0, and CAJIN models are recorded to be 19.45%, 29.83%, 65.16%, and 45.25%, respectively. The RMSE value of the Eigen‐Decomposition model is the smallest, when compared to JB2008, MSIS2.0, and CAJIN (Figure 3c), implying that the Eigen‐ Decomposition model predicts TMD better. Figures 3d and 3e show latitudinal and longitudinal behavior of the TMD deviation with error bars from Swarm‐C data for the Eigen‐Decomposition model (red), MSIS2.0 (gray), JB2008 (blue line) and CAJIN (black line) models. In Figures 3d and 3e, the MSIS2.0 model exhibits a 60% deviation at equatorial and low latitudes and 52% at middle and high latitudes. As for the performance of MSIS2.0, its deviation may stem from limitations in its composition and temperature parameterizations, particularly during disturbed conditions. Specifically, the fixed relationships used for neutral species (e.g., O, N 2 ) and the response of exospheric temperature to geomagnetic activity may not fully capture the dynamic variability of the thermosphere during solar maxima or storm‐time periods (e.g., Calabia & Jin, 2016b; Ruan et al., 2018). The Eigen‐Decomposition model and CAJIN model have near‐zero relative deviation in the latitudinal variation, better as compared to MSIS2.0 and JB2008 models. For longitudinal variation, the TMD relative deviation was about 55% and almost steady along longitudes. The MSIS2.0 model exhibits a relative deviation in TMD that is approximately twice as large as that observed in the JB2008 model and CAJIN has an interesting longitudinal variation centered at zero, suggesting that MSIS2.0 cannot accurately reproduce the longitudinal and latitudinal structures of TMD variability. Overall, Figure 3indicates that the Eigen‐Decomposition model reproduces the retrieved TMD variability better than the JB2008, MSIS2.0 and CAJIN models, consistent with previous studies Figure 3. (a) Orbit‐averaged thermospheric mass density (TMD) derived from Swarm‐C observations (black) normalized to 450 km from 2014 to 2020, compared with the Eigen‐Decomposition model (red), Jacchia‐Bowman 2008 (JB2008) (blue), NRLMSIS2.0 (gray), and CAJIN (black) models. (b) Seasonal mean ratios and standard deviations between model estimates and Swarm‐C observations. (c) Statistical distribution of relative TMD errors for the four models. (d–e) Latitudinal and longitudinal variability of TMD, with error bars, comparing Eigen‐Decomposition, NRLMSIS2.0, JB2008, and CAJIN models. Root mean square error =root‐mean‐square error. Space Weather 10.1029/2025SW004351 OWOLABI ET AL. 8 of 14 15427390, 2025, 7, Downloaded from https://agupubs.onlinelibrary.wiley.com/doi/10.1029/2025SW004351 by Spanish Cochrane National Provision (Ministerio de Sanidad), Wiley Online Library on [06/12/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 that demonstrated the improved performance of Eigen‐Decomposition models in capturing real TMD behavior (e.g., Calabia & Jin, 2016b; Lei et al., 2012). It is important to note that Figure 3serves as a benchmark within the model's training domain. To assess generalizability, we highlight the independent evaluation presented in the next section, which utilizes data from a non‐Swarm satellite. Moreover, in order to further assess the performance of the Eigen‐Decomposition model, Figures 4a–4d shows a scatter plot of daily TMD between models and Swarm‐C observations, while Figures (4e–4h) shows a histogram of TMD deviation between models and Swarm‐C observations. Each histogram of the relative deviation [(TMD data–Model value)/Model value)] was fitted with a normal distribution to estimate the average TMD deviation and the variability across the data set. The upper right corners of Figures (4e–4h) show the mean deviation (μ) and variance (σ2). The TMD shows a generally linear distribution with a positive slope (Figures 4a–4d) and the residuals follow a normal distribution (Figures 4e–4h). Furthermore, the majority of deviations (errors) for all the models are concentrated within the range of ±2 kg/m 3 . In this Figure, the Eigen‐Decomposition model has the lowest RMSE (0.084 kg/m 3 ), suggesting the lowest average prediction error. It also has the highest correlation coefficient (0.990), showing it better captures data patterns and variability. The JB2008 model has a slightly higher RMSE (0.106 kg/m 3 ) compared to the Eigen‐Decomposition model, indicating marginally greater prediction error. While its correlation coefficient (0.988) remains high, it is slightly lower than that of the Eigen‐ Decomposition model, suggesting a modestly reduced ability to capture data variability. The MSIS2.0 model has a strong correlation with the Swarm‐C observations (correlation coefficients can reach 0.962), which indicates that MSIS2.0 model is capable of representing the overall distributions and variability in the original TMD observations. In addition, the RMSE (0.181 kg/m 3 ) shows some accuracy of the MSIS2.0 model. Among the four models, CAJIN has the highest RMSE (0.382 kg/m 3 ), indicating the greatest average error. Its correlation coefficient (0.774) is lower than those of the Eigen‐Decomposition model, JB2008, and MSIS2.0 models, indicating a modest ability to capture data variability. This is consistent with the findings of Calabia et al. (2020), which showed that CAJIN underestimates TMD values by approximately 20% compared to Swarm‐C data. Conversely, the MSIS2.0 model tends to overestimate TMD. The JB2008 model shows slightly better correlation and accuracy than MSIS2.0, consistent with our findings. The Eigen‐Decomposition model was further validated using Figure 4. Scatter plots of the daily thermospheric mass density variability between models and Swarm‐C observations at 450 km altitude; (a) Eigen‐Decomposition model, (b) Jacchia‐Bowman 2008 (JB2008) (c) MSIS2.0 and (d) CAJIN models. Histograms of model‐observation deviation determined for individual models (blue bars) and Gaussian‐function regression that best fits the empirical distribution (red line) for (e) Eigen‐Decomposition model, (f) JB2008, (g) MSIS2.0, and (h) CAJIN models. The correlation coefficient (R), root‐mean‐square error (root mean square error), mean (μ) and variance (σ 2 ) are also given in the panels. Space Weather 10.1029/2025SW004351 OWOLABI ET AL. 9 of 14 15427390, 2025, 7, Downloaded from https://agupubs.onlinelibrary.wiley.com/doi/10.1029/2025SW004351 by Spanish Cochrane National Provision (Ministerio de Sanidad), Wiley Online Library on [06/12/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