Modelling H₃⁺ in planetary atmospheres: effects of vertical gradients on observed quantities

(a) Modelling overview

The majority of H₃⁺ ions in giant planet ionospheres are in photochemical equilibrium (PCE), as the H₃⁺ chemical lifetime is much shorter than the transport time scale at most altitudes, and therefore the ion continuity equation simplifies to equating local production and loss (i.e. P_s = L_s) [25,26]. While transport processes are still relevant—and highly so for H⁺, especially at high altitude—the dominance of chemical loss at low altitudes justifies the use of one-dimensional (1D) simulations over those regions, which offer the advantages of simplicity and computational freedom over more dynamically comprehensive three-dimensional simulations.

Ionization fractions at the giant planets are roughly of order 10⁻⁶ [27], indicating that ion chemistry and dynamics are largely inconsequential for the underlying neutral atmosphere. However, this generality is less true at auroral latitudes, where there appears to be a correspondence between auroral emission signatures with complex hydrocarbons and stratospheric hazes [28,29], and where H₃⁺ can act as a thermostat, helping to maintain a cooler thermosphere [30] and limit atmospheric escape.

Two models are adopted in the present work, which is focused on non-auroral latitudes at Jupiter and Uranus. The simulations are conducted in one dimension, owing to the prevalence of PCE for H₃⁺ distributions, and there are separate neutral and plasma modules in order to enable a more computationally efficient exploration of ion chemistry.

The neutral module is described in detail by Moses & Poppe [31], and is actually a combination of a meteoroid ablation code [32,33] with the Caltech/JPL 1-D KINETICS photochemical model [34,35]. KINETICS solves the coupled mass-continuity equations as a function of pressure, and (for the neutral species) includes molecular and eddy diffusion transport terms. It has proved to be effective and highly adaptable, having been applied to all of the giant planets, and currently treats 70 hydrocarbon and oxygen species that interact via approximately 500 recently updated chemical reactions [31,36]. Input for the meteoroid ablation code follows from revised constraints on interplanetary dust fluxes in the outer Solar System based on in situ spacecraft data [37]. The resulting oxygenated and hydrocarbon mixing ratios are in agreement with a wide range of observational constraints [31]. Therefore, after adjusting the KINETICS simulations for the solar and geometric conditions explored here, the resulting neutral atmospheres serve as an excellent background for exploring realistic ion-neutral photochemistry at the giant planets.

Plasma densities and temperatures follow from another 1D model called BU1DIM (the Boston University 1-D Ionosphere Model). BU1DIM was originally developed for Saturn [38,39], though has since been applied to Earth [40] and Mars [41]. Its most recent iteration has been generalized for application to any planetary atmosphere, and includes significantly expanded chemistry [25]. BU1DIM describes the time- and altitude-dependent structure of an ionosphere by solving the coupled continuity, momentum and energy equations for all ion species of interest. Jupiter's magnetic field is specified using results from the Juno spacecraft [42]. At Uranus, however, magnetic field measurements are limited to a single flyby [43]. The primary effect of magnetic fields on 1D ionospheric calculations is to constrain the plasma motion (e.g. introducing a sin² I term into the expression for vertical ion drift velocity, where I is the magnetic dip angle [38,44]). Therefore—partly due to incomplete knowledge of Uranus' magnetic field, and partly due to the predominance of PCE at H₃⁺ altitudes—magnetic field lines at Uranus are considered to be vertical here in order to focus investigations on the effect of vertical thermospheric gradients on derived H₃⁺ parameters. Modelled ion production rates follow from the attenuation of solar Extreme UltraViolet (EUV; 10–121 nm) and soft X-ray photons (combined, the XUV) [45], which are extrapolated to Jupiter and Uranus based on measurements from the Thermosphere Ionosphere Mesosphere Energetics and Dynamics Solar EUV Experiment (TIMED/SEE) [46]. In addition, secondary ionization and thermal electron heating rates are specified using parameterizations derived from coupled electron transport calculations at Saturn [47]. Aside from solar XUV radiation, no other sources of energy input are considered here (e.g. energetic particle precipitation).

