Improving the Thermosphere Ionosphere in a Whole Atmosphere Model by Assimilating GOLD Disk Temperatures

Global‐Scale Observations of Limb and Disk (GOLD) disk measurements of far ultraviolet molecular nitrogen band emissions are used to retrieve column integrated disk temperatures (Tdisk), which are representative of the lower‐and‐middle thermosphere. The present work develops a new approach to assimilate the Tdisk in the whole atmosphere community climate model with thermosphere‐ionosphere eXtension using the data assimilation research testbed ensemble adjustment Kalman filter. Nine days of data, 1–9 November 2018, are assimilated. Analysis state variables such as thermospheric effective temperature (Teff, airglow layer integrated temperature), ratio of atomic oxygen to molecular nitrogen column densities (O/N2), and column electron content are compared with a control simulation that is only constrained up to ∼50 km. It is observed that assimilation of the GOLD Tdisk improves the analysis states when compared with the control simulation. The analysis and model states, particularly, Teff, O/N2, and electron column density (ECD) are compared with their measurement counterparts for a validation of the assimilation. Teff and O/N2 are compared with GOLD Tdisk and O/N2. While, the ECD is compared with ground based total electron content measurements from global navigational satellite system receivers. Root mean square error (RMSE) improvements in Teff and O/N2 are about 10.8% and 22.6%, respectively. The RMSE improvement in analyses ECD is about 10% compared to the control simulation.

• A new approach has been developed to assimilate Global-Scale Observations of Limb and Disk (GOLD) T disk in whole atmosphere community climate model with thermosphere-ionosphere eXtensions which is validated using independent measurements • Analysis states of both the thermosphere and ionosphere show improved agreement with independent measurements • Results demonstrate a great potential of the GOLD T disk data to improve thermosphere-ionosphere data assimilation With time the whole atmosphere ionosphere thermosphere models are improving, and number of observations from the TI system and lower atmosphere are increasing. Therefore, we are in a great stage to do a whole atmosphere data assimilation by combing the models and the observations. There is a long-history of lower atmosphere data assimilation (Gelaro et al., 2017;Hersbach et al., 2020;Rienecker et al., 2011), but the whole atmosphere system data assimilation is relatively new. There have been significant developments in the assimilation of thermosphere-ionosphere observations such as, neutral density (Matsuo et al., 2013;Mehta et al., 2018;M. V. Codrescu et al., 2004;Ren & Lei, 2020;S. M. Codrescu et al., 2018;Sutton, 2018), thermospheric temperature (Laskar, Pedatella, et al., 2021), thermospheric airglow radiance (Cantrall et al., 2019), and electron content (Aa et al., 2016;Bust et al., 2004;Bust & Immel, 2020;Chen et al., 2016;Datta-Barua et al., 2013;Forsythe et al., 2021;He et al., 2020;Kodikara et al., 2021;Lee et al., 2012;Lin et al., 2015;Matsuo et al., 2013;Pedatella et al., 2020;Song et al., 2021). While these results were promising and showed that the assimilation of TI observations improves the model states, most were limited to using upper atmosphere only models or used limited thermospheric datasets from low-earth-orbit satellites or ionospheric only measurements or observing system simulation experiments. Furthermore, a majority of them have not combined lower, middle, and upper atmosphere data in the assimilation. Also, the spatial and temporal coverage of thermospheric data available earlier were limited.
Temperature is one of the basic parameters in whole atmosphere models. Neutral temperature retrieved from Global-scale Observations of Limb and Disk (GOLD) disk measurements have increased the number of thermospheric observation in the recent years, which enables scope for a better whole atmosphere data assimilation that can potentially improve the specification of the TI system. Laskar, Pedatella, et al. (2021) performed a set of Observing System Simulation Experiments (OSSEs) to evaluate the impact of assimilating GOLD disk temperature (T disk ) observations on thermospheric temperature and dynamics. They found that the OSSE that includes the GOLD T disk improved the model temperature root mean square error (RMSE) and bias by 5% and 71% when compared with the forecast state, and the improvements are 20% and 94% when compared with lower atmosphere only assimilation. Laskar, Pedatella, et al. (2021) also found that the migrating diurnal tide (DW1) and local diurnal tide over Americas improve by about 8% and 17%, respectively, upon assimilation of GOLD disk temperature (T disk ) observations. In the current study we assimilate actual GOLD T disk in a whole atmosphere data assimilation system and assess their impact on the thermosphere-ionosphere parameters by validating analysis states with their measurement counterparts.

