Polar heating in Saturn’s thermosphere

. A 3-D numerical global circulation model of the Kronian thermosphere has been used to investigate the inﬂuence of polar heating. The distributions of temperature and winds resulting from a general heat source in the polar regions are described. We show that both the total energy input and its vertical distribution are important to the resulting thermal structure. We ﬁnd that the form of the topside heating proﬁle is particularly important in determining exospheric temperatures. We compare our results to exospheric temperatures from Voyager occultation measurements (Smith et al., 1983; Festou and Atreya, 1982) and auroral H + 3 temperatures from ground-based spectroscopic observations (e.g. Miller et al., 2000). We ﬁnd that a polar heat source is consistent with both the Smith et al. determination of T ∞ ∼ 400 K at ∼ 30 ◦ N and auroral H + 3 temperatures. The required heat source is also consistent with recent estimates of the Joule heating rate at Saturn (Cowley et al., 2004). However, our results show that a polar heat source can probably not explain the Festou and Atreya determination of T ∞ ∼ 800 K at ∼ 4 ◦ N and the auroral H + 3 temperatures simultaneously.


Introduction
The temperature of Saturn's thermosphere is much higher than can be explained by solar EUV heating alone. Several studies have calculated the temperature contrast T =T ∞ −T m between the exospheric temperature T ∞ and the mesopause temperature T m . An early calculation by Strobel and Smith (1973) found T ∼10 K. Simple thermal balance calculations for Jupiter by Hunten and Dessler (1977) have recently been applied to Saturn by Yelle and Miller (2004), who found T ∼35 K. If T m =143 K (Yelle and Miller, 2004), then both these studies predict T ∞ <200 K.
The existing temperature measurements consist principally of solar and stellar occultations at low latitudes Correspondence to: C. G. A. Smith (chriss@apl.ucl.ac.uk) performed by the Voyager spacecraft (Smith et al., 1983;Festou and Atreya, 1982) and ground-based infrared spectroscopic measurements of H + 3 temperatures in the polar regions of the planet (Miller et al., 2000). The Voyager solar occultation experiment (Smith et al., 1983) returned an exospheric temperature of T ∞ ∼420 K at a latitude of 30 • N. In contrast, the stellar occultation experiment (Festou and Atreya, 1982) gave a value of T ∞ ∼800 K at 4 • N. In the decades since Voyager there has been considerable debate over the relative validity of these results. For a detailed analysis of the advantages and disadvantages of occultation measurements, and a discussion of the Smith et al. (1983) and Festou and Atreya (1982) measurements at Saturn, we refer the reader to the excellent review by Smith and Hunten (1990).
A recent reanalysis of the Voyager occultations (Vervack et al., 2005 1 ) suggests that T ∞ at 30 • latitude may be closer to ∼550 K. However, for the purposes of this study we will continue to consider the published results of Smith et al. (1983) and Festou and Atreya (1982). We will, however, discuss the implications of the reanalysis upon our results.
Infrared spectroscopic studies of the polar regions (Miller et al., 2000) broadly support the Voyager measurements. These initially indicated H + 3 temperatures of ∼600 K in the region from the pole to ∼15 • colatitude, encompassing the main auroral oval, but a more recent analysis of H + 3 spectra (Melin, 2005 2 ) gave an average polar cap temperature of 425±75 K. It is to be emphasised that these measurements represent temperatures at the H + 3 peak altitude, not exospheric temperatures.
Although there is a large range of values from the experimental determinations, the measured temperatures are all in the range 400-800 K. This is clearly much higher than the value of <200 K expected from solar EUV heating. There have been a number of attempts to explain this discrepancy. The two most commonly suggested sources for the extra energy are transport from the lower atmosphere, most probably by buoyancy waves or acoustic waves, or energy from the magnetosphere, which is expected to be deposited largely in the polar regions of the planet. The magnitude and spatial distribution of these supplementary heat sources is largely unknown.
The expected power input from the lower atmosphere has been studied in detail by a number of authors for the Jupiter case Matcheva and Strobel, 1999;Hickey et al., 2000;Schubert et al., 2003). The heating rates derived are in general not enough to explain the elevated thermospheric temperatures at Jupiter. Discussion of wave heating rates at Jupiter has focussed on the wave structures observed by the Galileo probe , and the heating rates these may imply. There are no equivalent measurements at Saturn, so the wave heating rate remains undetermined.
Magnetospheric energy sources have been studied in much more detail thanks to constraints provided both by in-situ measurements of the magnetosphere and by remote sensing of the auroral emissions. Most recently, Cowley et al. (2004) have presented a simple self-consistent model of the plasma flows across the whole of Saturn's polar ionosphere. They assume an "effective" Pedersen conductivity of 1 mho (Bunce et al., 2003) and then derive a total global particle heating power of the order of ∼0.1-0.2 TW and a total global Joule heating power of ∼10 TW distributed between the pole and 25 • colatitude. The particle heating power is consistent with previous estimates based on Voyager data (e.g Atreya, 1986) and on Hubble Space Telescope observations (Trauger et al., 1998). The Joule heating powers are considerably higher than earlier estimates: for example Atreya (1986) reported estimated Joule heating of 0.15 mWm −2 , which, integrated across the whole polar cap down to 25 • colatitude amounts to only ∼0.3 TW. The discrepancy is explained by the relatively large electric fields of the Cowley et al. model, since the Joule heating varies as the electric field squared. The open field line region, for example, subcorotates by 70% with respect to the planet, consistent with recent observations by Stallard et al. (2004); the Atreya study assumed subcorotation of only 10%, without specifying a detailed latitudinal structure.
It seems reasonable to conclude that particle heating is negligible compared to Joule heating. However, the observational evidence for particle heating powers is derived from UV observations only. Infrared observations (Stallard et al., 2004) indicate H + 3 emission across the entire polar cap. It is possible that this H + 3 is generated by soft precipitation that is not detected in the UV. We will discuss this possibility further below. Smith et al. (2005) have recently proposed two further sources of heating, also associated with the magnetosphere. Miller et al. (2000) first pointed out that the polar thermosphere of Jupiter may be spun into subcorotation with respect to the planet by the ion drag associated with the magnetospheric electric fields. Millward et al. (2005) have recently presented modelling results that support this conjecture. Similar behaviour is expected at Saturn. Smith et al. (2005) have quantified the kinetic energy input from this process in terms of the ion-neutral coupling parameter k, as defined by Millward et al. (2005). For plausible neutral wind configurations the magnitude of this K.E. input may be comparable to the Joule heating. In steady state, this K.E. entering the atmosphere through ion drag must be continually dissipated as heat. Although this may be an important energy input, the spatial distribution of the dissipation is unknown. Smith et al. (2005) also showed that the Joule heating will be increased if the electric field exhibits significant fluctuations about its mean value. This is known to be an important effect at Earth (Codrescu et al., 2000) and is likely to have an impact at Saturn. Again, its magnitude is unknown, but it is likely to be of a similar order to the mean-field Joule heating of 10 TW calculated by Cowley et al. (2004).
It is clear that some components of the polar heating are very well quantified -in particular heating associated with the main oval FUV aurora -but many uncertainties remain. The purpose of this paper is thus to study the effects of a general polar energy input, whose origin we do not specify.
In addition to the difficulties encountered in constraining the form of the polar heating, the question of whether the energy may be distributed globally by thermospheric winds remains unanswered. The principle barrier to redistribution of polar heating is thought to be the Coriolis force, which will turn westwards any equatorwards wind. The radius of Saturn is ∼10 times that of Earth, and its angular rotation speed ∼2 times as fast. Thus the Coriolis forces are expected to be ∼20 times those on Earth, and are likely to dominate thermospheric dynamics.
Given the lack of reliable measurements of the actual heating rates and wind velocities, we do not expect our model to represent a perfect physical match to the actual conditions at Saturn. Instead our intention is to understand in some detail the dynamics and temperature distributions produced by polar heating, and to assess the plausibility of polar heating as a global thermospheric energy source.
A preliminary study along these lines has already been carried out by Müller-Wodarg et al. (2005) 3 . They employed the same model that we describe below to investigate the power input from the lower atmosphere and from the auroral zones required to explain the observed temperatures. They were able to reproduce the Smith et al. measurement of ∼400 K at 30 • N by imposing either an arbitrary, globally distributed "wave heating" rate of ∼14 TW or a combination of a polar heat source of ∼10 TW of Joule heating (Cowley et al., 2004) with an arbitrary "wave heating" rate of ∼6 TW. Our study builds on this work by examining a full parameter space of simple, plausible polar heating distributions.
The plan of the paper is as follows. In Sect. 2 we describe our model and the basic heating distributions that we have employed. In Sect. 3 we describe the general response of the model to a typical polar heating input, and in Sect. 4 we go on to examine the sensitivity of this response to a full parameter space of heating distributions, and compare these results to observations. In Sect. 5 we discuss the limitations of our results and the effects of some simple modifications to our model. Finally in Sect. 6 we conclude.

