Correction of the Photoelectron Heating Efficiency Within the Global Ionosphere-Thermosphere Model Using Retrospective Cost Model Refinement

Many physics-based models are used to study and monitor the terrestrial upper atmosphere. Each of these models have internal parameterizations that introduce bias if they are not tuned for a specific set of run conditions. This study uses Retrospective Cost Model Refinement (RCMR) to remove internal model bias in the Global Ionosphere Thermosphere Model (GITM) through parameter estimation. RCMR is a low-cost method that uses the error between truth data and a biased estimate to improve the biased model. Neutral mass density measurements are used to estimate an appropriate photoelectron heating efficiency, which is shown to drive the modeled thermosphere closer to the real thermosphere. Observations from the Challenging Mini-Payload (CHAMP) satellite taken under active and quiet solar conditions show that RCMR successfully drives the GITM thermospheric mass density to the observed values, removing model bias and appropriately accounting for missing physical processes in the thermospheric heating through the photoelectron heating efficiency.


Introduction
The upper atmosphere is a complex system that responds to solar input, lower atmospheric driving, and internal interactions between the ionosphere and geomagnetic field as well as charged and neutral particles. Physics-based models utilize the current theoretical understanding of the geospace environment to specify the charged and neutral portions of the upper terrestrial atmosphere (Fuller-Rowell and Rees, 1980;Dickinson et al., 1984;Roble et al., 1988;Richmond et al., 1992;Roble and Ridley, 1994;Millward et al., 2001). These specifications can then be used to determine how the upper atmosphere is expected to react at times and places where measurements or observations may not be available, and to conduct theoretical studies (such as those of idealized impulse-response functions) (e.g., Ridley, 2008, 2011;Krall et al., 2009;Qian et al., 2010Qian et al., , 2012. Models, however, have internal biases that result in systematic errors in their output. Some of these biases are caused by inadequacies in the input drivers, others are due to physical approximations, and still more are created by observational uncertainties of internal physical parameters (Anderson et al., 1998;Shim et al., 2012;Schunk, 2014).
Data assimilation is an effective way to remove internal model biases. Ensemble Kalman filter methods have been successfully used in coupled thermosphereionosphere models to reduce the global bias in neutral mass density (Matsuo et al., 2012(Matsuo et al., , 2013 and ionospheric total electron content (Mandrake et al., 2005). Various data assimilation methods have been applied to the Global Ionosphere-Thermosphere Model (GITM) (Ridley et al., 2006;Kim et al., 2008;Morozov et al., 2013;D'Amato et al., 2013). Kim et al. (2008) used the Localized Unscented Kalman Filter to assimilate electron density, ion temperature, and ion velocity in a 1-dimensional version of GITM, and demonstrated that the full 3-dimensional GITM could be adapted to assimilate data using an Ensemble Kalman Filter. Morozov et al. (2013) further developed this method using the Data Assimilation Research Testbed. D' Amato et al. (2013) took another approach, showing how Retrospective Cost Model Refinement (RCMR) could be used within GITM to assimilate observations at a single location.
RCMR is a method for improving model performance by minimizing error over a trailing window using observed and modeled data (D'Amato et al., 2011a;Hoagg and Bernstein, 2012). RCMR is closely related to Retrospective Cost Adaptive Control (RCAC), which was developed for aerospace applications such as spacecraft attitude control (Cruz et al., 2012;D'Amato et al., 2011b;Fuentes et al., 2010). However, RCAC has proved useful in tackling related problems such as input estimation, state estimation, and model refinement for linear and nonlinear applications (D'Amato et al., 2011a;Yan and Bernstein, 2014).
RCMR has previously been used to reduce model bias and perform input and state estimation in the Global Ionosphere Thermosphere Model (GITM). D' Amato et al. (2013) demonstrated how model bias, introduced by neglecting atmospheric cooling processes within GITM, could be removed with RCMR using in situ neutral mass density observations. The present paper demonstrates how RCMR can be used to reduce model bias in the specification of the neutral mass density by improving the estimation of the photoelectron heating efficiency.
The absorption of solar radiation is the dominant process shaping the thermosphere and ionosphere. Solar Extreme Ultraviolet (EUV) and soft X-ray fluxes are responsible for photodissociation and photoionization, which heat the ion, electron, and neutral populations through direct and indirect, chemical, and transport processes. Although the total amount of solar energy input into the terrestrial atmosphere can be measured, the amount of energy that goes into each process is much harder to quantify. Complex models have been developed to handle the treatment of photoelectrons that result from the photoionization process in a much more rigorous manner than has been accounted for in any global ionosphere-thermosphere model (e.g., Torr et al., 1990;Richards, 1991Richards, , 2004Solomon and Qian, 2005, and references therein). Thus, heating efficiencies that define the percentage of the total energy that contributes to a particular process are typically defined (Schunk and Nagy, 2009).
Formally, GITM defines the photoelectron heating as the difference between the total absorption and the ionization rates, multiplied by an efficiency term. This photoelectron heating efficiency controls the fraction of the photoelectrons that deposit energy into the thermosphere, raising the neutral temperature along with heating from other sources, such as auroral electrons, Joule and frictional heating, exothermic chemical processes, conduction, and direct solar input. Calculated heating rates typically vary with altitude, local time, season, and solar cycle. Photoelectrons expend much of their energy during photoionization of neutral particles, though they also transfer heat through collisions with ambient electrons and neutral particles. A mid-latitude example of the importance of various heating sources is given in Figure 1. This figure shows the percent that each heating mechanism contributes to the total temperature as a function of altitude. Although photoelectron heating never clearly dominates the way that chemical and EUV heating do, it is a major contributor to the neutral temperature at high altitudes.
GITM does not currently capture spatiotemporal variations in the photoelectron heating efficiency. GITM models the EUV heating efficiency by setting it to 5%, as is done in other global ionosphere-thermosphere models do. This has changed since the publication of Ridley et al. (2006), which describes an altitude profile of the EUV heating efficiency. The altitude dependence of the EUV heating efficiency was eliminated shortly after the release of the model by completely separating chemical heating from EUV heating (as shown in Figure 1). The photoelectron heating efficiency in GITM is treated in a similar manner: a single value is specified upon input, and is used at all times, locations, and solar activity levels. This introduces an internal model bias that is corrected in the present study using RCMR.

Model
GITM is an ionosphere-thermosphere code that models the atmospheres of the Earth and several other solar planets and moons (e.g., Bell et al., 2011). The terrestrial GITM model solves the continuity, momentum, and energy equations for the major neutral and ion species in a self-consistent manner. GITM can be run at a single geographic location (1-dimensional version) or globally (3-dimensional version) using a spherical grid that stretches in latitude and altitude. Unlike other thermospheric models, GITM uses altitude as the vertical coordinate instead of pressure, facilitating the inclusion of non-hydrostatic solutions.
To refine the model of a physical system, RCMR reinterprets the signals (model inputs and parameters) and systems (physical processes) from an adaptive control perspective and then minimizes a performance variable by recursively estimating a single poorly modeled (or missing) subsystem. The goal of RCMR is to use the minimized performace variable, collected using retrospective optimization on recently obtained signals, to iteratively update the model. Figure 2 illustrates how RCMR is used to modify GITM in the context of this study. The three large grey boxes show the real and GITM upper atmospheres represented as systems of subsystems. The subsystems, or physical processes, are represented by the yellow boxes. External signals from the sun and the lower atmosphere are represented by the sun and earth, respectively. Note that in the initial GITM atmosphere, the physical processes related to heating are missing. This is because a particular parameter, u (the photoelectron heating), cannot be measured in the real upper atmosphere. RCMR (shown as the purple box) provides a way of correcting the GITM output without measuring u. Using accessible signals (the neutral density) from the real (y) and modeled (ŷ) upper atmosphere, the performance (z) is defined as the difference between y andŷ. RCMR then asks the question, "What value of the control inputs should have been used to get better performance?" A retrospective cost function, which relates z to u, is then minimized to obtain the reoptimized control inputs (û, or the photoelectron heating efficiency).û is then used in the next GITM timestep, providing an improved representation of the desired subsystem. These estimates cannot be applied retroactively, as the system has already evolved into a new state. However, the effort is not wasted as RCMR uses the optimized sequence of control inputs to modify the estimate, thereby improving future model performance. Essentially, the algorithm learns from its past "mistakes" to continually correct or refine the model.
A convenient feature of RCMR is the fact that it does not explicitly require a model of the main system. The only modeling information required is a collection of Markov parameters, which are coefficients that define the impulse response function of z to the change inû. In this case (as in most cases), only the first Markov parameter (H) is necessary. H can be obtained though a combination of numerical testing and studying the qualitative behavior of the system. Once the correct sign and order of this number are obtained, its magnitude can be adjusted to alter the convergence rate of the estimate. Figure 3 shows the estimated photoelectron heating efficiency as a function of model days for a series of RCMR-GITM runs using truth data. The truth data has a photoelectron heating efficiency of 0.08, but each RCMR-GITM run began with an initial guess of 0.12 and a value of H ranging from 1 to 1,000. RCMR successfully converges to the true photoelectron heating efficiency for all of these H values, though when H is closer to 50 the convergence is quicker and the variations are less extreme.
Seven runs are used in this study: two runs use simulated satellite data from 'truth' model runs, demonstrating that RCMR is capable of estimating a range of photoelectron heating efficiencies, two runs use GITM without RCMR, and three runs use real satellite data. The dates for these runs were chosen to lie near solar maximum (2002)  All of these GITM runs use realistic magnetospheric and solar drivers, an apex magnetic field, and the MSIS lower atmosphere and tides at the lower boundary. High latitude drivers are integrated into GITM using the Weimer electric potential patterns (Weimer, 1996) and the conductivity patterns determined by Fuller-Rowell and Evans (1987), which account for different auroral activity levels and provide GITM with appropriate electron precipitation. The runs use an H of 50, have a 5 • × 2 • (longitude by latitude) resolution, and have a variable altitude resolution of 50 grid points that span 100-600 km. The results of these model runs can be accessed at the Virtual Model Repository (VMR) at http://vmr.engin.umich.edu/Model/ gitm under the run names specified in Table 1.
The truth data used as input to RCMR for the two truth runs are produced by extracting neutral mass density from an artificial satellite track run through a GITM simulation. This truth satellite data follow the CHAMP orbit, but have a constant altitude of 400 km. The data are extracted at a 60 s resolution. F 10.7 is set to a constant value of 150 sfu, simplifying the solar EUV input to reduce temporal variations. The photoelectron heating efficiency is set to 0.08 and 0.04 for the first and second runs, respectively. These values flank the former standard GITM photoelectron heating efficiency of 0.06, which was chosen based on the values for photoelectron heating efficiency used in other global ionosphere-thermosphere models. GITM is allowed to run for one day  Chamberlin et al. (2007) and Chamberlin et al. (2008)). The 'realistic' 2005 run uses the same sources to describe the HP and the solar irradiance. The 'unrealistic' 2005 runs, however, have HP set to a constant low value and use the daily measured F 10.7 to represent the solar irradiance instead of the more accurate FISM model. These 'unrealistic' runs test the robustness RCMR, showing whether or not the model refinement method is able to handle incorrect energy inputs. In Table 1 and the rest of the paper the 2002 and 2005 runs performed with FISM and the NOAA HP are referred to as realistic, while the runs with F 10.7 and with constant low NOAA HP are referred to as unrealistic.

Observations
The neutral density for this study is measured by two satellite missions, the CHAllenging Minisatellite Payload (CHAMP, Reigber et al. (2002)) and the Gravity Recovery and Climate Experiment (GRACE, Tapley et al. (2004)). Both of these satellite missions carry accelerometers that allow accurate measurements of the neutral mass density. The CHAMP satellite flew between 15 July 2000 and 19 September 2010. It was initially launched into a near-circular orbit with an inclination of 87 • and an altitude of 454 km. With the aid of 2 orbital changes in 2002, the satellite altitude was approximately 385 km in 2002 and 366 km in 2005. The CHAMP neutral mass density is used as truth data in the second half of this study.
The GRACE mission consists of two satellites, which were both launched into orbit on 17 March 2002. These satellites have inclinations of 89.5 • and initial altitudes of 500 km. By 2005 the orbital altitude had decayed to about 482 km. The two satellites follow the same orbital path, and are separated by a distance of about 220 km. The GRACE neutral mass density is not assimilated in the  Figure 4).

Proof of Concept
Before attempting to correct internal model bias through a parameter using observations, it must first be proved that RCMR is capable of correcting a perfectly known bias. This is done using the two truth runs, which are described in Section 2. Figures 5 and 6 show the estimated photoelectron heating efficiency (a) and error in the estimated neutral mass density along the CHAMP (b) and GRACE (c) orbits. These plots show that, as RCMR drives GITM towards the truth neutral density near the CHAMP orbit, the estimated photoelectron heating efficiency converges to the value used in the production of the truth dataset. Comparing panel (c) with the other panels in Figures 5 and 6 shows that the GITM thermosphere is responding to the changes driven by RCMR at locations other than those of the estimation.
Comparing each panel in Figure 5 to the corresponding panel in Figure 6 shows that RCMR does not require that the initial guess for the photoelectron heating efficiency be particularly close to the truth value. In fact, RCMR converged to the truth photoelectron heating efficiency faster for the 0.04 run, where the truth and estimated photoelectron heating efficiencies differ by less than 1% after 4.83 days, than for the 0.08 run, which took 6.50 days to reach the same degree of agreement. Once the truth and estimated photoelectron heating efficiencies have reached this level of agreement, the neutral mass density error on average reduces by as much as two orders of magnitude, as specified in Table 2. These statistics are computed using either all of the data before the convergence of the photoelectron heating efficiency, or all of the data after convergence. Once the photoelectron heating efficiency has converged onto the truth value, the neutral mass densities from the RMCR run and the truth run agree equally well along the GRACE-B and CHAMP orbits.

Realistic Application
In this section, real observations are used as input to allow RCMR to adapt its internal model bias correction within GITM to realistic, changing geophysical conditions. Recall that the two time periods considered in this study, described in Section 2, encompass high and low solar activity levels as well as geophysically active and quiet times. Figure 7 shows the standard GITM and CHAMP neutral mass densities for the 2002 (a) and 2005 (b) periods, with the output averaged over an orbital period (about 92 minutes for CHAMP and about 94 minutes for GRACE-B) to smooth over the altitude-dependent and day-to-night variations in neutral density. Note that with realistic solar inputs (shown in Figure 7 (a), described in Section 2), GITM successfully produces the major features seen in the neutral mass density. However, in both cases GITM typically underestimates the neutral mass density, though not by a constant amount.
When GITM is run with RCMR using CHAMP observations as input, the modeled neutral mass density quickly begins to track the density specified by the observations. Figure 8 shows the results of the RCMR-GITM runs using   the same format as Figure 7. Now that there is no truth photoelectron heating efficiency for RCMR to converge to, the performance of the model must be evaluated on its ability to reproduce the satellite observations. To this end, an error criteria is established that considers both the difference between the RCMR-GITM and CHAMP neutral mass densities, and the rate of change of that difference with time. Examination of the runs showed that when the RCMR-GITM neutral density tracks the CHAMP observations, the difference is small and does not change rapidly, thus a maximum average neutral mass density difference of 5.0 × 10 −13 kg m −3 and a maximum rate of change of the average neutral mass density difference of 4.0 × 10 −16 kg m −3 s −1 was chosen for times not associated with gaps in the observational data. Using this difference criteria, the realistic runs take a similar amount of time as the truth runs to 'converge'. The 2002 run meets the difference criteria in 6.82 days, the unrealistic 2005 run meets the difference criteria in 4.21 days, and the realistic 2005 run meets the difference criteria in 4.19 days. Comparing the agreement of GITM and RCMR-GITM to the CHAMP neutral mass density measurements for the period of time where RCMR-GITM has met the difference criteria (see Table 3) shows that the percent difference between the model and measurements -as well as the variation about the difference -reduce significantly. Thus, not only has the offset between the model and the measurements been essentially eliminated, but the temporal variations are now being much more consistently reproduced.
Away from the CHAMP orbit, RCMR-GITM improves the comparison between the simulation and the measured variations in the neutral mass density. Figure 9 shows the GITM and GRACE-B neutral mass densities over the three RCMR runs. Note that, especially for the realistic runs, RCMR-GITM produces enhancements and depletions in the neutral mass density at the same time   Table 4 shows that RCMR-GITM consistently overestimates the neutral density measured by GRACE-B, even though using RCMR greatly reduces the variation of this offset. This implies that there could possibly be an offset between the CHAMP and GRACE neutral mass densities, which was also found by Matsuo et al. (2012), who had to scale the GRACE satellite data by a constant fraction in order to match the CHAMP data when conducting data assimilation experiments. These offsets that are observed in Figure 9 suggest that RCMR-GITM would be a useful tool for comparing such datasets. By using the CHAMP measurements to adjust the thermospheric temperature, RCMR-GITM has produced a thermosphere scaled to the CHAMP specifications that accounts for coupled processes. Thus CHAMP and GRACE thermospheric specifications can be compared even when the satellites are not co-located, though the comparison alone does not provide sufficient insight to allow investigation into the reason behind any differences in the satellite observations. The photoelectron heating efficiency estimates (shown in Figure 10) do not converge to a single value, as they did in the ideal cases. The variations reflect the combined changes in the photoelectron heating efficiency and other possibly missing physics or other incorrectly tuned parameterizations in the thermospheric temperature calculation. The difference in the degree of the estimated photoelectron heating efficiency variations between the realistic and unrealistic runs shows how much the internal bias is reduced by including the HP and using a more accurate solar EUV specification. It also shows that the estimated photoelectron heating efficiency does not need to settle to a constant or slowly varying value in order to compensate for the internal model error.
Comparing estimated photoelectron heating efficiencies after the data have met the difference criteria from the 2002 and 2005 realistic runs suggests that there may be a solar activity dependence in the photoelectron heating efficiency. The mean photoelectron heating efficiency estimate from 2002 is 0.164 ± 0.015, while the mean photoelectron heating efficiency estimate from 2005 is 0.134 ± 0.014. However, as highlighted by the unrealistic 2005 run, any variations in the photoelectron heating efficiency estimate will also include internal error corrections (such as those cause by modeling resolution, uncertainties in other physical parameters, and inaccuracies at the upper and lower boundaries and drivers). It is therefore impossible to examine the estimated photoelectron heating efficiency as if it were a corrected photoelectron heating efficiency analogous to that in the real thermosphere. Thus, while this study cannot state that there is a solar activity variation in the photoelectron heating efficiency, it can state that the solar activity sufficiently influences the ionosphere-thermosphere system (even using the most accurate specification of EUV irradiance currently available) that non-assimilative GITM runs should adjust the photoelectron heating efficiency for different levels of solar activity. It should also be emphasized that photoelectron heating efficiency estimates provided by RCMR-GITM would be most appropriately used in non-assimilative GITM runs with the same resolution, model version, and input types. Other global ionosphere-thermosphere models, however, can use the magnitude of the estimated photoelectron heating efficiency produced in the 2002 and 2005 realistic runs to guide their own parameterizations.

Conclusions
RCMR, a cost-effective method for model refinement, has been shown to successfully improve the performance of GITM, an upper atmospheric general circulation model. Using limited information about the internal physical systems and prior observations, RCMR was configured to use in situ measurements of neutral mass density to correct internal biases in the thermospheric temperature calculation. The configuration process involved choosing a single Markov parameter, which was found by numerical testing using a truth dataset, and an internal model parameter sensitive to both the neutral mass density and the thermospheric temperature. Once configured, RCMR-GITM proved capable of reproducing a specified thermospheric state at and away from the location of the assimilated data with initial photoelectron heating efficiency estimates that were close to and far from the truth value.
The results obtained using actual CHAMP measurements of the neutral mass density to drive the photoelectron heating efficiency with RCMR-GITM proved the usefulness and versatility of the method. RCMR-GITM performed well under both ideal (realistic solar and magnetospheric drivers during solar maximum) and unideal (parameterized solar input and missing magnetospheric drivers during solar minimum) conditions. Comparing the estimated photoelectron heating efficiencies produced by the realistic and unrealistic 2005 runs highlights the importance of using realistic solar and magnetospheric inputs. It also shows the robustness of the RCMR method in accounting for a wide variety of internal model biases even if these biases change rapidly in universal time. Comparing the estimated photoelectron heating efficiencies from realistic 2002 and 2005 runs suggested that solar activity should be considered when specifying this input parameter in non-assimilative GITM runs. This shows that RCMR is capable of reducing the number of model runs needed to perform a parameter adjustment study. Finally, the CHAMP-aligned RCMR-GITM thermosphere proved to be a useful tool for comparing the CHAMP and GRACE datasets.