Retrieving Neptune's aerosol properties from Keck OSIRIS observations. I. Dark regions

TL;DR

This study uses high-resolution spectroscopic data from Keck OSIRIS to analyze Neptune's dark atmospheric regions, revealing details about aerosol layers and methane distribution, with implications for understanding Neptune's atmospheric composition and structure.

Contribution

First detailed analysis of Neptune's aerosol and methane profiles in dark regions using 3D data cubes and advanced modeling techniques.

Findings

Neptune's cloud opacity is dominated by a compact, optically thick cloud near 3 bar.

A vertically extended haze of small particles exists above the deep cloud.

Aerosol optical depth decreases significantly from low to high southern latitudes.

Abstract

We present and analyze three-dimensional data cubes of Neptune from the OSIRIS integral-field spectrograph on the 10-m Keck telescope, from July 2009. These data have a spatial resolution of 0.035"/pixel and spectral resolution of R~3800 in the H and K broad bands. We focus our analysis on regions of Neptune's atmosphere that are near-infrared dark- that is, free of discrete bright cloud features. We use a forward model coupled to a Markov chain Monte Carlo algorithm to retrieve properties of Neptune's aerosol structure and methane profile above ~4 bar in these near-infrared dark regions. Using a set of high signal-to-noise spectra in a cloud-free band from 2-12N, we find that Neptune's cloud opacity is dominated by a compact, optically thick cloud layer with a base near 3 bar and composed of low albedo, forward scattering particles, with an assumed characteristic size of ~1$\mu$m. Above this cloud, we require a vertically extended haze of smaller (~0.1 $\mu$m) particles, which reaches from the upper troposphere (~0.6 bar) into the stratosphere. The particles in this haze are brighter and more isotropically scattering than those in the deep cloud. When we extend our analysis to 18 cloud-free locations from 20N to 87S, we observe that the optical depth in aerosols above 0.5 bar decreases by a factor of 2-3 or more at mid- and high-southern latitudes relative to low latitudes. We also consider Neptune's methane (CH$_4$) profile, and find that our retrievals indicate a strong preference for a low methane relative humidity at pressures where methane is expected to condense. Our preferred solution at most locations is for a methane relative humidity below 10% near the tropopause in addition to methane depletion down to 2.0-2.5 bar. We tentatively identify a trend of lower CH$_4$ columns above 2.5 bar at mid- and high-southern latitudes over low latitudes.

Tables

Table 1. Table 1: Overview of model parameters

	Parameter	Definition	Default value	Sections where investigated
Thermal profile
	$\delta$TP_trop	tropopause temperature offset (K)	no offset	5.3; 6.3
	$\delta$TP_strat	stratosphere temperature offset (K)	no offset	5.3
CH4 profile
	mCH4,t	deep troposphere CH4 mole fraction	0.04	5.2
	mCH4,s	upper stratosphere CH4 mole fraction	0.00035	5.2
	$R{H}_{ac}$	CH4 relative humidity (near tropopause)	1.0eeeAfter Section 5.2, we adopt a new nominal of $R{H}_{ac}=0.4$	5.2; 6.2
	$P_{dep}$	CH4 depletion pressure depth (bar)	no depletion	6.2
Aerosolsfffparameters may be varied in one or all aerosol layers
	Pmax	maximum (bottom) pressure (bar)	free	4; 5; 6
	Pmin	minimum (top) pressure (bar)	${10}^{-5}$	4.1
	hfrac	fractional scale height	free	4; 5; 6.1; 6.3
	$\tau$	optical depth at 1.6 $\mu$m	free	4; 5; 6
	$r_p$	particle radius ($\mu$m)gggradius of particles at the peak of the particle size distribution, which is described in Appendix A	0.1, 0.5 or 1.0	4.2
	$g$	phase function asymmetry parameter	free	4; 5.1
	${\omega }_0$	single scattering albedo	free	4; 5.1

Table 2. Table 2: Fit results, single aerosol layer ( 1L ) models

run name	$P_{max}$ (bar)	$P_{min}$ (bar)	$\tau$	$h_{frac}$	${\omega }_0$	$g$	$\log {\sigma }_t$	$DIC$
1L_nom	${10}_{-3}^{+5}$	$0.001$	$9_{-6}^{+20}$	${0.478}_{-0.007}^{+0.007}$	${0.41}_{-0.01}^{+0.01}$	${0.34}_{-0.03}^{+0.02}$	$-{6.48}_{-0.02}^{+0.02}$	$-25641$
1L_pmin1.4	${10}_{-3}^{+5}$	1.4hhhfixed at the Karkoschka and Tomasko, (2011) value.	${11}_{-7}^{+20}$	${0.471}_{-0.009}^{+0.009}$	${0.40}_{-0.01}^{+0.01}$	${0.33}_{-0.03}^{+0.03}$	$-{6.44}_{-0.01}^{+0.02}$	$-25430$
1L_pminfree	$9_{-3}^{+5}$	${0.007}_{-0.005}^{+0.01}$	$8_{-5}^{+20}$	${0.478}_{-0.007}^{+0.007}$	${0.41}_{-0.01}^{+0.01}$	${0.34}_{-0.03}^{+0.02}$	$-{6.48}_{-0.01}^{+0.02}$	$-25642$

Table 3. Table 3: Results of two aerosol layer retrievals from limb darkening dataset.

run name

$r_{p,\alpha }$ ($\mu$m)

layer $\alpha$ (top)

layer $\beta$ (bottom)

$\log {\sigma }_t$

$DIC$

$P_{max,\alpha }$

${\tau }_{\alpha }$

$h_{frac,\alpha }$

${\omega }_{0,\alpha }$

$g_{\alpha }$

$P_{max,\beta }$

${\tau }_{\beta }$

$h_{frac,\beta }$

${\omega }_{0,\beta }$

$g_{\beta }$

(bar)

(bar)

2L_hazeA

0.1

${0.54}_{-0.04}^{+0.03}$

${0.102}_{-0.005}^{+0.005}$

${0.41}_{-0.09}^{+0.07}$

${0.24}_{-0.03}^{+0.02}$

${0.37}_{-0.1}^{+0.06}$

${2.68}_{-0.06}^{+0.08}$

${0.61}_{-0.07}^{+0.08}$

${0.19}_{-0.01}^{+0.01}$

${0.28}_{-0.03}^{+0.05}$

$-{0.6}_{-0.2}^{+0.2}$

$-{6.94}_{-0.02}^{+0.02}$

$-27684$

2L_hazeB

0.5

${0.51}_{-0.01}^{+0.02}$

${0.138}_{-0.006}^{+0.005}$

${0.04}_{-0.02}^{+0.04}$

${0.16}_{-0.02}^{+0.02}$

${0.24}_{-0.09}^{+0.09}$

${2.94}_{-0.07}^{+0.1}$

${0.40}_{-0.04}^{+0.05}$

${0.26}_{-0.02}^{+0.02}$

${0.37}_{-0.03}^{+0.03}$

$-{0.8}_{-0.1}^{+0.2}$

$-{6.81}_{-0.02}^{+0.02}$

$-27142$

2L_hazeC

1.0

${0.51}_{-0.01}^{+0.03}$

${0.144}_{-0.006}^{+0.006}$

${0.04}_{-0.02}^{+0.04}$

${0.14}_{-0.02}^{+0.02}$

${0.1}_{-0.1}^{+0.1}$

${2.91}_{-0.07}^{+0.09}$

${0.33}_{-0.03}^{+0.04}$

${0.27}_{-0.02}^{+0.02}$

${0.45}_{-0.04}^{+0.04}$

$-{0.6}_{-0.2}^{+0.1}$

$-{6.82}_{-0.02}^{+0.02}$

$-27150$

2L_DISORT

0.1

${0.59}_{-0.03}^{+0.04}$

${0.019}_{-0.001}^{+0.002}$

${0.85}_{-0.06}^{+0.07}$

${0.91}_{-0.05}^{+0.06}$

${0.24}_{-0.03}^{+0.02}$

${3.3}_{-0.3}^{+0.4}$

$5_{-2}^{+4}$

${0.170}_{-0.009}^{+0.01}$

${0.45}_{-0.01}^{+0.01}$

${0.50}_{-0.02}^{+0.02}$

$-{7.03}_{-0.02}^{+0.02}$

$-28016$

Table 4. Table 4: Results of two aerosol layer retrievals from the high S/N μ = 0.8 𝜇 0.8 \mu=0.8 dataset.

run name	mCH4,t	mCH4,s	$R{H}_{ac}$	layer $\alpha$ (top)					layer $\beta$ (bottom)					$\log {\sigma }_t$	$DIC$	$\Delta DIC$iiichange in $DIC$ relative to the SN_nom model fit; other models of the SN dataset may also be compared directly using the $DIC$. A negative change of 10 or more is interpreted as an improvement in the fit quality. A positive change of $>10$ implies a worsening of the fit, whereas changes of less than $\pm 10$ are not conclusively better of worse.
				$P_{max,\alpha }$ (bar)	${\tau }_{\alpha }$	$h_{frac,\alpha }$	${\omega }_{0,\alpha }$jjj${\omega }_0$ and $g$ are retrieved in the SN_nom case and fixed near the best-fit SN_nom values in the remaining model cases. Two stream was used for all fits summarized in this table, and ${\omega }_0$ and $g$ are only appropriate for two stream at this particular viewing geometry.	$g_{\alpha }$	$P_{max,\beta }$ (bar)	${\tau }_{\beta }$	$h_{frac,\beta }$	${\omega }_{0,\beta }$	$g_{\beta }$
SN_nom	$0.040$	$3.5e-4$	$1.0$	${0.55}_{-0.07}^{+0.07}$	${0.050}_{-0.009}^{+0.01}$	${1.2}_{-0.2}^{+0.2}$	${0.23}_{-0.04}^{+0.04}$	$-{0.7}_{-0.2}^{+0.2}$	${3.1}_{-0.2}^{+0.3}$	${1.8}_{-0.5}^{+0.9}$	${0.18}_{-0.02}^{+0.02}$	${0.29}_{-0.03}^{+0.03}$	${0.0}_{-0.2}^{+0.1}$	$-{7.11}_{-0.04}^{+0.04}$	$-4800$
SN_CH4gridA	$0.022$	$3.5e-4$	$0.4$	${0.41}_{-0.07}^{+0.07}$	${0.053}_{-0.002}^{+0.002}$	${0.4}_{-0.2}^{+0.2}$	$0.2$	$-0.7$	${3.6}_{-0.2}^{+0.2}$	${1.07}_{-0.09}^{+0.1}$	${0.23}_{-0.02}^{+0.02}$	$0.3$	$0.0$	$-{7.28}_{-0.04}^{+0.04}$	$-4923$	$-123$
SN_CH4gridB	$0.022$	$3.5e-4$	$1.0$	${0.53}_{-0.06}^{+0.07}$	${0.052}_{-0.003}^{+0.003}$	${1.2}_{-0.2}^{+0.2}$	$0.2$	$-0.7$	${3.8}_{-0.3}^{+0.6}$	${1.2}_{-0.2}^{+0.5}$	${0.27}_{-0.03}^{+0.03}$	$0.3$	$0.0$	$-{7.14}_{-0.04}^{+0.04}$	$-4821$	$-21$
SN_CH4gridC	$0.022$	$1.5e-3$	$0.4$	${0.45}_{-0.07}^{+0.08}$	${0.057}_{-0.002}^{+0.002}$	${0.5}_{-0.2}^{+0.3}$	$0.2$	$-0.7$	${3.5}_{-0.2}^{+0.2}$	${1.04}_{-0.09}^{+0.1}$	${0.20}_{-0.02}^{+0.02}$	$0.3$	$0.0$	$-{7.20}_{-0.04}^{+0.04}$	$-4868$	$-68$
SN_CH4gridD	$0.022$	$1.5e-3$	$1.0$	${0.68}_{-0.08}^{+0.1}$	${0.053}_{-0.003}^{+0.004}$	${1.8}_{-0.2}^{+0.3}$	$0.2$	$-0.7$	${4.0}_{-0.4}^{+1.0}$	${1.3}_{-0.3}^{+2.0}$	${0.28}_{-0.04}^{+0.03}$	$0.3$	$0.0$	$-{7.08}_{-0.04}^{+0.04}$	$-4771$	$29$
SN_CH4gridE	$0.040$	$3.5e-4$	$0.4$	${0.41}_{-0.07}^{+0.06}$	${0.054}_{-0.001}^{+0.001}$	${0.4}_{-0.2}^{+0.2}$	$0.2$	$-0.7$	${2.94}_{-0.09}^{+0.1}$	${1.08}_{-0.09}^{+0.1}$	${0.16}_{-0.01}^{+0.01}$	$0.3$	$0.0$	$-{7.27}_{-0.04}^{+0.04}$	$-4921$	$-121$
SN_CH4gridF	$0.040$	$1.5e-3$	$0.4$	${0.45}_{-0.08}^{+0.08}$	${0.058}_{-0.002}^{+0.002}$	${0.5}_{-0.2}^{+0.3}$	$0.2$	$-0.7$	${2.9}_{-0.1}^{+0.1}$	${1.04}_{-0.09}^{+0.1}$	${0.14}_{-0.01}^{+0.01}$	$0.3$	$0.0$	$-{7.20}_{-0.04}^{+0.04}$	$-4867$	$-66$
SN_CH4gridG	$0.040$	$1.5e-3$	$1.0$	${0.67}_{-0.09}^{+0.1}$	${0.055}_{-0.003}^{+0.003}$	${1.7}_{-0.2}^{+0.2}$	$0.2$	$-0.7$	${3.1}_{-0.2}^{+0.8}$	${1.3}_{-0.2}^{+2.0}$	${0.20}_{-0.03}^{+0.02}$	$0.3$	$0.0$	$-{7.06}_{-0.04}^{+0.03}$	$-4753$	$47$
SN_CH4RHfreeA	$0.040$	$3.5e-4$	${0.016}_{-0.008}^{+0.01}$	${0.49}_{-0.07}^{+0.07}$	${0.046}_{-0.004}^{+0.002}$	${0.3}_{-0.2}^{+0.1}$	$0.2$	$-0.7$	${2.88}_{-0.07}^{+0.06}$	${0.96}_{-0.06}^{+0.06}$	${0.135}_{-0.009}^{+0.01}$	$0.3$	$0.0$	$-{7.47}_{-0.04}^{+0.04}$	$-5052$	$-252$
SN_CH4RHfreeB	$0.040$	$1.5e-3$	$8e-4_{-3e-4}^{+4e-4}$	${0.35}_{-0.05}^{+0.08}$	${0.033}_{-0.003}^{+0.003}$	${0.2}_{-0.1}^{+0.2}$	$0.2$	$-0.7$	${2.92}_{-0.05}^{+0.05}$	${0.97}_{-0.05}^{+0.05}$	${0.159}_{-0.008}^{+0.008}$	$0.3$	$0.0$	$-{7.47}_{-0.04}^{+0.04}$	$-5058$	$-258$
SN_CH4stratfree	$0.040$	$2.2e-4_{-2e-5}^{+4e-5}$	$0.4$	${0.42}_{-0.07}^{+0.06}$	${0.053}_{-0.002}^{+0.001}$	${0.4}_{-0.2}^{+0.2}$	$0.2$	$-0.7$	${3.0}_{-0.1}^{+0.1}$	${1.1}_{-0.1}^{+0.1}$	${0.16}_{-0.01}^{+0.01}$	$0.3$	$0.0$	$-{7.29}_{-0.04}^{+0.04}$	$-4928$	$-128$
SN_Ttrop	$0.040$	$3.5e-4$	$0.4$	${0.51}_{-0.05}^{+0.06}$	${0.054}_{-0.002}^{+0.002}$	${1.0}_{-0.1}^{+0.2}$	$0.2$	$-0.7$	${2.9}_{-0.1}^{+0.1}$	${1.1}_{-0.1}^{+0.1}$	${0.16}_{-0.01}^{+0.01}$	$0.3$	$0.0$	$-{7.21}_{-0.04}^{+0.04}$	$-4872$	$-72$
SN_Tstratp20	$0.040$	$3.5e-4$	$0.4$	${0.36}_{-0.04}^{+0.05}$	${0.052}_{-0.001}^{+0.001}$	${0.2}_{-0.1}^{+0.1}$	$0.2$	$-0.7$	${2.97}_{-0.09}^{+0.1}$	${1.09}_{-0.09}^{+0.1}$	${0.16}_{-0.01}^{+0.01}$	$0.3$	$0.0$	$-{7.31}_{-0.04}^{+0.04}$	$-4941$	$-141$
SN_Tstratm20	$0.040$	$3.5e-4$	$0.4$	${0.44}_{-0.07}^{+0.06}$	${0.054}_{-0.001}^{+0.001}$	${0.5}_{-0.2}^{+0.2}$	$0.2$	$-0.7$	${2.9}_{-0.1}^{+0.1}$	${1.1}_{-0.1}^{+0.1}$	${0.15}_{-0.01}^{+0.01}$	$0.3$	$0.0$	$-{7.27}_{-0.04}^{+0.04}$	$-4914$	$-114$

Table 5. Table 5: Results for initial, aerosol-only, fit to 18 cloud-free locations.

location	lat (∘)	${\mu }_H$	${\mu }_K$	layer $\alpha$ (top)			layer $\beta$ (bottom)			$\log {\sigma }_t$	$DIC$
				$P_{max,\alpha }$ (bar)	${\tau }_{\alpha }$	$h_{frac,\alpha }$	$P_{max,\beta }$ (bar)	${\tau }_{\beta }$	$h_{frac,\beta }$