Early theoretical models of giant planet ionospheres predicted electron densities that were up to an order of magnitude too large based on later spacecraft measurements, with Saturn exhibiting the most extreme discrepancy [27,48]. One commonly adopted mechanism for reducing modelled electron densities in order to better match observations was to convert H⁺ into a molecular ion via the reaction

Without the introduction of some form of ion-neutral charge-exchange reaction, such as (2.1), modelled H⁺—and hence electron density, n_e—is unrealistically large, as the radiative recombination rate coefficient for H⁺ is extremely slow (approx. 10⁻¹² cm³ s⁻¹ for typical giant planet thermospheric electron temperatures) [49]. The (2.1) reaction rate is thought to be near its maximum kinetic value [50,51], however, the fraction of molecular hydrogen in the fourth or higher vibrational state is not constrained by observations at present. For Jupiter, we adopt the vibrational density results from calculations by Majeed et al. [52], which lead to an effective H₂ vibrational rate coefficient in combination with [51]. For Uranus, as two of the dominant sources of vibrationally excited H₂ have been shown to be photon-induced fluorescence and dissociative recombination of H₃⁺ ions [52]—two solar-driven processes—we scale the fractional H₂ vibrational populations for Jupiter by (1/r²) to account for the diminution of solar photons with distance. Thus adjusted, the Majeed et al. results are then interpolated onto the appropriate Uranus pressure grid. Further model inputs and specific settings are discussed in relation to their corresponding results in §3.

(b) Observations and data reduction

While this is primarily a modelling study, there are two primary sources of data used to constrain the model results at Jupiter. First, the Galileo G0N radio occultation, obtained on 8 December 1995, sampled Jupiter's dusk ionosphere near 24° S latitude and 292° E longitude [53,54]. This measurement provides a representative ionospheric electron density profile suitable for demonstrating the effect of vertical temperature gradients on retrieved H₃⁺ parameters, at least when combined with model simulations that reproduce both the Galileo electron density profile and subsequent H₃⁺ column density observations.

Such H₃⁺ column densities are the second source of data in this study: ground-based, high-resolution spectroscopic observations of Jupiter in the L telluric window. Jupiter high-resolution (R∼25 000) spectroscopic data spanning 3.26–4 μm were obtained over four nights in April 2016 (14, 16, 20 and 23) using the Near InfraRed Spectrograph (NIRSPEC [55]) on the 10 m Keck II telescope, as described in Moore et al. [56]. Combined, these observations yielded global coverage of Jupiter, with overlapping data from 2+ nights for more than half of the planet. For the current work, a spectrum corresponding to the G0N latitude and local time was extracted and reduced as described in Moore et al. [56]. Briefly, using the line list of Neale et al. [57] and the partition function and total emission formulation of Miller et al. [24], assuming conditions of q-LTE, Gaussian line-fitting techniques [58] are used to fit the modelled spectra to observed H₃⁺ R- and Q-branch lines in order to retrieve column-averaged vibrational temperatures and column-integrated densities [59].

(a) Effect of vertical atmospheric gradients: H₃⁺ density

We first examine the effect of vertical gradients using a series of simplified, synthetic slab atmospheres. A column-integrated line-of-sight will transect the H₃⁺ layer in a planetary ionosphere at some oblique angle—or along the zenith for observations from directly overhead. Giant planet H₃⁺ emissions are optically thin [60], and lie in a spectral region with strong methane absorption (near 3.4 μm) [61]. The intensity of the H₃⁺ spectrum in each slab scales linearly with density and exponentially with temperature [30]. Furthermore, no self-absorption is considered between slabs. Thus, the observed spectrum follows from the net sum of each individual ‘slab’ emission element within the ionosphere, where the slab emission depends on the column density and temperature of a slab. The retrieved H₃⁺ ‘temperature’ will only represent the ‘true’ atmospheric temperature for an infinitely thin or isothermal H₃⁺ layer. In all other, more realistic cases, the derived temperature and density are thus necessarily a result of the convolution of the altitude structure in both density and temperature.