The Model
We use the Kronian thermosphere model described by Müller-Wodarg et al. (2005) 3 . The model is a 3-D General Circulation Model (GCM) which solves the coupled non-linear Navier-Stokes equations of energy, momentum and continuity on a latitude, longitude, pressure grid. The Appendix of Müller-Wodarg et al. (2005) 3 describes the equations in detail. The lower boundary is at the 100 nbar level (approximately 800 km above the 1bar level) and it has 28 levels with a vertical resolution of 0.4 pressure scale heights. The lower boundary condition is a fixed temperature of 143 K and a horizontal wind velocity of zero. We use a latitudinal resolution of 1 • in order to resolve a reasonably thin oval of auroral heating.
For this study, we have assumed zonal symmetry of the Kronian thermosphere. However, we calculate the EUV heating for a full range of longitudes, and then use the zonal average as the input to the model. The assumption of zonal symmetry makes a very small difference to our results. As discussed by Müller-Wodarg et al. (2005) 3 , solar EUV heating produces a diurnal temperature contrast of only 13 K. This temperature contrast is negligible compared to the high temperatures which we seek to generate with polar heating. For several of the runs described below we have performed validation runs using the full 3-D model, including diurnally variable solar heating. As expected, the resultant diurnal temperature differences are only of the order of ∼10-20 K.
The model also provides the facility to calculate transport of the major neutral gases (H, H 2 and He) by winds and by diffusion. We have chosen to omit these calculations for reasons of computational efficiency. The effect of turning on full composition calculations is discussed in Sect. 5. We have also set the eddy diffusion coefficient to zero in our standard calculations. The effect of this is that the increased thermal conductivity and viscosity associated with eddy diffusion processes are not included. This assumption will also be discussed briefly in Sect. 5.
Composition profiles were taken from  and the model was started with an initial uniform thermospheric temperature of 143 K, corresponding to the fixed lower boundary temperature. It was then run for the equivalent of 800 Kronian rotations, in an equinox configuration, with solar EUV heating as the sole energy input. This long initial run ensured it had reached near steady-state before the inclusion of a polar heat source. The "steady-state" exospheric temperature which this set-up run achieved was 180K at the equator and 160 K at the poles. This confirms the estimates of Yelle and Miller (2004)  We have then explored a parameter space of 25 runs, within which the total magnitude and peak altitude of the polar heating is varied. All runs were initialised using the "steady-state" run described above. For each simulation, an identical polar heat source is imposed in each hemisphere. Latitudinally, the heat source is distributed as a Gaussian with FWHM 3 • , centred at 14 • colatitude (see Fig. 1). We have chosen this form because we expect the polar heating to be concentrated in the auroral zones. However, we have experimented with more horizontally distributed forms, the effects of which are discussed in Sect. 5.