A	$20$	$0.66$	$0.54$	${0.7}_{-0.1}^{+0.2}$	${0.021}_{-0.001}^{+0.001}$	${2.1}_{-0.4}^{+0.4}$	${2.3}_{-0.1}^{+0.2}$	${0.58}_{-0.07}^{+0.1}$	${0.09}_{-0.04}^{+0.04}$	$-{6.74}_{-0.07}^{+0.06}$	$-2878$
B	$16$	$0.46$	$0.63$	${1.2}_{-0.1}^{+0.2}$	${0.026}_{-0.002}^{+0.002}$	${1.2}_{-0.1}^{+0.1}$	${2.2}_{-0.1}^{+0.2}$	${0.61}_{-0.07}^{+0.1}$	${0.06}_{-0.03}^{+0.04}$	$-{6.94}_{-0.08}^{+0.08}$	$-2904$
C	$10$	$0.26$	$0.52$	${0.9}_{-0.2}^{+0.2}$	${0.020}_{-0.001}^{+0.002}$	${0.9}_{-0.2}^{+0.2}$	${3.4}_{-0.6}^{+1.0}$	$4_{-3}^{+30}$	${0.15}_{-0.02}^{+0.02}$	$-{7.07}_{-0.08}^{+0.09}$	$-2994$
D	$5$	$0.80$	$0.84$	${0.6}_{-0.1}^{+0.2}$	${0.021}_{-0.001}^{+0.001}$	${0.8}_{-0.3}^{+0.3}$	${3.4}_{-0.6}^{+0.6}$	${10}_{-10}^{+50}$	${0.12}_{-0.02}^{+0.01}$	$-{6.90}_{-0.07}^{+0.07}$	$-2941$
E	$1$	$0.88$	$0.82$	${0.8}_{-0.1}^{+0.1}$	${0.0174}_{-0.0009}^{+0.001}$	${1.6}_{-0.2}^{+0.2}$	${3.7}_{-0.7}^{+0.9}$	${12}_{-9}^{+50}$	${0.14}_{-0.01}^{+0.01}$	$-{7.21}_{-0.08}^{+0.07}$	$-3042$
F	$0$	$0.51$	$0.74$	${0.8}_{-0.1}^{+0.1}$	${0.022}_{-0.002}^{+0.002}$	${0.7}_{-0.2}^{+0.2}$	${2.5}_{-0.2}^{+1.0}$	${0.8}_{-0.2}^{+5.0}$	${0.12}_{-0.04}^{+0.04}$	$-{6.96}_{-0.08}^{+0.08}$	$-2946$
G	$-7$	$0.93$	$0.89$	${0.3}_{-0.1}^{+0.3}$	${0.014}_{-0.001}^{+0.001}$	${0.7}_{-0.5}^{+0.7}$	${3.3}_{-0.5}^{+2.0}$	${2.2}_{-0.9}^{+20.0}$	${0.18}_{-0.03}^{+0.02}$	$-{6.62}_{-0.06}^{+0.06}$	$-2832$
H	$-9$	$0.57$	$0.80$	${0.7}_{-0.2}^{+0.2}$	${0.024}_{-0.001}^{+0.001}$	${0.8}_{-0.4}^{+0.3}$	${3.0}_{-0.4}^{+0.9}$	${1.8}_{-0.9}^{+7.0}$	${0.14}_{-0.02}^{+0.02}$	$-{6.77}_{-0.07}^{+0.07}$	$-2893$
I	$-17$	$0.95$	$0.81$	${0.7}_{-0.2}^{+0.3}$	${0.019}_{-0.001}^{+0.001}$	${1.4}_{-1.0}^{+0.7}$	${2.6}_{-0.1}^{+0.4}$	${12}_{-8}^{+50}$	${0.05}_{-0.02}^{+0.02}$	$-{6.44}_{-0.06}^{+0.06}$	$-2770$
J	$-21$	$0.99$	$0.94$	${0.41}_{-0.03}^{+0.04}$	${0.0265}_{-0.0008}^{+0.0007}$	${0.10}_{-0.06}^{+0.1}$	${2.4}_{-0.1}^{+0.2}$	${20}_{-10}^{+60}$	${0.03}_{-0.02}^{+0.01}$	$-{6.75}_{-0.06}^{+0.06}$	$-2882$
K	$-43$	$0.91$	$0.73$	${1.86}_{-0.1}^{+0.08}$	${0.038}_{-0.006}^{+0.005}$	${0.70}_{-0.07}^{+0.09}$	${2.5}_{-0.2}^{+0.2}$	${12}_{-8}^{+40}$	${0.04}_{-0.01}^{+0.02}$	$-{6.64}_{-0.06}^{+0.06}$	$-2831$
L	$-49$	$0.79$	$0.60$	${2.04}_{-0.1}^{+0.09}$	${0.065}_{-0.008}^{+0.01}$	${0.37}_{-0.03}^{+0.04}$	${2.6}_{-0.1}^{+0.2}$	${10}_{-8}^{+50}$	${0.023}_{-0.009}^{+0.01}$	$-{6.47}_{-0.06}^{+0.06}$	$-2741$
M	$-57$	$0.87$	$0.85$	${2.1}_{-0.2}^{+0.4}$	${0.08}_{-0.02}^{+0.1}$	${0.28}_{-0.05}^{+0.06}$	${2.7}_{-0.2}^{+0.7}$	$2_{-2}^{+20}$	${0.06}_{-0.03}^{+0.06}$	$-{6.15}_{-0.06}^{+0.05}$	$-2605$
N	$-64$	$0.53$	$0.34$	${1.78}_{-0.06}^{+0.08}$	${0.024}_{-0.003}^{+0.005}$	${0.7}_{-0.1}^{+0.1}$	${2.3}_{-0.1}^{+0.2}$	${11}_{-7}^{+40}$	${0.04}_{-0.02}^{+0.01}$	$-{6.63}_{-0.07}^{+0.06}$	$-2840$
O	$-75$	$0.68$	$0.67$	${2.13}_{-0.09}^{+0.3}$	${0.10}_{-0.02}^{+0.1}$	${0.24}_{-0.03}^{+0.04}$	${2.6}_{-0.2}^{+3.0}$	$2_{-2}^{+20}$	${0.03}_{-0.02}^{+0.1}$	$-{6.37}_{-0.05}^{+0.05}$	$-2706$
P	$-76$	$0.36$	$0.30$	${1.79}_{-0.07}^{+0.1}$	${0.025}_{-0.003}^{+0.004}$	${0.69}_{-0.09}^{+0.1}$	${2.5}_{-0.2}^{+0.2}$	$3_{-2}^{+10}$	${0.07}_{-0.01}^{+0.02}$	$-{6.84}_{-0.07}^{+0.08}$	$-2921$
Q	$-83$	$0.52$	$0.56$	${1.88}_{-0.08}^{+0.06}$	${0.039}_{-0.007}^{+0.009}$	${0.40}_{-0.06}^{+0.09}$	${2.3}_{-0.1}^{+0.2}$	$5_{-4}^{+30}$	${0.03}_{-0.01}^{+0.02}$	$-{6.42}_{-0.05}^{+0.05}$	$-2740$
R	$-87$	$0.48$	$0.47$	${1.91}_{-0.07}^{+0.05}$	${0.036}_{-0.004}^{+0.006}$	${0.50}_{-0.06}^{+0.08}$	${2.5}_{-0.2}^{+0.2}$	$4_{-2}^{+8}$	${0.05}_{-0.01}^{+0.02}$	$-{6.57}_{-0.05}^{+0.06}$	$-2812$

Table 6. Table 6: Results for variable-methane fits to 18 cloud-free locations. These retrievals include limited free parameters for aerosols and two free parameters describing the methane profile.

location	lat (∘)	${\mu }_H$	${\mu }_K$	layer $\alpha$ (top)		$R{H}_{ac}$	$P_{dep}$ (bar)	$\log {\sigma }_t$	$DIC$	$\Delta DIC$kkkChange in the $DIC$ relative to the initial aerosol-only fit for that location
				$P_{max,\alpha }$ (bar)	${\tau }_{\alpha }$
A	$20$	$0.66$	$0.54$	${0.56}_{-0.07}^{+0.07}$	${0.0208}_{-0.0009}^{+0.001}$	${0.01}_{-0.01}^{+0.02}$	${1.76}_{-0.04}^{+0.01}$	$-{6.71}_{-0.07}^{+0.06}$	$-2857$	$21$
B	$16$	$0.46$	$0.63$	${1.2}_{-0.1}^{+0.1}$	${0.0219}_{-0.0009}^{+0.001}$	${0.05}_{-0.02}^{+0.03}$	${1.98}_{-0.04}^{+0.04}$	$-{6.95}_{-0.08}^{+0.07}$	$-2923$	$-19$
C	$10$	$0.26$	$0.52$	${1.36}_{-0.09}^{+0.07}$	${0.0190}_{-0.0009}^{+0.001}$	${0.04}_{-0.01}^{+0.02}$	${2.02}_{-0.04}^{+0.03}$	$-{7.17}_{-0.08}^{+0.08}$	$-3029$	$-35$
D	$5$	$0.80$	$0.84$	${0.9}_{-0.1}^{+0.1}$	${0.0181}_{-0.0008}^{+0.0009}$	${0.02}_{-0.01}^{+0.02}$	${2.25}_{-0.02}^{+0.03}$	$-{7.06}_{-0.08}^{+0.08}$	$-2997$	$-57$
E	$1$	$0.88$	$0.82$	${0.73}_{-0.08}^{+0.09}$	${0.0165}_{-0.0007}^{+0.0008}$	${0.02}_{-0.01}^{+0.02}$	${2.15}_{-0.03}^{+0.03}$	$-{7.32}_{-0.08}^{+0.08}$	$-3084$	$-42$
F	$0$	$0.51$	$0.74$	${1.27}_{-0.1}^{+0.09}$	${0.022}_{-0.001}^{+0.001}$	${0.05}_{-0.02}^{+0.04}$	${2.11}_{-0.03}^{+0.04}$	$-{7.05}_{-0.08}^{+0.07}$	$-2979$	$-33$
G	$-7$	$0.93$	$0.89$	${0.4}_{-0.1}^{+0.2}$	${0.015}_{-0.001}^{+0.002}$	${0.05}_{-0.03}^{+0.09}$	${2.01}_{-0.05}^{+0.05}$	$-{6.58}_{-0.06}^{+0.05}$	$-2817$	$14$
H	$-9$	$0.57$	$0.80$	${1.3}_{-0.1}^{+0.1}$	${0.0222}_{-0.0007}^{+0.0007}$	${0.012}_{-0.004}^{+0.007}$	${2.02}_{-0.04}^{+0.03}$	$-{7.04}_{-0.09}^{+0.08}$	$-2982$	$-89$
I	$-17$	$0.95$	$0.81$	${1.62}_{-0.04}^{+0.04}$	${0.0118}_{-0.0007}^{+0.0006}$	${0.00004}_{-3e-05}^{+6e-05}$	${2.16}_{-0.03}^{+0.04}$	$-{6.68}_{-0.06}^{+0.06}$	$-2863$	$-92$
J	$-21$	$0.99$	$0.94$	${0.68}_{-0.06}^{+0.1}$	${0.0237}_{-0.002}^{+0.0009}$	${0.17}_{-0.1}^{+0.09}$	${2.7}_{-0.3}^{+0.3}$	$-{6.72}_{-0.06}^{+0.06}$	$-2870$	$12$
K	$-43$	$0.91$	$0.73$	${1.77}_{-0.06}^{+0.05}$	${0.029}_{-0.002}^{+0.002}$	${0.41}_{-0.08}^{+0.1}$	${3.9}_{-0.8}^{+0.9}$	$-{6.68}_{-0.06}^{+0.06}$	$-2846$	$-16$
L	$-49$	$0.79$	$0.60$	${2.05}_{-0.04}^{+0.04}$	${0.006}_{-0.002}^{+0.003}$	${0.016}_{-0.006}^{+0.009}$	${2.18}_{-0.06}^{+0.08}$	$-{6.48}_{-0.06}^{+0.06}$	$-2747$	$-6$
M	$-57$	$0.87$	$0.85$	${2.08}_{-0.06}^{+0.04}$	${0.006}_{-0.003}^{+0.004}$	${0.03}_{-0.02}^{+0.03}$	${2.38}_{-0.08}^{+0.08}$	$-{6.16}_{-0.05}^{+0.05}$	$-2605$	$0$
N	$-64$	$0.53$	$0.34$	${1.88}_{-0.04}^{+0.05}$	${0.021}_{-0.002}^{+0.002}$	${0.39}_{-0.07}^{+0.09}$	${4.0}_{-0.6}^{+0.9}$	$-{6.69}_{-0.06}^{+0.06}$	$-2862$	$-22$
O	$-75$	$0.68$	$0.67$	${2.30}_{-0.07}^{+0.05}$	${0.0012}_{-0.0003}^{+0.0009}$	${0.015}_{-0.004}^{+0.005}$	${2.23}_{-0.05}^{+0.04}$	$-{6.42}_{-0.05}^{+0.05}$	$-2735$	$-28$
P	$-76$	$0.36$	$0.30$	${1.94}_{-0.07}^{+0.1}$	${0.023}_{-0.003}^{+0.003}$	${0.3}_{-0.1}^{+0.1}$	$3_{-1}^{+1}$	$-{6.94}_{-0.07}^{+0.07}$	$-2934$	$-13$
Q	$-83$	$0.52$	$0.56$	${2.21}_{-0.05}^{+0.05}$	${0.0013}_{-0.0003}^{+0.0007}$	${0.006}_{-0.002}^{+0.003}$	${2.19}_{-0.05}^{+0.04}$	$-{6.54}_{-0.06}^{+0.06}$	$-2794$	$-54$
R	$-87$	$0.48$	$0.47$	${2.18}_{-0.03}^{+0.05}$	${0.006}_{-0.002}^{+0.002}$	${0.013}_{-0.006}^{+0.01}$	${2.05}_{-0.04}^{+0.04}$	$-{6.68}_{-0.06}^{+0.06}$	$-2850$	$-38$

Table 7. Table 7: Results for initial retrievals to the haze parameters (CH 4 fixed), including variations in the tropopause temperature as observed by Orton et al., ( 2007 ) .

location	lat	$\Delta T_{trop}\left(K\right)$	layer $\alpha$ (top)			layer $\beta$ (bottom)			$\log {\sigma }_t$	$DIC$	$\Delta DIC$lllrelative to initial model fit
			$P_{max,\alpha }$ (bar)	${\tau }_{\alpha }$	$h_{frac,\alpha }$	$P_{max,\beta }$ (bar)	${\tau }_{\beta }$	$h_{frac,\beta }$
M	$-57$	$-3.8$	${2.1}_{-0.2}^{+0.3}$	${0.09}_{-0.03}^{+0.08}$	${0.26}_{-0.03}^{+0.03}$	${2.8}_{-0.3}^{+0.7}$	$1_{-1}^{+8}$	${0.06}_{-0.04}^{+0.09}$	$-{6.15}_{-0.06}^{+0.06}$	$-2606$	$0$
P	$-76$	$+4.3$	${1.78}_{-0.06}^{+0.07}$	${0.025}_{-0.002}^{+0.003}$	${0.9}_{-0.1}^{+0.1}$	${2.4}_{-0.1}^{+0.2}$	$2_{-1}^{+4}$	${0.06}_{-0.01}^{+0.02}$	$-{6.80}_{-0.08}^{+0.07}$	$-2894$	$28$
R	$-87$	$+7.5$	${1.90}_{-0.07}^{+0.06}$	${0.034}_{-0.004}^{+0.006}$	${0.62}_{-0.09}^{+0.1}$	${2.5}_{-0.1}^{+0.2}$	$4_{-3}^{+9}$	${0.05}_{-0.02}^{+0.02}$	$-{6.55}_{-0.07}^{+0.07}$	$-2796$	$16$

Equations

$$\frac{I}{F}=\frac{r^2}{\Omega }\frac{F_N}{F_{\odot }}$$

$${\sigma }^2={\sigma }_n^2+{\sigma }_p^2+{\sigma }_t^2$$

$$N=N_0\ast \exp \left(\frac{z}{H\ast h_{frac}}\right)$$

$$n(r)\propto r^{6}exp\left(-6\ast \frac{r}{r_p}\right)$$

$$p(\mathbf{x}|\mathbf{d})\propto p(\mathbf{x})p(\mathbf{d}|\mathbf{x})$$

$$p(\mathbf{d}|\mathbf{x})=\prod _i{\left(\frac{1}{2\pi {\sigma }_i^2}\right)}^{1/2}\exp \left[-\frac{{\left(d_i-d_{\mathbf{x},i}\right)}^2}{2{\sigma }_i^2}\right]$$

$$\ln p(\mathbf{d}|\mathbf{x})=-\frac{1}{2}\sum _i\left[\frac{{\left(d_i-d_{\mathbf{x},i}\right)}^2}{{\sigma }_i^2}+\ln \left(2\pi {\sigma }_i^2\right)\right]$$

$${\sigma }^2={\sigma }_n^2+{\sigma }_p^2+{\sigma }_t^2$$

$$D(\mathbf{x})=-2\ln \left\{p(\mathbf{d}|\mathbf{x})\right\}+2\ln \left\{f(\mathbf{d})\right\}$$

$$DIC=\overline{D(\mathbf{x}})+p_D$$

$$p_D=\overline{D(\mathbf{x}})-D(\overline{\mathbf{x}})$$

$$p_D=\frac{1}{2}\mathrm{Var}(D(\mathbf{x}))$$

Full text

We construct a set of high signal-to-noise spectra spanning a range of viewing geometries to constrain the vertical structure of Neptune’s aerosols in a cloud-free latitude band from 2–12∘N. We find that Neptune’s cloud opacity at these wavelengths is dominated by a compact, optically thick cloud layer with a base near 3 bar. Using the pyDISORT algorithm for the radiative transfer and assuming a Henyey-Greenstein phase function, we observe this cloud to be composed of low albedo (single scattering albedo $=0.45^{+0.01}_{-0.01}$), forward scattering (asymmetry parameter $g=0.50^{+0.02}_{-0.02}$) particles, with an assumed characteristic size of $\sim 1\mu$m. Above this cloud, we require an aerosol layer of smaller ($\sim 0.1\mu$m) particles forming a vertically extended haze, which reaches from the upper troposphere ($0.59^{+0.04}_{-0.03}$ bar) into the stratosphere. The particles in this haze are brighter (single scattering albedo $=0.91^{+0.06}_{-0.05}$) and more isotropically scattering (asymmetry parameter $g=0.24^{+0.02}_{-0.03}$) than those in the deep cloud. When we extend our analysis to 18 cloud-free locations from 20∘N to 87∘S, we observe that the optical depth in aerosols above 0.5 bar decreases by a factor of 2–3 or more at mid- and high-southern latitudes relative to low latitudes.

We also consider Neptune’s methane (CH4) profile, and find that our retrievals indicate a strong preference for a low methane relative humidity at pressures where methane is expected to condense. When we include in our fits a parameter for methane depletion below the CH4 condensation pressure, our preferred solution at most locations is for a methane relative humidity below 10% near the tropopause in addition to methane depletion down to 2.0–2.5 bar. We tentatively identify a trend of lower CH4 columns above 2.5 bar at mid- and high-southern latitudes over low latitudes, qualitatively consistent with what is found by Karkoschka and Tomasko, (2011), and similar to, but weaker than, the trend observed for Uranus.

Keywords: Neptune, atmosphere; Radiative transfer; Spectroscopy

NOTICE: this is the authors’ version of a work that was accepted for publication in Icarus. Changes resulting from the publishing process may not be reflected in this document. Figures have been compressed; for higher resolution figures, please contact the corresponding author at shcook@amnh.org.

††slugcomment: Published in Icarus: 15 September 2016.