Based on previous modelling results [62,63], we begin by approximating Jupiter's H₃⁺ layer to have an equivalent slab width of approximately 1000 km in altitude, as the calculated H₃⁺ density is a nearly constant ∼5000 cm⁻³ over this altitude range [62]. This leads to an implied column density of ∼0.5 × 10¹⁶ m⁻². Over this same region, Jupiter's temperature increases from—very approximately—525 K to 775 K [12]. The temperature gradient is thus roughly 0.25 K km⁻¹, and the mean temperature of the H₃⁺ layer, T_mean, is approximately 650 K. We then divide this synthetic H₃⁺ layer into an arbitrary number n slabs, each with column density N_slab = 0.5 × 10¹⁶/n, a temperature T_slab based on the temperature gradient, and a corresponding slab emission. We choose n = 10 for this example, as fewer slabs are better-represented graphically, but discuss other variations below. This general atmospheric structure is illustrated in figure 1.

• Download figure
• Open in new tab
• Download PowerPoint

The composite H₃⁺ spectrum based on the structure from figure 1 will clearly be dominated by higher altitude (redder, hotter) slab spectra. By fitting this spectrum as described in §2b, we retrieve a column density and column-integrated temperature of 0.4 × 10¹⁶ m⁻² and 695 K, respectively. These values demonstrate an important complication with using observed H₃⁺ spectra to derive realistic column densities, as the retrieved density from the simulated spectral fit is only 80% of the true density. This discrepancy is due to the exponential relation between H₃⁺ temperature and emission, combined with the standard assumption that the ‘H₃⁺ layer’ is isothermal, and this weighting is also reflected in the retrieved temperature, 695 K, which is approximately 7% higher than the true mean temperature represented in figure 1 (650 K).

A more realistic treatment of the simplified situation depicted above would complicate matters further due to the problems common to all astronomical observations, such as contaminations due to other emissions and absorptions, detector noise, and so forth. Before stressing this initial (N_fit/N_true) = 0.8 result too strongly, however, it is important to investigate the various sensitivities that lead to such a discrepancy in the column density. For example, while the conditions in figure 1 have been chosen to roughly represent Jupiter's ionosphere, the mean temperature and the temperature gradient will be different at other locations at Jupiter and at other planets. Similarly, we might question what effect a more realistic density profile would have. The rest of this subsection is therefore devoted to outlining how different choices in representing the atmosphere and in generating the synthetic spectra affect the retrieved H₃⁺ column density.

First, we examine the choice of the number of slabs n on the derived column density. As demonstrated in figure 2, increasing the number of slabs leads to an asymptotic approach towards (N_fit/N_true) = 0.822, which is a approximately 3% increase over the 10 slab representation imagined in figure 1. There is a corresponding inverted trend in the total H₃⁺ emission with slab number, as representations with fewer slabs associate a wider range of the H₃⁺ layer with higher temperatures. In the example considered for figure 2, the net emission decreases by approximately 5% towards an asymptotic value of 21.35 μW m² sr⁻¹. For comparison, typical giant planet models blanket the ionosphere with fewer than 100 grid points. The derived column density will also depend on the mean temperature in the atmosphere, however. Figure 3 explores this effect while keeping the temperature gradient and the number of slabs fixed (i.e. the temperature gradient is as shown in figure 1, 0.25 K km⁻¹, and the mean temperature is varied).

• Download figure
• Open in new tab
• Download PowerPoint

• Download figure
• Open in new tab
• Download PowerPoint

Based on figures 2 and 3, derived H₃⁺ column density is highly dependent upon the atmospheric temperature profile in the H₃⁺ layer. Before exploring a wider range of temperature gradients, such as might be more widely representative of H₃⁺ in giant planet ionospheres, we investigate the additional impact of density gradients. (It should be noted that the N_fit/N_true ratio is not sensitive to the true slab column density, which is perhaps not surprising given the linear relation between density and emission and the fractional error introduced from temperature gradients.) Whereas the temperature profile is generally monotonically increasing in the lower thermosphere, there will be both positive and negative gradients in H₃⁺ density.