Data, Models, and Methodology
The primary data set used is the GOLD T disk , which has been assimilated in the whole atmosphere community climate model with thermosphere-ionosphere eXtension (WACCMX). In addition to T disk , lower and middle atmosphere data have also been assimilated. For validation of the analysis states from the assimilation system, independent measurements of GOLD O/N 2 and Global Navigation Satellite System Total Electron Content (GNSS-TEC) are also used. Further details of these data and models are given below.

GOLD T disk
GOLD observed the Earth's thermosphere in the far ultraviolet wavelengths for over 18.5 hr each day, from 0610 to 0040 Universal Time (UT) of the next day (Eastes et al., 2019McClintock et al., 2020;Laskar et al., 2020). The primary GOLD observations are emission intensities in the far ultraviolet (FUV) range of 134.5-166.5 nm. Data for one full disk scan are available at every 30 min from 6 to 23 UT (Eastes et al., 2019). The current investigation uses level 2 T disk data (version 3), that are retrieved from 2 × 2 binned level-1C data, which are available in the GOLD web-page, https://gold.cs.ucf.edu/ as 'Level 2-TDISK'. The retrieval algorithm is an improvement of the previously used methods for limb measurements (Aksnes et al., 2006;Krywonos et al., 2012).
The 2 × 2 binned data have a spatial resolution of 250 × 250 km near nadir and it gets slightly coarse at view angles higher than 45° from nadir. The GOLD daytime disk scans in N 2 Lyman-Birge-Hopfield (LBH) bands are used to retrieve T disk data. Effective altitude and contribution function (CF) of the T disk varies with solar zenith angle (SZA) and emission angle (EA). The SZA variation of the CF is well quantified (Laskar, Pedatella, et al., 2021) and thus is included in the present assimilation. However, the EA effects are not yet included in the assimilation. But, it has been observed that the EA does not have significant impact on the CF for EAs below 50°, so the T disk data having EA > 50° are not included in this assimilation and analysis. This limit also restricts the latitude and longitude coverage, as shown in Figure 1, to about ±50° in latitude and about −10°W to −90°W in longitude. Also, for high SZA observations the signal-to-noise ratio (SNR) is low, which for the current V3 T disk introduces a bias. Thus, the low SNR observations having SZA >65° are not considered in the analysis and assimilation.

GOLD O/N 2
GOLD disk measurements of OI-135.6 nm emission and N 2 -LBH bands in the ∼134-162 nm wavelength range are used to retrieve the ratio of atomic oxygen to molecular nitrogen column densities (ΣO/ΣN 2 ) (Correira et al., 2021). For simplicity we use the notation O/N 2 instead of ΣO/ΣN 2 . The disk O/N 2 has the same spatial and temporal coverage as T disk . O/N 2 data are used here only for the comparison and validation of the analyses O/N 2 . We use the 2 × 2 binned version 3 O/N 2 data, named as 'Level 2-ON2' in the GOLD data repository. Also, as the GOLD O/N 2 is not optimized for auroral latitudes (Correira et al., 2021), the latitudes above ±60° are not used in the current analysis. Typical random, systematic, and model uncertainties in the GOLD O/N 2 are about 5%, 5%, and 30%-40%, respectively. Note that the model uncertainty is a bias with an unknown sign (Correira et al., 2021).

