Assessment of ionospheric Joule heating by GUMICS-4 MHD simulation, AMIE, and satellite-based statistics: towards a synthesis

We investigate the Northern Hemisphere Joule heating from several observational and computational sources with the purpose of calibrating a previously identified functional dependence between solar wind parameters and ionospheric total energy consumption computed from a global magnetohydrodynamic (MHD) simulation (Grand Unified Magnetosphere Ionosphere Coupling Simulation, GUMICS-4). In this paper, the calibration focuses on determining the amount and temporal characteristics of Northern Hemisphere Joule heating. Joule heating during a substorm is estimated from global observations, including electric fields provided by Super Dual Auroral Network (SuperDARN) and Pedersen conductances given by the ultraviolet (UV) and X-ray imagers on board the Polar satellite. Furthermore, Joule heating is assessed from several activity index proxies, large statistical surveys, assimilative data methods (AMIE), and the global MHD simulation GUMICS-4. We show that the temporal and spatial variation of the Joule heating computed from the GUMICS-4 simulation is consistent with observational and statistical methods. However, the different observational methods do not give a consistent estimate for the magnitude of the global Joule heating. We suggest that multiplying the GUMICS-4 total Joule heating by a factor of 10 approximates the observed Joule heating reasonably well. The lesser amount of Joule heating in GUMICS-4 is essentially caused by weaker Region 2 currents and polar cap potentials. We also show by theoretical arguments that multiplying independent measurements of averaged electric fields and Pedersen conductances yields an overestimation of Joule heating. Correspondence to: M. Palmroth (Minna.Palmroth@fmi.fi)


Introduction
Joule heating, calculated as the scalar product of the current and electric field, is a term used to describe the Ohmic production of heat that occurs as the charged particles drifting in the direction of the electric field collide with the neutral particles of the resistive medium. In the ionosphere, the net charge due to the electron precipitation from the magnetosphere creates an electric field that pulls ions. Thus, the ions close the field-aligned currents (FACs) horizontally and in the process they undergo collisions with atmospheric neutral particles, which are usually in motion due to thermospheric winds. The electric field in the frame of the neutrals is E U =E+U×B, where E is electric field imposed on the ionosphere, U is the velocity of the thermospheric wind and B is the magnetic field. Therefore, the ionospheric Joule heating P J H is computed as where J is the electric current.
Since it is difficult to obtain global measurements of the neutral winds, the Joule heating estimates are carried out typically by assuming U=0 and measuring the electric field and the height-integrated Pedersen conductivity P , from which the J·E term in Eq. (1) can be calculated as P E 2 . P and E can be obtained either by using satellites (e.g., Foster et al., 1983;Heelis and Coley, 1988), radars (e.g. Fujii et al., 1999), or methods based on ground magnetometer and radar data (e.g. Ahn et al., 1983). Joule heating can also be assessed by field-aligned Poynting flux measurements (e.g. Gary et al., 1994Waters et al., 2004). For many purposes this is advantageous as the Poynting flux is the total electromagnetic energy that includes both the energy dissipation due to Joule heating and the mechanical energy consumed by the neutral winds. This can be seen in the Poynting theorem (µ −1 0 ∇·(E×B)=−J·E).
To date, the most comprehensive statistical study on Joule heating was published by Foster et al. (1983), who used data from about 25,000 Atmospheric Explorer (AE)-C satellite passes over the ionosphere and binned the data according to season and magnetic activity. They showed that the largest amount of Joule heating occurs at the oval near dawn and dusk. This was later confirmed with different measuring methods (e.g. Gary et al., 1995;Fujii et al., 1999;Waters et al., 2004;Olsson et al., 2004), indicating that the closure of Region 1 and Region 2 currents to each other plays a major role in the Joule heat production rate. Furthermore, both Foster et al. (1983) and Gary et al. (1995) emphasized the role of the dayside in producing Joule heat, since there the conductances are never small due to ionization caused by solar illumination. In fact, Foster et al. (1983) described that statistically the Joule heating pattern represents a horseshoe opening towards a minimum in the nightside.
Over the course of a few decades, the importance of Joule heating as a major sink of magnetospheric energy has become evident. Early studies of energy coupling between the solar wind and the magnetosphere and ionosphere (Akasofu, 1981) assumed that energy related to Joule heating is only a fraction of the ring current energy injection rate. When measurements of ionospheric properties yielded Joule heating estimates, as characterized by AE (e.g. Ahn et al., 1983) or K p (Foster et al., 1983) indices, it was realized that the ionospheric Joule heating has a more pronounced role in consuming energy associated with the solar wind-magnetosphere coupling. Presently, it is estimated that 50-60% of the estimated input energy is consumed in the Joule heating processes during isolated and storm-time substorms on both hemispheres (Østgaard and Tanskanen, 2003), whereas over 50% of total stormtime energy consumed by the ionosphere and ring current is consumed by Joule heating (e.g. Knipp et al., 1998;Pulkkinen et al., 2002). Estimates of the amount of Joule heating in the ionosphere are important, for example, for the space weather community, as the expansion of the atmosphere due to Joule heating increases satellite drag at low-altitude orbits. Furthermore, the amount of ionospheric Joule heating affects global thermospheric dynamics, and hence the quantitative assessment of high-latitude global Joule heating is of crucial importance to ionospheric physics, ranging from low to high latitudes.
Several attempts have been made to incorporate global models in the effort to estimate the amount of ionospheric Joule heating. Thayer et al. (1995) investigated the relationship between the electromagnetic energy, Joule heating, and the mechanical energy, using a model based on the National Center for Atmospheric Research's thermosphere-ionosphere general circulation model.  and Knipp et al. (1998) have estimated the Joule heating using the assimilative mapping of the ionospheric electrodynamics (AMIE) procedure (Richmond and Kamide, 1988;Richmond, 1992). Slinker et al. (1999) computed the Joule heating using a global MHD simulation and compared the result with the AMIE output. Palmroth et al. (2004a), also using a global MHD simulation, computed the amount of Joule heating along with the energy associated with the particle precipitation in the ionosphere. However, there are several difficulties in computing ionospheric Joule heating rates using global models. The amount of Joule heating depends quadratically on the electric field, and is thus highly sensitive to the accuracy of the polar cap potential structure. The limited amount of ionospheric processes implemented in the ionospheric descriptions of the global MHD simulations naturally affects the ability of the simulation to assess the total Joule heating rate. Although attempts have been made to include thermospheric physics to global MHD simulations (Raeder et al., 2001;Ridley et al., 2003;Wiltberger et al., 2004), the neutral winds are not usually considered in MHD simulations, and therefore the Joule heating associated with the ionospheric circulation set up by the neutral winds is not modeled. Along with the neutral winds, the discrete arc physics may not be properly modeled with global MHD simulations, as the scale size of discrete arcs is smaller than the spatial resolution of the ionospheric simulation, and hence the Joule heating associated with discrete arcs may not be correctly taken into account in the total Joule heating estimates. Palmroth et al. (2004a) used the GUMICS-4 MHD simulation  to estimate the total storm-time and substorm-time ionospheric energy consumption including the two largest ionospheric energy sinks: the Joule heating and particle precipitation. The temporal variation of Joule heating computed from the GUMICS-4 simulation compared well with an AE-based Joule heating proxy (Ahn et al., 1983), although the amount of energy in the GUMICS-4 ionosphere appeared to be less than that predicted by the Ahn et al. (1983) proxy; similar conclusions applied to the energy associated with the precipitating particles. Later, consistent results were obtained in simulations of other events (Palmroth et al., 2004b). Palmroth et al. (2004a, b) gave a mathematical expression describing the GUMICS-4 ionospheric power consumption, where P ionosphere is the summed power associated with the Joule heating and precipitation in the GUMICS-4 simulation, and ρ is solar wind density, v is velocity, B z,I MF is the IMF z component, and p dyn is the solar wind dynamic pressure. Scaling parameters ρ 0 =m p ·7.3·10 6 m −3 =1.22·10 −20 kgm −3 and v 0 =400 km/s are chosen as typical solar wind density and velocity, to obtain units of Watts for the constant C.
Since the components of P ionosphere are both independently in a temporal agreement with AE-proxies, and, on the other hand, P ionosphere correlates with the right-hand side of Eq. (2) with more than 80% of the correlation coefficients in all the simulated events, the right-hand side of Eq.
(2) can roughly predict the temporal behavior of ionospheric power dissipation as determined by the AE-proxies. It was also argued that if one scaled the GUMICS-computed total ionospheric power consumption to correspond with observational values, i.e. if Eq.
(2) was "calibrated" by increasing C, it could predict the ionospheric power consumption correctly, both temporally and magnitudewise, then the power law could be used even for space weather forecasting. Notice that Eq. (2) describes a directly-driven system, and therefore a time delay must be taken into account before it can be applied to real events. The purpose of this paper is to give a realistic scaling for Eq.
(2) as far as Joule heating is concerned; the calibration associated with the particle precipitation energy is the subject of further study. The calibration is broken downinto two aspects: 1) estimating the difference of the total Joule heating power in the GUMICS-4 ionosphere as compared to observational methods, and 2) determining if this difference stays constant in time. Hence, we aim to compare both the magnitude and the temporal variation of the total GUMICS-4 Joule heating with several observational methods. Furthermore, it is clear that any instantaneous model, such as the GUMICS-4 global MHD simulation, giving the Joule heating rate should be consistent with large statistical surveys. Therefore, in one case study, we compare the Joule heating in the GUMICS-4 global MHD simulation with the results of large statistical surveys (Foster et al., 1983;Olsson et al., 2004), commonly used proxies by Ahn et al. (1983), the AMIE procedure, and data from SuperDARN radars and and the UV and X-ray imagers on board the Polar satellite. Furthermore, we discuss the most important factors that affect the accuracy of the Joule heat estimates in the various methods.

