A Comparative Assessment of the Distribution of Joule Heating in Altitude as Estimated in TIE‐GCM and EISCAT Over One Solar Cycle

During geomagnetically active times, Joule (or frictional) heating in the Lower Thermosphere‐Ionosphere is a significant source of thermal energy, greatly affecting density, temperature, composition and circulation. At the same time, Joule heating and the associated Pedersen conductivity are amongst the least known parameters in the upper atmosphere in terms of their quantification and spatial distribution, and their parameterization by geomagnetic parameters shows large discrepancies between estimation methodologies, primarily due to a lack of comprehensive measurements in the region where they maximize. In this work we perform a long‐term statistical comparison of Joule heating as calculated by the NCAR Thermosphere‐Ionosphere‐Electrodynamics General Circulation Model (TIE‐GCM) and as obtained through radar measurements by the European Incoherent Scatter Scientific Association (EISCAT). Statistical estimates of Joule heating and Pedersen conductivity are obtained from a simulation run over the 11 year period spanning from 2009 until 2019 and from radar measurements over the same period, during times of radar measurements. The results are statistically compared in different Magnetic Local Time sectors and Kp level ranges in terms of median values and percentiles of altitude profiles. It is found that Joule heating and Pedersen conductivity are higher on average in TIE‐GCM than in EISCAT for low Kp and are lower than EISCAT for high Kp. It is also found that neutral winds cannot account for the discrepancies between TIE‐GCM and EISCAT. Comparisons point toward the need for a Kp‐dependent parameterization of Joule heating in TIE‐GCM to account for the contribution of small scale effects.

energy source for the system-as a source of thermal energy (Ching & Chiu, 1973;Knipp et al., 2004;Roble and Matsushita, 1975).For example, Knipp et al. (2004) used empirical models to estimate the relative contributions of solar EUV power, Joule heating power and precipitating particle power to the Earth's LTI energy budget over three solar cycles.They found the daily averages of solar EUV, Joule heating and precipitating particle power to be 464, 95, and 36 GW, respectively; however, they found over 60 extreme events where the Joule heating power greatly exceeded the solar EUV power, reaching over 1000 GW for 12 of these events.Another notable issue is the heating efficiency of each heating mechanism: whereas for particle precipitation and solar EUV the heating efficiencies are on the order of 50% (Rees et al., 1983;Torr et al., 1980), the efficiency of Joule power is nearly 100% according to Thayer and Semeter (2004), meaning that even slightly elevated Joule heating power can exceed solar EUV power as the dominant heating mechanism in the LTI.
Even though the physics of the process by which electromagnetic energy is converted into heat into the LTI is well understood (e.g., Aikio et al., 2012;Brekke & Kamide, 1996;Lu et al., 1995;Strangeway, 2012;Thayer & Semeter, 2004;Vasyliūnas & Song, 2005; X. Zhu et al., 2005), and is captured in physics-based Global Circulation Models (GCMs), the quantification and spatial distribution of Joule heating in altitude, latitude and longitude is still largely unknown (e.g., Palmroth et al., 2005Palmroth et al., , 2020; T. Sarris et al., 2023).A key reason is that Joule heating in the LTI is influenced simultaneously by neutral winds, plasma turbulence, variable electric fields, and modifications in conductivity by precipitating particles and by strong electric fields, all of which are largely variable during active times.In addition, the effects of small scale structures in the electric field, that are known to be present within the auroral oval, are not well quantified.It thus comes as no surprise that large discrepancies appear in estimates of Joule heating between different models, which often disagree by large factors (Perlongo et al., 2018; T. E. Sarris, 2019), making their cross-comparison and evaluation difficult (Scherliess et al., 2019).Furthermore, large discrepancies appear between proxies of Joule heating that are based on solar and geomagnetic activity indices, as shown in Figure 3 of Pirnaris et al. (2023).
A key unknown in the spatial distribution of Joule heating is the altitude where it maximizes and the dependence of this altitude on geomagnetic activity.For example, Griffis et al. (1981) used a few days of data (late December of 1974) from the Atmosphere Explorer C satellite as input to a computer program and derived the maximum Joule heating value at 115 km.Deng and Ridley (2007) evaluated model data from GITM, the Global Ionosphere Thermosphere Model (Ridley et al., 2006) in a case study and found the maximum at 120 km.Huang et al. (2012) used data from TIE-GCM, the NCAR Thermosphere Ionosphere Electrodynamics General Circulation Model (TIE-GCM) (Qian et al., 2014;Richmond et al., 1992), spanning about a month, and found the largest Joule heating deposition at 125 km.Banks (1977) used data from the Chatanika, Alaska auroral-zone incoherent-scatter radar to compute altitude profiles for the polar cleft region for the 4 August 1972, solar proton event, and found the maximum Joule heating rate at 128 km.The dependence of Joule heating in Magnetic Local Time (MLT) has been noted by Brekke and Rino (1978), whereas the dependence on neutral winds and geomagnetic activity on the altitude profiles of Joule heating has been highlighted by Thayer (1998) and Cai et al. (2013).
The exact quantification of Joule heating requires the simultaneous and co-located measurement of an extensive list of parameters, which is only available in situ.However, the in-situ observation of the Lower Thermosphere-Ionosphere at the 100-200 km transition region where Joule heating maximizes presents many challenges, as this altitude range is too high for balloons and too low for systematic measurements by orbiting satellites.Sounding rockets are able to sample this region (e.g., Sangalli et al., 2009), however their measurements are near-instantaneous, providing, essentially, snapshots of altitude profiles above the location of the launch site.It is for this reason that systematic measurements of global coverage do not exist, impeding the accurate representation of Joule heating in models of the upper atmosphere (Ruan et al., 2018).
In the absence of systematic in situ measurements of all relevant parameters needed to quantify Joule heating, one of the proposed methodologies that has been used extensively in the past to approximate Joule heating is via Incoherent Scatter Radar (IS) data; examples of such calculations are included in, for example, Thayer (1998), Aikio andSelkälä (2009), andAikio et al. (2012).It is noted that Joule heating and Pedersen conductivity cannot be directly obtained from the IS radar measurements, which provide electron densities and ion velocities to be used in the estimates; thus Joule heating is obtained via subsequent processing of these basic measurements with the aid of the empirical MSIS neutral atmosphere model, as described in the studies listed above, and references therein.The comparison between a commonly used physics-based GCM, such as TIE-GCM and GITM, and a ground-based radar such as EISCAT, can help in identifying corresponding discrepancies and missing physics or unknown assumptions, which is the scope of this paper.
Most of the previous studies have focused on quantifying Joule heating during single events or for periods of a few days.Instead, in this study we investigate the distribution of Joule heating statistically over long (multi-year) periods of time in terms of altitude, magnetic local time and geomagnetic activity, as estimated (a) via the NCAR TIE-GCM and (b) via measurements from the European Incoherent Scatter (EISCAT) Tromsø UHF radar (Rishbeth and Williams (1985)).Together with Joule heating, Pedersen conductivity is also evaluated.Both TIE-GCM model outputs and EISCAT statistics are gathered for the 11-year period from 2009 to 2019, corresponding approximately to solar cycle number 24.Statistics are inter-compared in terms of percentiles within altitude bins, local time sectors and geomagnetic activity levels.It is noted that the EISCAT Tromsø UHF radar does not operate continuously, but during campaigns, as discussed in further detail below.
In the following, in chapter 2 we present the methodology that is used in the estimation of Joule heating and Pedersen conductivity in TIE-GCM, including the parameterization of the TIE-GCM run, and the methodology by which Joule heating and Pedersen conductivity are obtained via EISCAT measurements.In chapter 3 we discuss the statistical sampling and the segmentation of data that was performed in terms of altitude, MLT and geomagnetic activity, and we present the results of the statistical analysis.In chapter 4 we discuss the results and potential reasons for the observed discrepancies, and in chapter 5 we summarize the key findings of this study, concluding with implications from the observed discrepancies and potential methodologies and measurements that are needed to resolve them.