discussed in
Vertically, the distributions ( Fig. 2) resemble Chapman profiles (Ratcliffe, 1972). Specifically, we specify the heating rate per unit mass q m as a function of the pressure p: Here q m∞ is the asymptotic value of the heating rate in the upper thermosphere. We fix q m , so the heating is not dependent on the thermal structure of the atmosphere. This means, however, that the heating rate per unit volume q v is dependent on the thermal structure. However, we can usefully specify the heating rate per unit area per unit pressure scale height q h -which is equivalent to q v if the atmosphere is isothermal -by multiplying q m by the mass density per unit area per unit pressure scale height p/g: where q h0 =q m∞ (p 0 /g)e −1 and g is the acceleration due to gravity. The peak heating rate is then q h0 at a pressure of p 0 . We use p 0 values of 60, 40, 20, 6, and 2 nbar (Fig. 2). We have quoted both of these formulae because it is useful to understand the behaviour of both q h and q m . The rate q h represents more accurately the location of the heating inputs: most of the heating enters the atmosphere at the peak of q h . However, q m is more representative of the atmospheric response to the heating, because it indicates the energy input relative to the local heat capacity of the atmosphere.
The physical origin of the Chapman profile is the description of atmospheric heating and ionisation rates resulting from monochromatic irradiation. However, in this context we employ the Chapman profile purely as a simple distribution which contains the main features we expect a polar heat source to exhibit. In particular, a sharp cut-off below the peak energy input at p 0 and an exponentially decaying tail on the topside. The Chapman profile is a very specific case in that q m tends asymptotically to a constant value with altitude. The consequences of deviations from this simple behaviour are discussed in Sect. 5.
The morphologies of the low and high altitude limits of the Chapman profile are both consistent with the general form of energy deposition produced by particle precipitation (e.g. Grodent et al., 2001). The Joule heating rate should also have a similar form, since its structure is largely controlled by that of the ionosphere. The ionosphere is believed to exhibit a flat bottomside (due to the sudden onset of hydrocarbon chemistry near the homopause, which leads to the rapid destruction of H + 3 ), and the topside ion densities are controlled by ambipolar diffusion, which produces an exponential decay with altitude (Moses and Bass, 2000). Both these features will also be present in the Joule heating rates; we discuss the case of Joule heating in more detail in Sect. 5.
Radiative cooling is not included in the model. The H + 3 molecular ion is a strong radiator in the near-infrared, and at Jupiter it is known to be an important coolant of the thermosphere. In the auroral regions of Jupiter it is an order of magnitude more important than thermal conduction, and in the Jovian equatorial regions it is of comparable importance. At Saturn, it is believed to be of some significance in the auroral regions (Miller et al., 2000). However, it is much less effective at lower latitudes where H + 3 densities and temperatures are lower. Thus it may affect the immediate energy balance of the polar regions, but not the effectiveness of midlatitude wind systems that may redistribute the polar energy. As such our polar heat source should be thought of as a net energy input, taking into account local radiative cooling in the auroral regions.
We calculate the total energy inputs by numerically integrating the normalised horizontal and vertical distributions across our model grid. We then set the parameter q h0 to give round values of the total energy input. Specifically, we use total polar energy inputs (integrated over both hemispheres) of 1, 2, 4, 6 and 10 TW. For brevity in the following discussion we will describe a run with an input of x TW peaking at y nb as x TW/y nb.
Each model was run for 400 Kronian rotations. To check that 400 rotations is sufficient to represent steady state, we have run the models at the "corners" of our parameter space -i.e. those representing inputs of 1 and 10 TW at pressures of 60 and 2 nbar -for 800 rotations. We find that by 800 rotations the model temperatures are exponentially approaching an asymptote. Figure 3 shows the time development of the exospheric temperature at 30 • N as a percentage of this asymptotic value. By 400 rotations the temperatures all lie well within 10% of the asymptote. Note that this represents a systematic underestimate of the temperature at 30 • N. We find 10% to be a perfectly acceptable error given the general nature of our modelling.