Figure 4 presents the combined effects of both density and temperature gradients in the H₃⁺ layer on simulated retrieved column densities, N_fit. These results follow from an atmospheric layer divided into n = 36 slabs, with a true slab column density N_true = 0.5 × 10¹⁶ m⁻² and mean temperatures T_mean of (a) 450 K, (b) 650 K and (c) 850 K. These calculations follow the approach outlined in figure 1, except the absolute H₃⁺ temperatures are allowed to range only between 150 and 1250 K, temperatures observed to be appropriate for giant planet upper atmospheres. This requirement means that, for the largest temperature gradients explored, there are regions of the H₃⁺ layer that are isothermal (i.e. at either 150 K or 1250 K). As expected due to the linear relationship between H₃⁺ emission and density, there is no discrepancy between modelled and ‘observed’ column density when there is no temperature gradient. The generally vertical contours in figure 4 indicate that (N_fit/N_true) is most sensitive to temperature gradients, though density gradients begin to play an important role when the temperature gradient is larger than approximately 0.5 K km⁻¹. For realistic temperature gradients at the giant planets, approximately 0.4–2 K km⁻¹ [12,18], the retrieved H₃⁺ column density ranges from approximately 20–90% of the true value, depending primarily on the mean temperature within the H₃⁺ layer. Observed uncertainties in N_fit are often comparable, so the effect illustrated in figure 4 might be absorbed into experimental errors for low S/N data. For instance, in a study of the anticorrelation between H₃⁺ density and temperature, Melin et al. [59] find that a S/N of 12 is required to achieve a 10% column density uncertainty for a column-averaged H₃⁺ temperature of 600 K.

• Download figure
• Open in new tab
• Download PowerPoint

A more realistic H₃⁺ layer would experience variations in both temperature and density gradients with altitude, and so could not be represented as simply as in figure 4. Careful examination of the evolution of dT/dz and dn/dz with altitude could give some idea of the net error induced in N_fit based on those gradients, though such an analysis would also rely on prior knowledge of the atmosphere, and so significantly reduce the value of added observation. Therefore, while the preceding figures establish that, unless the H₃⁺ layer is in an isothermal atmosphere, the H₃⁺ column density retrieved from observations will represent a lower limit, it is not practical at this stage to examine an all-inclusive range of possible atmospheric structures in order to assign definitive quantitative values to those lower limits. Instead, these results serve as further motivation for investigating other similar complications of interpreting H₃⁺ observations, and for outlining ‘toy’ model parameters that would be relevant for development of a full H₃⁺ retrieval model (e.g. such as in [64]).

(b) Effect of vertical atmospheric gradients: H₃⁺ temperature (Uranus)

While density and temperature structures within a planet's H₃⁺ layer affect the retrieved column density, they also affect the interpretation of the retrieved column-averaged temperature. In order to demonstrate this effect, we model the global distribution of H₃⁺ at Uranus. First, the background atmosphere is based on Moses & Poppe [31], appropriate for globally averaged conditions with dust-derived oxygen influxes of 1.2 × 10⁵ H₂O molecules cm⁻² s⁻¹, 2.5 × 10⁵ CO molecules cm⁻² s⁻¹ and 3.0 × 10³ CO₂ molecules cm⁻² s⁻¹, consistent with H₂O, CO and CO₂ observations [65–67]. This 1D atmosphere is applied uniformly at Uranus. While clearly not fully realistic, using a fixed neutral atmosphere, where ion and neutral chemistry is not fully coupled, allows for clearer elucidation of the effects of a varying H₃⁺ layer on retrieved temperatures, as will be demonstrated below. Next, a simulation date of 15 September 2017 is chosen. This choice is largely arbitrary for the purposes of demonstrating the effects of gradients on retrieved H₃⁺ temperatures; however, there do happen to be contemporaneous H₃⁺ observations from September 2017 [68], for which Uranus' sub-solar latitude was −38°. Finally, 1D ionospheric calculations are performed globally as described in §2a, with a 1° latitude resolution. Profiles of background neutral density, temperature and ion density at 30° S latitude are shown in figure 5. Electron temperatures (not shown) are calculated to diverge from the neutral temperature around 2000 km altitude and reach approximately 800 K at 30° S latitude, 12 solar local time (SLT). (Note that altitude levels throughout the text are referenced to the 1 bar pressure level.) Calculated H₃⁺ temperatures are found to be equal to the background neutral temperature.