Solar wind observations
The Wind spacecraft, located at (230, −22, −6) R E in geocentric solar ecliptic (GSE) coordinate system, recorded a southward turning of the interplanetary magnetic field (IMF) at 22:27 UT after several hours of northward orientation (Fig. 1a). The IMF B z rotated relatively smoothly from north to south, reaching about −8 nT at ∼23:00 UT, while the northward orientation was attained at 00:27 UT. Another southward turning at 02:15 UT was followed by a northward turning an hour later, at 03:14 UT. The IMF B y (Fig. 1b) fluctuated between 0 and −10 nT during the two southward IMF periods. The solar wind velocity (Fig. 1c) was between 450 and 500 km/s throughout the time interval of interest, whereas the solar wind density (Fig. 1d) remained in between 3 and 5 cm −3 . The solar wind dynamic pressure (Fig. 1e) fluctuated between 1 and 2 nPa. The bottom panel shows   the epsilon ( ) parameter (Akasofu, 1981) which is used as a proxy of the energy input to the magnetosphere in the geocentric solar magnetospheric (GSM) coordinate system. Concurrent with periods of southward IMF, the epsilon parameter shows two intervals of enhanced energy input, both of which exceed the energy input level that typically leads to a magnetospheric substorm (100 GW, Akasofu, 1981). The time delay between the Wind X GSE position and the magnetopause is ∼49 min using the average solar wind velocity in the X GSE direction (∼475 km/s).