General response
The general response of the model is illustrated for the 6 TW/40 nb run in Fig. 4. The colour contours represent temperature in Fig. 4a and zonal wind velocity in Fig. 4b. In both figures the vertical and meridional velocities are plotted as arrows. The vertical scale has been stretched for clarity, and the vertical wind speeds have been scaled accordingly to preserve the geometry of the flow. Thus it should be made clear that in the plot shown vertical velocities are typically <5 m/s, whereas peak meridional velocities are of the order of 300 m/s. The thickness of the arrows represents the combined vertical/meridional wind speed; this is dominated by the meridional component.
The model grid points are shown as dots. Note that the altitude of the upper pressure level of the model falls off strongly towards the equator, as a simple consequence of the thermal structure. The two crosses indicate the location of the energy inputs. The co-latitude of the crosses corresponds to the peak of the horizontal distribution shown in Fig. 1. The lower cross shows the location at which the heat input per scale height q h peaks (see Fig. 2). The upper cross shows the altitude at which the heat input per unit mass q m has risen to ∼90% of its asymptotic value q m∞ ; for brevity we will refer to this point as the "shoulder" of q m .
The general behaviour is as follows. The heated zone at 14 • colatitude exhibits strong upwelling driven by the heating. There is a temperature peak here at about 2000 km altitude. Below the peak the only significant energy input comes directly from the imposed heating profile. The atmosphere is cooled by a combination of downwards thermal conduction and adiabatic cooling of the rising gas. Above the peak, thermal conduction transfers heat upwards, and this is almost as important an energy source as the topside of the imposed heating profile.
The temperature peak does not correspond in altitude to the peak of q h , but rather to the "shoulder" of q m . Thus q m is a more useful measure of the heating rate than q h for determining the thermal structure. The vertical temperature profile follows that of q m quite closely: below the "shoulder" both T and q m increase with altitude; above the "shoulder" they are both approximately constant.
At low altitudes, pressure gradients drive meridional winds both poleward and equatorwards of the heated zone. The polewards flow is convergent as a consequence of the spherical geometry, and this forces a sinking cell to form inside the polar cap. The sinking gas undergoes adiabatic heating. This heating almost "fills in" the polar cap with the same temperatures observed in the main oval, although there is still a small residual temperature gradient driving the flow. The equatorwards flow is driven by a much steeper meridional pressure gradient. This flow is divergent, again as a consequence of the geometry. The divergence of the flow reduces the heating efficiency of the gas as it moves to lower latitudes, since the energy must be spread over a wider area. The gas is also cooled efficiently by thermal conduction as it moves to lower latitudes. The result is a persistent sharp pressure gradient equatorwards of the oval. At higher altitudes, the behaviour is essentially identical, except that the highest temperatures are found at the pole and one single equatorwards flow develops. The high-altitude hotspot at the pole results largely from upwards thermal conduction from the hot polar cap beneath.
The flows away from the oval are driven by pressure gradients, but are substantially modified by the Coriolis force to generate zonal winds (Fig. 4b). The polewards flow inside the main oval is forced into weak superrotation (eastwards flow). The much stronger equatorwards flow outside the oval is driven into a fast subcorotating (westwards) flow. We find that vertical viscous drag is the principle force limiting this zonal flow. In total, the Coriolis force is thus in balance with the combined meridional pressure gradient and vertical viscous drag forces. The importance of viscosity in controlling this behaviour will be discussed further in Sect. 5.
As the flow moves towards lower latitudes the pressure gradient, and thus the flow velocity, reduces significantly. As the gas slows it sinks, and this leads to some adiabatic heating of the low-latitude region. The adiabatic heating at the equator is significant, accounting for ∼0.04 mWm −2 compared to only ∼0.01 mWm −2 from solar heating.