• Download figure
• Open in new tab
• Download PowerPoint

Figure 6 presents global ionospheric density results at Uranus, plotted versus SLT and planetocentric latitude. Note that, as H₃⁺ temperatures were found to be identical to the neutral temperature, at least for the conditions in figure 5, these global simulations do not include plasma temperature calculations, and therefore any H₃⁺ temperature variations are a reflection of the vertical distribution of H₃⁺ and the background temperature profile. Following the preceding approach, N_fit and T_fit follow from generating and fitting synthetic H₃⁺ spectra based on modelled H₃⁺ distributions. These simulations use 138 H₃⁺ slabs (i.e. grid points in altitude), and therefore N_fit is within approximately 0.2% of its asymptotic value based on figure 2.

• Download figure
• Open in new tab
• Download PowerPoint

While direct reproduction wasn't the goal of these simulations—and would be at least partially coincidental anyway due to the homogeneous neutral atmosphere—calculated column densities are broadly consistent with H₃⁺ observations [6,68]. As expected from §3a, the true modelled H₃⁺ column density, figure 6a, is slightly larger than the retrieved value, figure 6b, with a global mean column density ratio, figure 6c, of 0.87. The column density ratio is nearest to 1.0 at dawn, and at high northern latitudes, indicating that more of the H₃⁺ layer there is in the isothermal region of the model atmosphere above 2800 km altitude, as would be expected based on solar zenith angle (SZA) effects and nightside recombination chemistry, which depletes lower-altitude layers near the electron density peak more rapidly [69]. This is more evident from figure 6d, which reveals that the column-averaged H₃⁺ temperature is higher in those regions, a direct response of the peak altitude of the H₃⁺ layer being shifted towards higher altitudes, as seen in figure 6e. This is primarily an SZA effect, as more oblique slant paths through the atmosphere generate higher-altitude photoionization; however, there is also a slight offset post-noon due partly to conversion of H⁺ to H₂⁺ (and thus, H₃⁺, following reaction with H₂) via reaction (2.1). The H₃⁺ lifetime at the H₃⁺ peak is on the order of ∼1 h for most of the Uranus day, as would be expected based on the modelled peak density (approx. 3000 cm⁻³) and dissociative recombination rate of H₃⁺ with electrons (approx. 10⁻⁷ cm³ s⁻¹ [70]). Column-averaged H₃⁺ lifetimes, weighted by H₃⁺ density, are slightly higher than at the peak altitude, but still typically less than 3 h (figure 6f ).

Owing to the simplified nature of the Uranus simulations, and owing to the sparse thermospheric constraints at present, a detailed model-data comparison is not warranted here. Instead, we emphasize the temperature variations shown in figure 6d. Despite the fact that the thermospheric temperature profile is identical at all latitudes, retrieved H₃⁺ temperatures vary by greater than 35 K, or roughly 5% of the mean temperature. This is simply understood as primarily a SZA effect: at low SZAs the H₃⁺ layer is lower in the ionosphere, probing the lower temperatures there, whereas the reverse is true for high SZAs. This result is, again, not surprising; however, it highlights two important points: (1) observed H₃⁺ temperature variations do not necessarily imply anything about the energetics of the thermosphere, and (2) these SZA effects should be accounted for when interpreting measured temperature variabilities. Finally, in reviewing the results of figure 6, it is important to also emphasize the elements that are missing from the simulations presented there. In particular, we have neglected energetic particle precipitation, which would be expected to lead to increased ionization and enhanced thermospheric temperatures, mainly at high magnetic latitudes. At Uranus, the magnetic polar regions are at mid- and low-latitude as a result of the tilted magnetic dipole axis [71], though their exact longitudes are unknown at present due to uncertainty in Uranus' rotation period. Thus, the chief value of figure 6 is in demonstrating the qualitative effects of H₃⁺ density and temperature gradients on retrieved parameters. More realistic global variations of ionization and heating at Uranus would be expected to lead to different quantitative structures.