Substorm evolution
The Nuuk (GHB) magnetometer station, located at a magnetic latitude of 70.5 • at the west coast of Greenland , recorded a negative bay in the H component at 23:45 UT (Fig. 2a), 29 min after the southward IMF had arrived at the subsolar magnetopause (Fig. 1). This corresponds to a magnetic local time (MLT) at GHB of ∼21:15. A second intensification occurred at 00:12 UT. Both deviations showed signatures of northward propagation. Similar characteristics were also observed by other magnetometer stations on the west coast of Greenland, as presented in Fig. 2a Figure 2b shows the AE index obtained from the WDC-2 in Kyoto. As characterized by the AE index, the first intensification reached almost 1000 nT, whereas the second was ∼500 nT in magnitude. Throughout the night after 00:00 UT the Kilpisjärvi (KIL) all-sky camera recorded auroral forms (not shown). Since the magnetic local time of KIL is dawnward of Greenland and the main activity location, a clear onset time of the auroral breakup could not be identified from the ground optical data.
The Los Alamos geostationary spacecraft 1990-095 traversed the pre-midnight region (∼21:00 MLT) at the time of the first activation onset. Figure 3 shows the electron data from the spacecraft 1990-095. The electron fluxes started to decrease at 23:20 UT, which coincides with the time of the southward IMF arrival at the magnetopause. At 23:48 UT a sharp increase in the electron fluxes was recorded, consistent with the timing of the negative bay onset in the ground magnetometer data (23:45 UT). The geostationary electron data support the global picture of a substorm sequence: At the time of the arrival of the southward IMF the tail field started to stretch in the nightside, as evidenced by the decrease in the electron flux at geostationary orbit near midnight. The injection marked the onset of the substorm expansion phase, consistent with the ground magnetometer recordings.
At 03:30 UT, the spacecraft 1990-095, now traversing near the midnight meridian (∼01:00 MLT), again recorded a decrease in the electron fluxes. At 03:57 UT, another sharp increase in the electron fluxes occurred. The data presented suggest that the 04:00 UT event is consistent with the picture of a magnetospheric substorm, as evidenced by a flux decrease that indicates tail stretching, and geostationary injection and northward propagation of the negative bay on the ground that indicates dipolarization of the tail magnetic field. Figure 4 presents the Northern Hemisphere instantaneous convection patterns as measured by the SuperDARN radars (Greenwald et al., 1985). The convection patterns are shown at the onset (23:50 UT) and maximum of the first substorm (00:40 UT), during the activity minimum between the substorms (03:00 UT), and close to the maximum of the second substorm (04:20 UT). The bottom panel of Fig. 4 gives the polar cap potential determined from SuperDARN measurements. At 23:50 UT the dusk convection cell is well-defined by the SuperDARN radar measurements, whereas the dawn cell has insufficient data coverage. The polar cap potential inferred from the radar measurements is 72 kV. At the maximum of the first substorm, 00:40 UT, both dusk and dawn cells have insufficient data coverage, although some convection features are fairly well-defined by measurements. Therefore, the polar cap potential (61 kV) at 00:40 UT is determined by a statistical model (e.g. Ruohoniemi and Baker, 1998 and references therein), which tends to underestimate the instantaneous potential (M. Ruohoniemi, private communication, 2004). At the minimum between the substorms (03:00 UT), and near the maximum of the second substorm (04:20 UT) the SuperDARN radars define the convection pattern again fairly well. The polar cap potentials were 42 kV and 53 kV, respectively. In conclusion, except for the maximum of the first substorm (00:40 UT), the convection pattern and the polar cap potential determined from the Super-DARN radars are quite reliable during the time period shown in Fig. 4. Figure 5 presents the Pedersen conductances for the same time periods as given in Fig. 4, i.e. at the onset and maximum of the first substorm, in between the two substorms, and near the maximum of the second substorm, respectively. The conductances, including the solar contribution and the precipitation component, are derived using a method described by Aksnes et al. (2002), which utilizes data from the UV and X-ray imagers on board the Polar satellite. At the onset of the first substorm ( Fig. 5a) the pre-midnight sector shows enhanced Pedersen conductance, verifying that the onset timing determined from the ground magnetometer data (Fig. 2) is correct, because the GHB station is located below the enhanced conductance region. At maximum, the Pedersen conductance, using the Aksnes et al. (2002) method, is more than 15 S, while on average the conductance is 4-10 S along the oval. Figure 5b indicates that at the maximum of the first substorm the enhanced conductance extends over the nightside, while the average Pedersen conductance is ∼10 S. In between the two substorms ( Fig. 5c) the nightside conductance has decreased to ∼4-6 S. The second substorm took place in the post-midnight sector, as indicated by the locations of the Pedersen conductance maxima (Fig. 5d). This confirms the timing of the second substorm using the ground magnetometers located in the dawn sector ( Fig. 2). At maximum activity, the Pedersen conductance is over 15 S, although variable between 6-10 S along the nightside oval. Figure 6 shows data from the onset and maximum of the first substorm (23:50 and 00:40 UT, first two rows), minimum in between the substorms (03:00 UT, third row), and near the maximum of the second substorm (04:20 UT, last row). In each plot, 12:00 MLT is up, 24:00 MLT is down and 18:00 (06:00) MLT is to the left (right); white grid lines indicate the magnetic latitude. The first column presents the global electric fields as measured by the SuperDARN radars. These electric fields may include large errors during insufficient radar data coverage. Furthermore, as the potentials are fitted to the radar data using a spherical harmonic fitting procedure, a residual potential sometimes appears equatorward of the zero potential line. To avoid overestimating the Joule heating due to the residual electric field, we have forced the electric field to be zero equatorward of the zero potential line. The second column gives maps of Joule heating, computed by taking the electric field from the SuperDARN radars and the Pedersen conductance from the Polar measurements ( Fig. 5), which are first interpolated to have the spatial resolution of the SuperDARN electric fields. Hereafter, this Joule heating is termed "the observed Joule heating". The last column gives the instantaneous Joule heating as computed from the AMIE procedure. The polar cap potential is overlaid on the AMIE Joule heating maps as white contours. At the onset of the first substorm, at 23:50 UT, the twocell potential pattern is developed and the polar cap potential difference is 69 kV in AMIE, while SuperDARN radars measure 72 kV (Fig. 4). While the AMIE results indicate that the maximum Joule heating appears in the nightside, the observed Joule heating is enhanced wherever the SuperDARN electric fields show large values, particularly in the dusk sector. While the maximum of the observed Joule heating is 25 mWm −2 at the onset, on average, the observed Joule heating is barely over 2 mWm −2 throughout the oval and is over 10 mWm −2 only in very limited regions. In AMIE, the maximum Joule heating at the onset is 38 mWm −2 , while dusk and dawn sectors produce ∼5-15 mWm −2 . Apart from the nightside maximum, AMIE produces Joule heating also at the equatorward edge of the convection cells, indicating that the closure of Region 1 and Region 2 field-aligned currents to each other produces significant Joule heating. This is consistent with a number of other research results (e.g. Foster et al., 1983;Gary et al., 1994;Olsson et al., 2004;Palmroth et al., 2004c). At the maximum of the first substorm, at 00:40 UT, the AMIE polar cap potential is 132 kV, while SuperDARN suggests 61 kV. The latter is probably an underestimation due to insufficient data coverage (see Fig. 4b). Statistically, for K p 4.7 at 00:40 UT the polar cap potential is ∼90 kV, as suggested by the empirical polar cap potential proxy of Boyle et al. (1997). The SuperDARN and Polar measurements show a maximum at ∼20:00 MLT (17 mWm −2 ). However, on average the observed Joule heating is ∼2 mWm −2 throughout the Northern Hemisphere. On the other hand, the AMIE results show that the Joule heating is largely concentrated on the dawnside Region 1 and Region 2 current closure region, although some Joule heating also appears on the nightside. The AMIE maximum is 106 mWm −2 .