Exploring parameter space
An important advantage of a GCM as compared to a 1-D model is that it allows the sparse and spatially separated observations of the Kronian thermosphere to be tested simultaneously. To compare the available observations with the model-derived temperatures it is, of course, necessary to associate each measurement with specific grid points within the model. In the case of the spacecraft occultation data we associate the measured exospheric temperatures of 400 K at 30 • N (Smith et al., 1983) and 800 K at 4 • N (Festou and Atreya, 1982) with the model-derived temperatures at the topmost pressure level and at the appropriate latitude. In the case of the ground-based spectroscopic measurements we assume that the observed emission is associated with the peak of the ionosphere, and that this in turn is, approximately, spatially coincident with the peak of q h . Thus we associate the observed ionospheric temperatures with the modelderived temperatures at the pressure p 0 where q h peaks vertically. The spectroscopic measurements represent an average temperature across the polar cap, so we average the temperature at the peak altitude across the whole polar cap within 14 • colatitude. Henceforth we will refer to this polar cap average as the "auroral" temperature.
The results of the complete grid of 25 model runs are summarised in Figs. 5 and 6. Figure 5 plots the auroral temperature for each model run along the x-axis and the corresponding 30 • N temperature along the y-axis. Figure 6 is identical, but with the 4 • N temperature plotted along the y-axis. Note the log scales. Points corresponding to the same height distribution are connected by lines (the line formats correspond to those in Fig. 2); the energy input is represented by the plot symbol.
The observational constraints are indicated by the horizontal and vertical grey bands. The horizontal band in Fig. 5 represents the Smith et al. (1983) temperature of 400 K at 30 • N, with a width of ±25 K. This range is intended not as an error bar, but as a guide to the scale of the plot and the closeness of the model-derived points to the constraint. The horizontal band in Fig. 6 represents the 800 K at 4 • N measurement of Festou and Atreya (1982) in the same way. The vertical bands on both of Figs. 5 and 6 represent the "auroral" measurement of 425 K (Melin et al., 2005) 2 with an error of ±75 K. The dark grey box is the "target" region within which both constraints are satisfied.  5. Comparison between auroral temperature (x-axis) and temperature at 30 • N (y-axis). The grey horizontal band represents the Voyager solar occultation measurement of ∼400 K (Smith, 1983), with a range of ±25 K. The vertical band represents the auroral temperature of 425±75 K (see text).
Focusing now on Fig. 5, we can see that our parameter space covers situations in which the atmosphere is both much colder and much hotter than the observations. For the 1 TW/60 nb run -which has the lowest power input at the lowest altitude -the auroral and mid-latitude temperatures are barely discernible from the case where there is no auroral input. At the other extreme the 10 TW/2 nb run exhibits an auroral temperature which approaches 2000 K and a midlatitude temperature greater than 1000 K.
It is also clear that the temperature distribution is sensitive both to the total power input and to the vertical distribution. Increasing either the total input power or the altitude of the input increases both temperatures. This is unsurprising, since in both cases the ratio of the energy input to the heat capacity of the gas being heated increases (since at higher altitude the gas is less dense and thus has a lower heat capacity). However, the loci of constant vertical distribution and varying energy input (connected by lines) generally have steeper gradients than the loci of constant energy input and varying vertical distribution (plot symbols). This indicates that the auroral temperature is most sensitive to the altitude of input, whereas the low latitude temperature is more sensitive to the total input power.
To reproduce either of the observational constraints individually presents no problem. However, the regions of parameter space that satisfy the 30 • N constraint tend to produce a rather overheated auroral zone, e.g. 4 TW/6 nb. Conversely, the regions that produce the correct temperatures in the auroral zone mostly underheat the low-latitude regions, e.g. 2 TW/20 nb. However, our parameter space does intersect the observations. The point at which both constraints are satisfied represents an input of ∼8 TW/60 nb. Indeed, energy inputs in the range 6-10 TW peaking between 60-40 nbar are all roughly consistent with both the auroral and mid-latitude temperatures. A significant factor in the ability of the model to reproduce both constraints simultaneously is the decision to associate the "auroral" constraint with the temperature at the point of peak energy input, rather than the exospheric temperature. As discussed in Sect. 3 the temperature structure responds more strongly to the per unit mass heating q m than to the per unit scale height heating q h . Thus the temperature peak lies above the peak of q h . This is also partly a consequence of the proximity of the peak of q h to the cold (143K), fixed temperature mesopause. This allows the region of peak energy input to be cooled efficiently by thermal conduction. Meanwhile, as described above, the upper altitude ranges are heated both directly by the topside of the heating profile, and indirectly through upwards thermal conduction from the temperature peak. Thus we find that the polar exospheric temperature may be considerably higher, by as much as hundreds of degrees Kelvin, as the auroral temperature (at the peak of q h ). It is thus, in principle, possible for the auroral temperature to be colder than the equatorial exospheric temperature while a pressure gradient still exists in the upper thermosphere to drive a redistributive wind system. This means that there is no a priori inconsistency between an auroral temperature of ∼425 K and an exospheric temperature at 4 • N of 800 K as reported by Festou and Atreya (1982).
However, looking now at Fig. 6, we see that when we try to reproduce the Festou and Atreya (1982) measurement at 4 • N, this situation does not arise. To reach 800 K at 4 • N we require an input of 10 TW/6 nb or, visually interpolating from the graph, ∼9 TW/2 nb, and this gives an auroral temperature three times the constraint. It is not possible to find a combination in our parameter space which gives both measurements simultaneously.
The recent reanalysis of the Voyager occultations by Vervack et al. (2005) 1 , suggests that the exospheric temperature at 30 • latitude may be closer to 550 K. This only represents an increase of 150 K over the constraint marked on Fig. 5, and it is clear that to achieve this temperature we require slightly more energy input at a slightly lower altitude: for example 10 TW/60 nb produces a good match between the new mid-latitude temperature and the 425 K auroral constraint.

Comparison to likely energy inputs
We have employed a general polar heating distribution rather than a specific physical heating mechanism because of the uncertainty in the exact nature of the polar heating. We now briefly assess the physical reasonableness of these inputs.

Joule heating
If the contribution of neutral winds to the electric fields is neglected, Joule heating is directly proportional to the ionospheric Pedersen conductivity. In the thermosphere the electron-neutral collision frequency is small and the Pedersen conductivity depends only on the ions. The Pedersen conductivity σ P i due to an ion i can be given in terms of the ratio r i =ν in / i where ν in is the ion-neutral collision frequency and i is the gyrofrequency: Here n i is the ion density, q i is the ion charge, and B is the magnetic field strength. The function f (r i ) is strongly peaked at r i =1; the altitude where this occurs depends upon the thermal structure. To calculate r i we take a nominal value of B 60 000 nT, appropriate for the polar regions of Saturn (Cowley et al., 2004), and use the expressions for ν in given by Banks and Kockarts (1973). We find that for the full parameter space of 25 runs r H + =1 in the range 50-80 nbar and r H + 3 =1 in the range 20-60 nbar. These ranges correspond well to the regions of parameter space in which we find best agreement between our model and the observations. However, if the ion densities are low at pressures greater than 10nbar, as suggested by ionospheric models (Moses and Bass, 2000), then the Pedersen conductivity will also be low in this region. The conductivity is more likely to peak close to the peak of ion density, above the 10nbar level. However, the ionospheric structure is complicated, and at high latitudes depends critically on the level of the homopause and the nature of any particle precipitation that may generate ionisation. A detailed treatment of these issues is beyond the scope of this work.
The total polar heating powers of 6-10 TW that give best agreement between our models and the data are also similar to the Joule heating powers predicted by Cowley et al. (2004). However, these powers depend critically upon the assumption of a height-integrated Pedersen conductivity of 1mho. The value of 1mho has its origin in magnetic field measurements of regions of the magnetosphere that map to the thermosphere equatorwards of the main oval (Bunce et al., 2003). It is not clear whether the same values would pertain to the polar cap region, especially given the H + 3 emission across the polar cap observed by Stallard et al. (2004), which would suggest an enhanced conductivity in this region. If the conductivity were greater than 1mho within the polar cap, the total Joule heating powers would increase correspondingly. Thus it is possible that a total power of 10TW is an underestimate. Again, a quantification of the true conductivity and heating rates is beyond the scope of this work. The issue is also complicated by the effect of strong neutral winds on the Joule heating rate (e.g Smith et al., 2005).