(c) Vertical structure of H₃⁺ density and temperature (Jupiter)

One of the most sought-after observables for planetary atmospheres is the thermal structure, as so much of the rest of planetary dynamics and energetics depend upon it. Furthermore, as established in §3a,b, atmospheric gradients in the H₃⁺ layer can confuse measurements of H₃⁺ density and temperature, meaning the very parameter we want to constrain, T(z), is itself limiting observational insight. We now introduce a potential method of unravelling this confusion by combining column-integrated H₃⁺ measurements with forward modelling.

For this proof-of-concept, we use the Galileo G0N radio occultation, obtained on 8 December 1995, which sampled Jupiter's dusk ionosphere near 24° S planetocentric latitude and 292° E longitude [54]. Figure 7 presents a model reproduction of the G0N occultation, along with corresponding background neutral atmospheric parameters and modelled ion densities. The model simulation is for 24° S latitude, with a forced vertical drift W_D of 50 cm s⁻¹, and plasma density comparisons are extracted for 18 SLT in accordance with the dusk terminator measurement of G0N. As described in §2a, the solar flux is specified using extrapolated TIMED/SEE measurements, the secondary ionization and photoelectron heating rate are parameterized [47], and the effective (2.1) reaction rate comes from Majeed et al. [52], with the rest of the chemistry as specified in Moore et al. [25].

• Download figure
• Open in new tab
• Download PowerPoint

The good model-data agreement in figure 7 gives some confidence that the broad electron density features are well represented (that is, aside from the noisy densities in the topside and the narrow structures below the peak, which might be attributed to gravity waves [72]). As an additional test of the model simulation, we turn to H₃⁺ observations obtained in April 2016 using Keck/NIRSPEC (§2b). First, a series of spectra obtained with the NIRSPEC slit oriented E-W near Jupiter's Great Red Spot are combined and reduced in order to obtain a calibrated H₃⁺ spectrum at 24° S latitude, 17 SLT. Second, we derive a spectral fit to the data, and compare the fit parameters to the H₃⁺ column density from the ionospheric simulation shown in figure 7. An extracted portion of the H₃⁺ spectrum, its corresponding spectral fit, and the modelled column densities are shown in figure 8. Retrieved parameters are T_fit = 774 ± 64 K and N_fit = (2.31 ± 0.88) × 10¹⁵ m⁻², and the latter is over plotted on the modelled column densities with horizontal error bars that identify the range of local times that contributed to the observed spectrum.

• Download figure
• Open in new tab
• Download PowerPoint

There is good model-data agreement in figure 8b, which lends additional confidence that the model is an accurate representation of both the observed electron density (figure 7) and H₃⁺ column density (figure 8b). Before moving on, however, there are a couple of important caveats to emphasize. First, and perhaps most important, the ionospheric simulation is for 20 April 2016 during solar minimum, in accordance with the Keck H₃⁺ observations, whereas the Galileo radio occultation was obtained on 8 December 1995, nearly 21 years prior, also during solar minimum. It would be far more surprising if Jupiter's ionosphere had not changed in those intervening years than if it had, especially given the variability present in ionospheric radio occultations at giant planets [53,73]. Therefore, the fact that the model agreement is good in both figures 7 and 8 is most likely a coincidence rather than an indication of atmospheric stability, though the fact that the approximately 21 year separation between the two datasets is nearly 2 full solar cycles allows some minimum of hope for a happy coincidence to be maintained. Second, the Galileo radio occultation was at 68° W (System III) longitude, whereas the Keck observations were centred at 308° W longitude, a slightly different magnetic environment. The majority of the modelled H₃⁺ layer is still in photochemical equilibrium, meaning that this variation should have minimal effect, at least if solar photons are the main ionospheric driver as assumed here, though high altitude ion drifts would be altered. Nevertheless, given that the model is able to reproduce both available datasets, and given that there are no better combinations of ionospheric constraints available, we shall progress forward under the assumption that the ionospheric model simulation provides as accurate a representation of the H₃⁺ density structure as possible at present.