TIE-GCM
The NCAR TIE-GCM is a first-principles model of the coupled thermosphere and ionosphere system that solves the three-dimensional momentum, energy and continuity equations for neutral and ion species at each time step (e.g., Qian et al., 2014).Solutions are performed on 29 constant-pressure levels, extending from approximately 97-500 km in intervals of one-half scale height.Main assumptions of TIE-GCM include hydrostatic equilibrium, incompressibility on constant pressure surfaces, constant gravity, and steady-state ion and electron energy equations.
The electric field in TIE-GCM is imposed based on an externally driven geopotential field, which, in the case of the model runs of this study, is introduced based on the Weimer electric potential specification (Weimer, 2005).Furthermore, TIE-GCM assumes that the electric field is always perpendicular to the magnetic field of the Earth.An IGRF model is used for the Earth's magnetic field (Thébault et al., 2015).The high latitude energy input associated with auroral particle precipitation is represented by an analytical auroral model (Emery et al., 2012;Roble & Ridley, 1987).Further to the above, additional external drivers and parameterizations that are used in the TIE-GCM include: an empirical model that is used to specify photoelectron heating; an empirical model that is derived from two-stream calculations to specify the production of secondary electrons; empirical formulations that are used to specify the upper boundary conditions for electron heat transfer and electron number flux; an eddy diffusion formulation to include the effects of mixing caused by gravity waves; and the Global Scale Wave Model to specify atmospheric tides at the lower boundary.
The Joule heating rate in TIE-GCM is calculated based on the following equation (e.g., Lu et al., 1995, Equation 3): where    is the neutral wind vector, σ P is the Pedersen conductivity,  ⃗  is the local geomagnetic field vector and  ⃗  is the electric field component that is perpendicular to the geomagnetic field.In TIE-GCM it is assumed that magnetic field lines are equipotentials, and thus that  ⃗ ‖ = 0 and therefor  ⃗  =  ⃗ ⟂ .
Pedersen conductivity in TIE-GCM is calculated by the following equation (e.g., Schunk & Nagy, 2009, Equation 5.117):  Schunk and Nagy (2009).In a recent study, Ieda (2020) discussed the effects of including the non-resonant interactions in the calculation of the ionospheric conductivities.For the O 2 −   + 2 case they concluded that non-resonant collisions are possible for temperatures below 600 K (in contrast with the traditionally assumed 800 K limit) thus lowering the transition altitude to about 145 km, inside the E region.For the O−O + case, the transition temperature is quite low (about 231 K) and this pair is mostly located at higher altitudes.The current study focuses on altitudes lower than 150 km and omits the non-resonant interactions.As future work, it would be possible to include the Equation A1 of Ieda (2020) in the TIE-GCM source code and study its effect on the calculation of the conductivities and, subsequently, its effect on Joule heating.TIE-GCM can be executed with various time steps and grid resolutions.For the current study a grid size of 2.5° in geographic latitude and 2.5° in geographic longitude was used.This corresponds to a horizontal distance of ∼284 km in latitude and ∼95 km in longitude at the location of the EISCAT radar.In the vertical dimension, TIE-GCM uses log pressure Z = ln(p 0 /p) as the vertical coordinate with the reference pressure set to p 0 = 5 × 10 7 hPa.Four grid points per scale height were used, leading to 57 pressure levels.For each pressure level the corresponding altitude was also calculated.For the purpose of this study, TIE-GCM was run over the period from 2009 until 2019, with a run-time resolution of 30 s and an output time resolution of 2 hr.In total, Joule heating was estimated at ∼10 8 points.