Spatial variation
In between the two substorms, at 03:00 UT, the AMIE polar cap potential has decreased to 27 kV, while SuperDARN radars suggest 42 kV. The observed Joule heating is now 10 mWm −2 at maximum (dayside poleward edge), while on average the Joule heating is ∼1 mWm −2 throughout the oval. The observed Joule heating is enhanced on the dusk oval. On the other hand, the AMIE procedure yields virtually no Joule heating (maximum ∼5 mWm −2 ). Near the maximum of the second substorm, at 04:20 UT, the polar cap potential is 54 kV in AMIE, while SuperDARN radars obtain 53 kV. Generally, the spatial variation of the observed Joule heating agrees with the AMIE results, as both show enhanced Joule heating at the vicinity of 70 • on virtually all MLTs. Additionally, the observed Joule heating shows a maximum at ∼20:00 MLT. The maximum of the observed Joule heating is now 11 mWm −2 , but on the average ∼ 1 mWm −2 . The AMIE maximum during that period is 12 mWm −2 .
3.2 Total integrated Joule heating Figure 7 presents the total Joule heating integrated over the Northern Hemisphere computed using various methods during the 28-29 March 1998 event. In Fig. 7a, Joule heating is computed by taking the global electric field from Super-DARN radars (Fig. 6), while the Pedersen conductance is taken from the Polar global imagers (Fig. 5). Thus, Fig. 7a presents the integration of the observed Joule heating in Fig. 6; the 10-min spatial resolution is due to the integration of UVI and PIXIE instruments used in deriving the Pedersen conductance. Figure 7b gives estimates of the integrated Joule heating rate based on AE (Ahn et al., 1983) and K p indices (Foster et al., 1983), while Fig. 7c gives the total Joule heating as computed by the AMIE procedure. Ahn et al. (1983) used radar and ground magnetic field measurements to develop global conductance distributions, which were then used to compute the electric field distributions with the method introduced by Kamide et al. (1981). Separate equations were given for the standard AE index (based on 12 magnetometers) and for an AE index derived from a larger number of magnetometers. Figure 7 shows the Ahn et al. (1983) proxy for the standard AE index. Regardless of the poor temporal resolution, the Foster et al. (1983) proxy P J H (GW )=4+20K p at equinoxes gives an estimate of the level of Joule heating based on the most comprehensive statistical study published up to date. The vertical lines correspond to Fig. 4, Fig. 5, and to instantaneous global maps presented in Fig. 6. As indicated in Fig. 4, the SuperDARN radars have insufficient data coverage during 00:30-02:30 UT, and therefore both the temporal variation and the magnitude of the Joule heating in Fig. 7a during this period include large uncertainties. Compared to the AE-proxy in Fig. 7b, the onset of the second substorm occurs slightly earlier (later) in the observed Joule heating (AMIE), indicating that the various methods do not have a one-to-one correspondence in terms of temporal evolution. In terms of magnitude of the total Joule heating, the best correspondence is obtained during the second substorm, when the observed Joule heating, the AE-proxy, and the AMIE results are in quantitative agreement, given that the error estimate of the Ahn et al. (1983) proxy is ±50%. During the first substorm, the Joule heating from the different methods are far from a quantitative agreement: The AMIE results show twice as large values as compared to the AEproxy, which shows twice as large values as compared to the K p -proxy and the observed Joule heating.