https://doi.org/10.1016/j.icarus.2016.04.032

1 Introduction

Nearly three decades after the flyby of Voyager 2, Neptune’s upper troposphere and stratosphere remain enigmatic in terms of composition, circulation, and spatial and temporal variability. Aerosols represent a key probe of the atmospheric physics at these altitudes: aerosols trace dynamical motions associated with large-scale circulation patterns and localized storms, and influence the structure of the atmosphere by contributing opacity and by modifying the composition and thermal profile.

Numerous photometric and spectroscopic studies have targeted the distribution and properties of Neptune’s aerosols; the majority of these efforts have focused on the discrete, bright clouds, which dominate Neptune’s near-infrared (NIR) appearance (e.g., Roe et al., 2001b, ; Sromovsky et al., 2001b, ; Gibbard et al., , 2003; Max et al., , 2003; Luszcz-Cook et al., 2010b, ; Irwin et al., , 2011, 2014; de Pater et al., , 2014). These clouds tend to be concentrated in latitude bands (Smith et al., , 1989, Fig. 1), and are typically found to be at altitudes near or below the tropopause, which may suggest a composition of methane ice (Roe et al., 2001b, ; Sromovsky et al., 2001b, ; Max et al., , 2003; Gibbard et al., , 2003; Luszcz-Cook et al., 2010b, ; Karkoschka and Tomasko, , 2011). Methane has an equilibrium condensation temperature of $\sim$ 80K, which corresponds roughly to the 1.5-bar level in Neptune’s atmosphere. Methane ice clouds should therefore form at or above this level. In some cases, a fraction of the observed clouds appear to be well (more than a scale height) above the tropopause (Gibbard et al., , 2003; Irwin et al., , 2011; de Pater et al., , 2014). It is possible that these clouds originate from the settling of hydrocarbons produced by photochemistry (Gibbard et al., , 2003); from localized upwelling of methane within an anticyclone (de Pater et al., , 2014); or from a mix of methane ice and settling hydrocarbon haze particles (Irwin et al., , 2011). Latitudinal trends in the vertical location of clouds have been observed: Gibbard et al., (2003) note an increase in the altitudes of clouds from the south pole to northern midlatitudes. Higher resolution studies indicate more complex patterns in cloud altitudes (Irwin et al., , 2011; Fitzpatrick et al., , 2014; de Pater et al., , 2014).

In this paper we focus on Neptune’s three-dimensional aerosol structure within NIR “dark regions” – regions free of discrete bright clouds, rather than on the bright clouds themselves. In particular, we seek to characterize the other layers of aerosols that are both expected and observed in Neptune’s atmosphere, including global cloud decks, tropospheric haze, and stratospheric haze. The properties and distribution of aerosols in dark regions have implications for Neptune’s global circulation pattern and the mechanisms of cloud and haze production, and serve as a baseline atmospheric structure that can be used when interpreting the spectra of discrete bright clouds.

Thermochemical equilibrium models indicate that Neptune’s condensible gases should form a series of global cloud layers in the troposphere (de Pater et al., , 1991). Based on these models and radio observations, we expect the presence of an H2S ice cloud with a base pressure of 6–7 bar (equilibrium condensation temperature of $\sim 140$K). Initially, photometric and spectroscopic data suggested the presence of an optically thick cloud with a top around 3–5 bar (Hammel et al., , 1989; Baines and Smith, , 1990; Baines et al., 1995a, ; Baines et al., 1995b, ), which could be the putative H2S ice cloud. Since then, several analyses have presupposed a cloud with these properties (e.g. Sromovsky et al., 2001b, ; Gibbard et al., , 2002; Luszcz-Cook et al., 2010b, ). However, other results indicate this cloud may be optically thin, or may not exist at all, with the tropospheric opacity instead arising from vertically extended haze (Sromovsky et al., 2001b, ; Karkoschka and Tomasko, , 2011; Karkoschka, 2011a, ). If an optically thick cloud deck is indeed present at these depths, it must have a very low single scattering albedo in the NIR (less than 0.1–0.2 at 1.6 microns – Sromovsky et al., 2001b ; Roe et al., 2001a ).

Above the H2S cloud deck, thermochemical models predict CH4 condensation (de Pater et al., , 1991). Some atmosphere models posit a transparent, horizontally homogenous layer of condensed methane in addition to any discrete, NIR-bright methane ice clouds (Baines and Smith, , 1990; Baines and Hammel, , 1994; Baines et al., 1995a, ; Sromovsky et al., 2001a, ). The base of this methane haze is presumed to be near the methane condensation level. Sromovsky et al., 2001a determined that the global average optical depth of this layer is 0.11 at 0.75 microns, higher than the values of 0.03 and 0.05 inferred by Baines and Smith, (1990) and Baines and Hammel, (1994), respectively for the global average, but less than the values found for latitudes of 22–30∘S by Pryor et al., (1992); Baines and Hammel, (1994), assuming isotropic scattering. Assuming the methane haze is composed of $\gtrapprox 1\mu$m Mie scatterers (Conrath et al., , 1991; Baines and Hammel, , 1994; Burgdorf et al., , 2003), the NIR optical depth of this haze is predicted to be comparable to that at 0.75 $\mu$m. Other atmospheric models have a very different aerosol structure in the upper troposphere: Gibbard et al., (2002) found that a haze at an altitude of 0.3 bar provides a good fit to their NIR imaging data of cloud-free regions, although more complex, multi-layer aerosol distributions are also permitted. Irwin et al., (2011) treat NIR-dark regions in the same ways as bright regions– with a two-layer aerosol structure consisting of a cloud deck at 2 bar and a horizontally varying upper cloud deck at 0.02–0.2 bar. In their 2014 reanalysis of these data, they model a dark belt near 10∘S and placed the upper cloud at 0.3 bar. The extended haze model of Karkoschka and Tomasko, (2011) is free of aerosols above the 1.4 bar pressure level, aside from discrete clouds, although Karkoschka, 2011a consider the possibility of an extended upper tropospheric haze between the tropopause and 1.4 bar, which is 10 times thinner than their main haze layer and is latitudinally variable.

In the stratosphere, methane is photodissociated by solar ultraviolet photons, which results in the production of a number of hydrocarbons– primarily ethane and acetylene (Romani et al., , 1993). As these hydrocarbons settle, they condense into optically thin haze layers. The observations of Pryor et al., (1992) and Baines and Hammel, (1994) indicate that the global average optical depth for the stratospheric aerosols is of order 0.1 at 0.75 $\mu$m, with a preferred mean particle radius of 0.2 $\mu$m. Using Mie theory to extrapolate to NIR wavelengths, we find the predicted opacity of such particles is $<0.02$ at 1.6 $\mu$m. Karkoschka and Tomasko, (2011) find the stratospheric haze optical depth to be an order of magnitude lower than Pryor et al., (1992); they note that the discrepancy is primarily due to different assumptions for the single scattering albedo by the two authors. Meridional variability – specifically a deficit in stratospheric haze near the equator – was observed by Baines and Hammel, (1994).

With the objective of improving the characterization of Neptune’s aerosols, we have obtained H- and K- broadband (1.47 – 2.38 $\mu$m) spectroscopic observations of Neptune with the OSIRIS AO-assisted integral field spectrograph on the Keck II telescope. The NIR is well suited to studying the distribution of aerosols in the upper troposphere and lower stratosphere, since wavelength variations in methane opacity result in sensitivity to a broad range of atmospheric depths (Fig. 2). Using these data, we retrieve Neptune’s aerosol properties in dark regions at latitudes from 20∘N to 87∘S, considering the effects of signal to noise (S/N) and viewing geometry changes (center-to-limb variations) on our results. Our observations offer several advantages over previous studies for characterizing Neptune’s dark regions: whereas traditional slit spectroscopy does not facilitate the clean separation of the quiescent background from the bright features that typically dominate the observed NIR intensity, OSIRIS produces three-dimensional data cubes ($x$, $y$, and wavelength); these data offer moderate (R$\sim 3800$) spectral resolution at very good spatial resolution of 0.06–0.08” across more than 90% of Neptune’s visible hemisphere, and are therefore well-suited for extracting spectra of dark regions that are relatively free of contamination from nearby bright features. Furthermore, in contrast to slit spectroscopy, the OSIRIS cubes are ideal for selecting spectra at desired latitudes and viewing geometries. These data represent more than a factor of two increase in spatial resolution over the NIFS integral field spectrograph observations previously reported by Irwin et al., (2011, 2014) and the first published integral field spectrograph observations of Neptune in K band. As demonstrated in Fig. 2, the addition of the K-band data increases our sensitivity to higher altitudes (see also Supplementary Fig. 29).

The structure of this paper is as follows: in Section 2 we describe the observations and data reduction. Section 3 presents an overview of the radiative transfer code and atmospheric retrieval procedure. In Section 4, the latitude band from 2–12∘N is analyzed using binned spectra spanning a range of viewing geometries. In Section 5, we evaluate the effects of the vertical temperature and CH4 profiles on our derived aerosol structure for a high signal-to-noise (S/N), binned spectrum. Section 6 broadens our investigation to additional locations on Neptune. A summary and discussion are presented in Sections 7 and 8.

2 Observations and data reduction

2.1 Observations

We observed Neptune in the H (1.47 – 1.80 $\mu$m) and K (1.97–2.38 $\mu$m) broadband filters (referred to as Hbb and Kbb, respectively) on 25 and 26 July 2009 with the near-infrared imaging spectrometer OSIRIS on the 10-meter W. M. Keck II telescope. The intention was to obtain data from both hemispheres of Neptune, which are observable on consecutive nights due to Neptune’s rotation period of $\sim 16$ hours (Karkoschka, 2011b, ). However, the seeing conditions on 25 July 2009 were unstable, resulting in poor quality data on that night, particularly in H band. Therefore we present only the data from 26 July 2009.

In broad band mode, $\sim 16\times 64$ lenslets in the OSIRIS integral field spectrograph image sampler are illuminated. We selected an instrumental plate scale of 0.035”; at this resolution, a 0.56”$\times$2.24” patch of the sky is sampled with each exposure. The light passing through each individual lenslet is diffracted to produce a moderate-resolution (R$\sim 3800$) spectrum at each of 1019 locations.

Neptune’s angular diameter on 26 July 2009 was 2.35” and its distance from Earth was 29.1 AU. At this distance, each 0.035” pixel corresponds to 740 km at disk center. We oriented the exposures so that the long direction of each frame pointed east-west (parallel to the equator) along the planet. Since Neptune’s diameter was slightly larger than the length of the field of view, the western limb of the planet was not observed. A set of eight 300-second exposures were taken in each of the Hbb and Kbb filters. The first and last frames of each set were ‘sky’ frames, for subtracting the infrared background. The remaining exposures were positioned to obtain a series of contiguous horizontal slices starting at Neptune’s south pole and stepping across the disk until the northern limb was reached. During the night, several A0-type stars with known J, H, and K magnitudes were also observed, to serve as telluric and photometric standards. Two exposures were taken of each star, dithering between the two halves of the field of view so that the exposures could serve as sky frames for one another. The airmass of the Hbb observations of Neptune was 1.3–1.4; the airmass of the Kbb observations was 1.2.

2.2 Initial processing

The initial processing of the scientific and calibration data is accomplished using version 2.3 of the OSIRIS data reduction pipeline111http://irlab.astro.ucla.edu/osiris/pipeline.html, which performs sky subtraction, cleans the data of cosmic rays and common instrumental issues, extracts spectra from the raw frames, and combines the individual science exposures into mosaicked Hbb and Kbb spectral cubes that cover more than 90% of Neptune’s visible hemisphere. This last step is performed using the telescope right ascension and declination, which has an estimated astrometric accuracy of 40 mas222https://www2.keck.hawaii.edu/inst/osiris/OSIRIS_Manual_v2.3.pdf. The pixel values of locations covered by multiple pointings are computed by averaging the pixel values of the individual exposures. The data reduction pipeline also produces mask cubes that flag bad pixels identified by the pipeline. We flag 30 additional pixels near the edge of each of the unmosaicked frames, which appear anomalously bright due to lenslet masking errors; and 5 edge columns of the mosaicked Kbb cube, which are also anomalously bright. The initial Kbb mosaic exhibited brightness fluctuations between the individual exposures used to create the mosaic, indicative of an imperfect sky subtraction due to short-timescale fluctuations in the sky brightness. To improve the sky subtraction, at each wavelength we take the median of all unflagged background (off-disk) pixels in each individual exposure, to create a median background spectrum for that exposure. These background spectra are subtracted from the appropriate pixels in the mosaicked Kbb cube. A similar correction in Hbb is not possible, because Neptune itself is much brighter at these wavelengths, so that pixels near to the disk are dominated by flux from the nearby planet rather than by the sky.

2.3 Flux calibration

In order to calibrate our target data we use the calibration star observed under the most similar atmospheric conditions (in time and airmass) to those of the Neptune observations (SAO164840, airmass=1.2). After initial processing of the pair of star exposures in each filter, stellar spectra are extracted using aperture photometry, then divided by a model of an A0 star spectrum that is calibrated to the 2MASS J, H and K magnitudes (Castelli and Kurucz, , 2004) and adjusted for reddening (Fitzpatrick and Massa, , 1999). The two calibration spectra are averaged and corrected for the difference in airmass between the science target and calibrator observations. Then, each pixel in the science data cubes is divided by this stellar spectrum to correct for telluric lines and convert from units of counts into units of flux density. Finally, we convert the data from observed flux density into units of $I/F$, defined as (Hammel et al., , 1989):

[TABLE]

where $r$ is Neptune’s heliocentric distance in AU, $\pi F_{\odot}$ is the sun’s flux density at Earth’s orbit (Colina et al., , 1996), $F_{N}$ is Neptune’s observed flux density, and $\Omega$ is the solid angle subtended by a pixel on the detector. By this definition, $I/F=1$ for uniformly diffuse scattering from a Lambert surface when viewed at normal incidence. The reduced and calibrated data, averaged over wavelength, are shown in Fig. 1.

2.4 Spatial and spectral resolution

The spatial resolution of the data is estimated from the stellar spectra. We find that the full width at half maximum (FWHM) of the point spread function (PSF) core is 2.0–2.2 pixels in Hbb, and 1.7–1.8 pixels in Kbb. We note that the PSF is non-Gaussian, having in addition to the narrow core a broad component (halo, FWHM $\sim$6–10 pixels) which complicates the interpretation of discrete features and the nearby background- this will be discussed further in Paper II. In this work, we avoid whenever possible locations that are within 10 pixels of identified bright features, to minimize their influence on our dark-region spectra.

Prior to analysis, the data are smoothed to a spectral resolution of R$\sim$600 to increase the signal to noise (S/N) and match the resolution of our methane absorption coefficients (Appendix A.1.3). Additional binning of the Kbb data, to increase S/N, is performed for the analysis in Section 6.

2.5 Data uncertainties

To measure the wavelength-dependent noise in the data ($\sigma_{n}$), we use the unflagged pixels in each cube (Hbb and Kbb) that lie outside of Neptune’s disk as determined by image navigation (Section 2.6). An uncertainty spectrum is calculated as the standard deviation of all background pixels at each wavelength, after rejecting outlier pixels (pixels with a value more than 30 and 15 standard deviations away from the mean background in the Hbb and Kbb filters, respectively). This rejection typically discards 0–3 pixels per wavelength, or $\leq 0.03$% of the background pixels at a given wavelength. The mean uncertainty in I/F in a single pixel due to random noise is 0.001. The uncertainties are highest at wavelengths where the atmospheric transmission is low: near the edges of the bands and in the CO2 features at 2.0 and 2.06 microns (Fig. 3).

A second contribution to the uncertainty is systematic, introduced by the photometric calibration and telluric correction. We consider two components to the uncertainty in the calibrator spectrum, which we call “overall photometry” and “relative photometry” uncertainties. The overall photometry uncertainty refers to gross errors in the photometry due to, e.g., uncertainties in the H and K magnitudes of the calibrator star or loss of stellar flux off the edges of the detector. Based on comparisons with other calibrator stars observed on the same night and previous experience our conservative estimate for these uncertainties is 15–20%. These errors would systematically offset all wavelengths in a given band, and we choose not to incorporate this term into the data uncertainty estimate, but do test the effects of such errors on our model fits. This is done by selecting a spectrum from the data and generating a grid of test spectra with offsets of $\pm 20\%$ in the Hbb and/or Kbb photometry from the nominal calibration. We then perform a sample analysis, akin to those discussed in Section 6, for each of the test spectra and compare the results. We find that variations in the model results between test spectra are generally less than the parameter uncertainties from the model retrievals (Section 3.2). The one exception is that the single scattering albedos of the aerosols increase/decrease by 15–20% when the Hbb photometry is adjusted by $\pm 20\%$, respectively.

The relative photometry error ($\sigma_{p}$) refers to the random noise and small slopes (in wavelength) introduced by the calibrator spectrum, which we infer and characterize based on observed differences between the pair of spectra extracted for the calibrator star. We estimate $\sigma_{p}$ by determining the difference between our two stellar calibration spectra at each wavelength. We find that the differences between the two calibration spectra are correlated in wavelength – that is, data from neighboring wavelengths are more likely to show similar relative offsets than data from widely separated wavelengths. However, there is not enough information in our calibration spectra to remove these systematic effects. Fortunately, differences between our calibration spectra are small, averaging 1.7% of the flux in Hbb and 1.8% in Kbb. We therefore include in our uncertainty estimate the term $\sigma_{p,\lambda}=0.017\times(I/F)_{\lambda}$ in H band and $\sigma_{p,\lambda}=0.018\times(I/F)_{\lambda}$ in K band. In the spectrum of a single pixel, the relative photometry error is small compared to the random error; however, when spectra are averaged, the relative photometry errors can dominate, especially at the brightest wavelengths (Fig. 3). The calculation of the overall uncertainty term from these contributions is discussed in Section 3.2.

2.6 Image navigation

The brightness and variability of Neptune’s clouds present a challenge for automated methods of image navigation. Therefore, to determine planetocentric latitude, longitude and emission angle ($\theta$, defined as the angle between the line of sight and the local normal to the planetary surface, specified throughout by $\mu\equiv\cos\theta$ ) for each pixel in our data, we fit a circle to the limb of planet by eye using the wavelength-collapsed Hbb and Kbb cubes. By testing small perturbations to the disk center positions, we estimate the accuracy of our centering to be $\lesssim 1$ pixel – roughly half the FWHM of the PSF core. The JPL Horizons333http://ssd.jpl.nasa.gov/?horizons ephemerides provide the sub-observer latitude (-28.7∘) and longitude (152∘ and 126∘ at the midpoint of the Hbb and Kbb observations, respectively). The phase angle of the observations was 0.7∘.

3 Modeling