Particle heating
As discussed in the introduction, there is good evidence that particle heating associated with the UV aurora is negligible. However, Stallard et al. (2004) report observations of the IR aurora which indicates strong emission by the H + 3 molecular ion across the whole polar cap. The presence of significant H + 3 emission inside the polar cap indicates either that there is a local ionisation source, or that the H + 3 is transported inwards from the main oval.
For transport to explain the presence of H + 3 there would need to be strong polewards thermospheric winds that were capable of transporting H + 3 across the polar cap within its recombination timescale. This timescale is given by τ ∼[10 −13 n e ] −1 , such that if the electron density is n e ∼10 9 m −3 then τ ∼10 000 s. The radius of Saturn's polar cap is approximately 15 000 km. Thus wind speeds of the order of 1.5 km/s would be required to fill the polar cap with H + 3 . It is not clear whether these wind speeds are plausible, nor how they would be generated.
Furthermore, the ions would need to be fully coupled to the neutrals for this transport to be effective. This requires the ratio r i discussed above to be greater than 1, so that ionneutral collisions dominate over electromagnetic forces. This occurs below the level where r i =1, i.e. below the 20 nbar level (see above). In this region the hydrocarbon densities are high enough to rapidly destroy both H + 3 and H + (Moses and Bass, 2000). Thus the only region where significant ion transport is likely to occur lies in the region of the thermosphere where ion lifetimes are extremely short. For this reason we believe the H + 3 densities are unlikely to be explained by transport.
The H + 3 must therefore be produced by local ionisation. If this exists it must be soft precipitation that does not produce a bright UV aurora. Soft precipitation is consistent with the generation of H + 3 since it would result in ionisation above the homopause, where hydrocarbon densities are low and the H + However, there is no reason why it should not deposit some energy at lower altitudes also. It is difficult to place any sound limits on the distribution or quantity of this energy source.

Other heat sources
In the introduction we also discussed the possibility of heating due to the dissipation of kinetic energy associated with ion drag, and increased levels of Joule heating caused by fluctuations about the mean electric field. The possible magnitudes of these energy sources has been discussed by Smith et al. (2005). However, their spatial distribution is unknown, and we defer detailed discussion of these energy sources to a future study.

Sensitivity to vertical energy distribution
As described in Sect. 2, the functions that we have chosen to describe the vertical distribution of our polar heating source have a very specific form. In particular the shape of the topside heating, which tends asymptotically to a constant per unit mass heating rate in the upper thermosphere, represents a very specific case.
In order to test the importance of the topside heating rate, we have repeated the 6 TW/40 nb run with modified vertical heating profiles (Fig. 7). We modify the topside by altering the exponential decay constant of the heating rates above the peak such that the total heating rate is either increased or decreased by 1 TW (dotted lines); we modify the peak by adding or subtracting a Gaussian from the peak, again to increase or decrease the total rate by 1 TW (dashed lines).
The results are shown in Fig. 8, which plots the same parameters as Fig. 5. The solid line is equivalent to the solid line in Fig. 5, i.e. it represents the parameter space of 40 nb runs. The dotted lines join the runs with modified topside heating, and the dashed lines join the runs with modified peak heating.
The most striking feature of Fig. 8 is that the dotted lines are almost vertical. Thus the topside modification has a negligible effect on the auroral temperature, and a large effect on the exospheric temperature at low latitudes: the exospheric temperature is, to an extent, decoupled from the auroral temperature. Furthermore, adding or subtracting 1TW from the topside changes the exospheric temperature at mid-latitudes by ∼100 K. This represents a modification of less than 20% to the total heating rate, but a change of ∼40% in the temperature contrast with the mesopause. The principal conclusion to draw from this observation is that it is essential to correctly quantify topside heating when assessing the effectiveness of a polar heat source. It is common in the literature to quote values for the total or height-integrated energy input: however these are clearly unhelpful if the structure of the topside -which represents a relatively small percentage of the total energy input -is improperly constrained.
The peak modification has a strong effect on the auroral temperature: this is unsurprising as we have defined the auroral temperature as the temperature at the heating peak. However, it has a similar effect on the exospheric temperature. Thus changes to the peak heating rate have a similar influence throughout the thermosphere. This is unsurprising, since the peak heating is at low altitude, and the rest of the thermosphere thus "sits" on top of it.