EISCAT Radar Measurements and Data Analysis
Incoherent Scatter Radars (ISRs) have often been used to provide detailed information on both the vertical structure and the temporal evolution of the LTI, including estimates of Joule heating (e.g., Aikio & Selkälä, 2009;Aikio et al., 2012;Banks, 1977;Cai et al., 2013;Fujii et al., 1999;Thayer, 1998Thayer, , 2000;;Wickwar, 1974;Vickrey et al., 1982).In the present study, Joule heating was estimated statistically from measurements provided by the European Incoherent Scatter (EISCAT) radar facility in Ramfjordmoen, near Tromsø, Norway (geographic: 69.59°N, 19.22°E, cgm: 66.58°, 102.94°).The site hosts a UHF radar, which operates in the 930 MHz band with transmitter peak power 2.0 MW, a 12.5% duty cycle and 1 µs-10 ms pulse length with frequency and phase modulation capability.The antenna is a 32 m mechanically fully steerable parabolic dish used for transmission and reception.EISCAT measurements used in this study were collected from an ensemble of campaigns conducted in the period from 2009 to 2019.EISCAT measurements correspond to a total of 29,806 Pedersen conductivity and 23,938 Joule heating samples that are used in this study.These samples come from 134 different days.The 3-hr Kp index was between 0 and 2 during 118 days, between 2 and 4 during 80 days and between 4 and 9 during 9 days.It is noted that the radar did not operate continuously, and thus the total number of data points from all the campaigns are not evenly distributed during the above-mentioned 11 year period; however measurements that were obtained from the model run cover all relevant ranges of Kp, MLT and altitude.
For the electric field estimation from the EISCAT Tromsø UHF radar the beam-swing method was utilized, where the antenna points in different directions in a short cycle to capture three non-coplanar ion velocity components; these experiments are referred to as cp2 or ip2.The EISCAT radar observes the autocorrelation function (ACF) of the ionospheric scattering process in many range gates along the radar beam, and plasma parameters are fitted to the observed ACFs by means of least-squares fits using the Grand Unified Incoherent Scatter Design and Analysis (GUISDAP) software (Lehtinen & Huuskonen, 1996).The data analysis is run by the EISCAT association, and the results were downloaded from the Madrigal database.The electric field analysis was carried out as follows: F-region line-of-sight ion velocity measurements from at least three different beam pointing directions from 200 to 350 km altitudes were combined.The ion velocity vector was solved by means of a least-squares fit to line-of-sight velocities using a statistical inversion in the same manner as described in Nygrén et al. (2011).It is noted that a fit can fail if data from some of the beam directions is missing, for example, due to technical problems.It is also noted that unreliable results have been removed if the following criteria were met: (a) the chi-square of the velocity vector fit was larger than 50; (b) the standard deviation of the electric field was larger than 200 mV/m; (c) the error correlation (i.e., the error covariance divided by the standard deviations) of the two electric field components was larger than 0.95, which usually occurs when data from one of the beam directions is noisy.The electric field was obtained by assuming that ions follow the  ⃗  × ⃗  drift in the F region and the magnetic field vector was obtained from the IGRF model.The median largest horizontal distance between the three beams at 300 km altitude is 85 km.
Altitude profiles of Pedersen conductivity were calculated via the electron density profiles measured by the EISCAT UHF radar looking in the field-aligned direction, using the methodology described in, for example, Brekke and Hall (1988), Moen et al. (1990), and Aikio et al. (2012).The equation used for Pedersen conductivity is basically the same as Equation 2, but the densities of molecular ions NO + and   + 2 are combined, because they cannot be separated in the EISCAT analysis.The ion-neutral and electron-neutral collision frequencies needed in the Pedersen conductivity were calculated according to Brekke and Hall (1988), and the required neutral parameters were obtained from the NRLMSISE-00 empirical model (Picone et al., 2002).
The radar scan cycle duration is typically 4 or 6 min, and 6 min resolution was used for the electric field analysis for consistency.As a result, one or two field-aligned beams that produce Pedersen conductivity profiles are included in each 6-min period.Electric field values were calculated for all Pedersen conductivity profiles by means of linear interpolation.As the electric field analysis sometimes fails even if the field-aligned electron density profile was successfully measured, and the electric field over failed fits is not interpolated, the number of Pedersen conductivity profiles in the final data set is somewhat larger than the number of Joule heating profiles.
The Joule heating rate was subsequently calculated according to: where σ P is the Pedersen conductivity and  ⃗  is the electric field.It is noted that the electric field in EISCAT is measured in a coordinate system that is fixed to the Earth.