To analyze our OSIRIS spectra, we use a Python-based atmospheric retrieval code which has two main components: a forward model, which performs the radiative transfer, and a retrieval algorithm to determine the posterior probability distribution of the model parameters. In this section we present a brief overview of our retrieval code and a summary of the parameters most relevant to the analysis described in Sections 4 – 6. A more detailed description of the model may be found in Appendix A.

3.1 Forward model

We construct a $\sim$100-layer model atmosphere that extends from pressures of 20 bar up to $10^{-5}$ bar from an input thermal structure (temperature-pressure profile); atmospheric composition as a function of depth; gas opacities as a function of temperature and pressure; and a description of the aerosols. The nominal thermal and composition profiles for this study (Appendix A) are illustrated in Fig. 4. Our model accepts any number of aerosol layers, which can be placed at any depth within the model atmosphere and can overlap one another – a detailed description of the aerosol layer parameterization is also available in Appendix A. The key model parameters for this study, and their default values, are summarized in Table LABEL:table:model; these parameters involve perturbations to the nominal thermal and methane profiles and variations on the aerosol properties and distribution.

Given a model atmosphere, we solve the radiative transfer equation using either a two-stream approximation or a Python implementation of the discrete ordinate method for radiative transfer (pyDISORT). The two stream radiative transfer code is a two-point quadrature method, following Toon et al., (1989); this code has been previously used in NIR studies of Jupiter, Uranus, and Neptune (de Pater et al., 2010b, ; de Pater et al., 2010a, ; de Kleer et al., , 2015; Luszcz-Cook et al., 2010b, ). The pyDISORT software was introduced by Ádámkovics et al., (2015) for use with Titan, and is an implementation of CDISORT444www.libradtran.org. The NIR scattering phase functions of Neptune’s aerosols are poorly known; for simplicity, we assume single Henyey-Greenstein phase functions in our pyDISORT calculations, which are expanded into infinite series of Legendre polynomials. The radiative transfer is performed using four moments (and streams), which is the number we find adequately represents this simple phase function.

The two-stream approximation sacrifices computational accuracy for the sake of speed, an essential compromise for this analysis: using two stream, a typical retrieval for a single model and single spectrum takes of order one day to complete on one of our 12-core machines. An equivalent pyDISORT run with four streams is a factor of $\sim 3$ slower, with computation time increasing roughly as the third power of the number of streams. In the interest of time, our initial retrievals (Sections 4.1 – 4.3) utilize the two stream algorithm. The effects and limitations of using this approximation are examined in Appendix B and briefly reviewed in Section 4.4; in each subsequent retrieval, the adopted radiative transfer algorithm is clearly indicated.

3.2 Retrieval

As in de Kleer et al., (2015), we pair our forward model with a Python implementation of the Goodman and Weare, (2010) affine-invariant Markov chain Monte Carlo (MCMC) ensemble sampler called emcee (Foreman-Mackey et al., , 2013), as a means of estimating model parameter values and uncertainties. Like any MCMC algorithm, emcee generates a sampling approximation to the posterior probability function by constructing one – or in this case, an ensemble of – chain(s) sampled from the desired probability distribution by random walk. In contrast to the simpler and more common Metropolis-Hastings MCMC sampling, this algorithm is more efficient, insensitive to the aspect ratio of the distribution, requires much less hand-tuning, and is easily parallelized. The details of our retrieval method are described in Appendix A. One modification of the retrievals in de Kleer et al., (2015) is that our uncertainty estimate now includes a free parameter, $\sigma_{t}$, which is a model “tolerance”, representing any unknown uncertainties, particularly deficiencies in the model (for example, in the model opacities or in assumptions for the composition or thermal structure) that prevent the model from matching the data. The full model uncertainty is thus defined as:

[TABLE]

where $\sigma_{n}$ is the random noise component and $\sigma_{p}$ is the relative photometry uncertainty, which is a scale factor multiplied by the model reflectivities (Section 2.5). In practice, we find that including $\sigma_{t}$ does not appreciably change the estimated values of other model parameters, but improves convergence and causes the parameter uncertainties to be more realistic.

Once we have used emcee to approximate the posterior probability distribution for a given model and dataset, we report the 16th, 50th, and 84th percentiles from the marginalized parameter distributions (e.g., Table LABEL:table:1L). In some cases, we also provide plots of the one- and two-dimensional projections of the posterior probability distributions (e.g., Fig. 5). We frequently refer to the 50th percentile of the marginalized distributions as the model “best-fit” for convenience, and refer to the 16th and 84th percentiles of the marginalized distributions as the “1$\sigma$ parameter uncertainties.” We note that these best-fit and uncertainty values quantify the spread of parameter values for a given model, and do not necessarily encompass the range of values a parameter might take for different atmospheric models (different model assumptions and/or free parameters). To compare different models, we calculate the Deviance Information Criterion ($DIC$, Spiegelhalter et al., , 2002). The $DIC$ includes a term for goodness-of-fit with a penalty term for model complexity. The $DIC$ for a single model is not meaningful, and may be a positive or negative number. However, differences in the $DIC$ provide a metric for evaluating the relative success of different models for the same dataset, with a difference in the $DIC$ ($\Delta DIC$) of 10 or more indicating a preference for the model with the lower $DIC$. See Appendix A for more information on the $DIC$. We note that a large value of $\sigma_{t}$ may be another indicator that a model is not a good match to the data, as it implies that the differences between the data and model are not well captured by the known data uncertainties ($\sigma_{n}$ and $\sigma_{p}$).

4 Aerosol properties and structure: 2–12∘N

In this section, we focus our analysis on constraining the properties and structure of Neptune’s aerosols. We take advantage of the spatial resolution of our data to select out only spectra in the latitude range of 2–12∘N. Across this range, discrete clouds are generally not observed n Neptune (and none are observed in our data). By selecting out a relatively thin latitude band, we hope to minimize inherent variations in the properties of the atmosphere (temperature, composition, and aerosol structure). We separate the data into bins at intervals of 0.1 in the cosine of the emission angle, centered at $\mu=$ 0.3, 0.4, 0.5, 0.6, 0.7, and 0.8. The bins contain seven ($\mu=0.3$) to 106 ($\mu=0.8$) individual spectra (Fig. 6). Within each bin, we calculate the mean spectrum (considering only unflagged pixels) and the standard deviation as a function of wavelength, which is divided by the square root of the number of input spectra to estimate the noise in the averaged spectrum. The uncertainties estimated in this way are somewhat larger than ($1.1$–$2.2$ times) those found from the background pixels (Section 2.5); this may be due to real variations in the atmospheric properties within the latitude and $\mu$ ranges included in each bin. We define $\sigma_{n}$ for the six mean spectra using these higher uncertainty estimates.

We then perform a series of retrievals, considering the six spectra simultaneously and assuming they all reflect the same atmosphere in terms of aerosol structure, temperature, and composition. This dataset is referred to as the ‘limb-darkening’ dataset, since it contains information on the aerosol structure that is revealed due to changes in viewing geometry. We use the default thermal structure (Appendix A.1.1); set mCH${}_{4,t}=0.04$, mCH${}_{4,s}=3.5\times 10^{-4}$, and$RH_{ac}=100\%$ (refer to Table LABEL:table:model for definitions); and assume no shallow tropospheric methane depletion (Appendix A.1.2). In Sections 4.1 – 4.3 we use the two-stream algorithm to solve the radiative transfer equations; in Section 4.4 we evaluate the limitations of this approximation and present a retrieval using pyDISORT.

4.1 Single aerosol layer (1L) models

Considering the uncertainty in Neptune’s aerosol structure, our approach is to determine the simplest structure that can reproduce our data. Therefore, we begin by testing models with only a single layer of aerosols, defined by its maximum pressure ($P_{max}$), total optical depth ($\tau$), single scattering albedo ($\omega_{0}$), phase function asymmetry factor ($g$), and fractional scale height ($h_{frac}$). Previous models of Neptune’s cloud-free regions or disk-averaged aerosol structure (e.g. Baines and Hammel, , 1994; Karkoschka and Tomasko, , 2011; Karkoschka, 2011a, ; Irwin et al., , 2011, 2014) find that Neptune’s aerosol opacity is greatest in the troposphere, at pressures $>1$ bar, where particle sizes $\gtrapprox 1\mu$m are indicated (Conrath et al., , 1991); therefore we assume that the extinction cross section is given by a particle distribution with $r_{p}=1\mu$m (Appendix A.1.4).

4.1.1 1L_nom model

We first perform a retrieval in which the aerosol layer top pressure $P_{min}$ is set to 1 mbar. This parameterization permits either a single, compact cloud deck (for $h_{frac}<1$) or a haze that extends from the base pressure to the top of the observable atmosphere (for $h_{frac}\sim 1$). The results from this retrieval are summarized in Table LABEL:table:1L and Fig. 7, and the aerosol distribution corresponding to these results is illustrated in Fig. 8. We retrieve an optically thick, tropospheric aerosol layer with a deep base at 4–30 bar and a moderate scale height of $\sim 0.5$ times the gas scale height. As illustrated by Fig. 5, $P_{max}$ and $\tau$ are highly correlated, with a deeper cloud base pressure corresponding to a higher total optical depth, such that the optical depth per bar at pressures that the observations are sensitive to (Fig. 2) does not vary appreciably between models (Fig. 8). Figure 5 also shows that the posterior probability distributions of $P_{max}$ and $\tau$ are roughly flat across a broad range of values. This makes sense, since adding additional cloud opacity below an already optically thick cloud should not influence the resulting model spectrum.

The retrieved value of the asymmetry factor $g$ is $0.34^{+0.02}_{-0.03}$, in contrast to a 1.6 $\mu$m value of 0.76 predicted by Mie theory for 1-$\mu$m particles (see Fig. 24). The retrieved cloud is also fairly dark, with $\omega_{0}\approx 0.4$ for the best fit. However, as we will show in Section 4.4, $\omega_{0}$ and $g$ are strongly influenced by the use of two stream for the radiative transfer, and the results for these parameters should be interpreted with great caution for all models which rely on the two-stream algorithm.

4.1.2 1L_pmin1.4 and 1L_pminfree models

As shown in Fig. 8, the 1L_nom model is not sensitive to the precise value of the parameter $P_{min}$, since the aerosol opacity drops off below the altitude of the tropopause. However, this model setup does not capture a scenario in which the aerosols are vertically uniform below some cutoff depth in the troposphere, as inferred by Karkoschka and Tomasko, (2011); Karkoschka, 2011a . These authors found that Neptune’s aerosols are uniformly distributed between 1.4 bar and a pressure of at least 10 bar. We now consider whether this type of aerosol structure is compatible with our observations. In the 1L_pmin1.4 retrieval, we perform a run with the same free parameters as in the 1L_nom retrieval, but set $P_{min}$ to the Karkoschka and Tomasko, (2011) haze upper bound of 1.4 bar. As shown in Table LABEL:table:1L, the retrieved aerosol properties for the 1L_pmin1.4 retrieval are consistent with the 1L_nom case; however, the fit quality as measured by the $DIC$ is significantly worse when the aerosols are restricted to $P>1.4$ bar. Finally, we consider the case where $P_{min}$ is a free parameter in the retrieval, to determine whether our data indicate a preference for a different value of $P_{min}$, and whether this parameter influences the probability distributions of other model parameters. This run is referred to as the 1L_pminfree retrieval. We find that, as in the 1L_pmin1.4 case, the retrieved aerosol properties are consistent with the 1L_nom retrieval. The fit quality does not improve when $P_{min}$ is allowed to be a free parameter in the retrieval.

4.2 Two aerosol layer (2L) models

We next consider models containing two layers of aerosols. As in the $1L$ fits, we allow $P_{max}$, $\tau$, $h_{frac}$, $\omega_{0}$, and $g$ to be free parameters with flat (for $h_{frac}$, $\omega_{0}$, $g$) or log-flat (for $P_{max}$, $\tau$) priors. For simplicity, we restrict $\log{P_{max,\alpha}}(bar)=[-6.8,0.9]$ ($P_{max,\alpha}$ between 1 mbar and 2.5 bar) and $\log{P_{max,\beta}}(bar)=[-2.3,2.3]$ ($P_{max,\beta}$ between 0.1 and 10 bar). For the bottom layer (layer $\beta$) we assume the extinction cross section is given by the particle distribution with $r_{p}=1\mu$m. For the top layer (layer $\alpha$), we consider three particle distributions, with characteristic particle sizes of $r_{p}=0.1,0.5,$ and 1.0 $\mu$m (Appendix A.1.4).

This functional form for the aerosol distribution encompasses the parameterizations of several previous efforts: for example, Sromovsky et al., 2001a find that the background aerosol structure consists of a homogeneous reflecting layer at 3.8 bar and a transparent reflecting layer near 1.3 bar. Gibbard et al., (2002) assume an optically thick cloud layer at $P\geq 3.8$ bar and find that the addition of a thin haze layer at 0.3 bar matches their cloud-free region data well, although more complex aerosol structures are also deemed possible. Irwin et al., (2011) find that a dark region near the equator (their Pixel 2, Fig. 7) in their Gemini data from September 2009 are well-matched by a moderate optical thickness cloud at $\sim 2$ bar, and a thinner cloud between 50 and 100 mbar. In their 2014 reanalysis, these authors revise the depth of the upper cloud to near 300 mbar.

The retrieved parameters for our initial 2L fits are summarized in the top of Table LABEL:table:2L: the $r_{p}=0.1,0.5,$ and 1.0$\mu$m cases for the top aerosol layer are labeled 2L_hazeA, B and C, respectively. We find a preference for small particles in the top aerosol layer (2L_hazeA). All three cases result in a similar aerosol structure, consisting of an optically thin ($\tau_{\alpha}\approx 0.1$) upper aerosol layer with a base around 0.5 bar, and a more optically thick ($\tau_{\beta}\sim$ 0.3–0.6) layer near 3 bar. In all cases the bottom layer is found to be compact. The 2L_hazeA models are illustrated relative to the data in Fig. 9, and the posterior probability distributions for this retrieval are shown in Fig. 10. The model fits show a clear improvement over the single aerosol layer (1L) models, which is reflected in the substantial decrease in the $DIC$ ($\Delta DIC<-2000$). The improvements are concentrated at the shortest wavelengths and the $\sim$1.65–1.8 $\mu$m spectral region.

Figure 8 permits a visual comparison of the aerosol distributions for the best one-layer retrieval (1L_nom) and the 2L_hazeA retrieval. We observe that aerosol layer $\beta$ in the 2L_hazeA model is qualitatively similar to the single cloud layer in the 1L_nom retrieval. Fig. 9 shows, in addition to spectra corresponding to the best fit and randomly selected 2L_hazeA models, two models that correspond to only layer $\alpha$ or layer $\beta$ of the best-fit 2L_hazeA model. These plots highlight that layer $\alpha$ is primarily responsible for the improved fit at $\lambda\lessapprox 1.5$ $\mu$m and $\lambda\sim$1.65–1.8 $\mu$m.

4.3 Atmospheres with more than two aerosol layers

Residuals in our two-layer model fits are greater than our estimated data uncertainties; one potential contribution to this discrepancy is that a two-layer aerosol model is insufficient to describe the true distribution of aerosols. We attempted to model the atmosphere with three aerosol layers using a number of prior distributions and initial values for the parameters. We found that the inclusion of a third aerosol layer caused the retrieval to become unstable: parameters had multimodal posterior probability distributions that depended substantially on the parameter initialization. Aerosol layers frequently overlapped substantially and appeared less physically compelling than our two-layer solutions. Therefore, we restrict ourselves to two aerosol layers for the remainder of the analysis.

4.4 Limb-darkening analysis with pyDISORT

In Appendix B, we show that two stream is accurate only when modeling a single spectrum at a single viewing geometry, and when the two stream asymmetry factor $g$ is allowed to be a free parameter in the retrieval. The tests described in Appendix B strongly suggest that while two stream is sufficient for estimating the vertical profile of aerosols, it is not well suited to fully take advantage of the information content of our binned, limb-darkening dataset – that is, to place the best possible constraints on the optical properties of the aerosols from the changes in the spectrum due to variations in viewing geometry. Despite the added computational cost, we require a more accurate radiative transfer solver to achieve this goal.

We therefore perform a new retrieval of the aerosol properties for the limb-darkening dataset using the pyDISORT algorithm with four streams (section 3.1, Appendix A). This run, referred to as the 2L_DISORT model, corresponds exactly to the 2L_hazeA run, except for the use of the more accurate radiative transfer solver. The main purpose of the 2L_DISORT retrieval is to accurately constrain Neptune’s aerosol properties (including $\omega_{0}$ and $g$) from 2–12∘N. Secondary goals are to determine 1) whether the use of the two stream approximation contributes to the residuals in our earlier model fits, and 2) if and how the use of the two stream approximation influences the retrieved values of model parameters.

The results of the 2L_DISORT run are presented in Table LABEL:table:SN and Fig. 11; the models are shown relative to the data in Fig. 12. We find that the more accurate radiative transfer solver produces a better fit to the data than two stream with $\Delta DIC=-332$ for the 2L_DISORT retrieval relative to the 2L_hazeA case. While the parameter probability distributions retrieved by the two algorithms are generally different with statistical significance, we find that the aerosol structure derived by the two algorithms is qualitatively similar, consisting of an optically thin ($\tau_{\alpha}\lessapprox 0.1$) and vertically extended upper layer with a base at 0.5–0.6 bar, and a higher optical depth, more vertically compact bottom aerosol layer near 3 bar. The 2L_DISORT retrieval favors a totally optically thick ($\tau_{\beta}>1$) bottom aerosol layer. As observed in the 1L retrievals, the base of an optically thick cloud is not precisely constrained and is strongly correlated with layer optical depth. The upper aerosol layer is observed to be more vertically extended in the 2L_DISORT solution, with a factor of two larger aerosol scale height than in the best-fit two stream model. A visual comparison of the retrieved aerosol distributions for the two models can be found in Fig. 8.

Consistent with our expectations from Appendix B, we find that the retrieved optical properties of the aerosols are most affected by the choice of radiative transfer solver. In the 2L_DISORT case, a moderately forward scattering phase function ($g_{\beta}=0.50^{+0.02}_{-0.02}$) is preferred for the deeper aerosol layer, as compared to the backscattering phase function suggested by two stream ($g_{\beta}=-0.6^{+0.2}_{-0.2}$). Using pyDISORT, we retrieve a more reflective, but still relatively dark, bottom cloud ($\omega_{0,\beta}=0.45^{+0.01}_{-0.01}$). For the upper aerosol layer, the 2L_DISORT retrieval indicates nearly perfectly reflecting aerosols, with $\omega_{0,\alpha}=0.91^{+0.06}_{-0.05}$ – much brighter than what was found by the 2L_hazeA retrieval.