Sensitivity to horizontal energy distribution
Latitudinally, our polar heat source is concentrated in a thin oval close to the observed UV auroral oval at ∼14 • colatitude. However, it is possible that the polar heating may be distributed more widely across the polar cap. For example, the Joule heating calculated by Cowley et al. (2004) is distributed between the pole and ∼25 • colatitude. Similarly, the possible soft precipitation within the polar cap, on which we have speculated above, would presumably generate heating across the whole polar region.
We have investigated the effect of a more distributed polar heat source by repeating the 6 TW/40 nb run using a latitudinal distribution representing a polar cap "filled in" with uniform heating. This represents the opposite extreme to The plus signs show the temperatures produced by modifying the heating distributions, as marked. Dotted lines join together runs with modified topside heating. Dashed lines join together runs with modified peak heating. The dot-dash line is joined to the run using a "filled in" polar cap distribution. concentrating the heating in a narrow oval. The distribution is shown in Fig. 1 by the dashed line. In the runs using this distribution, we rescaled our peak heating rate q h0 to produce identical total energy inputs. The effect of this modification upon the auroral and mid-latitude temperatures is marked on Fig. 8.
We find that the only significant effect of this modification is to alter the thermal structure and dynamics of the polar cap. Instead of the temperature peak just within the main oval, the temperature peak develops at the pole; and the two circulation cells redistributing heat into and away from the oval are replaced by a single circulation cell which redistributes heat equatorwards. The auroral temperature, which we defined as an average across the polar cap, increases slightly (Fig. 8) because the hotspot generated by the heating lies fully within the region over which we average. This is thus largely a consequence of our definition of the auroral temperature.
The effect on the low latitude temperature is more interesting, in that it stays almost exactly the same. There is a slight increase but it is only of the order of 10-20 K. This insensitivity to the polar cap energy distribution appears to be a consequence of the behaviour of the winds within the polar cap, which, as described in Sect. 3, almost "fill in" the polar cap with a uniform temperature. The roughly uniform polar cap temperature then "feeds" the global circulation. The exact polar cap energy distribution has a small effect on the resultant mid-latitude temperatures, but it is effectively negligible.
Thus we conclude that while the mid-latitude thermal structure is very sensitive to the vertical distribution of energy, it is almost entirely decoupled from the horizontal distribution of energy in the polar cap.

Importance of thermal conduction and viscosity
In the thermosphere mean free paths are large and molecular diffusion processes dominate the physics. Of particular importance are molecular thermal conduction (diffusion of heat energy) and molecular viscosity (diffusion of momentum and kinetic energy).
Thermal conduction is important because, in the absence of significant radiative cooling, it is the principle coolant of the thermosphere. The thermal conductive flux F (positive downwards) in the thermosphere is given by: where κ is the thermal conductivity, T is temperature, and z is altitude. The altitude, z, depends on T because if the atmosphere is heated it will expand. Thus heating the thermosphere does not significantly change dT /dz. If κ was independent of temperature then the ability of thermal conduction to remove energy from the thermosphere would remain approximately constant as temperatures increased. However, κ is given by an expression of the form κ=AT s where A is a constant and s ∼0.7. Thus the thermal conductivity increases with temperature, but the non-linearity of the dependence means that the thermal conductivity increases more slowly with temperature than the internal energy of the gas U =c P T . Due to this non-linearity, heating the thermosphere has a weak insulating effect, because higher temperatures lead to less efficient thermal conduction relative to the energy stored in the gas. The consequence of this is that at higher temperatures the redistribution of energy by winds is more efficient, because conduction cools the gas less before it reaches mid-latitudes.
This causes the lines in Figs. 5 and 6 to curve upwards as greater amounts of energy are input -the hotter models are closer to a situation in which the global temperature distribution is uniform.
The viscosity is also of importance to the effectiveness of the redistribution. The zonal force balance at high altitudes between vertical viscous drag and Coriolis forces is important in determining the speed of meridional winds. If the viscous drag is reduced a faster westward flow develops, increasing the poleward Coriolis forces, which slow the equatorwards flow. Thus a lower viscosity leads, perhaps counterintuitively, to less effective redistribution.
It is worth asking whether we may have underestimated the viscosity. The molecular viscosity µ is given by expressions of the form µ=λT s , where s∼0.7. For molecular hydrogen, the major gas in our model, λ=1.46×10 −7 and for atomic hydrogen, the only other gas important at high altitudes, λ=2.07×10 −7 . (The units of µ are kg m −1 s −1 .) Thus atomic hydrogen is approximately 40% more viscous than molecular hydrogen. In our startup profiles atomic hydrogen represents no more than ∼10% of the gas in the upper altitude ranges. This does not change because we run the models with fixed composition profiles. To significantly increase the viscosity would require a greater proportion of atomic hydrogen by tens of percent. However, the strong vertical winds in the auroral regions are likely to mix greater quantities of molecular hydrogen into the high altitude regions, thus reducing the viscosity. To test this hypothesis we repeated the 6 TW/40 nb run with full composition transport calculations switched on. As expected, molecular hydrogen is transported into the high altitude regions. However, the effect on the thermal structure is negligible: temperature differences are of a few degrees Kelvin or less throughout the thermosphere. This justifies our decision to run the models with composition transport switched off.
To investigate the effect of different viscosities and thermal conductivities, we have repeated the 6TW/40nb run with the viscosity and thermal conductivities artificially modified. We have increased and decreased both parameters by a factor of 2. The effect of this is illustrated in Fig. 9. The increasing (decreasing) thermal conduction simply results in globally lower (higher) temperatures, as expected. The effect of altering the viscosity is more interesting. We expect a higher (lower) viscosity to improve (inhibit) the redistribution of energy. This is the case: the temperature at mid-latitudes increases (decreases) by ∼10-20 K. However, the effect on the auroral temperature is much stronger. It decreases (increases) by about 50 K. Thus the strongest effect that the viscosity has on the redistribution is to mediate the degree to which winds remove energy from the auroral regions. It has a weaker effect on the temperatures at mid-latitudes.
Finally we have investigated the effect of including eddy diffusion processes in our calculations of thermal conductivity and viscosity, applying a uniform eddy coefficient K=1.0×10 3 m 2 s −1 as discussed by Müller-Wodarg et al. (2005) 3 . The result of this is shown in Fig. 9 by the point marked with the dot-dash line. The resultant temperature distribution is uniformly cooler, but remains close to the locus of 40 nb runs (the solid line). Thus the effect of including eddy processes seems to be an increase in the total energy required, but not in the distribution of that energy. The simple reason for this is that the eddy contribution to thermal conduction and viscosity is only important at very low altitude. The effect of the eddy viscosity is negligible, while the eddy conduction only has an influence very close to the mesopause. The result is that the eddy conduction simply makes the base of the thermosphere colder, and the temperatures throughout the rest of the thermosphere are uniformly shifted accordingly.