Statistical Distributions
For the purpose of this study, and in order to compare TIE-GCM and EISCAT data, TIE-GCM outputs were used from the geographic latitudes that are closest to the geographic latitude of the EISCAT Tromsø radar, which is located at 69.6°.These include the TIE-GCM grid points with geographic latitudes 68.75° and 71.25°.In altitude, TIE-GCM outputs were averaged in ranges of 4 km each, spanning altitudes from 100 to 150 km.In MLT, outputs were averaged in four MLT sectors, ranging from 15:00 to 21:00, 21:00 to 03:00, 03:00 to 09:00, and 09:00 to 15:00.In terms of geomagnetic activity, outputs were binned in three ranges in terms of the Kp index, with ranges of 0 ≤ Kp < 2, 2 ≤ Kp < 4, 4 ≤ Kp < 9, which in the following are referred to as low, medium and high activity, respectively.Within each of the above bins, the 10th, 25th, 50th, 75th, and 90th percentile values were calculated.EISCAT outputs were also binned similarly in terms of MLT and geomagnetic activity.Within each of the above bins, similarly to the calculations of TIE-GCM gridded data, the 10th, 25th, 50th, 75th, and 90th percentile values were calculated.As it is expected, the 11 year period is dominated by quiet days, with the ratio of low to medium activity periods being 2.2 to 1 and the ratio of low to high activity periods being 5.6 to 1.

Results
The statistical distribution of Joule heating as a function of altitude, MLT and Kp as obtained in TIE-GCM and EISCAT are shown in the top and bottom panels of Figure 1 respectively.The results of TIE-GCM are colored in blue hues and of EISCAT in green hues, and the same coloring scheme is followed throughout this paper to differentiate between TIE-GCM and EISCAT measurements.The three rows in each panel correspond to the three ranges in Kp: low (0-2), medium (2-4) and high (4-9) levels of geomagnetic activity.The four columns in each figure correspond to the four MLT sectors, divided in the afternoon (15-21), midnight (21-03), morning (03-09) and noon (09-15) regions.Color shades under the altitude distribution curves in these figures indicate the percentiles of Joule heating, with the darker shaded area under the thicker black line corresponding to the median values of Joule heating as a function of altitude.Lighter-shaded areas under the thin gray curves show progressively higher and lower percentiles to the right and left of the median, respectively.
For TIE-GCM statistics, the number marked in blue at the lower right corner of each panel indicates the height-integrated value of the median of Joule heating, in units of mW/m 2 .This corresponds to the area at the left of the median (thick line) and was computed by first calculating the median Joule heating at each altitude, then constructing a median Joule heating altitude profile, and finally integrating this profile over altitude.The number in black at the top right corner indicates the total number of grid points that were used in the statistical calculations.Similarly, for EISCAT statistics, the number in green indicates the height-integrated value of the median of Joule heating, whereas the number in black indicates the total number of radar profiles that were used in the statistical calculations.
The comparison between the upper and lower panels of Figure 1 shows that Joule heating estimates are lower in TIE-GCM than in EISCAT for higher levels of geomagnetic activity, and, inversely, that they are generally higher in TIE-GCM than in EISCAT for lower levels of geomagnetic activity.An exception is observed for the noon sector (MLT 09-15), where EISCAT measurements show very low levels of Joule heating for all levels of geomagnetic activity.It is also observed that EISCAT observations show considerably larger deviations from the median, in particular for higher Joule heating values.
With respect to the altitude of the peak of Joule heating, in TIEGCM the maximum is found within the 119-123 km altitude bin, whereas in EISCAT it is found within the 120-121 km altitude range for low kp, 119-121 km for medium Kp and 118-121 km for high kp levels.Thus, on average, the altitude of the peak of Joule heating is found to decrease slightly on average for higher levels of geomagnetic activity.It is also noted, however, that the distribution shows a large deviation from the median, due to the smaller number of profiles for high kp, and thus this dependence of the peak altitude on geomagnetic activity could be a statistical artifact.
Figure 2 displays the Pedersen Conductivity as a function of altitude, MLT and Kp, as calculated from TIE-GCM (top panel) and EISCAT radar measurements (bottom panel), in a similar format as that of Figure 1.
From Figure 2 it can be seen that, similarly to Joule heating estimates, Pedersen Conductivity is again overestimated in TIE-GCM for low Kp, whereas it is underestimated for high Kp.Furthermore, also similarly to Figure 1, values from TIE-GCM display a more uniform distribution than those from EISCAT.Finally, the peak of Pedersen conductivity is found in the altitude range from 119 to 123 for TIE-GCM and within the 120-121 km altitude range for low Kp, 119-121 km for medium Kp and 118-120 km for high kp levels for EISCAT, demonstrating a dependence on geomagnetic activity, with higher Kp leading to lower altitudes of the peak of Pedersen conductivity.