5 Atmospheric composition and temperature, 2–12∘N

We now investigate the effects of varying our assumed methane and thermal profiles. The full limb-darkening dataset described in the previous section, while ideally suited for placing center-to-limb constraints on the scattering properties of Neptune’s aerosols, is expensive to fit in terms of CPU hours: each retrieval to the limb-darkening dataset involves simultaneously fitting six spectra, and the variation in viewing geometry necessitates the use of the slower pyDISORT radiative transfer algorithm. Therefore, for this next set of retrievals, we use only the $\mu=0.8$ spectrum from the limb-darkening dataset, which has the highest S/N of the six binned spectra from the previous section. We refer to this set of retrievals as the $SN$ retrievals, to reflect that they are performed on the ‘high S/N’ dataset.

5.1 Nominal model, high S/N dataset

Since the $SN$ dataset consists of only a single spectrum at a single viewing geometry, we expect (from Appendix B) that we should be able to use the more computationally-efficient two stream algorithm to perform the radiative transfer, as long as we are interested primarily in the aerosol structure (as opposed to the aerosol scattering properties). Before exploring the role of composition and temperatures in our model fits, we first establish a nominal model to the $SN$ dataset. This SN_nom retrieval has the same model parameters as the 2L_hazeA and 2L_DISORT runs, but is performed using two stream on the single spectrum of the $SN$ dataset only. The results of this retrieval are presented in Table LABEL:table:SN and illustrated in Fig. 13. Encouragingly, we find that, with the exception of the optical properties, the parameter values retrieved using the SN_nom model are in good agreement with those from the more rigorous 2L_DISORT retrieval described in the previous section. In particular, the retrieved base pressures of both aerosol layers agree to within 1$\sigma$ for the two runs. The SN_nom retrieval finds aerosol layer $\alpha$ to be slightly more optically thick ($\tau_{\alpha}=0.050^{+0.01}_{-0.009}$ vs. $\tau_{\alpha}=0.019^{+0.002}_{-0.001}$ for 2L_DISORT) and vertically extended ($h_{frac,\alpha}=1.2^{+0.2}_{-0.2}$ vs. $h_{frac,\alpha}=0.85^{+0.07}_{-0.06}$ for 2L_DISORT). Both models find that aerosol layer $\beta$ is optically thick and vertically compact. Figure 8 illustrates the aerosol distributions for the SN_nom model relative to the 2L_DISORT and 2L_hazeA solutions.

Consistent with the results of the tests described in Appendix B, the values of the optical parameters ($\omega_{0},g$) retrieved by the SN_nom model differ from both the 2L_DISORT and 2L_hazeA values. The SN_nom optical parameters should be interpreted as the two-stream equivalent to the pyDISORT optical parameters for $\mu=0.8$ only. As illustrated by Fig. 13, the SN_nom fit quality is excellent for the $\mu=0.8$ data for which the retrieval was performed, but degrades towards higher emission angles. While the optical parameters from the SN_nom model should not be extended to more accurate radiative transfer formulations (or other viewing geometries), the parameters describing the aerosol structure should be robust. Furthermore, the SN_nom retrieval is expected to be more accurate than the 2L_hazeA retrieval, since the 2L_hazeA solution represents a forced compromise in which two stream attempts to accommodate the range of viewing geometries, for which a single set of two-stream $\omega_{0},g$ values is likely not appropriate.

5.2 Varying the CH4 profile

We now use the SN dataset to explore the relationship between our model fits and our choice of vertical methane profile. Specifically, we wish to determine whether our data indicate a preference for particular values of mCH4,t, mCH4,s, and $RH_{ac}$, as defined in Appendix A.1.2; and how the values of these parameters influence the retrieved aerosol structure. For these retrievals, we utilize two stream for the radiative transfer. We assume a two-layer aerosol structure, allowing $P_{max}$, $\tau$, and $h_{frac}$ for each layer to vary, while fixing $\omega_{0}$ and $g$ to the best-fit values from the SN_nom retrieval. We perform a set of retrievals equivalent to the SN_nom fit, but fix the methane parameters to different combinations of values from the following set: mCH${}_{4,t}=0.22,0.04$; mCH${}_{4,s}=3.5\times 10^{-4},15\times 10^{-4}$; and $RH_{ac}=40\%,100\%$. Since the Karkoschka and Tomasko, (2011) results do not indicate any methane depletion at these latitudes, we do not consider methane depletion here. The case where mCH${}_{4,t}=0.04$; mCH${}_{4,s}=3.5\times 10^{-4}$; $RH_{ac}=100\%$ corresponds to the SN_nom case and is not rerun. The seven remaining cases are labeled SN_CH4gridA-F and the results are reported in Table LABEL:table:SN.

We find that, for all eight of the CH4 profiles considered, the derived aerosol structure is qualitatively consistent: in all cases, we find aerosol layer $\alpha$ to have a base pressure of 0.4–0.7 bar and a 1.6-$\mu$m optical depth of 0.05–0.06. Aerosol layer $\beta$ is vertically compact ($h_{frac,\beta}<0.3$), has an optical depth near unity, and has a base pressure of 3–4 bar.

We also observe that the choice of CH4 profile systematically affects the retrieved probability distributions of the properties describing the distribution of aerosols, with all three of the parameters describing the CH4 profile having some effect on our retrievals. The parameters that show the strongest dependences on the adopted CH4 profile are the base pressures of both aerosol layers and the scale height of the upper aerosol layer. The retrieved base pressure of the deeper aerosol layer ($\beta$) depends on the assumed methane abundance in the troposphere, with lower values of mCH4,t resulting in higher retrieved base pressures of layer $\beta$. This result is unsurprising, since a lower methane mixing ratio implies that one must traverse to greater depths to encounter the same total methane column. A secondary trend of our retrievals is that, for a given choice of mCH4,t and mCH4,s, the retrieval using the lower methane relative humidity near the tropopause favors a deeper, but slightly more extended and optically thick aerosol layer $\beta$.

The retrieved properties of the upper aerosol layer are independent of the deep tropospheric methane mole fraction mCH4,t but do depend on the other two methane parameters: higher values of $RH_{ac}$ result in retrievals with a more vertically extended upper haze layer ($h_{frac,\alpha}>1$ when $RH_{ac}=100\%$ vs. $h_{frac,\alpha}\sim$ 0.4–0.5 when $RH_{ac}=40\%$). For a given value of $RH_{ac}$ and mCH4,t, the best-fit model of layer $\alpha$ with the higher value of mCH4,s has a deeper base, a higher total optical depth, and a higher fractional scale height. However, the differences are small ($<1.5\sigma$).

Using the $DIC$, we next evaluate whether our analysis favors a particular choice of methane profile. We find that we are not able to constrain the deep methane abundance, as parameterized by mCH4,t, at all: pairs of models which differ only in this methane parameter – for example, SN_CH4gridA and SN_CH4gridE – have mutually consistent values of the DIC. This is unsurprising, considering that the transmission and contribution functions for the wavelengths of our data (Fig. 2) indicate that the pressures influenced by mCH4,t contribute little at the wavelengths of our observations. Our poor constraint on mCH4,t translates into an uncertainty in $P_{max,\beta}$ beyond what is captured by the statistical uncertainty reported for an individual retrieval.

The methane parameter that we most strongly constrain is $RH_{ac}$: we find that all four models with $RH_{ac}=40\%$ have significantly lower (better) values of the $DIC$ than any of the four models with $RH_{ac}=100\%$. As shown in Fig. 14, the fit improvement is most evident in the reflectivity peak near 1.6 $\mu$m. To further investigate Neptune’s relative humidity, we perform two additional SN retrievals in which we allow $RH_{ac}$ to be a free parameter (referred to as SN_CH4RHfreeA,B). For these retrievals, mCH4,t remains fixed at 0.04, and we assume mCH4,s is either $3.5\times 10^{-4}$ (case A) or $15\times 10^{-4}$ (case B). The results are summarized in Table LABEL:table:SN. We find that both of the RHfree retrievals show an improvement over the SN_CH4grid models, with $\Delta DIC=-252$ and $-258$ from the SN_nom retrieval, and $\Delta DIC<-100$ from the best SN_CH4grid retrieval. The retrieved relative humidities are near zero ($\lessapprox$2–3%). Such a low methane relative humidity is inconsistent with, e.g, Karkoschka and Tomasko, (2011); however, we note that our observations are most sensitive to the CH4 abundance at altitudes above the CH4 condensation level ($\sim 1.5$ bar), whereas Karkoschka and Tomasko, (2011) are most sensitive to the methane abundance at deeper pressures – this is discussed further in Section 8. By comparing the SN_CH4gridE,F and SN_CH4RHfreeA,B retrievals (Table LABEL:table:SN, Fig. 8), we note that a decrease of $RH_{ac}$ from 40% to $\lessapprox$2–3% has little influence on the retrieval of the aerosol distribution: with the exception of $\tau_{\alpha}$, the aerosol parameter distributions for these four retrievals agree at the $\sim 1.5\sigma$ level.

Our SN_CH4grid retrievals appear to indicate a dependence of the DIC on the value of mCH4,s, with models having a lower stratospheric methane abundance resulting in better fits for a given tropospheric methane abundance and relative humidity. However, the two SN_CH4RHfree runs, differing only in the value of mCH4,s, have equivalent DIC values ($|\Delta DIC|<10$). This leads us to suspect that the preference for a smaller value of mCH4,s in the CH4grid models may not actually be an indication of a lower value of the methane mixing ratio throughout the stratosphere, but could instead reflect the preference for a lower methane abundance near the tropopause. We try one additional retrieval (SN_CH4stratfree) in which mCH4,t and $RH_{ac}$ are fixed (at 0.04 and 0.4, respectively) but mCH4,s is added as a free parameter. The retrieved value of mCH4,s is only $2.2e-4$, lower than our nominal choice of $3.5e-4$ which itself is at the low end of the values determined by previous authors. However, the DIC indicates no improvement over the SN_CH4gridE retrieval, and the aerosol structure is effectively unchanged by decreasing mCH4,s from $3.5e-4$ to $2.2e-4$.

In the retrievals that follow, we adopt the SN_CH4gridE methane profile as our new nominal model. This profile has a CH4 relative humidity of 0.4 near the tropopause; as discussed above, our RHfree retrievals indicate a humidity that is lower by more than an order of magnitude. While we have no reason to reject this finding, we adopt the more moderate value of RHfree for consistency with previous studies (Karkoschka and Tomasko, , 2011; Irwin et al., , 2011), to simplify comparison with these works. As noted previously, this choice has little effect on the retrieval of the aerosol distribution, with the exception that models with the higher value of $RH_{ac}$ have higher – but qualitatively similar – values of $\tau_{\alpha}$. We revisit Neptune’s methane profile in Section 6.2.

5.3 Varying the thermal profile

We perform a final set of model runs on the SN data set, to investigate the influence of possible errors in the adopted thermal profile. In our model, the thermal profile influences the gas opacity coefficients as well as the CH4 profile, through its effect on the saturation vapor pressure curve of methane. Using the SN_CH4gridE methane profile, we consider three perturbed thermal profiles. In the SN_Ttrop case, we adjust the tropopause temperature from its nominal value of 54.9 K to the value of 56.4 K, as found by Orton et al., (2007) (‘case B’). In SN_Tstratp20 and SN_Tstratm20 we increased/decreased, respectively, the temperature in the stratosphere by 20 K at 1 mbar. The details of these adjustments are described in Appendix A.1.1. The results of these model runs are shown in Table LABEL:table:SN. We find that the posterior probability distributions for all parameters are consistent with (within $1\sigma$ of) those found for the unperturbed thermal profile, with the exception that increasing the temperature of the tropopause results in a higher value of $h_{frac,\alpha}$. The 2L_Tstratp20 retrieval has a $\Delta DIC=-20$ relative to the SN_CH4gridE retrieval; the other two modified thermal profiles result in a decreased fit quality relative to the SN_CH4gridE retrieval.

6 Individual cloud-free regions

We follow our intensive analysis of the limb darkening data set with a broader study of many cloud-free regions across Neptune, spanning latitudes from 20∘N to 87∘S and viewing geometries from $\mu=0.26$ to $\mu=0.99$. This analysis permits the consideration of spatial variations in the distribution of aerosols, temperature, and methane abundance. It also allows us to evaluate the relative information content of a single spectrum at a single viewing geometry relative to the extensive limb-darkening dataset.

Our approach for selecting individual locations for study is as follows: for every 10∘ in latitude, we seek two locations that are 1) dark (cloud-free), 2) similar in latitude but with different viewing geometries, 3) at longitudes for which we have both Hbb and Kbb data, and 4) at least 10 pixels from any bright clouds, to avoid contamination due to the broad PSF halo (Section 2.4). At some latitudes (e.g., 20∘N and 60∘S), nearby bright bands lead us to relax the final criterion; at other latitudes (e.g., 30∘S) these bright bands prevent the identification of any suitably dark locations. Further, at high southern latitudes, the prevalence of bright clouds has resulted in a more uneven latitude spacing of selected suitable dark locations. Our final sample consists of 18 specific locations, shown and labeled A–R in Fig. 1. Two of the locations (C and D) are within the original limb-darkening latitude band from 2–12∘N.

To extract a spectrum for each of these locations, we first extract the Hbb spectrum from a single pixel and record the latitude and cosine of the emission angle, $\mu_{H}$, for that pixel. We then extract a Kbb spectrum corresponding to a similar physical location on Neptune in the following way: we determine the time interval between the Hbb and Kbb observations of the relevant latitude. Next, we estimate the expected rotation for a patch of atmosphere at the given latitude using the planetary rotation rate and the mean zonal wind profile from Sromovsky et al., (1993). Finally, we estimate the Kbb X,Y location using the Hbb X,Y location and latitude, the centers of the Hbb and Kbb cubes, and the expected amount of rotation. The data from the two bands are considered a single ‘location’ for the purposes of atmospheric retrieval, and are modeled simultaneously, with the caveat that we use the $\mu_{H}$ and $\mu_{K}$ value appropriate for each band. We note that the feature coordinates are imprecise due to the $\lesssim 1$ pixel uncertainty in the image centering. A one-pixel error in position would correspond to roughly 3∘ error in latitude/longitude for a feature near disk center, with larger errors possible at other viewing geometries (see Fig. 4 of Martin et al., , 2012). Large-scale meridional trends retrieved by our analysis should be robust to possible navigation errors. Our retrieval procedure also ignores the finite width of the PSF for our observations: models are produced for the singular $\mu_{H}$ and $\mu_{K}$ corresponding to the expected viewing geometry of the extracted data pixels, neglecting any changes in the measured intensity due to contributions from the wings of the PSF. To estimate the errors introduced by this approximation, we compared a synthetic, unsmoothed model cube to one convolved with the PSF, as estimated from stellar spectra. We find that for the three locations at the highest Hbb viewing angles (B,C,P), the mean error introduced by neglecting the finite width of the PSF is comparable to the mean uncertainty due to random noise ($\bar{\sigma_{n}}\approx 0.001$). Further, at locations with moderate Hbb viewing angles (F,H,N,Q,R), the maximum error introduced by this approximation is greater than $\bar{\sigma_{n}}$. We do not incorporate an estimate of this error – which is wavelength, location, and model-dependent – into our analysis: doing so would decrease the unknown component of the uncertainty ($\sigma_{T}$) for high viewing angle models, and would in some cases alter the relative weighting of the data points in the retrievals: errors from ignoring the width of the PSF are greatest in the 1.6-$\mu$m I/F peak.

The spectra corresponding to the 18 selected locations are shown in Fig. 15. Visual inspection of these plots indicates a north-south trend in the Hbb spectral shape and peak intensity: spectra at northern latitudes tend to have higher 1.5 $\mu$m intensities and lower peak intensities than spectra from southern latitudes.

6.1 Local aerosol distribution

Using the nominal thermal profile, the SN_CH4gridE methane profile, and the particle size distributions from Section 4.2, we retrieve the aerosol properties for each location A–R. Since the 18 cloud-free locations span a broad range of viewing geometries, we cannot utilize two stream for the radiative transfer; we instead use pyDISORT and adopt the best-fit values of $g$ and $\omega_{0}$ from the 2L_DISORT retrieval (Section 4.4). The free parameters of these retrievals are the base pressure, optical depth, and fractional scale height of the two aerosol layers. The results are summarized in Table LABEL:table:DR and shown in Figs. 15 – 20.

While the majority of the 18 retrievals result in well-behaved posterior probability distributions (e.g. Fig. 16), three retrievals (locations I, M, and O) have posterior probability distributions that are double-peaked in the parameters that describe aerosol layer $\alpha$ while otherwise meeting our criteria for convergence (see Fig. 17 for an example). We cannot rule out that with additional iterations, emcee may have converged on a single solution. The values quoted in Table LABEL:table:DR are, as always, the 16th, 50th, and 84th percentiles for the marginalized distributions, which, for double peaked distributions, do not clearly identify the highest probability values. For two of these locations (M, O), one of the two posterior probability distribution peaks corresponds to an aerosol layer based near the maximum allowed $P_{max,\alpha}$ value of 2.3 bar, such that layer $\alpha$ substantially overlaps with aerosol layer $\beta$ (Fig. 17). Interpretation of a solution with such strongly overlapping aerosol layers is not straightforward and caution is recommended for these locations. Figs. 16 and 17 also highlight the strong correlations between parameters.

The plots of optical depth per bar (Fig. 18) and cumulative optical depth (Fig. 19) illustrate the range of likely aerosol distributions for each of the 18 locations. In these plots, we also show the best-fit solution from the 2L_DISORT retrieval for reference. It is clear from these plots that a single spectrum provides a much weaker constraint on the aerosol distribution, relative to the more extensive limb-darkening dataset considered in the 2L_DISORT retrieval – in particular the fractional scale height of the upper aerosol layer is poorly constrained at many locations. However, we do observe a number of interesting trends in the single-spectrum retrievals: at every location, the deep aerosol layer ($\beta$) dominates in terms of total optical depth, and is generally observed to be physically compact ($h_{frac,\beta}<0.2$) and optically thick. The exceptions are locations A and B in the north, and perhaps location F at the equator, which have best-fit values of the optical depth $\tau_{\beta}<1$. The base pressure of this bottom cloud layer is always below 2.2 bar; since the cloud is generally found to be optically thick, solutions with higher base pressures (and correspondingly larger total optical depths) are allowed in many locations.