GUMICS-global MHD simulation
GUMICS-4 (Janhunen, 1996) is a computer simulation designed specifically for solving the coupled solar wind-magnetosphere-ionosphere system. The code consists of two computational domains: The MHD domain including the solar wind, and the magnetosphere and the electrostatic domain including the ionosphere. The fully conservative MHD equations are solved in a simulation box extending from 32 R E to −224 R E in the X GSE direction and ±64 R E in the Y GSE direction and Z GSE direction. Near the Earth the MHD simulation box reaches a spherical shell with a radius of 3.7 R E , which maps along the dipole field to approximately 60 • in magnetic latitude. The grid in the MHD simulation box is a Cartesian octogrid, and it is adaptive, meaning that whenever the code detects large spatial gradients, the grid is refined. Furthermore, in order to save computation time, the code uses subcycling , in which the time step varies with the local travel time of the fast magnetosonic wave across a grid cell. Solar wind density ρ, temperature T , velocity v and magnetic field B are treated as boundary conditions along the sunward wall of the simulation box; outflow conditions are applied on the other walls of the simulation box.
The MHD magnetosphere is coupled to a high-resolution electrostatic ionosphere. The ionosphere is a spherical shell at an altitude of 110 km, which is mapped to the magnetosphere using dipole field lines. The region between the ionosphere and the 3.7 R E shell is a passive medium, which only transmits electromagnetic effects, and where no currents flow perpendicular to the magnetic field. A triangular finite element grid of the ionosphere is fixed in time, although refined to ∼100 km spacing in the auroral oval region. The ionosphere-magnetosphere feedback loop includes the mapping of the field-aligned currents and the electron precipitation from the magnetosphere to the ionosphere. Precipitation from the magnetosphere is assumed to originate from a Maxwellian source population having the plasma sheet characteristics, and which has a finite probability to fall into the loss cone. The electron precipitation affects the ionospheric electron densities, which are calculated from ionization and recombination rates in 20 nonuniform height levels ranging from 64 km to 194 km. The electron densities are used in the calculation of the Pedersen and Hall conductivities, which are height-integrated to give the conductances at 110 km altitude. The dayside conductances are computed from F10.7 flux according to the empirical formulas of Moen and Brekke (1993), and 0.5 S is assumed to be the background conductance originating from ion precipitation, galactic UV and Xrays, energetic neutral atoms, and UV and X-rays reflected from the Moon. The MSIS (Hedin, 1991) model is used for the thermospheric model, giving for example, the collision frequencies needed in the computation of the conductances. The ionospheric potential is solved by using the field-aligned currents as a source for the horizontal ionospheric currents, which in turn are defined by the height-integrated Ohm's law in the ionosphere. The ionospheric potential is then mapped to the inner shell of the magnetosphere, where it is used as a boundary condition for the MHD equations. More information on the ionospheric computation can be found in Janhunen and Huuskonen (1993).
The 28-29 March 1998 event was simulated using Wind observations as an input to the code. The IMF B x was set to zero to ensure a divergenceless input magnetic field. The smallest grid size in the simulation was set to 0.25 R E . As the neutral winds are not modeled in the simulation, the GUMICS-4 Joule heating P J H G is calculated as where the electric field and the Pedersen conductances are determined in the ionosphere; dS is the area element of the ionosphere and the integration is carried out over the Northern Hemisphere (between 0 • and 90 • in magnetic latitude), although usually the Joule heating in GUMICS-4 is concentrated poleward of ∼50 • . Figure 8 presents global maps of the GUMICS-4 Pedersen conductance, field-aligned currents, and Joule heating at 23:50 UT, 00:40 UT, 03:00 UT, and 04:20 UT; in each plot the orientation is as in Fig. 6, and the latitude grid is indicated as black contours. The onset is characterized as enhanced Pedersen conductance in the 21:00 MLT sector, consistent with the measurements made by the Polar global imagers (cf. Fig. 5). However, the magnitude of the GUMICS-4 Pedersen conductance at the oval (∼5 S) is lower than what is observed (∼12 S, Fig. 5). Region 1 currents are enhanced, whereas the Region 2 currents are weak in the dusk and dawn sectors. The polar cap potential pattern is similar to the observed pattern in Fig. 4, as the convection cells are tilted towards the dawn. Furthermore, the GUMICS-4 Joule heating pattern has similarities with the observed Joule heating in Fig. 6, as the oval ranging from the dusk to postmidnight, and the polar cap in between the convection cells show enhanced heating. The magnitude of GUMICS-4 Joule heating is much lower than the observations show in Fig. 6.
Near the maximum of the first substorm, at 00:40 UT, the GUMICS-4 Pedersen conductance has remained essentially  similar as during the onset, whereas the Region 1 currents have intensified. Joule heating is also enhanced compared to the onset conditions; however, the heating occurs essentially in the same regions (oval and polar cap in between the convection cells) as during the onset. Also, the observed Joule heating in Fig. 6 shows this behavior, as both the oval and the polar cap show enhanced heating. At the minimum between the substorms, at 03:00 UT, the Pedersen conductance, field-aligned currents and Joule heating have all decreased. Although at 03:00 UT the IMF z component is positive, the NBZ current system is not fully developed, possibly because of a simultaneous considerable IMF y component that moves the reconnection regions towards lower latitudes (Kallio and Koskinen, 2000). A fully developed NBZ current system, where a negative field-aligned current is generated in the dawnside high latitudes, is observed later in the GUMICS-4 results as the IMF clock angle is more purely oriented northward. Joule heating, however, continues to show the same morphology as before: the faint maxima of heating occur at the oval and in between the convection cells, consistent with the observed Joule heating in Fig. 6. Near the maximum of the second substorm, at 04:20 UT, the observations show that the Pedersen conductance is enhanced in the post-midnight region (Fig. 5), whereas the GUMICS-4 results show the enhanced Pedersen conductance in the premidnight region. Nevertheless, the field-aligned currents and the Joule heating are again enhanced, and again consistent with the observed Joule heating; the heating occurs at the oval and in between the convection cells where the electric field is increased. Figure 9 presents the total integration of GUMICS-4 Joule heating in the Northern Hemisphere as a function of time. The total integrated Joule heating computed using the Su-perDARN electric fields and Polar measurements of Pedersen conductance (from Fig. 7) are shown with scaling on the right axis. Figure 9 shows clearly that the temporal variation of GUMICS-4 Joule heating is quite similar with the observed Joule heating. Compared to the AE-proxy and the observed Joule heating (Fig. 7), the GUMICS-4 results show a somewhat simultaneous increase and decrease during the first substorm. The second substorm onset in the GUMICS-4 results occurs slightly later than in the observed Joule heating; however, the second onset is simultaneous in the GUMICS-4, AE-proxy and AMIE results, although the GUMICS-4 (AMIE) Joule heating increases faster (slower) than the AE-proxy (Fig. 7). Figure 7 shows that the first substorm has two intensifications; this double-peak property of the total Joule heating is also captured by the GUMICS-4 simulation. The relative amplitude of the GUMICS-4 intensifications is similar to those in the AE-proxy: the Joule heating level during the second substorm is about two-thirds of the first substorm. Figure 9 clearly suggests that the magnitude of GUMICS-4 Joule heating is quite consistently a factor of 10 smaller than the observed Joule heating. This is also in accordance with the magnitude of the Foster et al. (1983) proxy. Table 1 presents the summary of Joule heating estimates from the Ahn et al. (1983) proxy, AMIE, SuperDARN and Polar measurements, and GUMICS-4, respectively. The observed Joule heating at 00:40 UT is most probably underestimated since the statistical method used to derive the convection pattern underestimates the polar cap potential during insufficient data coverage of the SuperDARN radars. The most distinctive feature of Table 1 is that the different Joule heating estimates are relatively closer together towards the end of the substorm period, and that the relative differences among the various estimates during the maximum of the first substorm is large.