GNSS-TEC
The GNSS-TEC data used in this study are obtained from the madrigal database (https://cedar.openmadrigal.org). Madrigal TEC maps are derived from worldwide Global Navigational Satellite System (GNSS) ground-based receivers. The vertical TEC data are available at 5 min temporal and 1° by 1° spatial bins. Details on the TEC retrieval algorithm can be found in Rideout and Coster (2006) and Vierinen et al. (2016). In the current study the TEC maps are averaged over 20 min centered at every UT hour to compare them with the analysis ECD from assimilation. The 20 min averaging is chosen to get enough satellite passes over a particular spatial grid.

WACCMX
The WACCMX version 2.1 is a whole atmosphere general circulation model extending from the surface to the upper thermosphere (∼500-700 km depending on solar activity) . Whole Atmosphere Community Climate Model with thermosphere-ionosphere eXtension includes the chemical, dynamical, and physical processes that are necessary to model the lower, middle, and upper atmospheres. The thermosphere and ionosphere processes in WACCMX are similar to those in the NCAR Thermosphere-Ionosphere-Electrodynamics General Circulation Model (TIE-GCM), including the transport of O + and self-consistent electrodynamics as well as realistic solar and geomagnetic forcing. The model horizontal resolution is 1.9° × 2.5° in latitude and longitude, and the vertical resolution is 0.25 scale height above ∼50 km.

SD-WACCMX
In this simulation the WACCMX horizontal winds and temperature are relaxed toward Modern-Era Retrospective analysis for Research and Applications, Version 2 (MERRA2) (Gelaro et al., 2017;Rienecker et al., 2011), so we name it as Specified Dynamics WACCMX (SD-WACCMX). The relaxation or nudging to MERRA2 is up to 50 km altitude, and the model is free-running above this altitude (Marsh, 2011). The SD-WACCMX is used in this study as a control simulation. In addition to MERRA2, SD-WACCMX simulations (often referred here as SD) also use operational solar F10.7 cm flux and geomagnetic Kp index for forcing and thus they can be used as a control simulation for the assessment of the data assimilation states.
In the present work we assimilate lower and middle atmosphere as well as thermosphere observations in the WACCMX + DART. The lower atmosphere measurements include conventional meteorological observations (i.e., temperatures and winds from aircraft, radiosonde measurements, etc.), as well as GNSS radio occultation refractivity. Assimilation of these observations improves specifications of the troposphere-stratosphere globally, which is important for the studies of the vertical coupling of waves from lower-atmosphere to the thermosphere (McCormack et al., 2017;Pedatella et al., 2014Pedatella et al., , 2018Wang et al., 2011).
In addition to lower altitude observations, middle atmosphere temperatures from Sounding of the Atmosphere using Broadband Emission Radiometry (SABER) instrument on the Thermosphere Ionosphere Mesosphere Energetics and Dynamics (TIMED) satellite and Aura Microwave Limb Sounder (Aura-MLS) are also used. Altitude coverage of temperature profiles extends from stratosphere to mesosphere-lower-thermosphere (MLT) altitudes (∼15-105 km for TIMED-SABER and ∼15-90 km for Aura-MLS). The latitude coverage of TIMED-SABER retrieved temperature alternates between 83°S and 52°N (south viewing mode) and 83°N-52°S (north-viewing mode) (Remsberg et al., 2008). We performed 9 days (1-9 November 2018) of data assimilation, during which TIMED-SABER was in the north-viewing mode on 1 November only. From 2 to 9 November it was in the south viewing mode. While for the Aura-MLS it varies from 82°S to 82°N (Schwartz et al., 2008). Though Aura-MLS and TIMED-SABER temperatures are middle atmospheric observations, for simplicity we refer to them here as part of lower atmosphere observations. Assimilation of these data has previously been demonstrated to improve specification of the MLT state and dynamics (Laskar et al., 2019;McCormack et al., 2017;Pedatella et al., 2014).
In addition to lower atmosphere observations, GOLD T disk are used in the whole atmosphere assimilation. As the thermospheric dynamics can quickly change in response to changes in forcing conditions, we use a 1 hr assimilation frequency. Additionally, Pedatella et al. (2020) have shown that using a 1 hr data assimilation cycle and removal of second-order divergence damping in WACCMX + DART significantly improves tidal amplitudes, which were previously found to be too small (Pedatella et al., 2018). As full disk images are available at 30 min intervals during sunlit hours, a 1 hr interval will have sufficient data in the thermosphere. Also, the lower atmosphere analysis states in WACCMX + DART agree well with other lower atmospheric assimilations, for example, MERRA2 (McCormack et al., 2021). Figure 1 shows the locations (in a) and altitude or pressure versus number of observations (in b) that are assimilated successfully during a representative hour on a particular day. The red points show the GOLD observations and blue points are the rest of the observations, which we term as lower atmosphere observations, including TIMED-SABER and Aura-MLS. Note that the peak altitude of the N 2 -LBH emission is shown here as a representative altitude of about 150 km, but in the assimilation the impact of T disk is distributed over altitudes based on the SZA dependent CF (Laskar, Pedatella, et al., 2021). One can see that about 70,000 observations per hour are assimilated. On average about 1.5 million observations per day are assimilated. The simulations used in this study are listed in Table 1. The SD-WACCMX is used in this study as the control simulation.
We have performed two WACCMX + DART assimilations. One that assimilates lower atmosphere and GOLD T disk observations, but the direct impact of T disk has been restricted only to the model temperature, referred to as DA1 in Table 1. The second experiment assimilates the same observations as the first experiment, but the T disk observations directly impact the model T, O, O 2 , and O + , referred to as DA2 in Table 1. We used 40 ensemble members in the assimilation. In order to achieve sufficient spread in the ensemble members, we used Gaussian distributions of solar and geomagnetic forcing parameters with mean as the actual value and standard deviations of 15 sfu for F10.7 cm flux and 1 for Kp index (i.e., 10.7 ∼ 10.7, 15 2 and ∼ , 1 2 ). We reset any F10.7 value less than 60 sfu to 60 sfu and any negative Kp to 0. The forcing perturbation for each ensemble member remains the same for all the days. To avoid artifacts arising from initial ensemble members, the spinup duration for the two assimilation runs are about 2 weeks that is, each assimilation run starts from 15 October 2018.

Results
In order to assess and validate the performance of the assimilation we compare the ensemble averaged analysis states to their measurement counterparts. For example, effective temperature (T eff ) from model simulation is compared with GOLD T disk ; O/N 2 is compared with GOLD O/N 2 ; and Electron Column Density (ECD) is compared with the GNSS-TEC. Note that T eff here refers to the vertically integrated GOLD equivalent temperature that is calculated by integrating the model temperature profile weighted by the SZA dependent CFs. Also, the ECD is similar to TEC, but the column integration is only to the topmost layer of WACCMX, which is about 480 km for the current cases. Figure 2 shows a comparison of the local time and latitude variation of the GOLD T disk with T eff from ensemble averaged states of the DA1 (DA1 T eff ) and SD-WACCMX (SD T eff ) for 2 different days. The latitudes and local times are restricted to only those locations and times where GOLD T disk is being assimilated. Beyond those local time and latitudes GOLD data are available, but we are not using them in the assimilation as explained in Sections 2.1 and 2.2. Note that in this figure only a representative longitude of 48°W is shown, which is close to the sub-satellite point of GOLD.
It can be noted from Figure 2 that the broad variations between T disk and DA1 T eff are similar on both the days. On 5 November 2018 there was a moderate geomagnetic storm for which the average temperature is more than 100 K higher than 3 November 2018. Moreover, the morning temperatures are relatively warmer, particularly between 40° and 50°S. These variations of the GOLD T disk during geomagnetic events have been reported and discussed in . These results suggest that the data assimilation is driving the model temperature in the right direction that is, closer to those observed. A quantitative estimate of the differences between them are given later. Note that since both the assimilation experiments updated temperature directly at every assimilation step, the T eff are almost the same for both DA cases. So, the T eff for only the DA1 is shown here. Note. U, V, T, N/A, SD, and DA stands for zonal wind, meridional wind, temperature, not applicable, specified dynamics, and data assimilation, respectively. Also, O, O 2 , and O + refers to the mass mixing ratio of atomic oxygen, molecular oxygen, oxygen ion, respectively.

Table 1 Whole Atmosphere Community Climate Model With Thermosphere-Lonosphere Extension Simulation and Data Assimilation Experiments Used in This Study are Listed
A change in temperature also impacts other states by altering the model dynamics. Therefore, assimilation of T disk can also impact the O/N 2 ratio, which is another primary data set from the GOLD mission. Figure 3 shows  For a quantitative estimation of the above observed differences between actual measurements and their data assimilation equivalents we calculate the Root Mean Square Error (RMSE). The RMSE in SD T eff , DA1 T eff , and DA2 T eff with respect to GOLD T disk are shown in Figure 4a for all 9 days. The RMSE for each day is calculated over the whole disk and local time range as shown in Figure 2 for temperature and Figure 3 for O/N 2 . Note that the temperature RMSEs in the two data assimilation runs are clearly smaller than the SD. Also, the temperature RMSE for the two assimilation runs are almost the same, which is expected as both the assimilations updated model temperature directly. The RMSEs in O/N 2 are shown in Figure 4b. The average O/N 2 RMSEs are better for the two assimilation runs compared with the SD, and DA2 has the best RMSE. The pre-storm RMSEs are smaller compared with storm onset and recovery phase. Average RMSE improvements in effective temperature and O/N 2 compared to the SD are about 10.8% and 22.6%, respectively. The improvements of pre-storm RMSE in T disk and O/N 2 are about 6.4% and 27.9% while during the storm they were about 15.5% and 17.4%, respectively. These results suggest that even though the storm times RMSEs are larger, the percentage improvements are larger too.
For a more robust diagnosis of the relationship between SD T eff , DA1 T eff , and DA2 T eff with respect to GOLD T disk for all the available latitudes and longitudes in the disk scans between 10 and 20 UT during 1-9 November 2018 we make scatter diagrams as shown in Figure 5, where the red color represents high density points. Red (solid) and blue (dashed) lines represent least square fitted straight line and one-to-one (45° slope or gradient equal to one line) relationship. Correlation coefficients and fitted linear equations are also given. From these scatter plots it can be seen that the majority of the T disk versus DA2 T eff points (in a) fall on the one-to-one line. But, for the T disk versus SD T eff (in e) comparison, the highest density observations (red points) deviate away from the one-to-one linear relationship. Also, the correlation coefficient and gradient of the fitted lines are better for the assimilation runs. Note that here also, only those observations are shown that fall within the 50° EA and 65° SZA limits. As the GOLD T disk has higher spread compared to DA2 T eff , DA1 T eff , and SD T eff the shape of the scatter plot is elongated toward the T disk axis (in a, c, and e). Similar to temperature, the O/N 2 scatter diagrams are shown in Figures 5b, 5d, and 5f but the EA and SZA restrictions are not applied here. The correlation coefficients for O/N 2 are small, though they are statistically significant as p-values (probability that the correlation arises from noise) are zero, suggesting a weak linear relationship. As the high density (red) points are mostly located around a circle for the two assimilation cases, the linear correlation would not be a great measure of the relationship between them. Therefore, we calculated the RMSE for the two assimilations and SD with respect to GOLD O/N2. The RMSEs for the DA1, DA2, and SD with respect to GOLD are 0.20, 0.17, and 0.23, respectively, suggesting . Note that the temperature RMSEs in the two DA runs, are clearly smaller than the SD. Also, the average O/N 2 RMSEs are better for the two assimilation runs compared to the SD, and DA2 has the best RMSE.
10.1029/2021JA030045 8 of 13 that the two DA runs perform better compared to SD. The distribution of points in the GOLD T disk versus DA1 T eff and GOLD T disk versus DA2 T eff is nearly identical because the temperature was updated directly in both the assimilations. However, the distributions in O/N 2 in Figures 5b and 5d are significantly different.
The 23% improvement in DA2 O/N 2 , as seen in Figure 4, motivated us to analyze the electron content derived from the assimilations and compare them with independent TEC measurements. Figure 6 shows a latitude versus day-to-day variation of ECD in SD (a), DA1 (b), DA2 (c), and GNSS-TEC (d) centered at around 60°W (±5°) longitude. This spatial bin has been chosen due to the greater availability of GNSS data in this region. As mentioned in Section 2.3, the GNSS-TEC data are averaged over 20 min duration centered at every hour. Note that even with the 20 min averaging, there are missing data, specifically between 20° and 40°N. This figure shows that the magnitudes of electron densities and some of the shape and temporal variabilities of Equatorial Ionization Anomaly (EIA) in DA2 has better agreement with GNSS-TEC compared to the DA1 and SD. Particularly, the northern mid-latitude enhanced DA2 ECDs are in better agreement with GNSS-TEC. A quantification of the improvements is given at the end of this section. Though there are improvements in DA2, the EIA latitude extent and hemispheric asymmetries are not yet well reproduced in the assimilations. This could be due to the fact that the temperature variability cannot fully reflect the changes in the ionosphere as the ionosphere is also influenced by E-region winds in addition to neutral and ionospheric composition changes. We expect to have better agreement in the future when the GOLD O/N 2 and other ionospheric data set are assimilated in addition to the T disk .
For a qualitative assessment of the improvements seen in the ionospheric electron content, a comparison between SD-ECD (green), DA1 ECD (cyan), DA2 ECD (red), and GNSS-TEC (blue) for a limited spatial region is shown in Figure 7a. The RMSE (in Figure 7b) and bias (in Figure 7c) with respect to GNSS-TEC are also shown. Except for November 1st and second and the night hours of each day (shaded regions, when GOLD data are not assimilated), the other days' DA2 ECD has better agreement with GNSS-TEC as can be inferred from the smaller values of the RMSE and bias. Some of the local time variabilities also have better agreement with DA2. For example, the two-peak structures in daytime GNSS-TEC on days 3 and 5 are better reproduced in the DA2 ECD, while that on eighth has not been reproduced. The two peak structure is particularly strong on 3rd November as indicated by downward arrows. Note the dates are in local time at 60°W. Also, the broader shape of the local time variability in GNSS-TEC match better with DA2 ECD as can be seen on most days in Figure 7a. Except for 1st and 2nd November, the night sector (shaded regions) has higher RMSEs, in general during the last 4 hr of each day and particularly at the end of 6th November. This is expected because the GOLD temperature are assimilated only during daylight sector and therefore they are not able to constrain the night-time dynamics. Including ionospheric and O/N 2 observations in the assimilation would improve the results. The purpose of this comparison is to demonstrate that the ionosphere is also improved upon assimilation of GOLD T disk , though there are still large RMSEs and biases. Quantitative estimates of the differences, that vary with latitude and time, are given in Figure 8 and its discussion as given below.
In Figure 3 we show that the GOLD O/N 2 has latitudinal differences from the DA O/N 2 . Also, we have seen in Figure 6 that the agreement between DA2 ECD and GNSS-TEC varies with latitude. To investigate these latitudinal differences in TEC we have calculated the RMSE and mean bias (GNSS-TEC minus DA-ECD) at every 10° latitude bin during the 9 days. The RMSE and mean bias in electron contents from SD, DA1, and DA2 with respect to GNSS TEC are shown in Figure 8. One can note that the lowest values of the RMSE and bias are observed for the DA2, the red lines marked with stars. The RMSE and bias at every latitude bin is calculated from all the 24 × 9 = 216 hr of data. The percentage improvements in RMSE for DA1 ECD and DA2 ECD with respect to GNSS-TEC are about 3% and 10%, respectively. The 9 days average mean biases with respect to GNSS-TEC for the SD, DA1, and DA2 are about 1.9, 0.5, and 0.2 TECu, respectively. Also, the latitudinal average of absolute-biases are 1.92, 1.37, and 1.44 for SD, DA1, and DA2, respectively. Though the latitudinal average of the mean biases is slightly smaller for the DA2 compared to DA1, it is clear, from the absolute values, that the biases are smaller for both the assimilations compared to SD. Also, the mean bias is positive at higher latitudes (>30°) as seen in Figure 8b. Since O/N 2 and TEC vary in proportion, to a large extent, the smaller O/N 2 (from GOLD as shown in Figure 3b at the higher latitudes compared to SD and DA) may produce the positive mean biases in TEC. Negative bias and high RMSE between 0 and 20°S for the DA2 also imply that the equatorial elctrodynamics, which is controlled by ionospheric E-region winds and composition, are not well constrained in the assimilations. Also, the night-time (when GOLD data are not assimilated) electrodynamics, particularly pre-reversal enhancement that is highly variable, contributes to poorer low-latitude results. But, overall these results further emphasize that the DA2-Where in addition to temperature the O, O + , and O 2 mixing ratios are updated directly -has the most improved thermosphere and ionosphere. Overall, it can be observed that the RMSEs are lower in the Northern hemisphere compared to the Southern hemisphere, which suggests that the Northern hemispheric variabilities are better reproduced in the assimilation.

Conclusions
An investigation of the impact of GOLD T disk assimilation on thermosphere-ionosphere states is carried out using WACCMX + DART analysis states, GOLD measurements, and GNSS-TEC. The salient results of this investigation are:  Clearly, for the DA2 the RMSE is smaller compared to other two cases and bias is closer to zero.