The upper aerosol layer ($\alpha$) exhibits a more pronounced latitudinal variability than the deeper cloud. At location A, at 20∘N, the aerosols increase in concentration with altitude ($h_{frac,\alpha}\sim 2$), indicating the presence of aerosols above the tropopause. Locations B–I, between 16∘N and 17∘S have aerosol scale heights similar to the gas scale height ($h_{frac,\alpha}\sim 1$) and the aerosol distribution is consistent with the 2L_DISORT best-fit solution. Location J is an outlier, with a preference for a vertically compact aerosol layer near 0.4 bar, rather than a more extended haze. This location may be influenced by reflectivity from a nearby bright cloud – see Fig. 1. Further south (locations K–R, 43–87∘S), the optical depth of aerosols in the stratosphere and upper tropopause drops below the reference 2L_DISORT solution (Fig. 19), and plots of $\tau$/bar (Fig. 18) show no gap between aerosol layer $\alpha$ and the deep cloud. The details of the retrieved solution may depend somewhat on viewing geometry: locations M and O, for example, are at relatively low emission angles ($\mu_{H}>0.6$), and have more compact solutions for the upper aerosol layer than neighboring locations N and P, which are nearer the limb.

Figure 20 presents another visualization of the latitudinal trend in the aerosol structure. In this set of plots, we show the cumulative optical depth of aerosols at 0.1, 0.5, 1, and 2 bar, as a function of latitude. Above a pressure of 0.1 bar, we observe that the cumulative optical depth in aerosols south of 20∘S is very low. At low latitudes we observe a large scatter in the retrieved optical depths; however, much higher cumulative optical depths are permitted. Above 0.5 and 1 bar, there is a clear difference between cumulative optical depths observed for locations north and south of 20∘S. By a pressure depth of 2 bar, the cumulative optical depth in aerosols is similar near the equator and mid- and high southern latitudes.

6.2 Local methane profile

Next, we revisit Neptune’s methane profile, to search for latitudinal variations. Using Hubble STIS spectroscopy, Karkoschka and Tomasko, (2011) found evidence of a shallow depletion of CH4 below $\sim 1.2$ bar at some latitudes. To investigate whether our data exhibit signatures of similar behavior, we added to the parameterization of the previous section a methane depletion parameter, $P_{dep}$, as defined in Appendix A.1.2 and Table LABEL:table:model, which controls the CH4 profile below the methane condensation pressure of $\sim 1.7$ bar. The results of this set of retrievals were difficult to interpret: retrieved values of the methane depletion parameter ranged from 2–3 bar to more than 50 bar – well below our optically thick cloud $\beta$. There was no observable latitudinal trend in the retrieved parameters. Furthermore, we found very different depletion depths for spectra from similar latitudes, and high values of depletion near the equator, where it was not observed by Karkoschka and Tomasko, (2011). We concluded that this model parameterization was too poorly constrained by our data to yield useful information.

Recalling the results from Section 5.2, we hypothesize that we are less sensitive to the relatively deep CH4 depletion than we are to the methane relative humidity near the tropopause. To allow for more freedom in the CH4 profile, we perform a set of retrievals in which both $RH_{ac}$ (which controls the CH4 mole fraction near the tropopause) and $P_{dep}$ (which controls the CH4 mole fraction below $\sim 1.7$ bar) are free parameters. For this test, we fix the base pressure, optical depth, and fractional scale height of layer $\beta$ to be $P_{max,\beta}=3.3$ bar, $\tau_{\beta}=5$, and $h_{frac,\beta}=0.1$. In other words, we assume that the properties of the bottom cloud are spatially constant, and that any observed variations in the 1.6 $\mu$m spectral peak are due to spatial variations in the methane column above that cloud rather than due to variations in aerosol layer $\beta$. We also fix $h_{frac,\alpha}$ to be $0.8$; while there is no physical motivation for forcing the scale height of the top aerosol layer to be constant, doing so makes it easier to interpret the remaining free parameters in our model. As in the previous set of retrievals, $\omega_{0}$ and $g$ for both layers are fixed to the 2L_DISORT best-fit values. The results of these retrievals are summarized in Table LABEL:table:CH4fit.

Of our 18 locations, we find that three locations (A, G, and J) are less well modeled by this parameterization relative to the original aerosol-only parameterization. Location A, demonstrates the greatest decrease in fit quality ($\Delta DIC=21$); the retrieved models for this location are in particularly poor agreement with the data near 1.6 $\mu$m. The likely cause of the decrease in fit quality is that location A is not well-matched by the assumption that cloud $\beta$ is optically thick. For location J, the initial aerosol-only retrieval finds layer $\alpha$ to be vertically compact; the decreased fit quality in the latter retrieval may be due to forcing $h_{frac,\alpha}=0.8$. We do not identify an obvious reason for the decrease in fit quality at location G. Of the remaining 15 locations, 13 are significantly better matched by the retrievals that allow the methane to vary, whereas the $DIC$ at two locations are comparable ($|DIC|\leq 10$) for the aerosol-only and variable-methane parameterizations. For the locations demonstrating an improvement in fit quality, this improvement is greatest near 1.5–1.6 $\mu$m (Fig. 21).

At every location in our analysis, the variable-methane retrieval favors a tropospheric CH4 mole fraction below our nominal SN_CH4gridE methane profile (see Supplementary Fig. 30). This is achieved through a combination of subsaturation near the tropopause and depletion below the CH4 condensation pressure. In the majority of cases, we find a best-fit value of $RH_{ac}$ below 10% coupled with CH4 depletion to a depth of $P_{dep}\approx$2.0–2.5 bar. For four locations, an alternative solution is observed: the relative humidity remains higher (best-fit values of 20–40%) but the retrieved depletion depth is higher as well. For three of these four cases (locations J, K, and P), the posterior probability distributions are actually double peaked, with one local minimum at the higher relative humidity and higher depletion pressure and a second at lower relative humidity and lower depletion pressure. As before when we had double-peaked posterior probability distributions, we note that the MCMC chains may have converged with additional steps. In Fig. 22 (left) we plot depletion pressure as a function of latitude. While latitudinal variations in the depletion depth are not ruled out, differences between the two families of solutions (low $RH_{ac}$ and $P_{dep}$ vs. higher $RH_{ac}$ and $P_{dep}$) dominate the scatter in $P_{dep}$. Latitudinal trends are more compelling when both relative humidity and depletion below the condensation level are considered: in the righthand panel of Fig. 22, we show the CH4 abundance above 2.5 bar for the retrieved CH4 profiles. This plot shows that both families of solutions represent similar CH4 columns above 2.5 bar; at pressures greater than 2.5 bar, the deep cloud generally becomes optically thick, and we are not sensitive to variations in the CH4 abundance. We also indicate in this plot the CH4 columns for the nominal SN_CH4gridE methane profile ($RH_{ac}$=0.4) and for the SN_CH4RHfreeA retrieval from Section 5.2, in which $RH_{ac}$ (but not $P_{dep}$) was allowed to vary ($RH_{ac}\approx 0.02$). The retrieved CH4 column at many locations near the equator is roughly consistent with the SN_CH4RHfreeA CH4 column. While not conclusive, this plot may suggest a local minimum in the upper tropospheric CH4 abundance just north of the equator, and a decreased abundance at high southern latitudes relative to low latitudes.

Finally, we observe that the variable-methane parameterization does not qualitatively influence the trend in the upper aerosol layer $\alpha$ observed in the previous set of retrievals: this is evident by comparing Fig. 23 to Fig. 18. However, the optical depth of the upper aerosol layer is, in general, lower in the variable-methane retrieval than in the aerosol-only retrieval, with best-fit values of $\tau_{\alpha}<0.01$ at a number of locations in the south for the variable-methane case.

6.3 Thermal variations

In Section 5.3 we found that perturbations to the stratospheric temperature did not cause appreciable changes to the posterior probability distribution of the aerosol model parameters. We also found that the 1.5 K adjustment to the tropopause temperature from its nominal value of 54.9 K to the Orton case B value (56.4 K) had a negligible effect on the retrieved aerosol structure. We now revisit the question of possible effects due to variations in the tropopause temperature, for three specific locations at which the local temperature is expected, from Orton et al., (2007), to exhibit large deviations from the nominal value. Two of the locations (P, R) are at high southern latitudes and should have elevated temperatures relative to the nominal value.The third location (M) is at mid southern latitudes and should have a low tropopause temperature, relative to the nominal value. We repeat the initial (aerosol-only, no methane depletion) fits for these three locations, and present the retrieved parameters for these revised thermal profiles in Table LABEL:table:DR_TP. None of the three locations exhibit an improvement in the DIC for retrievals with the revised thermal profile. The posterior probability distributions from the nominal retrievals and the temperature-perturbed retrievals are effectively the same for all three locations.

7 Summary of findings

We present a comprehensive analysis of OSIRIS Hbb and Kbb integral field spectrograph observations of the planet Neptune from 26 July 2009, focusing on regions free of discrete NIR-bright clouds. Our analysis involves a series of atmospheric retrievals using an atmospheric model and radiative transfer code, coupled to an MCMC algorithm that provides an estimate of the posterior probability distribution of the model parameters. The Deviance Information Criterion (DIC) is used to compare different atmospheric models for the same set of data. Taking advantage of the high spatial resolution of these data, we first construct a set of six high S/N spectra, representing the atmosphere in a latitude range of 2–12∘N at each of six values of the cosine of the emission angle. Our atmospheric retrievals using this high-quality limb darkening dataset reveal the following:

• •

Neptune’s cloud opacity at these wavelengths and latitudes is dominated by a compact cloud layer with a base near 3 bar. Using pyDISORT for the radiative transfer and assuming a Henyey-Greenstein phase function, we find that the base of this layer is located at $P_{max,\beta}=3.3^{+0.4}_{-0.3}$ bar, and the layer is optically thick ($\tau_{\beta}=5^{+4}_{-2}$ at 1.6 $\mu$m for the 2L_DISORT retrieval). Using pyDISORT, we observe this cloud to be composed of low albedo ($\omega_{0,\beta}=0.45^{+0.01}_{-0.01}$), moderately forward scattering ($g_{\beta}=0.50^{+0.02}_{-0.02}$) particles, compared to the Mie theory prediction of $\omega_{0}=1.0$ and $g=0.76$ at 1.6 $\mu$m for $r_{p}=1.0$ $\mu$m.

• •

A second aerosol layer, at lower pressures but also with a pressure base in the troposphere, is required in our models to match the data. The 2L_DISORT model retrieves $P_{max,\alpha}=0.59^{+0.04}_{-0.03}$ bar, $\tau_{\alpha}=0.019^{+0.002}_{-0.001}$, and $h_{frac}=0.85^{+0.07}_{-0.06}$ for the upper aerosol layer, indicating nearly uniform mixing of the aerosols with the gas, up into the stratosphere. For the wavelength-dependent cross section of upper aerosol particles, a particle size distribution of $r_{p}=0.1$ $\mu$m is preferred over one dominated by larger ($r_{p}=0.5$ or 1.0 $\mu$m) particles. Using pyDISORT and a Henyey-Greenstein phase function, we retrieve $\omega_{0,\alpha}=0.91^{+0.06}_{-0.05}$ and $g_{\alpha}=0.24^{+0.02}_{-0.03}$.

Focusing on the highest S/N spectrum from this set of six, we investigate Neptune’s methane and thermal profiles. We find:

• •

Our retrievals indicate a strong preference for a methane relative humidity of $RH_{ac}=$40% over 100% relative humidity. When $RH_{ac}$ is allowed to be a free parameter, the retrieved relative humidity is less than 2% ($0.016^{+0.01}_{-0.008}$ for the SN_CH4RHfreeA retrieval). Our retrievals indicate a possible preference for a stratospheric CH4 mixing ratio at the low end of the range found by previous authors ($3.5\times 10^{-4}$ or less); however, based upon our RHfree retrievals, we interpret this result as a preference for models with a decreased CH4 abundance near the tropopause, rather than a sensitivity to the CH4 abundance in the upper stratosphere. We are not able to constrain the value of the tropospheric CH4 abundance, $mCH_{4,t}$, with our data. We caution against over-interpreting the meaning of the values of the methane parameters retrieved by this study, in light of the very simple parameterization used for the methane profile (Appendix A.1.2).

• •

Overall, our retrieved aerosol structure is qualitatively insensitive to the assumed methane profile, being comprised in all cases of a vertically compact bottom cloud with a base pressure 3–4 bar, and a vertically extended, optically thin upper aerosol layer with a base pressure of 0.4–0.7 bar. The parameters most affected by the assumed methane profile are the base pressures of both aerosol layers and the scale height of the upper aerosol layer. The base pressure of the deeper aerosol layer depends on the tropospheric methane mixing ratio, with the highest base pressures corresponding to the lowest tropospheric methane mixing ratios. This result is not surprising, since a lower methane mixing ratio implies that greater distances must be traveled in order to encounter the same total methane column. The base pressure of the upper aerosol layer is slightly but systematically influenced by the value of the stratospheric methane mole fraction, with a higher values of mCH4,s resulting in a higher value of $P_{max,\alpha}$. The compactness of the upper aerosol layer, as parameterized by $h_{frac,\alpha}$, depends most strongly on the relative humidity of CH4 near the tropopause, with higher values of $RH_{ac}$ resulting in a more vertically extended upper haze layer.

• •

For the same methane profile, a model with a higher stratospheric temperature results in an improved fit, with no change to the retrieved aerosol parameters. An increase to the tropopause temperature or decrease to the stratospheric temperature result in worse fits. As with CH4, we consider only a simple parameterization to the thermal profile and only a few different test cases.

We follow this analysis of the strip from 2–12∘N with a study of 18 cloud-free locations across Neptune, spanning latitudes from 20∘N to 87∘S and viewing geometries from $\mu=0.26$ to $\mu=0.99$. We use these data to evaluate the information content of a single spectrum at a single viewing geometry relative to the more complete and higher S/N limb-darkening dataset, and to explore latitude variations in atmospheric properties. We find:

• •

While we do not attempt to fit the optical properties of the aerosols using single spectra, retrievals of the other aerosol parameters ($P_{max}$, $\tau$, and $h_{frac}$) for locations C and D from within the limb-darkening study region (latitudes of 10∘N and 5∘N, respectively) are remarkably consistent with – though not strictly within 1$\sigma$ of – the values found by the full 2L_DISORT retrieval. The posterior probability distributions for the layer $\alpha$ base pressure and optical scale height have greater uncertainties in the fits to individual locations than observed for the limb darkening dataset retrievals. This is likely due to the fact that most of the information about the upper aerosol layer is derived from low S/N portions of the spectrum.

• •

The properties of the deep aerosol layer $\beta$ are qualitatively consistent across all locations and with the full 2L_DISORT retrieval. At latitudes of 10∘N–20∘S, retrieved values of $P_{max,\beta}$ are within $1\sigma$ of the $3.3^{+0.4}_{-0.3}$ bar found by 2L_DISORT; outside of this latitude range, the cloud base is typically at a slightly lower pressure, with best-fit values of 2.2–2.7 bar. As in the 2L_DISORT case, we observe this cloud to be optically thick, with the exceptions of locations A and B at 20 and 16∘N, respectively, and possibly location F at the equator, for which this bottom cloud may be optically thin. At all locations, the retrieval favors a highly compact layer $\beta$, with $h_{frac,\beta}$ values similar to or less than the value retrieved in the full 2L_DISORT case.

• •

The structure of aerosol layer $\alpha$ varies with latitude: at low latitudes (20∘N – 20∘S), we observe a vertically extended haze with a base in the troposphere above layer $\beta$, with $P_{max,\alpha}$ typically 0.3–1.4 bar. A gap is observed between the top of layer $\beta$ and the base of layer $\alpha$, and the upper layer typically extends into the stratosphere. At mid- and high-southern latitudes, the base of aerosol layer $\alpha$ is deeper (from 1.7 to $>2$ bar), such that this aerosol layer overlaps in most cases with the top of layer $\beta$. We retrieve much lower aerosol optical depths in the upper troposphere for mid- and high-southern latitudes relative to low latitudes.

• •

Fits in which we fix all aerosol parameters except for $P_{max,\alpha}$ and $\tau_{\alpha}$, but allow two parameters describing the CH4 profile – $RH_{ac}$ and $P_{dep}$ – to vary, result in improved fits to the spectra from most locations. The predominant solution has $RH_{ac}<10\%$ and $P_{dep}\approx$ 2.0–2.5 bar. At four locations, our retrievals favor a solution with a higher methane relative humidity but a deeper methane depletion. These results indicate that the CH4 abundance in the upper troposphere is lower than our nominal CH4 profile at all locations, with tentative evidence of meridional variations in the methane column above 2.5 bar. Our results may also suggest that our simple parameterization does not fully characterize the true shape of Neptune’s CH4 profile and variations within.

• •

Thermal variations near the tropopause, reaching 7.5K near the south pole (Orton et al., , 2007) do not appear to have a dominant effect on the retrieved aerosol properties.

8 Discussion

8.1 Comparison with previous work

To our knowledge, this work represents the most thorough NIR characterization of Neptune’s cloud-free regions. Resolving the aerosol structure and composition in these regions is important for characterizing Neptune’s large-scale meridional circulation patterns and determining the processes responsible for aerosol production. These efforts also serve as context for characterizing localized storm activity that occurs across the planet. Comparison of this analysis with previous studies is complicated by differences in the information content of the data (due to different wavelength coverage and spatial extent) and differences in the modeling assumptions. Furthermore, temporal variations in Neptune’s aerosol structure and atmospheric composition are likely: Voyager, along with complementary Hubble and ground-based observations since, have shown that Neptune’s visible and near-infrared (NIR) appearance exhibit both rapid, localized change (Limaye and Sromovsky, , 1991) and decadal-scale brightness trends (Lockwood and Jerzykiewicz, , 2006; Hammel and Lockwood, , 2007; Karkoschka, 2011a, ). As a result, discrepancies between studies separated by several years may reflect real changes in atmospheric properties over time.

Despite these complications, we note a number of key similarities and differences between our results and previous efforts. In contrast to Karkoschka and Tomasko, (2011), who fit their 2003 HST-STIS (300–1000 nm wavelength) observations with a semi-infinite haze below 1.4 bar, our analysis strongly favors a compact and generally optically thick cloud layer at $\sim 3$ bar. This result is consistent with earlier work (e.g., Hammel et al., , 1989; Baines and Smith, , 1990; Baines et al., 1995b, ; Baines et al., 1995a, ; Roe et al., 2001a, ), and more recently, with an analysis of Gemini NIFS (1.477–1.803 $\mu$m wavelength) observations by Irwin et al., (2011, 2014). As noted in the literature (e.g., by Sromovsky et al., 2001b, ; Roe et al., 2001a, ), an optically thick 3-bar cloud must be dark in the NIR to match observations: we retrieve a NIR single scattering albedo of $0.45^{+0.01}_{-0.01}$ for this layer assuming a Henyey-Greenstein phase function and using pyDISORT for the radiative transfer. This is intermediate between the values of 0.1 – 0.2 found by Sromovsky et al., 2001b ; Roe et al., 2001a and the value of $\sim 0.75$ assumed by Irwin et al., (2011, 2014). We note that the retrieval of $\omega_{0,\beta}$ is correlated with the asymmetry parameter of the scattering phase function (Fig. 11), which may partly explain the discrepancy: for example, Irwin et al., (2011) fix the asymmetry parameter $g$ of their Henyey-Greenstein phase function to be 0.7, whereas Sromovsky et al., 2001b assume an isotropic phase function. We retrieve a value of $0.50^{+0.02}_{-0.02}$ for the Henyey-Greenstein asymmetry parameter.