Discussion
In comparing the long-term (solar cycle) averages of Joule heating and Pedersen conductivity between TIE-GCM and EISCAT, it can be seen in Figures 1 and 2 that, in general, TIE-GCM tends to under-estimate (over-estimate) Joule heating and Pedersen conductivity for high (low) levels of geomagnetic activity compared to EISCAT measurements, with the exception of the noon sector (09-15), where EISCAT measurements are generally low.This could be attributed to the fact that, in the dayside, the auroral oval and the plasma convection cells are typically pole-ward of Tromsø, and therefore the electric field as well as auroral conductances (and hence Joule heating) are small.
A key aspect to note when assessing the causes of the observed discrepancies between EISCAT and TIE-GCM is that EISCAT calculations of Joule heating do not take into account the effects of neutral winds, as the latter cannot be measured directly from ISRs.
In some cases, the E-region neutral wind has also been derived from the IS radar observations (see, e.g., Aikio et al., 2012;Cai et al., 2013;Nygrén et al., 2011;Thayer, 1998); we note, however, that the inversion works more reliably (with less uncertainty) for the electric field than the neutral wind, which, unlike the electric field, depends on altitude (Nygrén et al., 2012).For the large statistical EISCAT data set used in this study, the neutral wind analysis was not applied, since that would have seriously diminished the data coverage.Thus, Joule heating in EISCAT is calculated through Equation 3 by using only measurements of the electric field, which are performed in a coordinate system that is fixed to the earth, whereas Joule heating in TIE-GCM is calculated through Equation 1, which is in the reference frame of the neutrals.If the neutral wind is not known or observed, then it is a source of uncertainty or error in the Joule heating estimates, which becomes significant during geomagnetic storms (Mannucci et al., 2022).
In order to assess the extent to which the discrepancies between Joule heating in TIE-GCM and EISCAT depend on the effect of the neutral winds on a long-term average, we estimated the Joule heating in TIE-GCM after subtracting the effect of neutral winds.This was done according to Equation 12 of T. E. Sarris et al. (2023), where, by expanding Equation 1for Joule heating, we obtain: In the above equation, the term marked as q c is commonly referred to as convection heating (e.g., Billett et al., 2018;Lu et al., 1995), and corresponds to the Joule heating rate in the absence of neutral winds, whereas the second term, marked as q w , allows for the quantification of the contribution of the neutral winds to Joule heating.
In order to better assess the Kp levels and the areas in local time and altitude where the effects of neutral winds are more evident, in Figure 3 we plot with a solid blue line the median altitude profiles of Joule heating from TIE-GCM with the effects of the neutral winds (q c + q w ) and with a dashed blue line Joule heating in TIE-GCM without including the effect of the neutral winds (q c only).In this plot, the visual comparison between the median altitude profiles is further quantified arithmetically with the calculation of the percentage difference between the two height-integrated Joule heating estimates that correspond to the areas under the solid and dashed blue lines, respectively.The percentage difference is defined here as the height-integrated value of Joule heating with the effects of neutral winds minus Joule heating without the effects of neutral winds, divided by Joule heating with the effects of neutral winds.From Figure 3 it can be seen that Joule heating in TIE-GCM without taking into account the effect of the neutral winds is, in most cases, higher than Joule heating when including the effect of neutral winds, except for the morning and noon sections for low Kp values, where Joule heating has the lowest values among all MLT and Kp bins.The effects of the neutral winds maximize in the afternoon sector (MLT 15-21) for all Kp levels, with an average percentage difference of 26.3% in these sectors; whereas the percentage difference when taking into account all local time sectors and all Kp ranges is on average 16.2%.
From the above estimates and the calculated percentage differences, it appears that the effect of the neutral winds cannot account for the differences between TIE-GCM and EISCAT, shown in the percentages that are listed in the upper right corner of each panel.In comparison, Thayer (1998) used the Sondrestrom ISR to sample currents, conductivities, electric fields, and neutral winds in the E region, and evaluated height-integrated E region Joule heating rates, investigating in particular the influence of the neutral wind on these estimates.They found that the E region height-integrated Joule heating rate for the particular time period they investigated (which corresponded to solar minimum, daytime conditions with periods of moderate to strong geomagnetic activity) was 40% lower when neutral wind estimates were included compared to a fixed neutral atmosphere without winds.This is higher than the percentage difference found above, but within order of magnitude.Our result of lower heating rate with neutral winds compared to zero neutral winds is consistent with ion drag causing a partial neutral wind acceleration into the ion drift's direction.As a consequence, this would reduce the Joule heating rate estimate compared to a fixed neutral atmosphere without winds.Furthermore, Thayer (1998) found that, whereas neutral winds reduce the local Joule heating rate in the upper E region, they enhance the local Joule heating rate in the lower E region.
Looking in closer detail at the altitude profiles in Figure 3, we can see that a similar behavior is observed in the altitude profiles as obtained from TIE-GCM: for example, the dashed line (no neutral wind) in the MLT 15-21 sector for the Kp 4-9 range exceeds the values of the solid line (with neutral wind) at altitudes above 120 km, but has lower values than the solid line at altitudes below 120 km.Sangalli et al. (2009) used sounding rocket data and compared Joule heating rates with observed and omitted winds.For the event that they studied they found similar differences between the two cases as those presented herein.Aikio et al. (2012), using EISCAT data, found that the wind effect was to decrease Joule heating rates in the evening sector, but to increase them in the morning sector.For medium to high activity levels, the positive contribution of winds during the morning maximum was about 20%-30%.During geomagnetically quiet times, the main contributions to Joule heating were typically from the neutral winds.
It is noted that, if the true winds in the calculations presented above are more variable, then the impact on EISCAT estimates could be considerably higher.As noted by, for example, Lin et al. (2005) and Lu et al. (2008), neutral winds can be a significant ionospheric driver, and their variability can indeed have impacts on the EISCAT altitude profiles, especially when there are few measurements, as is the case for high Kp values.Furthermore, given that the EISCAT data are not uniformly sampled, there could be seasonal effects that are not being considered.As summarized in Dhadly et al. (2023), there are significant gaps in the quantification and understanding of neutral winds in the thermosphere, mostly due to a lack of measurements.They also highlight that, whereas there are considerable advances in physics-based and empirical models, these advances are still to be tested, especially at high latitudes, where the statistical behaviors diverge significantly from simulations (e.g., Liuzzo et al., 2015;Wu et al., 2012Wu et al., , 2015)).New measurements could conclusively address the issue of the effect of neutral winds in the LTI.
Another source of uncertainty in the calculations of Joule heating involves the spatio-temporal resolution of the estimation methodology.It is well known that small scale variations can contribute significantly to the total Joule heating.This is primarily due to the fact that variations of the electric field around a mean do not average to zero, as the square of the electric field contributes to Joule heating in Equation 1.For example, electric fields below 10 km are known to be significant at times, with variability at these scales often exceeding the average.In more detail, Codrescu et al. (1995) used theoretical calculations and found that E-field variability could change the distribution of Joule heating significantly, and can introduce inter-hemispheric asymmetries.Emery et al. (1999) multiply Joule heating by 2.5 in order to parameterize the effects of small scale structure in the northern hemispheres.Matsuo and Richmond (2008) incorporated characteristics of residual electric fields derived from observations into the forcing of the TIE-GCM and found that the additional Joule heating due to the subgrid-scale electric field becomes relatively more important during geomagnetically quiet periods.Deng et al. (2009) coupled TIE-GCM with a quantitative empirical model of the high-latitude forcing of the thermosphere based on observations from Dynamic Explorer 2 and reported that including electric field variability into the energy calculation increases the Joule heating by more than 100%; Cosgrove and Codrescu (2009) compared Joule heating as measured from Sondrestrom ISR and as modeled by the Assimilative Mapping of Ionospheric Electrodynamics procedure and found that small-spatial-scale electric field variability can possibly be responsible for the large disagreement, as only the radar is sensitive to 10-km-scale variations of the electric field.Q. Zhu et al. (2018) used GITM simulations and found that Joule heating can be significantly elevated by the small-scale and mesoscale electric field variability (∼27% globally).In the simulations presented herein, the Weimer model has been used for the electric field specification.It is noted that the resulting Joule heating is dependent on the prescribed electric field model, as has been investigated in great detail in a recent study by Pirnaris et al. (2023).
The sensitivity of the results presented herein to the selection of the electric field model is an important topic that needs to be investigated in further detail.
Since Joule heating in TIE-GCM is calculated in a grid of 2.5° in both latitude and longitude, the corresponding spatial resolution is on the order of ∼280 km in the meridional direction and ∼80 km in the azimuthal direction at auroral latitudes; it is thus expected that TIE-GCM cannot resolve small-scale variability.In comparison, the EISCAT beam-swing method utilized in this study for electric field estimates results in horizontal spatial resolutions on the order of several tens of km.It is noted that Kavanagh et al. (2022) compared ion velocity estimates from the technique used in that paper with the EISCAT tristatic method, showing that even small scale changes impacted both the variability and the Joule heating estimates (strongly dependent on local time and geomagnetic activity) by a significant amount.They also found that averaging the tristatic measurement to match the temporal resolution of the beam swing method did not account for the variability.
In order for TIE-GCM to obtain neutral temperatures that are in agreement with statistical observations, such as are obtained, for example, via NRLMSISE-00 (Picone et al., 2002), TIE-GCM in the estimate of Joule heating includes an empirically-derived multiplication factor, termed JOULEFAC, which increases Joule heating by a fixed factor of 1.5 by default, to account for sub-grid-scale and related effects (see, e.g., NCAR, 2016).As noted in, for example, Release 1.9 of TIE-GCM, the value of JOULEFAC is a matter of debate, and its value has been replaced/adjusted in several studies to other fixed levels; for example, Emery et al. (1999) needed to multiply the calculated Joule heating by 2.5 in the northern hemisphere and by 1.5 in the southern hemisphere during the 2-11 November 1993 storm event in order to reproduce observed thermospheric responses.The results from the comparison presented above indicate that a Kp-dependent, variable Joule heating factor, JOULEFAC(Kp) that results in lower Joule heating for lower levels of geomagnetic activity and higher Joule heating for higher levels of geomagnetic activity would lead to a closer agreement between TIE-GCM and EISCAT.
In the retrieval of the height-integrated Pedersen conductivity (conductance) and Joule heating from EISCAT measurements, knowledge of the neutral density and neutral temperature is required.In this study, neutral density and temperature are obtained from the NRLMSISE-00 empirical atmosphere model (Picone et al., 2002).In order to investigate the sensitivity of the resulting conductivity on the inputs from the empirical atmosphere, the neutral temperature and density values from NRLMSISE-00 were varied by a factor of two and conductivity was re-calculated.It was found that the sensitivity of the resulting height-integrated conductivity on neutral temperature was very small, whereas an increase in neutral densities by a factor of 2 increased the height-integrated conductivity by ∼15% on average and between the 10th and 90th percentiles the change varied from 4% to 25%.
In the results shown above it is noted that the selection of the Kp ranges and the selection of the Kp index itself can affect the comparison between TIE-GCM and EISCAT and the definition of the variable Joule heating factor, JOULEFAC(Kp), as described above.This is because, whereas the geomagnetic Kp index is extensively used as an indicator of geomagnetic activity, both for scientific and operational purposes (see, e.g., Matzka et al., 2021), it also has limitations.These limitations are due in particular to its 3-hr cadence, which means that rapidly changing processes, such as those on substorm scales, can be missed.To this direction, subsequent studies that are based on other indices can yield the sensitivity of these results to the selected index.
Another limitation in the comparisons that were performed is related to the fact that EISCAT is generally not sampling uniformly in time, either through each year or through the solar cycle, which means that there may be some biases introduced into the statistics that have not been considered.These effects could be quantified by restricting the TIE-GCM data set that is used to obtain the altitude profiles of Figures 1-3 to only those times when EISCAT observations are available and investigate the sensitivity of the results to the selected times.This sensitivity analysis has not been performed herein, and is the topic of a future investigation.
Finally, it is noted that the temporal scales over which Joule heating is resolved in TIE-GCM and EISCAT could also contribute to inaccuracies in Joule heating and Pedersen conductivity and resulting discrepancies between the two estimation methodologies.For example, as discussed in Rodger et al. (2001), a low time resolution in EISCAT will generally lead to an underestimation of Joule heating.In this study, the time resolution for Pedersen conductivity is ∼1 min and for the electric field it is ∼6 min, leading to a time resolution for Joule heating of ∼1 min.This is within the range of temporal resolutions used in past studies; for example, Aikio and Selkälä (2009) used a resolution of ∼2 min for both the electric field and conductivity in the tri-static EISCAT data analysis.The effect of different time resolutions in TIE-GCM and EISCAT on the resulting Joule heating needs to be further explored through parametric studies, and is beyond the scope of this study.