Theoretical aspects on estimating Joule heating
In general, conductance and electric field are spatially anticorrelated, such that electric fields tend to circumvent regions of high conductance (e.g. Evans et al., 1977;Mallinckrodt and Carlson, 1985). In the case of perfect spatial anticorrelation of electric fields and conductance, no Joule heating would be produced, since no charge carriers would carry electric current along the electric field. In the ionosphere, however, the global convection electric field and ionization due to precipitation and EUV radiation ensure that generally electric fields and regions of high conductance overlap spatially. Therefore, in the large scale the electric fields and conductance are correlated in the sense that both exist globally in the ionosphere. In smaller scales, the average of the global electric field and conductance may overestimate the true Joule heating in regions where electric fields and conductances are spatially anticorrelated. In fact, it can be shown (see Appendix A) that under the anticorrelation assumption where the right-hand-side of Eq. (4) denotes the true Joule heating in area A i , and the left-hand-side denotes the Joule heating obtained from independent averaged measurements of the electric field E i and conductance σ i . The relation has also been shown to hold observationally, as Foster et al. (1983) showed that excluding small regions in the dayside where electric fields and conductances correlate positively, the left-hand-side of Eq. (4) is typically 0.5-4 mW/m 2 larger than the right-hand-side of Eq. (4). The spatial anticorrelation of the electric field and Pedersen conductance can also be seen in Fig. 10, where GUMICS-4 Joule heating, electric field, and Pedersen conductance are plotted along the terminator at the first substorm onset (23:50 UT). Clearly, Fig. 10 indicates that the peaks of the electric field do not appear at the same location with the peaks of the Pedersen conductance. Furthermore, near the pole, where the electric field is enhanced, the Pedersen conductance is lower on average than equatorward of the oval. To prove that the spatial anticorrelation of the electric field and the Pedersen conductance in the simulation is not only a property of the one time instant plotted in Fig. 10, Fig. 11 presents the electric field and Pedersen conductance data pairs taken at random locations over the Northern Hemisphere poleward of 60 • during 23:30-06:00 UT, such that each data point corresponds to a location and time in the simulation ionosphere. Since the whole data set contains over 277 000 values, we plot only every 100th datapoint. Clearly, Fig. 11 suggests that in the simulation the Pedersen conductance is on average low on locations where the electric field is high, and vice versa. The correlation coefficient is −0.66 in both the plotted subset and the large data set (277 000 data points). Notice that in Fig. 11 the lowest conductance values (0.5) result from the background conductance. Theoretically, the overestimation of Joule heating due to independent measurements of the electric field and conductance may be infinite. For example, element A i may include piecewise defined asynchronous step functions of the electric field and conductance such that, for example, E i varies between 0 and 2 mV/m, while conductance varies between 2 and 0 S. Therefore the averages of the electric field and conductance yield 1 mW/m 2 for Joule heating in A i , although the true Joule heating would be 0 mW/m 2 . Assessment of the overestimation of Joule heating (P corr ) due to spatial anticorrelation of E and P can be computed using the correlation coefficient and standard deviations of E 2 and P (see Appendix A) as where Corr i ∈ [−1, 0] and ST D σ and ST D E 2 are standard deviations of Pedersen conductance and the square of the electric field, respectively, in an area element A i . A global estimate of the correction to Joule heating due to spatial anticorrelation would therefore require measurements of E and P with high spatial resolution over the Northern Hemisphere.
One possible situation where the above arguments (Joule heating is overestimated due to spatial anticorrelation of E and P ) are not valid is related to the so-called Cowling channel (e.g. Boström, 1974) of the polarization electric field. Namely, if a primary (convection) electric field is aligned with a slab of higher conductance (such as an auroral arc), a secondary electric field driving secondary currents appears at the ends of the slab. These secondary currents carry Joule heating, and the higher the conductance in the slab, the higher the secondary current and hence the higher the Joule heating. Thus, in this case the electric field and conductance are correlated, which increases the Joule heating. However, the convection electric field is parallel to the auroral arc generally only in the nightside. On the other hand, for example, Marklund et al. (2001a) found that in the auroral bulge the secondary electric field is efficiently closed by local field-aligned currents, thus making the Cowling channel model inadequate to describe the phenomena related to the substorm current wedge.