Previous aerosol models of Neptune’s NIR-dark regions typically include a second, optically thin tropospheric aerosol layer above the $\sim 3$ bar cloud deck, consistent with our findings (Pryor et al., , 1992; Baines and Hammel, , 1994; Sromovsky et al., 2001a, ; Gibbard et al., , 2002; Irwin et al., , 2011, 2014). One exception to this is Karkoschka and Tomasko, (2011); these authors find the troposphere to be clear above 1.4 bar. However, their visible-wavelength observations are most sensitive to altitudes below methane condensation. In the Karkoschka, 2011a analysis of HST/WFPC2 images from 1994 to 2008, the author allows for a possible upper tropospheric haze between the CH4 condensation level and the tropopause, which is 10 times thinner than the main tropospheric aerosol layer – similar to the model proposed by Baines and Hammel, (1994). Our best-fit model for the limb-darkening dataset from 2–12∘N likewise suggests the presence of a vertically extended upper tropospheric aerosol layer, although we find a somewhat higher bottom pressure of $P_{max,\alpha}=0.59^{+0.04}_{-0.03}$ (for the 2L_DISORT retrieval), roughly a scale height above the CH4 condensation level. Other models (e.g., Sromovsky et al., 2001a, ; Gibbard et al., , 2002; Irwin et al., , 2011) have presupposed that the upper tropospheric aerosol layer is vertically compact; our retrievals do not favor this possibility. The reanalysis of the Irwin et al., (2011) data by Irwin et al., (2014) shows that a vertically extended upper haze is equally allowed by their data. Our analysis suggests two likely contributors to the variations in solutions across previous studies: first, the results of our methane profile tests show that the retrieved base pressure and scale height of the upper aerosol layer depend strongly on the adopted CH4 profile. Secondly, many studies consider either a global average of Neptune’s aerosol structure, or analyze only one specific latitude band. Our analysis demonstrates that this upper aerosol layer is spatially variable. Whereas Irwin et al., (2011, 2014) use a single scattering albedo of 0.45 and asymmetry parameter of 0.7 for the upper aerosol layer, we retrieve values of $\omega_{0,\alpha}=0.91^{+0.06}_{-0.05}$ and $g_{\alpha}=0.24^{+0.02}_{-0.03}$. These numbers are similar to what one would expect from Mie theory for small particles ($r_{p}\approx 0.1\mu$m). The single scattering albedo that we retrieve for the upper haze is also more consistent with what is observed for the particles in Neptune’s discrete clouds (e.g., Irwin et al., , 2011).

Despite some discrepancies between the Karkoschka and Tomasko, (2011); Karkoschka, 2011a vertical aerosol profile and our results, these previous efforts provide a valuable point of comparison for our analysis of latitudinal variations in Neptune’s discrete-cloud free atmosphere. For one, the Karkoschka and Tomasko, (2011); Karkoschka, 2011a analyses represent the most comprehensive previous investigations of Neptune’s dark regions. Secondly, these studies, which utilize 300–1000 nm HST data, provide complementary constraints on Neptune’s upper troposphere and stratosphere to our Hbb and Kbb NIR data. Whereas our observations are most sensitive to the aerosol structure between mbar pressures and the top of the $\sim 3$ bar cloud deck, the peak sensitivity of the Karkoschka and Tomasko, (2011) data is below the CH4 condensation level. Likewise, Karkoschka and Tomasko, (2011) are most sensitive to Neptune’s CH4 distribution at deeper pressures.

Karkoschka and Tomasko, (2011) identify several prominent features in Neptune’s reflectivities: the first is an increase in reflectivity from low latitudes to high southern latitudes, with the transition occurring around $-30^{\circ}$ latitude. A similar trend is clearly evident in our own dark-region spectra (Fig. 15). Karkoschka and Tomasko, (2011) explain this feature as the result of latitude variations in the methane mixing ratio around the 2-bar layer. Our analysis shows a signature, at low significance, of a decreased CH4 column above the 2.5-bar pressure level at mid- and high-southern latitudes, in qualitative agreement with the Karkoschka and Tomasko, (2011) findings. However, our results suggest that the dominant cause of the variation from low latitudes to mid- and high-southern latitudes is variation in Neptune’s upper tropospheric haze: we find that the optical depth of this haze is higher at low latitudes ($-20$ to $+20^{\circ}$) than at more southern latitudes. The second feature observed by Karkoschka and Tomasko, (2011) is a narrow dark belt south of Neptune’s equator, which these authors attribute to a dip in the optical depth of Neptune’s tropospheric haze at this latitude. While our sensitivity to variations in the opacity of the deep tropospheric cloud deck is limited, we note that the three dark locations in our data for which we retrieve a bottom cloud opacity $\tau_{\beta}<1$ are at 20∘S, $16^{\circ}$S, and the equator – these latitudes roughly correspond to the minima in tropospheric haze optical depth/bar found by Karkoschka and Tomasko, (2011). Further, we tentatively retrieve a local maximum in the CH4 column near 10∘S, which might correspond to the dark band identified by Karkoschka and Tomasko, (2011). In Karkoschka, 2011a , the author finds evidence for either an increase in methane relative humidity or decrease in upper tropospheric aerosols in a band from 5–$20^{\circ}$N; this may correspond to the peak in CH4 column above 2.5 bar that we find at 10–20∘N. Finally, Karkoschka and Tomasko, (2011) observe a dark belt near 60∘S in Neptune’s visible continuum, which they ascribe to a decrease in the single scattering albedo of tropospheric aerosols at long wavelengths. We do not identify this feature in our data.

Taken together, our results and the Karkoschka and Tomasko, (2011); Karkoschka, 2011a analyses highlight the complementary nature of data from the optical and NIR wavelength regimes, and also the need to consider both the CH4 and aerosol vertical profiles in order to understand meridional variations in Neptune’s feature-free atmosphere. For CH4, both studies are consistent in finding latitude variations in the methane abundance around the 2-bar layer, with the Karkoschka and Tomasko, (2011) data demonstrating this result at higher significance. Further, the Karkoschka and Tomasko, (2011) data are critical for showing that the depletion is shallow, reaching a maximum depletion depth of 3.3 bar at high southern latitudes, whereas our data are most robust at higher altitudes, showing that the CH4 abundance near the tropopause, as parameterized by relative humidity, is lower than expected at all latitudes. For the aerosols, resolving discrepancies between the Karkoschka and Tomasko, (2011); Karkoschka, 2011a studies and our own work, particularly regarding the compactness of the deep aerosol layer, is more challenging. Since the two wavelength regimes are optimally sensitive to different particle sizes and atmospheric depths, it is perhaps unsurprising that these discrepancies exist. Simultaneous observations and coincident analysis of data from both wavelength regimes would allow us to better elucidate Neptune’s true tropospheric aerosol structure, and to rule out differences caused by temporal variability.

8.2 Comparison with Uranus

Studies of this nature, which involve extracting spectra far from discrete bright clouds across a range of viewing angles and latitudes, are far easier for the relatively cloud-free Uranus. Increasingly complex models of the aerosol structure – some with 5 aerosol layers or more – have been used to describe the available Uranus data. As discussed by de Kleer et al., (2015), these models tend towards the same vertical profile of aerosols regardless of the parameterization. Furthermore, a two-layer model is sufficient to match the Uranus Hbb (and Kbb) data considered by de Kleer et al., (2015). Therefore, we focus on comparisons with a two-cloud parameterization of the Uranus aerosol structure.

Uranus is observed to have a deep, compact cloud in the 2–3 bar range (Karkoschka and Tomasko, , 2009; Irwin et al., , 2012; Tice et al., , 2013; de Kleer et al., , 2015), topped by a vertically diffuse haze with a base near 1 bar, increasing in concentration with altitude ($h_{frac}>1$). Both aerosol layers on Uranus decrease in optical depth from equator to pole, with the lower cloud exhibiting a sharp peak in optical depth near the equator and the upper haze trending more gradually (e.g., Irwin et al., , 2007; Sromovsky and Fry, , 2007; Karkoschka and Tomasko, , 2009). The latitudinal trends in aerosol structure are mirrored by an increase in methane depletion from no depletion near the equator to $P_{dep}>10$ bar near the poles (e.g., Karkoschka and Tomasko, , 2009; Sromovsky et al., 2011a, ). Considering the marked differences in (discrete) cloud activity between these two planets, the similarities between the background aerosol structures – in particular, in terms of base pressures and aerosol scale heights – are notable. The total optical depth of Neptune’s upper aerosol layer is greater than observed on Uranus, consistent with Neptune’s higher abundance of CH4 and other hydrocarbons in the upper atmosphere. Like for Uranus, we observe a dramatic decrease in the optical depth of Neptune’s upper haze from the equator towards the south pole. However, we do not observe strong signatures of decreased opacity in Neptune’s bottom cloud like are observed for Uranus (Sromovsky et al., 2011a, ; de Kleer et al., , 2015). While CH4 depletion towards the poles may be present for Neptune, we do not find evidence for the relatively deep tropospheric methane depletion observed for Uranus. Further, the “proportionally descended gas” model of CH4 depletion (Sromovsky et al., 2011a, ), which matches the Uranus data well, may not be appropriate for Neptune: our data suggests that Neptune’s atmosphere is relatively dry at altitudes above the CH4 condensation pressure, and models including both variations in relative humidity (above the CH4 condensation level) and CH4 depletion (below the condensation level) are better at reproducing the observations than models in which only depletion is considered.

8.3 Implications for global circulation pattern

A primary goal of studies of Neptune’s aerosol structure and composition is to understand the global circulation patterns. Enhanced cloud activity at midlatitudes, patterns in temperature, ortho/para H2 ratios, and mid-infrared and radio brightness temperature measurements have led to a picture in which air rises above southern and northern midlatitudes, and dry air sinks over Neptune’s equator and poles (Bezard et al., , 1991; Conrath et al., , 1991; de Pater et al., , 2014). de Pater et al., (2014) suggest that this circulation extends all the way from the deep troposphere up into the stratosphere. This circulation pattern contrasts dramatically with the circulation proposed for Uranus, in which air is expected to rise at the equator and sink at the poles, and may be structured in three vertically stacked layers rather than extending across many scale heights as a single circulation cell (Sromovsky et al., , 2014).

Our retrieved latitude trends in Neptune’s background aerosol structure and methane profile are hard to reconcile with the current model of its circulation. An upper tropospheric/stratospheric haze at low latitudes, which is either confined to deeper pressures or not present at all at mid-and high southern latitudes, does not seem to be a natural outcome of a two-cell-per-hemisphere circulation pattern that extends from the deep troposphere to upper stratosphere. It is also surprising, in the context of this presumed circulation pattern, that the deep aerosol layer on Neptune appears to be so consistent in its properties as a function of latitude, with the exception of a decrease in its optical depth at the highest observed northern latitudes. While we do tentatively observe a meridional trend in Neptune’s 1–3 bar methane abundance, our data suggest that southern midlatitudes are as dry as Neptune’s south pole. Other observations have hinted that this relatively simple circulation model does not adequately describe Neptune’s circulation: for example, 1-cm maps of Neptune show an enhancement in the radio brightness temperature at mid-southern latitudes, indicative of relatively dry air in a region where moist, rising air is predicted (de Pater et al., , 2014). Additionally, clouds are often observed near the south pole and equator – regions of apparent subsidence (e.g., Luszcz-Cook et al., 2010a, ). This is similar to the situation for Uranus, in which models have difficulty reproducing all of the detected aerosol layers (Sromovsky et al., , 2014). In summary, a coherent dynamical picture for either Neptune or Uranus remains elusive; our observations reinforce this point. The vertical and latitudinal distribution of discrete, NIR-bright storms provides additional constraints on circulation models; in a followup paper, we will analyze these features in our OSIRIS data, and combine our findings for both cloud-free and cloudy regions to place improved constraints on Neptune’s dynamics.

Acknowledgements

The data presented were obtained at the W.M. Keck Observatory, which is operated as a scientific partnership among the California Institute of Technology, the University of California and the National Aeronautics and Space Administration. The Observatory was made possible by the generous financial support of the W.M. Keck Foundation. The authors wish to recognize and acknowledge the very significant cultural role and reverence that the summit of Mauna Kea has always had within the indigenous Hawaiian community. We are most fortunate to conduct observations from this mountain. This work was supported by NSF Grant AST-0908575 to UC Berkeley; SHLC was also supported by NASA Headquarters under the NASA Earth and Space Science Fellowship program - Grant NNX10AT17H; and a Kalbfleisch Postdoctoral Fellowship at the American Museum of Natural History. SHLC would like to thank A. Veicht and D. Zurek for computing assistance and support. The authors would like to thank L. Fletcher for providing his temperature profiles, and two anonymous reviewers for their insightful and valuable comments, which improved this paper.

Appendix A Atmospheric retrieval code

In this appendix, we provide additional details of our atmospheric retrieval code, introduced in Section 3.

A.1 Inputs to the forward model

As described in Section 3.1, we solve the radiative transfer equation for a model atmosphere using either a two-stream approximation or a Python implementation of the discrete ordinate method for radiative transfer (pyDISORT). The fundamental input parameters to the forward model include a thermal structure (temperature-pressure profile); the atmospheric composition as a function of depth; gas opacities as a function of temperature and pressure; and a description of the aerosols.

A.1.1 Thermal profile

A range of thermal profiles have been derived for Neptune; many of these are summarized in Figs. 5 and 6 of Luszcz-Cook and de Pater, (2013). For this analysis, we assume the Lindal, (1992) value of 71.5 K for the temperature at 1 bar. Below 1 bar, we extrapolate adiabatically to higher pressures assuming a dry adiabat with a deep methane abundance of 4%; H2O and H2S enrichments equivalent to 50 times solar O and S; and a solar enrichment in N as NH3 (Romani et al., , 1989; de Pater et al., , 1991; Deboer and Steffes, , 1996). We find that adjustments to the deep abundance of CH4 (Appendix A.1.2) have a negligible effect on the adiabatic thermal profile at the pressures considered here. Above the 1-bar level, we adopt the thermal profile described by Fletcher et al., (2014): the tropospheric profile is consistent with the model of Moses et al., (2005). This default thermal profile is shown in Fig. 4.

Latitudinal temperature variations near the tropopause have been observed by Conrath et al., (1998); Orton et al., (2007); Fletcher et al., (2014). These authors find a clear temperature minimum around 40–60∘S and maxima around the equator and at the south pole. In order to investigate the possible influence of thermal variations of this type, the atmospheric model accepts offsets to the tropopause temperature, which are applied as a sinusoidal perturbation to the thermal profile, such that the profile matches the nominal profile at 0.01 and 1.0 bar, but is higher or lower at the tropopause by the specified tropopause temperature offset (Fig. 4). Such variations are considered in Sections 5.3 and 6.3.

In Neptune’s stratosphere, disk-averaged thermal profiles differ from one another by as much as $\sim 20$K (see Luszcz-Cook and de Pater, , 2013), and stratospheric temperatures also may vary spatially. For example, Bezard et al., (1991) measured a broad local maximum in stratospheric emission at low latitudes. This was not seen in the 2005 ground-based data of Hammel and Lockwood, (2007); rather, these authors detected evidence for a 4–5 K temperature enhancement near the south pole. Orton et al., (2007) found evidence of a localized stratospheric warm region near 70∘ S, and Hammel et al., (2006) noted short-term variations in the disk-integrated ethane emission which could indicate discrete features. In contrast, Fletcher et al., (2014) found that, with the exception of the south pole, the stratosphere is latitudinally uniform, and temporal variability, if present, is limited to $<\pm$5 K.

In this analysis, we adopt the mean stratospheric temperature structure from fits to Keck/LWS (2003) data (Fletcher et al., , 2014). The model will accept an offset to the stratospheric temperatures, which is specified in terms of two values: a stratosphere temperature offset (K), and a stratosphere offset pressure (bar, set to $10^{-3}$ in this analysis). The full temperature offset is applied at pressures lower than the offset pressure. To ensure a smooth output thermal profile, smaller offsets, linear with log(pressure) are applied at pressures between the offset pressure and a pressure of 100 mbar. See Fig. 4 for examples. The effects of varying the stratospheric temperature on our fits are briefly investigated in Section 5.3.

A.1.2 Composition

Neptune’s upper atmosphere is composed primarily of H2 and He. We maintain a He/H2 ratio of 0.15/0.847 and an N2/H2 ratio of 0.003/0.847 by number throughout our model atmosphere. These values are consistent with Conrath et al., (1993) reanalysis of Voyager occultation measurements (Conrath et al., , 1991), which placed constraints on the atmospheric mean molecular weight, and the Burgdorf et al., (2003) analysis of spectra from the Infrared Space Observatory Long-Wavelength Spectrometer (ISO-LWS). We assume ‘intermediate’ hydrogen in the troposphere, in which the ortho and para states of hydrogen are in equilibrium at the local temperature, but the specific heat is near that of ‘normal’ hydrogen. This situation is described in Massie and Hunten, (1982). Fast vertical mixing from the interior could bring the ortho/para ratio closer to the 3:1 ratio expected for normal hydrogen, but previous results indicate that that the ortho/para ratio is in fact close to that of equilibrium hydrogen (Orton et al., , 1986; Baines et al., 1995b, ; Burgdorf et al., , 2003). We do not consider the effects of varying the ortho/para ratio. However, previous studies have found evidence for spatial variations in the ortho/para ratio: Conrath et al., (1998), in their analysis of Voyager 2/IRIS data, found sub equilibrium para ratios from the equator down to 50∘S, and higher para fractions relative to equilibrium from 50-90∘S and in the northern hemisphere. Fletcher et al., (2014) performed a reanalysis of the Voyager data using updated absorption spectra; they found a more symmetric para-H2 distribution, with super-equilibrium para ratios from 20∘S to 20∘N and near the poles; and sub-equilibrium para ratios at midlatitudes.