To re-state the problem: the observed H₃⁺ spectrum is a function of the integrated density N(z) and temperature T(z) profiles from the emitting H₃⁺ layer. Thus, the IR spectrum I(λ) also contains information about both of them. Essentially, there are three unknowns, so if either N(z) or T(z) can be convincingly constrained then the other can in principle be derived when combined with the observed spectrum. The spectrum in this example is known from Keck/NIRSPEC observations. Based on the good model-data agreement for the electron density profile and the H₃⁺ column density, we proceed under the assumption that N(z) is appropriately constrained. Therefore, we should also be able to reconstruct at least some limited representation of T(z) from the above two inputs.

First, as in §3a, the modelled H₃⁺ layer is idealized as n slabs, each of column density N_slab = N_true/n. Second, a temperature is assigned to each slab, randomly selected from 150 to 1250 K, and a synthetic H₃⁺ spectrum (i.e. the sum of the slab spectra) is computed. In principle, the specified temperatures of each slab could be independent of each other, but for the present case a monotonically increasing (or isothermal) temperature profile is enforced, consistent with observation [53,74]. This step is repeated until the modelled and observed spectra converge to within some pre-defined tolerance. In this case, the tolerance is set to a maximum of 0.03% disagreement between the slab-modelled spectrum and the spectrum obtained from the H₃⁺ fit. Once converged, it is useful to also provide some estimate of the sensitivity of the result. For this purpose, n-1 slab temperatures are held fixed to their converged values while the other is freely varied until the modelled column-averaged temperature exceeds the measured temperature uncertainty. The derived temperature uncertainties thus represent ∼8% errors in this case, as T_fit = 774 ± 64 K. Finally, based on the converged slab temperatures weighted by their uncertainties, an analytical temperature profile is derived following [75]:

where z is the altitude element, T_exo the exospheric temperature and A_i are constants of the fit.

Results from the above process with n = 4 are shown in figure 9. Each slab is represented by a shaded column, with varying vertical extent due to enforced equal slab column densities and variable slab number densities. The original modelled H₃⁺ density is given by the red curve. These variable slab widths enable higher altitude resolution near the H₃⁺ density peak and allow for more reasonable error bars than significantly smaller slab widths would. Corresponding slab temperatures are shown as black circles on the right, along with the derived temperature profile in red. Horizontal grey error bars indicate the 1σ uncertainty in slab temperature, and vertical grey error bars demarcate the vertical extent of each slab. The best-fit parameters for the temperature profile in figure 9 are T_exo = 788 K, A₁ = 136 K, A₂ = 159 km, and A₃ = 121 km² K⁻¹, and are only representative of the altitude range with significant H₃⁺ density (e.g. between ∼300 and 2500 km).

• Download figure
• Open in new tab
• Download PowerPoint

Combining the results of figure 9 with those from §3a, it appears that the Keck/NIRSPEC observations were minimally affected by gradients in the H₃⁺ layer. The derived thermal gradient is effectively zero above 1050 km (i.e. less than or equal to 0.05 K km⁻¹), less than 0.2 K km⁻¹ for altitudes between 750 and 1050 km, approximately 0.3 − 1.4 between 550 and 750 km, and approximately 1.8 K km⁻¹ at the bottom side. In situ measurements at Jupiter [12] and Saturn [18] find approximately 2 K km⁻¹ near 400 km and 0.4 K km⁻¹ at the base of the thermosphere, respectively. Based on figure 4, this implies that retrieved H₃⁺ column densities are less than 80% of the modelled values only for altitudes less than 750 km, a range that encompasses 12% of the total column density. Meanwhile, absolute modelled density gradients are less than 8 cm⁻³ km⁻¹ everywhere, and only approximately 3 cm⁻³ km⁻¹ where there is also a substantial temperature gradient (near 400 km altitude). Progression of the impact of the derived density and temperature gradients on retrieved column density is shown by the dashed gray curve in figure 4b. In total, the calculated error in column density, based on the combination of figure 4b and the density structure from figure 9, is approximately 10%, all associated with gradients at altitudes below 800 km. This error is comparable to the observational uncertainties (figure 8).