Discussion and conclusions
In this paper we have estimated the Northern Hemisphere Joule heating with various techniques during a substorm event of 28-29 March 1998, aiming to quantitatively estimate the ability of GUMICS-4 global MHD simulation to predict ionospheric Joule heating. The purpose of this work is to calibrate the total ionospheric power consumption equation (Palmroth et al., 2004a) given by Eq. (2). As was mentioned earlier, Palmroth et al. (2004a) showed that Eq. (2) was able to predict the ionospheric total power consumption in the simulation from solar wind parameters with more than an 80% correlation coefficient. Since in the simulation both the Joule heating and precipitation were in temporal agreement with the AE index, it was argued that the right-hand side of Eq. (2) could roughly predict the temporal behavior of ionospheric power dissipation as determined by the AE-proxies. Although the components of P ionosphere in Eq. (2) yielded lower ionospheric power dissipation as compared to the AE proxies, Palmroth et al. (2004a) speculated that Eq. (2) could still be used to predict ionospheric power dissipation if the constant C in Eq. (2) was scaled to account for the "missing" ionospheric power in the GUMICS-4 simulation. Palmroth et al. (2004a) showed that in all simulated events the fitting yields similar values for the exponents a, b, and d in Eq. (2), given by 0.8, 2.8, and −2, respectively. Therefore, approximately Eq. (2) is proportional to solar wind kinetic energy flux ρv 3 multiplied by a term characterizing the IMF z effect, i.e. reconnection. The approximate ρv 3 behavior of Eq. (2) is supported by a statistical survey carried out by Olsson et al. (2004), who found that total ionospheric Joule heating is highly correlated with solar wind kinetic energy flux ρv 3 . This agreement between the instantaneous GUMICS-4 results and the statistical results by Olsson et al. (2004) suggests that GUMICS-4 Joule heating is in accordance with the system behavior over longer time periods. We now focus on the scaling of C in the following way: First, we evaluate whether GUMICS-4 Joule heating follows temporally the Joule heating occurring in the ionosphere during this event. Second, we estimate roughly the correct magnitude of Joule heating during the event. Third, we discuss the physical problems associated with the scaling and whether the scaling of GUMICS-4 results leads to a reasonable estimate of global Joule heating. Finally, we discuss whether the calibration of Eq.
(2) is successful as far as Joule heating is concerned. A similar calibration of the precipitation power will be carried out in a future investigation.
Although at global scales the Pedersen conductance is almost constant in time due to the large contribution of the dayside EUV radiation, Fig. 5 shows that at auroral latitudes the precipitation governs the spatial distribution of Pedersen conductance. However, Joule heating depends quadratically on the electric field. Therefore, the temporal variation of the total hemispheric Joule heating most probably follows the temporal variation of the global electric field, and consequently the polar cap potential. The polar cap potential difference has been found to be linearly proportional to the AE index (Ahn et al., 1984;Weimer et al., 1990), and therefore the temporal variation of the polar cap potential should follow that of the AE index (Fig. 2). Hence, the temporal variation of the AE proxy presented in Fig. 7 gives most probably also the temporal variation of the total integrated Joule heating. As was seen from Fig. 9, the temporal evolution of the GUMICS-4 Joule heating is well-correlated with that of the AE index. This is not a new result: in all events simulated with GUMICS-4 to date, the total integrated Joule heating follows the temporal variation of the AE index. It is thus likely that the temporal evolution of the global integrated Joule heating in the Northern Hemisphere is well reproduced by the GUMICS-4 global MHD simulation. However, the level of Joule heating in all the simulated events, the present event included, appears to be less than suggested by, for example, the AE-proxies.
Although the Foster et al. (1983) statistics were binned using the K p index, which has a poor temporal resolution, the results are based on the most comprehensive statistical survey published to date, in which the data are gathered from direct measurements on board the AE-C spacecraft. Furthermore, to our knowledge, the Foster et al. (1983) paper is the only statistical study using measurements of P and E in which the spatial anticorrelation of the electric field and Pedersen conductance is taken into account, as Foster et al. (1983) measure the two parameters simultaneously using the same satellite recordings. Therefore, the magnitude of global Joule heating may not be overestimated in the original database, although their K p -proxy may sometimes underestimate the Joule heating due to the poor temporal resolution of the K p . Olsson et al. (2004) carried out another statistical survey on the Northern Hemisphere field-aligned Poynting flux using the Astrid-2 satellite electric and magnetic recordings on ∼3000 orbits. They found a rather good agreement with the Foster et al. (1983) results, although strictly speaking the Poynting flux includes the energy driving both Joule heating and neutral winds, whereas the Foster et al. (1983) results deal with P E 2 only. Nevertheless, the agreement is worth noticing, particularly because the direct Poynting flux measurements naturally do not include any overestimations due to multiplication of the averaged electric fields and Pedersen conductances from different sources. Furthermore, the two statistical surveys are carried out using satellite recordings utilizing different types of instruments (double probes on board Astrid-2, and electric field and ion drift meter on board AE-C), and the time resolutions are different (1 s for Astrid data and 15 s for AE-C data). Therefore, the Olsson et al. (2004) results are to be taken as strong evidence that the much larger data set of Foster et al. (1983) gives the P E 2 more or less correctly. On the other hand, the Olsson et al. (2004) results are also in accordance with the ρv 3 behavior of Eq. (2).
For the present event, the K p proxy of Foster et al. (1983) suggest Joule heating values below 100 GW for the Northern Hemisphere. The global observations during this event suggest ∼70 GW of Joule heating during the first onset, whereas the value during the first maximum (60 GW) is most likely to be an underestimation due to the lower radar data coverage. The Ahn et al. (1983) Joule heating proxy suggests ∼200 GW at the peak of the first substorm, while AMIE suggests ∼400 GW. Evidently, the large scatter among the different Joule heating estimates adds confusion as to what value should be used when quantitatively calibrating the power law of Eq. (2). However, Figure 7 and Table 1 suggested that the different estimates of Joule heating agree quantitatively towards the end of the substorm sequence. Figure 9 suggests that during the second substorm the GUMICS-4 result multiplied by 10 is in quantitative agreement with the observed Joule heating, which does not suffer from poor radar data coverage at that time period. With a factor of 10 lower, the GUMICS-4 Joule heating result is also consistent with the AMIE Joule heating from 02:30 UT onwards, although during the maximum of the first substorm the AMIE results show quite large Joule heating rates. Given that the error estimate of the Ahn et al. (1983) proxy is ±50%, the lower limit of the AE-proxy is also consistent with the lower than 100-GW estimate. Hence, three of the four estimates for Joule heating presented in this paper essentially agree with each other during the course of the event, while all estimates agree with each other during the latter part of the simulated period. This leads us to conclude that the "true" hemispheric Joule heating during the 28-29 March 1998 substorm is below 100 GW at maximum, and the temporal variation follows that of the AE index.
The scaling of results from a global MHD simulation may seem disturbing prima facie, particularly when considering that the global MHD simulations do not generally include the neutral winds, and GUMICS-4 does not include the discrete arc physics as the parallel potential drop is set to zero in the current version. However, for example, Lu et al. (1995) estimate that the neutral winds take up to 6% of the electromagnetic energy. This is further supported by the quantitative agreement of the Foster et al. (1983) and Olsson et al. (2004) results, where the energy of neutral winds are excluded in the former but included in the latter. The fact that the two statistical surveys are in quantitative agreement indicates that the neutral winds do not consume large portions of the electromagnetic energy statistically, although they may be important during individual events (Thayer, 1998). Furthermore, Wiltberger et al. (2004) showed that the Joule heating stays essentially the same whether or not it was calculated from a global MHD simulation coupled to a simulation modeling the thermospheric neutral winds. Therefore, it is likely that the lack of neutral winds in GUMICS-4 is not a severe drawback as far as the total integrated Joule heating is concerned. It is not clear, however, how the discrete arc physics would affect the global Joule heating. Namely, if the current associated with discrete arcs closes primarily locally, the potential difference over which the current closes would be small. In this case the global Joule heating may remain unaffected or decrease. On the other hand, if the return current region adjacent to the discrete arc has a very small conductivity due to escaping electrons (Marklund et al., 2001b), the Joule heating associated with the global horizontal current would increase. Therefore, since the effect of the discrete arcs to global Joule heating is not resolved, it is unclear what their effect is on the GUMICS-4 results. Essentially, the scaling of GUMICS-4 results by a constant factor can be justified by the good reproduction of the global electric field. This can be seen in Fig. 12, which shows that the polar cap potential is well-reproduced temporally in the GUMICS-4 simulation, given that there are also uncertainties in the temporal variation of the SuperDARN measurements due to the occasinal lower data coverage. As the magnitude difference between the GUMICS-4 polar cap potential with the Super-DARN measurements during good data coverage is always quite close to what is observed in Fig. 12, we have reason to believe that the scaling factor is also similar in other events. We conclude that the GUMICS-4 Joule heating multiplied by a factor 10 to achieve an agreement with observational methods also asserts that the observational methods may overestimate Joule heating during the maximum of the first substorm. As can be seen in Table 1, the scaling of the GUMICS-4 results by a factor of 10 yields an estimate that agrees with the observational methods only towards the end of the simulated period. During the onset and the first maximum, even the scaled GUMICS-4 results are lower than the observational methods suggest. The reason for this may be that the observational methods do not take into account the anticorrelation effect, which is intrinsically taken into account in the ionospheric computation of GUMICS-4 (Figs. 10 and 11). Equation (5)  the overestimation made by estimating the Joule heating using independent average measurements of P and E depends on the standard deviations of these two variables. Sugino et al. (2002) presented 10 years worth of European Incoherent Scatter Radar (EISCAT) data for the electric fields and conductances. Among other issues, the Sugino et al. (2002) data set shows clearly that the larger the intensities of the electric fields and conductances are, the larger their standard deviations. It is clear that the magnitudes of the electric fields and Pedersen conductances are positively correlated with increasing magnetic activity. Therefore, Eq. (5) suggests that the overestimation of Joule heating due to spatial anticorrelation of E and P increases with increasing magnetic activity, both in absolute and relative terms. In practice, this would indicate that Joule heating during intense substorms and storms would be overestimated, while during less intense times the observational methods would give better predictions of the Joule heating. Naturally this does not apply to the Foster et al. (1983) results that do take the anticorrelation effect into account (to the extent that their 15-s temporal resolution allows); however the poor temporal resolution of the K p index may sometimes lead to situations where the Foster et al. (1983) proxy yields ambiguous Joule heating values.
The overestimation made by taking independent averages of electric field and Pedersen conductances over regions where the two parameters are spatially anticorrelated can be assessed quantitatively. Matsuo et al. (2003) estimated the electric field variability statistically at high latitudes using DE-2 satellite measurements. They found that near equinox the standard deviation of the electric field over high latitudes is rather steadily ∼30 mV/m above 60 • MLAT, while the maximum of the standard deviation of the electric field can reach as high as 78 mV/m. Therefore, the standard deviation of the electric field is of the order of the measured electric field; for the standard deviation of the Pedersen conductance we can assume a modest 2-5 S, although at the oval boundary this might be an underestimation. If, for the purpose of assessing the Joule heating overestimation, one assumes that ST D E 2 =(ST D E ) 2 , over the area covered in Matsuo et al. (2003), one obtains approximately 2 S·(0.03 V/m) 2 ·10 12 m 2 ≈ 1.8-4.5 GW in the case of perfect anticorrelation (Corr i =−1). If this value remains steady over the high latitudes, over the ionosphere above 60 • MLAT (containing ∼35 area elements of the size 10 12 m 2 ), the total value sums up to 63-153 GW, depending on the choice of the standard deviation of the Pedersen conductance (2 S or 5 S, respectively). For example, correlation coefficient of −0.7 the above estimate would yield 44-107 GW. Therefore, the error made in estimating the total Joule heating by using independent averages of electric fields and conductances can become as large as the estimate itself.
Naturally the overestimation of total Joule heating due to spatial anticorrelation of P and E applies only over those regions in the ionosphere where the electric field and conductance are anticorrelated. As shown by Foster et al. (1983), this condition holds generally in the nightside, while there are regions adjacent to the cusp, where the electric field and conductance are spatially correlated. In the case of spatial correlation of the electric field and conductance, the Joule heating would actually be underestimated while ignoring the standard deviations of the two quantities. The total value of the overestimation or underestimation is affected not only by the values of the standard deviations, but also the total areas over which the quantities are spatially anticorrelated or correlated, respectively. According to Foster et al. (1983), the region of spatial correlation of P and E is generally smaller than the area of spatial anticorrelation. Therefore, while integrating over the whole Northern Hemisphere, ignoring the standard deviation would most probably still yield an overestimation of the Joule heating. Note that this also applies to our estimate of the Joule heating based on global observations: we took averages of the Pedersen conductance and the electric field from different sources, and therefore by definition this method yields an overestimation of Joule heating elsewhere except in the cusp region.
In conclusion we find that at any given time the GUMICS-4 estimate of the Joule heating during the 28-29 March 1999 event is 10 times less than the "true" Joule heating. The lesser amount of Joule heating in the GUMICS-4 global MHD simulation is essentially caused by lower polar cap potentials and weaker Region 2 currents. As the closure of Region 1 and Region 2 currents to each other at the oval is the strongest source of Joule heating in the polar ionosphere (e.g. Gary et al., 1994;Olsson et al., 2004), the weaker Region 2 currents in the MHD simulations lead to a situation in which smaller amounts of current close over the oval, decreasing the Joule heating associated with the Region 1 and Region 2 field-aligned current closure. The temporal evolution of Joule heating in GUMICS-4 is, on the other hand, consistent with the "true" Joule heating. To calibrate the power formula of Eq. (2), we find that as far as Joule heating is concerned, the constant C in Eq. (2) should be multiplied by a factor of 10. From the physics point of view, multiplying the simulated Joule heating by a constant is of course not satisfactory, but it helps us to obtain an estimate of the Joule heating in cases where observational data are not available, for space weather prediction purposes, or in other cases where an estimate of the Joule heating is required.