Summary and Conclusions
In this work, Joule heating rates and Pedersen conductivity were calculated statistically as a function of altitude, MLT and geomagnetic activity (Kp), using a long-term (11 year) simulation of the LTI based on TIE-GCM and a corresponding period of EISCAT radar observations.Model and radar data were obtained from 2009 until 2019 for both methodologies.Through this comparison, it was found that both the TIE-GCM run and EISCAT estimates agree on the shape of the altitude distribution.With respect to the approximate altitude of the peak in Joule heating and Pedersen Conductivity distributions, the maximum value of Joule heating and Pedersen conductivity are found at ∼120 km of altitude in both TIE-GCM run and EISCAT.
By subtracting the effect of neutral winds from the estimates of Joule heating in TIE-GCM, it is found that neutral winds can account, on average, for ∼16.2% of Joule heating in TIE-GCM.An MLT dependence is found in the effect of the neutral winds, with strongest effects appearing in the afternoon (∼26.3%),followed by the other sectors (∼12%).
It is concluded that this difference cannot account for the discrepancies in Joule heating between TIE-GCM and EISCAT, and thus, that the largest effect leading to the discrepancies between the Joule heating estimations in TIE-GCM and EISCAT are not due to the lack of neutral wind velocities when estimating Joule heating in EISCAT.Instead, these discrepancies should probably be attributed to small-scale effects that amount to sub-grid variability that cannot be resolved in TIE-GCM, and which are currently parameterized/ adjusted by a multiplication factor, or to the lack of small temporal scales in TIE-GCM.It was also found that TIE-GCM tends to under-estimate both Joule heating and Pedersen conductivity for high levels of geomagnetic activity and, inversely, to over-estimate their values for low levels of geomagnetic activity, compared to EISCAT measurements.This suggests that replacing the constant multiplication factor, JOULEFAC, that is currently used in TIE-GCM by a Kp-dependent factor would yield better results for both low and high activity levels.
The results presented herein underline the need for more accurate knowledge of the contribution of all relevant spatial and temporal scales, which is paramount to accurately parameterize Joule heating in models such as TIE-GCM.Revealing and accurately resolving the contribution of all scales of interest to Joule heating, while including the contributions of all parameters involved, is best achieved by in situ measurements of the lower thermosphere-ionosphere, as it has been emphasized in many recent studies (e.g., Heelis and Maute, 2020;Palmroth et al., 2021;T. E. Sarris et al., 2020;T. Sarris et al., 2023).Further to the improvements in the unambiguous characterization of Joule heating that will be enabled through comprehensive in situ observations in the LTI, improvement of models will also emerge from the upcoming developments in radars: for example, the enhanced spatial and temporal resolution that will be offered by EISCAT_3D, currently under development, will greatly improve estimates of Joule heating in small scales (Stamm et al., 2021).In addition, the planned scanning Fabry-Perot Interferometer facility that will be implemented around EISCAT_3D will be a key step in incorporating the effects of neutral winds in the resulting Joule heating.The combination of EISCAT_3D measurements with in situ measurements from a spacecraft mission quantifying the variability of the electric fields and all involved parameters in small scales, such as described in, for example, T. Sarris et al. (2023), will enable resolving conclusively the key missing pieces of Joule heating and Pedersen conductivity in the LTI.