Angular momentum transport
Another perspective on the importance of viscosity is apparent if we consider the conservation of angular momentum. In our rotating frame, the axial angular momentum of a parcel of gas consists of two components: its angular momentum due to the rotation of the planet, and that due to its own velocity relative to the planet, i.e. that due to the zonal winds. In the absence of viscosity, the sum of these components must be conserved (see Read (1986) for a detailed discussion). To move gas from the polar regions (where the angular momentum due to the rotation of the planet is small) to the equatorial regions (where the angular momentum due to the rotation of the planet is much greater) would thus require a reduction in the angular momentum due to the zonal winds. This implies the generation of fast westwards winds -comparable in magnitude to the rotation velocity of the planet -to keep the total angular momentum constant. In practice there is considerable viscosity in the thermosphere, and this supplies extra angular momentum from the lower atmosphere, allowing the gas to move equatorwards without the generation of unexpectedly fast zonal winds. Thus a higher viscosity leads to improved energy redistribution because it is able to supply more angular momentum to the thermosphere.
This analysis motivates us to consider two other possible mechanisms for supplying angular momentum to the thermosphere: buoyancy waves from the lower atmosphere and ion drag. The influence of buoyancy waves on the upper atmospheric flows are complicated at Earth and, given the very different nature of the deep atmosphere at Saturn, they are likely to have a different, and unknown, effect upon the upper atmosphere. It is possible that although supplying little energy to the upper atmosphere, buoyancy waves are able to provide a zonal drag that improves redistribution from the polar regions.
The effects of ion drag are slightly better understood. It is established that the inner magnetosphere almost corotates with Saturn (Cowley et al., 2004). It is driven into corotation by forces exerted by magnetosphere-ionosphere coupling currents. These currents transfer angular momentum from the conducting layer of the ionosphere to the magnetosphere. One would expect the same transfer of angular momentum to occur to high altitude regions of the thermosphere that also subcorotate. Thus the existence of a significant Pedersen conductivity at high altitudes may allow the corotating lower thermosphere, where the conductivity is high, to transfer angular momentum to higher altitudes, where the conductivity is lower, but non-negligible. We hope to address this effect in more detail in the future.

Electrodynamics
One possibility for explaining the high low-latitude temperatures may lie in the electrodynamics of the equatorial region. The Earth's equatorial region ionosphere is complex, due to the juxtaposition of a horizontal magnetic field and neutral wind-induced electric fields. The electron density is a minimum at the magnetic equator, and maximises each side of it due to the "fountain effect" (Ratcliffe, 1972). Also there is a greatly enhanced conductivity in the E-region leading to the equatorial electrojet. Solving for these effects selfconsistently is complex due to the global nature of the electrodynamical coupling (see e.g. Millward et al., 2001). The extent of these effects on Saturn is not known, although, like the Earth, Saturn has a strong equatorial magnetic field. It is possible there is substantial equatorial heating from electrojet currents. However, there is no a priori reason to think this might be the case, as the wind-induced electric fields on Saturn are not well-constrained. In the long-term, this will be investigated in a later version of our model.

Conclusions
A series of simulations using a basic GCM for the thermosphere of Saturn has looked at the effect of a general polar heating input. We find that for our chosen heating distributions an input of ∼6 TW peaking at 40 nb can explain both the Smith et al. (1983) determination of temperature at 30 • N and auroral H + 3 temperatures (Miller et al., 2000;Melin et al., 2005 2 ). However, it is not possible to find a reasonable energy input distribution which is consistent with both the auroral measurements and the Festou and Atreya (1982) determination for 4 • N. By seeking to match both auroral and midlatitude constraints we have been able to constrain the possible heating distributions more narrowly than would have been possible with one constraint only.
We have also indicated the sensitivity of the thermal structure to the topside heating profile, and shown that the latitudinal structure of the polar heating is less important than the total heating powers and vertical heating profile for determining the mid-latitude thermal structure.
It will be necessary in future to study the effects of specific heating sources in more detail, in particular the effects of ion drag on the dynamics and energetics at high latitudes.