Appendix A Correlation and Joule heating in the averaging process
Let σ (x) and E(x) 2 be the conductance and square of the electric field in some domain A, respectively, and let P = A d 2 xσ (x)E(x) 2 denote the corresponding Joule heating power. Consider idealised measurements of σ and E 2 that amount to spatial averaging over small subdomains (grid cells) A i of A, yielding the discrete quantities σ i ≡(1/A i ) A i σ (x)d 2 x and E 2 i ≡ (1/A i ) A i E(x) 2 d 2 x, i=1...N. We want to show that the Joule heating computed from the measurement data {σ i , E 2 i },P ≡ i A i σ i E 2 i always overestimates the true Joule heating P if σ and E 2 are anticorrelated in a subgrid scale, i.e. if from which it follows that Summing over i and using the definitions of P andP we ob-tainP >P , i.e. that the Joule heatingP computed from averaged data yields an overestimation of the true Joule heating P if σ and E 2 are anticorrelated in a subgrid scale in the sense of Eq. (A1). We can solve the magnitude of the Joule heating overestimation due to spatial anticorrelation of E and P . Using Eq. (A2), Eq. (A1) can be written as As the square root terms on the right-hand side are essentially standard deviations of σ and E 2 , the magnitude of Joule heating overestimation P corr in area A i is therefore calculated as where Corr i ∈ [−1, 0] and ST D σ and ST D E 2 are standard deviations of σ and E 2 , respectively, in area A i .