Data Availability Statement
EISCAT data availability: The EISCAT data are available in the Madrigal database: https://madrigal.eiscat.se/madrigal/ and the corresponding data of the altitude profiles can be found at: https://zenodo.org/doi/10.5281/zenodo.10017765.
Data processing source code: The processing of TIE-GCM output data and EISCAT data was done using the code that can be found at: https://github.com/baldimitris/TIEGCM_Statistics(Baloukidis et al., 2023).

Figure 1 .
Figure 1.Altitude profiles of Joule heating in Thermosphere-Ionosphere-Electrodynamics General Circulation Model (TIE-GCM) (top panel) and European Incoherent Scatter Scientific Association (EISCAT) (bottom panel), in units of 10 −8 W/m 3 .The curves represent from left to right the 10th, 25th, 50th, 75th, and 90th percentile of Joule heating and their area is colored with light, medium, dark, medium and light shade respectively.The central bolder curve denotes the median value.Marked in blue (green) bold numbers are the height integrated values of the median of Joule heating in TIE-GCM (EISCAT), in units of mW/m 2 ; marked in smaller black characters are the number N of grid points (top panel) and samples (bottom panel) that were used to produce the altitude profiles.

Figure 2 .
Figure 2. Altitude profiles of Pedersen Conductivity in Thermosphere-Ionosphere-Electrodynamics General Circulation Model (TIE-GCM) (top panel) and European Incoherent Scatter Scientific Association (EISCAT) (bottom panel), in units of mS/m.The curves represent from left to right the 10th, 25th, 50th, 75th, and 90th percentile of Pedersen conductivity and their area is colored with light, medium, dark, medium and light shade respectively.The central bolder curve denotes the median value.Marked in blue (green) bold numbers are the height integrated values of the median of Pedersen conductivity in TIE-GCM (EISCAT), in units of Siemens; marked in smaller black characters are the number N of grid points (top panel) and samples (bottom panel) that were used to produce the altitude profiles.

Figure 3 .
Figure 3.Comparison of the median of Joule heating rates estimated with European Incoherent Scatter Scientific Association (green line), with Thermosphere-Ionosphere-Electrodynamics General Circulation Model (TIE-GCM) when the effects of the neutral winds are included (blue solid line) and with TIE-GCM when the effects of the neutral winds are not included (blue dashed line).The number at each sub-figure marks the percentage difference between the height-integrated value of Joule heating in TIE-GCM with the effects of neutral winds and without the effects of neutral winds, divided by Joule heating with the effects of neutral winds.

,
+ and r e are the collision frequency to gyrofrequency ratios of O + ,   + 2 , NO + and e respectively, as obtained using the collision frequencies ofSchunk and Nagy (2009), and

,
+and N e are the number densities of these species, in units of m −3 .Collision frequencies of the above species are calculated for collisions with the following neutral species: O, O 2 , N 2 .It is noted that TIE-GCM considers only resonant interactions (i.e., interactions between an ion species and its parental neutral) for the O 2 −