Methane absorption has a dominant effect on Neptune’s NIR spectrum. A number of estimates have been made of Neptune’s CH4 abundance (mole fraction) at various atmospheric depths: based on Voyager occultation measurements, Lindal, (1992) determine a methane mole fraction of $0.02\pm 0.02$ below the methane condensation level at the occultation ingress latitude of $\sim 60^{\circ}$N. The Baines et al., 1995b analysis of hydrogen quadrupole and methane lines yields a disk-averaged tropospheric CH4 mole fraction of $0.022^{+0.005}_{-0.006}$, in good agreement with the Lindal, (1992) result. A CH4 mole fraction of 0.022 has since been adopted by many studies of Neptune’s upper atmosphere (e.g. Roe et al., 2001b, ; Sromovsky et al., 2001c, ; Gibbard et al., , 2003; Luszcz-Cook et al., 2010b, ; Irwin et al., , 2011); this value implies an enrichment factor of $\sim 50$ over the protosolar C/H ratio (Asplund et al., , 2009). More recently, Karkoschka and Tomasko, (2011) find a significantly higher best-fit value of $0.04\pm 0.01$ for the CH4 mole fraction in the troposphere at low latitudes; they also find that at high southern latitudes, the methane abundance is depressed between 1.2 and 3.3 bar. A similar shallow high-latitude depletion has been observed on Uranus (Karkoschka and Tomasko, , 2009; Sromovsky et al., 2011b, ; Sromovsky et al., , 2014; de Kleer et al., , 2015).

At the cooler temperatures of the upper troposphere, methane will condense out of the atmosphere; thermochemical models predict that the cloud base will be at the 1–2 bar level. This is consistent with the Lindal, (1992) interpretation of the Voyager occultation data, which placed the cloud base at 1.9 bar for their nominal model. The relative humidity above the methane condensation layer in the Lindal, (1992) nominal model is 20%; Karkoschka and Tomasko, (2011) found that humidities of 40-100% are required for good fits to their data.

The stratospheric CH4 abundance should be set by the saturation value at the temperature minimum of the tropopause– $\sim 2\times 10^{-4}$ for our thermochemical equilibrium calculations. However, observations indicate the stratospheric CH4 mole fraction is much higher than this: Orton et al., (1987, 1990, 1992); Orton and Yanamandra-Fisher, (2005) found values of 0.75–1.5$\times 10^{-3}$ from analyses of Neptune’s thermal spectrum; Yelle et al., (1993) found a similar range (0.6–5.0 $\times 10^{-3}$ in the lower stratosphere). Baines and Hammel, (1994) favored a somewhat lower mole fraction of 3.5$\times 10^{-4}$, whereas the more recent studies of Fletcher et al., (2010) and Lellouch et al., (2010) derived abundances of $(9\pm 3)\times 10^{-4}$ and $(1.5\pm 0.2)\times 10^{-3}$, respectively, for the middle of the stratosphere.

Like in the troposphere, the stratospheric methane abundance may vary with latitude: Conrath et al., (1998); Orton et al., (2007); Fletcher et al., (2014) showed that the tropopause temperature varies, with warmer temperatures near the south pole; this could allow methane to leak into the stratosphere at high southern latitudes, and then be redistributed across the planet (Orton et al., , 2007). However, Greathouse et al., (2011) suggest that stratospheric CH4 abundances may actually peak near the equator; consistent with the distribution found by Karkoschka and Tomasko, (2011) for the troposphere.

We parameterize the vertical methane profile in the following way: we define the methane mole fraction below the condensation level, mCH4t. At the methane condensation pressure, the methane mole fraction drops, and we adopt a constant relative humidity $RH_{ac}$. We allow the methane abundance to increase again above the tropopause, holding the relative humidity constant, until the stratospheric methane mole fraction mCH4s, is reached. Our initial CH4 profile assumes mCH4t=0.04, $RH_{ac}$ = 1.0, and mCH4s=0.00035. These choices are reevaluated in Section 5.2.

The model also includes an additional parameter, $P_{dep}$, the pressure depth down to which the methane abundance is smoothly depleted. $P_{dep}$ is defined using the proportionally descended gas model described in Sromovsky et al., 2011b ; Sromovsky et al., (2014) and de Kleer et al., (2015), with a shape parameter $vx=2$. We note that this depletion model was developed for Uranus and may not accurately represent the physical situation on Neptune. Methane depletion is considered in Section 6.2. Examples of the parameterized CH4 profile and CH4 depletion are shown in Fig. 4.

A.1.3 Gas opacity

The gas opacity in the NIR is dominated by H2 collision-induced absorption (CIA) and CH4 opacity. For CIA, we use the coefficients for hydrogen, helium and methane from Borysow et al., (1985, 1988); Borysow, (1991, 1992, 1993). For CH4 we use tabulated k-coefficients based on the improved line list and recommendations for outer planet NIR spectra described in Sromovsky et al., (2012): shortwards of 1.80 $\mu$m we adopt the line shape between that of Hartmann et al., (2002) and de Bergh et al., (2012) as described by these authors. Between 1.80 and 2.08 $\mu$m, we use the Karkoschka and Tomasko, (2010) k-coefficients, and at longer wavelengths we adopt the Sromovsky et al., (2012) line list with the Hartmann et al., (2002) line shape. All tabulated coefficients have a 10 cm*-1* resolution and 5 cm*-1* spacing.

A.1.4 Aerosols

A wide range of haze and cloud structures have been used in relation to Neptune. In our model, we aim to parameterize the aerosols in a simple, straightforward way that is nonetheless capable of representing realistic aerosol structures and reproducing previous models. Our model accepts any number of aerosol layers, which can be placed at any depth within the model atmosphere and can overlap one another. The parameters that define each aerosol layer are summarized in Table LABEL:table:model. They are:

• •

Pmax: Maximum (bottom) pressure of the aerosol layer.

• •

Pmin: Minimum (top) pressure of the aerosol layer. Unless indicated, Pmin is set to the top of the model atmosphere at $10^{-5}$ bar.

• •

hfrac: Fractional scale height of the aerosol layer. The number density of haze particles therefore varies across a layer of thickness $z$ (cm) as:

[TABLE]

where $N_{0}$ is the number density at the bottom of the atmospheric layer being considered, and $H$ is the pressure scale height of the gas. A very thin aerosol layer (“cloud”) might have $h_{frac}<0.1$, while a “haze” might have $h_{frac}\sim 1$.

• •

$\tau$: Total optical depth for the haze layer at 1.6 $\mu$m.

• •

$\sigma$: mean extinction cross section of the aerosols. For this analysis, we consider ensembles of particles with distributions of the form:

[TABLE]

where $n(r)$ is the number density of particles of radius $r$ and $r_{p}$ is the peak in the particle distribution (Hansen and Pollack, , 1970)131313This particle size distribution can be written in terms of an equivalent gamma distribution. For $r_{p}=1.0$, the equivalent gamma distribution has a mean of $7/6$ and a variance of $7/36$; for $r_{p}=0.1$, the equivalent gamma distribution has a mean of $0.7/6$ and variance of $0.07/36$.. A range of values have been adopted for the characteristic size of Neptune’s aerosols in the past; these values generally range from 0.1 to 2.5 $\mu$m (e.g. Pryor et al., , 1992; Baines and Hammel, , 1994; Moses et al., , 1995; Irwin et al., , 2014). We test three characteristic particle sizes in this study: $r_{p}=0.1,0.5$ and 1.0 $\mu$m (Fig. 24). The mean extinction cross section is calculated from the ensemble of particles using Mie theory; since the compositions of the condensates in Neptune’s atmosphere are not well known, and the refractive indices of the proposed condensates are also not well known, we use a ‘typical’ hydrocarbon refractive index of $1.4+0i$ for all particles.

• •

g: the asymmetry factor, which describes the shape of the phase function. Values of -1 (strongly backscattering) to 1 (strongly forward scattering) are accepted. In this study, $g$ is assumed to be wavelength-independent, and is generally a free parameter in the retrievals (or fixed according to the results of earlier retrievals). Figure 24 illustrates the wavelength-dependent value of $g$ derived by Mie theory for our three test particle size distributions. In cases where the two stream approximation is used, g fixes the fraction of forward and backscattered light according to the quadrature method as described in Meador and Weaver, (1980). For the pyDISORT algorithm, we interpret g as the asymmetry factor defining the Henyey-Greenstein phase function; the first four moments of the Legendre polynomial expansion are used in the radiative transfer calculation.

• •

$\omega_{0}$: single scattering albedo. Values between 0 (all extinction due to absorption) and 1 (all extinction due to scattering) are permitted. Mie theory finds $\omega_{0}=1.0$ for the particle ensembles considered here; however, this value is inconsistent with limb darkening data (see Irwin et al., (2011) and Section 4). Therefore, we allow $\omega_{0}$ to be a free parameter (constant with wavelength) in the retrievals.

Once the properties of each aerosol layer are specified, the total optical depth and the weighted mean asymmetry parameter and asymmetry factor are calculated for each depth in the atmospheric model.

A.2 Retrieval

To estimate model parameters and uncertainties, we pair our forward model with an affine-invariant Markov chain Monte Carlo (MCMC) ensemble sampler called emcee (Foreman-Mackey et al., , 2013). Emcee generates a sampling approximation to the posterior probability function $p(\boldsymbol{x}|\boldsymbol{d})$, where $\boldsymbol{x}$ is the vector of model parameters and $\boldsymbol{d}$ are the data, by constructing an ensemble of chain(s) sampled from the desired probability distribution by random walk. As noted in Section 3.2, this algorithm possesses several advantages over the simpler and more common Metropolis-Hastings MCMC sampling. The method involves stepping an ensemble of K “walkers” according to a proposal distribution that is based on the current position of the other K-1 walkers; a more complete description of the algorithm can be found in Foreman-Mackey et al., (2013) and references within; here we explain how emcee is employed in our retrievals. A similar implementation is described in de Kleer et al., (2015).

A.2.1 Implementation of the MCMC algorithm

Implementation of emcee requires specification of the posterior probability function, which is, up to a constant:

[TABLE]

where $p(\boldsymbol{x})$ is the prior probability distribution function and $p(\boldsymbol{d}|\boldsymbol{x})$ is the probability of drawing the observed values given the model parameters, or likelihood function. Any preexisting knowledge about the model parameters is encoded in $p(\boldsymbol{x})$; for the parameters in this analysis, we use either uniform or log-uniform priors, within an acceptable parameter range. For example, the prior on $\omega_{0}$ is uniform between 0 and 1; the prior on $\ln{\tau}$ is defined to be uniform between $-7$ and 5.

We assume a gaussian likelihood function:

[TABLE]

where $d_{i}$ is the reflectivity of the $i$th data point, $d_{\boldsymbol{x},i}$ is the value returned by the model for point $i$ given parameters $\boldsymbol{x}$, and $\sigma_{i}$ is the uncertainty for point $i$. In terms of log-likelihood (which is the quantity accepted by emcee), this becomes:

[TABLE]

As noted in Section 3.2, we define the uncertainty as:

[TABLE]

where $\sigma_{n}$ is the random noise component, $\sigma_{p}$ is the relative photometry uncertainty, and $\sigma_{t}$ is a model “tolerance”, representing any unknown uncertainties that prevent the model from matching the data. The primary effect of including $\sigma_{t}$ is to cause the parameter uncertainties to be more realistic.

Each atmospheric retrieval proceeds in the following way: we select a dataset of interest $\boldsymbol{d}$. We then set up the model atmosphere as described in Section 3.1, and identify the set of parameters $\boldsymbol{x}$ that may vary in the model, specifying whether each parameter will be considered in linear or log space. We also specify an allowed range and initial guess value for each parameter. The starting position of each of $>100$ walkers (30–75 walkers per free parameter) in the emcee ensemble is initialized to have a small random offset from that initial guess. We generally run emcee for 500–750 steps, and then check whether the “burn-in” phase – in which the walkers lose memory of the initial conditions and begin to reasonably sample the posterior probability distribution – is complete, 2) the chains have run for at least 10 autocorrelation times after burn-in, and 3) acceptance fraction is reasonably high (0.2–0.5) for all walkers. If these conditions are not yet met- for example, if the walkers still appear to be in the transient “burn-in” phase, we continue to run the chains. In some cases, a low acceptance fraction is observed for one or more walkers, which have gotten stuck in a local maximum in the posterior probability. In these cases, we either continue to run emcee until all walkers appear to have found the global maximum, or we initialize a new ensemble of walkers and begin a new run.

Once a run is complete, we produce a final set of samples by discarding the initial steps (burn-in) and combining the remaining samples from all walkers (see Fig. 25 for an example). Our results are represented by plots of the one- and two- dimensional projections of the posterior probability distributions (e.g. Fig. 5), and by the 16th, 50th and 84th percentiles from the marginalized parameter distributions.

A.2.2 Model comparison

While MCMC algorithms are a powerful tool for estimating model parameters and uncertainties, they do not provide a direct metric for evaluating the efficacy of a model, or for comparing different models (different sets of assumptions and free parameters). Large values of $\sigma_{t}$ imply the differences between the data and model are not well captured by the known uncertainties ($\sigma_{n}$ and $\sigma_{p}$), and therefore serve as one indicator that a model is not a good match to the data. A more robust way to evaluate the relative success of a model fit is to calculate the Deviance Information Criterion ($DIC$) (Spiegelhalter et al., , 2002). The $DIC$ includes a term for goodness-of-fit with a penalty term for model complexity. The $DIC$ is based on the deviance, $D(\boldsymbol{x})$, defined as the difference in log-likelihoods between the model fit and a perfect fit to the data. The deviance is given by:

[TABLE]

where $f(\boldsymbol{d})$ is a function of the data alone, and therefore cancels when considering the difference in deviance between models for the same dataset. The log-likelihood $\ln{p(\boldsymbol{d}|\boldsymbol{x})}$ is defined in Eq. 5. $DIC$ is then defined as

[TABLE]

$\overline{D(\boldsymbol{x})}$ is the mean deviance, which is easily calculated from the log-likelihoods found by emcee; $p_{D}$ is an effective number of parameters. Spiegelhalter et al., (2002) defined $p_{D}$ in the following way:

[TABLE]

where $\overline{\boldsymbol{x}}$ is the posterior mean of the parameters (the 50th percentiles from the marginalized parameter distributions). We find that, for cases where the posterior probability distribution is not multivariate normal, this formulation of $p_{D}$ does not produce an accurate estimate of the degrees of freedom of a model (and is often negative). Therefore, we adopt the alternative formulation of $p_{D}$ following Gelman et al., (2003):

[TABLE]

which we find to be substantially more robust.

The goodness of fit, as estimated by $\overline{D(\boldsymbol{x}})$, will improve for better models, but also with the addition of more parameters; $p_{D}$ compensates for this by favoring models with fewer parameters. Therefore, a smaller value of the $DIC$ (which may be a positive or negative number) indicates a better model. We consider a difference in the $DIC$ ($\Delta DIC$) of 10 or more to indicate a preference for the model with the lower $DIC$.

We note that the $DIC$ for a single model cannot be interpreted in terms of the success or failure of that model in reproducing the data; the value of the $DIC$ is in comparing two different models. Furthermore, the models being compared must be for the same dataset: values of the $DIC$ cannot be compared for models of different datasets (in this case, different spectra or sets of spectra) since $f(\boldsymbol{d})$ in the equation for $D(\boldsymbol{x})$ varies between datasets, and is not known.

Appendix B Comparison and analysis with pyDISORT

Ideally, the errors introduced by our radiative transfer method will be less than, or at worst comparable to, other sources of uncertainty in the retrievals, including errors in the gas opacity coefficients, uncertainties and approximations in the model atmosphere properties, and data uncertainties. Using this same software and comparable data for Uranus, de Kleer et al., (2015) found that the differences between two stream and the more accurate pyDISORT algorithm varied with wavelength, from $<5$% to $\sim$10%. Here we perform a similar comparison of the two stream and pyDISORT models using the best-fit 2L_hazeA model parameters.

Since two stream and pyDISORT treat the asymmetry factor $g$ differently, we do not expect $g_{\alpha}$ and $g_{\beta}$ to be the same for the two algorithms; therefore, we perform a simple, unweighted $\chi^{2}$ fit to calculate $g_{\alpha,D}$ and $g_{\beta,D}$: the pyDISORT asymmetry factors that provide the best match between the $\mu=0.8$ models from the two algorithms. We find that $g_{\alpha,D}=0.5$ and $g_{\beta,D}=-0.3$ (vs. $g_{\alpha}=0.4$ and $g_{\beta}=-0.6$ for two stream). The absolute difference between the two models, $|\delta m|\equiv|I/F_{DISORT}-I/F_{twostream}|$, is shown in Fig. 3. It can be seen in this figure that the errors introduced by the two stream approximation are wavelength dependent, and comparable in magnitude to the other sources of known error in our retrievals, $\sigma_{n}$ and $\sigma_{p}$. In Hbb, where the model differences are greatest, $|\overline{\delta m}|=1.7e-4$, compared to $\overline{\sigma_{n}}=1.5e-4$, for the $\mu=0.8$ spectrum.

More generally, with the exception of the largest emission angles, for which we were not able to match the pyDISORT and two stream models, we find that the two stream solution at a single viewing geometry is a good approximation to the more accurate pyDISORT model, allowing for different values of $g$ for the two algorithms. As emission angle varies, so do the values of $g_{\alpha,D}$ and $g_{\beta,D}$ that generate the best agreement between the two models. In other words, a single two stream phase function does not match the pyDISORT Henyey-Greenstein phase function for all viewing geometries. This is perhaps not surprising, since the two algorithms treat the phase function differently. To illustrate the divergence of the two stream and pyDISORT solutions, we plot (Fig. 26) the 2L_hazeA two stream and pyDISORT models for all six $\mu$ values in the binned dataset, with $g_{\alpha}=0.4$ and $g_{\beta}=-0.6$ for two stream and $g_{\alpha,D}=0.5$ and $g_{\beta,D}=-0.3$ for pyDISORT. We observe that the match between the two models diverges towards high emission angles, such that by $\mu=0.3$ the fractional difference between models is of order unity. This divergence appears to be much more extreme for the parameter values in this particular (2L_hazeA) model atmosphere than for models containing scatterers with more forward scattering phase functions: when we modify the two stream models such that $g=0.7$ in both aerosol layers, we find a good match to the $\mu=0.8$ model spectra for $g_{\alpha,D}=0.7$ and $g_{\beta,D}=0.6$. Furthermore, the resultant two stream and pyDISORT models are in very good agreement at all values of $\mu$ (Fig. 27).

In light of these results, we conclude that while two stream is sufficient for modeling a single spectrum, particularly for locations far from the limb, a more accurate treatment of the radiative transfer is warranted when multiple viewing geometries are considered simultaneously. We emphasize that the two stream asymmetry factor $g$ does not apply for more than a single viewing geometry or to other radiative transfer algorithms, such as pyDISORT.

Appendix C Supplementary figures