ARTICLE IN PRESS

JID: YICAR[m5G;January 2, 2017;11:46]

Icarus 0 0 0 (2017) 1–18

Contents lists available at ScienceDirect

Icarusjournal homepage: www.elsevier.com/locate/icarus

A post-new horizons global climate model of Pluto including the N2 ,CH4 and CO cycles
F. Forget a,∗, T. Bertrand a, M. Vangvichith a, J. Leconte a,1, E. Millour a, E. Lellouch baLaboratoire de Météorologie Dynamique (LMD/IPSL), Sorbonne Universités, UPMC Univ Paris 06, PSL Research University, Ecole Normale Supérieure,
Université Paris-Saclay, Ecole Polytechnique, CNRS, Paris, Franceb
Laboratoire d’Etudes Spatiales et d’Instrumentation en Astrophysique (LESIA), Observatoire de Paris, CNRS, UPMC, Université Paris Diderot, F-92195
Meudon, Francea r t i c l e

i n f oArticle history:
Received 29 March 2016
Revised 3 October 2016
Accepted 29 November 2016
Available online xxx
Keywords:
Pluto
Pluto, atmosphere
Atmospheres, composition
Atmospheres, dynamics

a b s t r a c t
We have built a new 3D Global Climate Model (GCM) to simulate Pluto as observed by New Horizonsin 2015. All key processes are parametrized on the basis of theoretical equations, including atmospheric
dynamics and transport, turbulence, radiative transfer, molecular conduction, as well as phases changes
for N2 , CH2 and CO. Pluto’s climate and ice cycles are found to be very sensitive to model parameters
and initial states. Nevertheless, a reference simulation is designed by running a fast, reduced version of
the GCM with simpliﬁed atmospheric transport for 40,0 0 0 Earth years to initialize the surface ice distribution and sub-surface temperatures, from which a 28-Earth-year full GCM simulation is performed.
Assuming a topographic depression in a Sputnik-planum (SP)-like crater on the anti-Charon hemisphere,a realistic Pluto is obtained, with most N2 and CO ices accumulated in the crater, methane frost covering both hemispheres except for the equatorial regions, and a surface pressure near 1.1 Pa in 2015 withan increase between 1988 and 2015, as reported from stellar occultations. Temperature proﬁles are inqualitative agreement with the observations. In particular, a cold atmospheric layer is obtained in thelowest kilometers above Sputnik Planum, as observed by New Horizons’s REX experiment. It is shown toresult from the combined effect of the topographic depression and N2 daytime sublimation. In the reference simulation with surface N2 ice exclusively present in Sputnik Planum, the global circulation is only
forced by radiative heating gradients and remains relatively weak. Surface winds are locally induced bytopography slopes and by N2 condensation and sublimation around Sputnik Planum. However, the circulation can be more intense depending on the exact distribution of surface N2 frost. This is illustrated in
an alternative simulation with N2 condensing in the South Polar regions and N2 frost covering latitudes
between 35°N and 48°N. A global condensation ﬂow is then created, inducing strong surface winds everywhere, a prograde jet in the southern high latitudes, and an equatorial superrotation likely forced bybarotropic instabilities in the southern jet. Using realistic parameters, the GCM predict atmospheric concentrations of CO and CH4 in good agreement with the observations. N2 and CO do not condense in theatmosphere, but CH4 ice clouds can form during daytime at low altitude near the regions covered by N2
ice (assuming that nucleation is eﬃcient enough). This global climate model can be used to study many
aspects of the Pluto environment. For instance, organic hazes are included in the GCM and analysed in acompanion paper (Bertrand and Forget, Icarus, this issue).
© 2016 Elsevier Inc. All rights reserved.

1. Introduction
Only six terrestrial bodies in our solar system (Venus, Earth,Mars, Titan, Triton, Pluto) possess an atmosphere thick enough to
be governed by the same equations of meteorology as on Earth,or able to support clouds or hazes. Among them, Pluto presents∗

Corresponding author.
E-mail address: forget@lmd.jussieu.fr (F. Forget).
1
Laboratoire d’Astrophysique de Bordeaux, Univ. Bordeaux, F-33270 Floirac,
France.

a unique case, with an atmosphere signiﬁcantly warmer than theunderlying surface, long radiative timescales, and a circulation
dominated by condensation/sublimation process of the main atmospheric component. Studying this exotic case can provide new
insight into the physics of terrestrial atmosphere.
The observations made by the New Horizons spacecraft have revealed the nature of the surface of Pluto and have provided unprecedented constraints on the state of the atmosphere in 2015
(Stern et al., 2015; Gladstone et al., 2016; Grundy et al., 2016;
Moore et al., 2016). Within that context, it is now interesting to

http://dx.doi.org/10.1016/j.icarus.2016.11.038
0019-1035/© 2016 Elsevier Inc. All rights reserved.

Please cite this article as: F. Forget et al., A post-new horizons global climate model of Pluto including the N2 , CH4 and CO cycles, Icarus(2017), http://dx.doi.org/10.1016/j.icarus.2016.11.038

JID: YICAR

ARTICLE IN PRESS

2[m5G;January 2, 2017;11:46]

F. Forget et al. / Icarus 000 (2017) 1–18

test our ability to create a 3D numerical simulator of the Pluto
climate system, analogous to the climate models already used on
the Earth as well as on Mars, Venus and Titan. Conversely, the output of such a Global Climate Model is useful to interpret the available atmospheric measurements, and can even shed light on somegeological observations.
1.1. Pluto’s ices and atmosphere observations
The presence of a signiﬁcant atmosphere on Pluto was demonstrated in 1988 by observing a stellar occultation by Pluto(Hubbard et al., 1988; Elliot et al., 1989). This atmosphere was predicted to be mainly composed of molecular nitrogen in vapourpressure equilibrium with N2 ice deposits observed on the surface.
In 2015, New Horizons determined a pressure of about 10 μbar, i.e.
1 Pa (Hinson et al., 2015b; Gladstone et al., 2016). However, a series of stellar occultations conducted since 1988 have shown thatthe pressure at a speciﬁc reference level (e.g. 1275 km from Pluto’s
center) has increased by a factor of three during that period (Elliot
et al., 20 03a; 20 03b; Olkin et al., 2015), suggesting a similar rise of
the surface pressure.
Prior to the New Horizons ﬂyby, spectroscopic observations had
demonstrated that, in addition to N2 ice, Pluto’s surface is coveredby patches of CH4 and CO ices (Owen et al., 1993; Douté et al.,
1999; Grundy et al., 2013). New Horizons was able to map theseices and revealed that the main reservoir was a thick ice cap informally named Sputnik Planum, with thinner N2 frost covering the
mid-northen latitudes and CH4 frost possibly everywhere except
in the equatorial dark regions (Grundy et al., 2016).
Accordingly, CH4 and CO gas were observed from Earth to bepresent in present-day Pluto’s atmosphere (Young et al., 1997; Lellouch et al., 2009; 2011a; 2015; 2016), with, during the 2008–2012
period, volume mixing ratios near 0.05% for CO and 0.5% for CH4 .
(Lellouch et al., 2011a; 2015). CH4 has also been oberved by New
Horizons Alice spectrograph. in 2015 and estimated to range between 0.6 and 0.84% in the lower atmosphere (Gladstone et al.,2016). Using the hydrostatic equation, atmospheric temperature
proﬁles have been derived from vertical density proﬁles retrieved
from Earth-based stellar occultations (Elliot et al., 1989; 2003b;
20 07; Person et al., 20 08; Young et al., 20 08; Dias-Oliveira et al.,2015; Sicardy et al., 2016), and from radio and solar occultations
performed by New Horizons (Hinson et al., 2015b; Gladstone et al.,2016). The latest stellar occultations and New Horizon’s data consistently show that the temperature proﬁle is characterized by asteep temperature gradient in the lower atmosphere, with temperature increasing from surface values (38–55 K) at 0 km to about
110 K at 20 km. This has been interpreted as resulting from the absorption of near-infrared solar radiation by gaseous methane (Yelleand Lunine, 1989; Strobel et al., 1996). Above 30 km, the temperature appears to decrease with altitude to reach about 70–80 K
around 200 km (Dias-Oliveira et al., 2015; Gurwell et al., 2015;
Gladstone et al., 2016). Such a structure requires infrared-coolingspecies acting only at a speciﬁc altitude range. C2 H2 and HCN –respectively detected by New Horizons (Gladstone et al., 2016) andfrom the ground (Lellouch et al., 2016) – have been proposed, but
the details of exactly how Plutos upper atmosphere is being cooled
remains poorly understood (Gladstone et al., 2016). Finally, stellar occultation observations suggest that the temperature proﬁles
are affected by oscillation that can be related to gravity waves orthermal tides (Elliot et al., 2003b; Person et al., 2008; Toigo et al.,2010).
1.2. 3D modelling of the pluto surface-atmosphere system
To improve our understanding of the complex Pluto surfaceatmosphere system, we have built a new Global Climate Model(GCM) including a full simulation of the nitrogen, methane, andcarbon monoxide cycles. This GCM computes the temporal evolution of the variables which control the meteorology and the climate of the planet in different points of a 3D grid covering the entire atmosphere. On the Earth, GCMs have been applied to weather
forecasting and climate change projections. Because these models
are almost entirely built on physical equations (rather than empirical parameters), several teams around the world have been able
to succesfully adapt them to the other terrestrial planets or satellites that have a solid surface and a thick enough atmosphere. The
Pluto GCM presented in this paper is derived from the LMD Global
Climate Model of planet Mars (Forget et al., 1999) which has beenused for numerous applications including simulating CO2 ice caps
analogous to Pluto’s N2 ice caps (Forget et al., 1998), the thermosphere (González-Galindo et al., 2009), photochemistry (Lefèvre
et al., 2008) or pal eoclimatology (e.g. Forget et al., 2006). The LMD
GCM has been adapted to Venus (Lebonnois et al., 2010) and Titan(Hourdin et al., 1995; Lebonnois et al., 2012). All these GCMs havebeen able to predict or accurately reproduce the observed thermal
structure and circulation, giving us some conﬁdence in its ability
to predict the characteristics of the Pluto atmosphere in spite of
the scarcity of observations.
For Pluto, after the simpliﬁed General Circulation Model (without phase changes) presented by Zalucha and Michaels (2013) for
Pluto and Triton, a realistic model was developed by Toigo et al.
(2015) a few months before the New Horizons encounter. Thismodel includes a “robust treatment of nitrogen volatile transport”, and initializes the full GCM using a two dimensional surface volatile exchange model and a one-dimensional radiativeconductive-convective model. In this paper we present a newmodel with a different origin and which beneﬁts from the New
Horizons observations. We include an improved N2 condensationsublimation scheme, the full CO and CH4 cycles, and explore theeffect of topography. Nevertheless, we use an analogous strategy
for the initialization.
In Sections 2, we provide a detailed description of the different components of our Pluto Global Climate Model, and in
Section 3 we discuss how the different model parameters werechosen and how the 3D GCM is initialized for our two baseline
simulations. The model results for temperature and winds and for
the CH4 and CO cycles are then presented in Sections 4 and 5, before the conclusion.
2. Model description
2.1. Generalities
As mentioned above, our Pluto Global Climate Model is derived
from the LMD Mars GCM (Forget et al., 1999) , with several newparameterizations. Its core is a hydrodynamical code dedicated tothe temporal and spatial integration of the equations of hydrodynamics, used to compute the large scale atmospheric motions andthe transport. The equations are solved using a ﬁnite difference
scheme on an “Arakawa C” grid (Arakawa and Lamb, 1977). Sucha scheme is equally valid for the Earth, Mars or Pluto. Therefore
the hydrodynamical core has not been modiﬁed for Pluto. Whilethe estimated surface pressure on Pluto (10 μbar or 1 Pa) is much
lower than on the Earth or even on Mars, the atmosphere is thickenough to be modeled with the primitive equations of meteorology
used in the model. In fact, it is generally found that GCMs dynamical cores are valid almost up to the exobase. For instance, on Mars
our dynamical core has been used sucessfully up to the thermosphere at pressures lower than 10−7 Pa (González-Galindo et al.,2009).
In this paper, we present simulations with a horizontal gridof 32 × 24 points to cover the planet, that is a grid-point spac-

Please cite this article as: F. Forget et al., A post-new horizons global climate model of Pluto including the N2 , CH4 and CO cycles, Icarus(2017), http://dx.doi.org/10.1016/j.icarus.2016.11.038

JID: YICAR

ARTICLE IN PRESS[m5G;January 2, 2017;11:46]

F. Forget et al. / Icarus 000 (2017) 1–18

ing of 7.5° latitude by 11.25° longitude. The corresponding spatial
resolution is of about 150 km, which is equal or better to the typical resolution used in planetary GCMs, and which is suﬃcient toresolve possible planetary waves. We also performed simulations
with a doubled resolution (64 × 48) and even an experimental
run with a 360 × 180 grid, and did not ﬁnd any fundamental differences in the results that could change the conclusions of thispaper. Their analysis is out of the scope of this paper and will
be presented in a future article, in which we will take into account a more realistic topography. In the vertical, the model uses
the terrain-following “sigma” coordinate system in ﬁnite difference
form (i.e. each layer is deﬁned by a constant value of the ratiopressure devided by surface pressure). 25 levels are typically used.
In the baseline model, most of the levels are located in the ﬁrst 15
km to obtain a good resolution close to the surface, in the boundary layer. The altitude of the ﬁrst mid-layers are 7 m, 15 m, 25 m,
40 m, 80 m etc.. Above 10 km, the resolution is about one scaleheight, with the upper pressure level equal to 0.007 times surface pressure, i.e. up to 250 km. (Note that in a companion paper
dedicated to the study of atmospheric hazes Bertrand and Forget(2016a), the top of the model is extended to about 600 km to include the altitudes of methane photolysis).
2.2. Radiative transfer
The incident insolation upon each modeled atmospheric column is calculated at each timestep, taking into account the variation of the Pluto-Sun distance throughout its orbit, the seasonal
inclination and the diurnal cycle.
While N2 is the major constituent of the atmosphere of Pluto,
its radiative effects are neglected in the lower atmosphere since N2
is transparent at solar and infrared wavelengths. Nevertheless weaccount for (1) the radiative heating and cooling by CH4 , which canvary in space and time depending of the results of the methane
cycle model (see Section 2.8) (2) cooling by the thermal infrared
rotational lines of CO, which volume mixing ratio is prescribed at
0.05 % everywhere (Lellouch et al., 2011a; 2016) and 3) the effectof other infrared emitting species in altitude.
2.2.1. Radiative transfer through CH4 and CO
For CH4 and CO, we use a correlated k-distribution radiativetransfer model, with 17 spectral bands in the thermal infraredand 23 for solar wavelengths. The bands are designed to well
represent the 1.6, 2.3 and 3.3 μm CH4 vibrational bands in thenear infrared as well as the 7.6 μm CH4 emission band in thethermal infrared. To calculate the k absorption coeﬃcients ineach bands, high resolution line-by-line spectra combining CO
and CH4 were computed from the HITRAN 2012 database using
the open-source “kspectrum” tool. Spectra and k coeﬃcientswere calculated to ﬁll a look up table matrix (from which the kcoeﬃcients are interpolated by the GCM in each spectral band)
comprising 8 temperatures × 7 log-pressure × 7 CH4 volume
mixing ratio grid, with T = {30, 40, 50, 70, 90, 110, 150, 200} K, p ={10−4 , 10−3 , 10−2 , 10−1 , 1, 10, 100} Pa, and [CH4 ] = {10−4 , 10−3 , 5 ×
10−3 , 10−2 , 5 × 10−2 , 10−1 , 5 × 10−1 } kg/kg. We found that no lessthan 33 points were needed for the g-space integration to get
accurate results throughout the matrix space (g is the cumulated
distribution function of the absorption data for each band).
2.2.2. Non local thermal equilibrium processes
In the low-pressure, low temperature Pluto environment, a major diﬃculty (and therefore uncertainty) in the radiative transfer
calculations results from the fact that the methane lines can be far
from Local Thermal Equilibrium (LTE). It is not the case of CO rotational lines which are assumed to remain in LTE for the pressurelevels that we model in this paper.

3

To account for non-LTE effects for the 7.6 μm CH4 band, we
modify the LTE cooling rates obtained with the correlated kdistribution radiative transfer model as in Strobel et al. (1996).
However, the total CH4 cooling rates we obtain are found to bemuch lower than shown in Strobel et al. (1996), and signiﬁcantly
smaller than the CO cooling rates. This is also found in recent models from the same authors (D. Strobel, personnal communication)
The difference is thought to result from the updated spectroscopic
database (HITRAN 2012 vs HITRAN 1986) and the fact that thetemperatures used in Strobel et al. (1996) are larger than here. The
uncertainties on the NLTE calculations for the 7.6 μm CH4 band
have thus a limited effect on our results.
For the near-infrared solar bands, we ﬁrst reproduced the calculations from Strobel et al. (1996) updated by Zalucha et al.
(2011) for each of the 2.3 and 3.3 μm bands. We had no information on the 1.6 μm band. Within that context, and given the overall uncertainty in the NLTE calculations (Boursier et al., 2003), we
authorized some empirical modiﬁcations of the theoretical NLTE
variations with atmospheric density (while keeping the theoretical shape) to adjust the heating rates in order to get temperatures
closer to the thermal structure observed by New Horizons. Therefore, the ability of our GCM to roughly reproduce the observed
mean thermal structure should not be regarded as a success of
our radiative transfer model. In practice, we multiply the total CH4
heating rate provided by the LTE radiative transfer code by a vertically varying non-LTE eﬃciency coeﬃcient ε NLTE .

εNLTE = 0.1 +0.9,
1 + ρ.55 /ρ

(1)with ρ the atmospheric density (kg m−3 ), and ρ .55 the reference density for which εNLTE = 0.55. After tuning, we set ρ.55 =2 × 10−6 kg m−3 .
2.2.3. Additional radiative coolers
As mentioned in the introduction, the presence of radiatively
cooling species at a speciﬁc altitude has been suggested to explain the decrease of temperature above 30 km (Dias-Oliveira et al.,2015; Gladstone et al., 2016). Using the cooling-to-space approximation, we phenomenologically represent this effect with the following cooling rate for pressures below 0.12 Pa:∂T= −5 × 10−11 B(λ0 , T )∂t(2)

with T the atmospheric temperature (K) and B(λ0 , T) the Planck
function (in W m−2 μm−1 sr−1 ) at wavelength λ0 . We use λ0 =14 μm since the main emission bands of the most likely coolingspecies C2 H2 and HCN (Gladstone et al., 2016) are respectively centered at 13.7 and 14.05 μm (we here neglect the rotational bandsof HCN at submillimeter wavelengths). The value −5 × 10−11 waschosen to simulate a moderate cooling yielding temperatures below 90 K in our reference simulation.
2.3. Atmospheric molecular thermal conduction and viscosity
We account for the effect of molecular conduction on temperature and molecular viscosity on winds. Both processes are governed by similar equations. Assuming the plane-parallel approximation, for thermal conduction we get:



∂T
1 ∂
∂T
=
k
∂t
ρ cp ∂ z ∂ z(3)where T is the temperature (K), ρ the density (kg m−3 ) and k
the thermal conduction coeﬃcient (J m−1 s−1 K−1 ), expressed ask = k0 T s , with k0 = 5.63 × 10−5 J m−1 s−1 K−(1+s ) and s = 1.12
(Hubbard et al., 1990).

Please cite this article as: F. Forget et al., A post-new horizons global climate model of Pluto including the N2 , CH4 and CO cycles, Icarus(2017), http://dx.doi.org/10.1016/j.icarus.2016.11.038

ARTICLE IN PRESS

JID: YICAR
4[m5G;January 2, 2017;11:46]

F. Forget et al. / Icarus 000 (2017) 1–18

For molecular viscosity:∂S 1 ∂∂S=μ∂t ρ ∂ z∂z(4)

where S stands for the components of the horizontal wind (m s−1 )and μ is the coeﬃcient of molecular viscosity (kg m−1 s−1 ), that is
related to the thermal conduction coeﬃcient by k = 14 [9c p − 5(c p −R )]μ. Given its similarity, both equations are discretized and solved
using the same implicit numerical schemes.
2.4. Surface temperatures and thermal conduction in the subsurface
Surface temperature evolution Ts is governed by the balance between solar insolation, thermal emisssion in the infrared, latent
heat exchanges (see Section 2.6), sensible heat ﬂux from the atmosphere (usually negligible on Pluto, but taken into account in
the model) and thermal conduction in the soil. On a weakly irradiated body like Pluto, the radiative ﬂuxes are small compared tothe internal heat stored in the ground. In particular, the subsurface
heat stored during one season can play a major role in the controlof the surface temperature at the opposite season.
The heat ﬂux from and to the subsurface is computed using aclassical model of the evolution of the subsurface temperatures T
as a function of time t and depth z. It satisﬁes the following equation:



∂T
∂
∂T
λ
C
=
∂t
∂z ∂z

(5)

where λ is the heat conductivity of the ground, (J s−1 m−1 K−1 )and C the ground volumetric speciﬁc heat (J m−3 K−1 ). This equation is solved using a ﬁnite differences approach and an implicit
Euler scheme. The key parameter which controls the inﬂuence ofthe subsurface heat storage and √conduction on the surface temperature is the thermal Inertia I = λC . In practice, we thus use I as
the key model parameter, assuming a constant value for C = 106 J
m−3 K−1 and making λ vary accordingly.
On Pluto the discretization requires a special attention compared to the Earth or Mars because one need to simultaneously
capture 1) the short period diurnal thermal waves in the nearsurface, low thermal inertia terrain and 2) the much longer seasonal thermal waves which can penetrate deep in the high thermal inertia substrate. In this paper, we assumed a relatively low
diurnal thermal inertia Iday = 50 J s−1/2 m−2 K−1 , slightly higher
than the 20–30 SI range reported by Lellouch et al. (2011b) fromtheir Spitzer data analysis. For the seasonal thermal inertia, we set
Iyear = 800 J s−1/2 m−2 K−1 , which corresponds to a low porosity
ice/rock-like substrate.
The skin depth of a thermal wave of period P (s) is:

I
δP =
C



P

π

(6)

The modeled diurnal and annual skin depths are thus 0.02 m
and 40 m respectively. To represent this accurately, the subsurface is divided into N = 22 discrete layers, with a geometricaly
stretched distribution of layers with higher resolution near the surface and a coarser grid at depth:zk = z1 2k−1

(7)where z1 = 1.414 × 10−4 m is the depth of the ﬁrst layer. The deepest layer depth is thus near 300 m.
2.5. Mixing in the boundary layer
Turbulent mixing and convection are parameterized as in Forgetet al. (1999). In practice, the boundary layer dynamics is accounted

for by a Mellor and Yamada (1982) unstationary 2.5-level closure
scheme, used to compute turbulent mixing coeﬃcients induced bywind shears depending on the temperature proﬁle stability and theevolution of turbulent kinetic energy. It is completed by a “convective adjustment” scheme which rapidly mixes the atmosphere in
the case of unstable temperature proﬁles (rare on Pluto).
Turbulence and convection mix energy (potential temperature),
momentum (wind), and tracers (gases and aerosols). In the surface
layer, the turbulent surface ﬂux is given by

F = ρCdU1 (q1 − q0 ),(8)

where q1 and q0 are the variable values in the ﬁrst atmosphericlayer and at the surface (q0 = 0 for winds), U1 is the horizontal
wind velocity in the ﬁrst layer, and Cd is the drag coeﬃcient. Because of the small depth of the ﬁrst layer z1 , we assume that thewind proﬁle in the ﬁrst meters above the surface is logarithmic
and not inﬂuenced by stability, and simply use


Cd =

κ

2
(9)

ln zz10where κ is the von Karman constant (κ = 0.4) and z0 is the roughness coeﬃcient, set to z0 = 0.01 m everywhere like in the Mars
GCM (Forget et al., 1999).
Turbulent mixing is negligible outside the boundary layer(which is often shallow on Pluto because of the positive lapse rate
above the surface). In our GCM there is no other vertical “eddy diffusion” process. In particular, species are only transported upwards
by the large scale circulation.
2.6. N2 condensation and sublimation
The condensation and sublimation of nitrogen ice must be
carefully computed in the Pluto environment. The amount of energy and the relative mass of the atmosphere involved in phaseschanges at each timestep can be very signiﬁcant. Locally, it not
only changes the surface temperature and pressure, but it also
modiﬁes the structure of the boundary layer by “pumping” the airwhen condensation occurs on the surface, and by releasing large
amount of cold, pure nitrogen (with no horizontal velocity) whenN2 sublimes. Our scheme is adapted from Forget et al. (1998).
However, we found it necessary to make several changes in theequations to better represent the intense condensation and sublimation at the surface of Pluto.
The variation of the condensation temperature Tc with nitrogen partial pressure PN2 is derived from the thermodynamic relations computed by Fray and Schmitt (2009), taking into accountthe transition from the α to the β crystalline form near 35.61 K(corresponding to PN2 = 0.53 Pa):

if PN2 < 0.53 Pa : Tc =if PN2 > 0.53 Pa : Tc =







296.9251
PN2−
ln35.600
1.09LN2
0.508059

296.9251
PN2−
ln
63.147
0.98LN212557.

−1

(10)

−1

(11)

with LN2 = 2.5.105 J kg−1 the latent heat of condensation for nitrogen.
2.6.1. Surface condensation and sublimation
The condensation and sublimation of nitrogen on the ground isprimarily controlled by energy and mass conservation. At a given
timestep, if the surface temperature predicted by radiative andconductive balance T0∗ falls below the condensation temperature atsurface pressure Tc 0, an amount δ m0 (kg m−2 ) of N2 condenses,

Please cite this article as: F. Forget et al., A post-new horizons global climate model of Pluto including the N2 , CH4 and CO cycles, Icarus(2017), http://dx.doi.org/10.1016/j.icarus.2016.11.038

ARTICLE IN PRESS

JID: YICAR[m5G;January 2, 2017;11:46]

F. Forget et al. / Icarus 000 (2017) 1–18

releasing the latent heat required to keep the solid-gas interface at
the condensation temperature (T0 = Tc0 ):δ m0 =

cs(LN2 + c p (T1 − Tc0 ))(Tc0 − T0∗ )(12)

cs is the surface heat capacity (in J m−2 K−1 ), cp the air speciﬁcheat at constant pressure (set to 1040 J kg−1 K−1 for N2 ) and LN2
the latent heat of N2 (2.5 105 J kg−1 ).
The term c p (T1 − Tc0 ) (J kg−1 ) corresponds to the extra heatbrought by the atmosphere (assumed to be at temperature T1 in
the ﬁrst model layer) when cooled to the condensation temperature Tc0 just above the surface. Because Pluto’s lower atmosphereis a warm stratosphere lying just above a surface, we found that
this term can be signiﬁcant. With T1 typically 10 K above Tc0 whenN2 condenses in the topics, it reaches 4% of the latent heat. Conversely, when surface N2 ice predicted temperature T0∗ is above thefrost point Tc0 , N2 sublimes and δ m0 is negative:

δ m0 =

cs(Tca 0 − T0∗ )
LN2

(13)

We set T0 = Tc0 , unless all the local ground ice of mass m0(kg m−2 ) completely sublimes. We then set δ m0 = −m0 and the
new surface temperature is expressed as: T0 = T0∗ − LN2 m0 /cs The
formation or disapearance of nitrogen ice on the substrate is taken
into account in the calculations of the surface albedo and emissivity.
2.6.2. Atmospheric condensation and sublimation
In the atmosphere, things are, in theory, more complex. The
condensation of a gas involves various microphysical processes: supersaturation, nucleation, crystal growth, sedimentation, etc... In
our model, we have kept the detailed Mars GCM CO2 ice sheme
described in the appendix of Forget et al. (1998) and directly
adapted it to N2 ice. Supersaturation is neglected and atmosphericcondensation and sublimation are computed using energy conservation principles as above. We do not simulate the growth andtransport of nitrogen ice particles. Instead, after condensing at agiven model level, we assume that N2 ice falls through the atmospheric layers located below it (where it can sublimate), possibly
down to the ground within a model timestep.
Because the atmosphere is warmer than the surface most of thetime, we have found that atmospheric condensation is a processesof little importance on Pluto as we model it with a 150 km resolution. In reality, ascending motions induced by local slopes or gravity waves could trigger condensation in N2 ice covered regions. We
will explore that in future versions of the model.
2.6.3. Computing mass, momentum and heat vertical ﬂuxes inducedby N2 condensation and sublimation
The condensation and sublimation of nitrogen induce signiﬁcant transport of air (mass, heat, momentum, tracers) through themodel layers as well as to and from the surface. These processes
must be taken into account on Pluto where an atmospheric layerof several tens of meters thick can undergo a phase change at eachtimestep. The numerical resolution of these processes in the “σ ”
vertical coordinates used in the GCM (see Section 2.1) is given in
the appendix.
2.7. Organic hazes
New Horizons revealed the presence of extensive hazes thought
to be primarily composed of organic particles indirectly producedby methane photolysis. Our GCM includes a model of the formation and transport of these particles. This model and its outputs
are described and analyzed in a companion paper by Bertrand andForget (2016a), and not detailed here.

5

2.8. Methane cycle and CH4 ice clouds
The 3D evolution of CH4 on the surface and in gaseous andsolid phase in the atmosphere is simulated taking into account1) the condensation and sublimation at the surface and in theatmosphere (see below), 2) the transport by the general circulation
using the “Van-Leer I” ﬁnite volume scheme from Hourdin andArmengaud (1999), 3) the mixing in the atmosphere by turbulent
diffusion and possibly convection (see Section 2.5), and 4) thegravitational sedimentation of CH4 ice particles (see below).
2.8.1. Surface condensation and sublimation
The mass ﬂuxes of methane to and from the atmosphere arecomputed using Eq. (8), with q0 and q1 the mass mixing ratios (kg/kg) just above the surface and in the middle of the atmospheric ﬁrst layer, respectively. Note that an important consequence of Eq. (8) is that the sublimation rate of methane is proportional to the horizontal wind velocity in the lower atmosphere.
When pure methane is on the surface, q0 is set equal to thesaturation vapour pressure mass mixing ratio of methane q sat CH4 ,calculated as a function ot temperature T (K) and pressure p using
the following expression derived from Fray and Schmitt (2009):

qsatCH4 = 0.117 × 105 e

6.12×105R(1/90.7−1/T ) × MCH4 × 1

Mair

p

(14)

Here MCH4 /Mair is the ratio of molecular masses use to convert
volume mixing ratio into mass mixing ratio and R = 8.314/MCH4 =519 m2 s−2 K−1 the methane gas constant. When both methaneand nitrogen ices are present on the surface and methane is subliming, we assume that methane is diluted in a solid solution
N2 :CH4 with 0.3% of methane (Merlin, 2015). Applying Raoult’slaw, we thus set q0 = 0.005q sat CH4 If the total amount of methaneon the surface is sublimed within a model timestep, the ﬂux tothe atmosphere is limited accordingly. If no methane ice is present
on the surface, then q0 = q1 if q1 < qsat CH4 (no condensation) and
q0 = qsat CH4 if q1 > qsat CH4 (direct condensation onto the surface).
The latent heat released by methane surface condensation andsublimation is taken into account in the surface energy budget assuming a latent heat LCH4 = 5.867 × 105 J kg−1 (Fray and Schmitt,2009).
2.8.2. Atmospheric condensation and CH4 cloud formation
Methane can also condense (and then sublimate) in the atmosphere when the CH4 mixing ratio exceeds the saturation mixing
ratio provided by Eq. (14). We do not know if CH4 can easily nucleate or if large super-saturation is required. Organic particles resulting from the photochemistry in the upper atmosphere probably
offer condensation nuclei suitable for heterogeneous condensation.
In the GCM we assume that all atmospheric methane in excess of
saturation condenses to form ice cloud particles.
The amount of latent heat released by methane condensation or
sublimation is far from being negligible. We ﬁnd that it can locallychange the atmospheric temperature by more than 10 K. Moreover, latent heating actually limits the amount of methane thatcondenses when the atmosphere is supersaturated. If CH4 condensation is calculated without simultaneously taking into accountlatent heat release, or using an explicit numerical scheme, themodel predicts very unrealistic temperatures (e.g. changes larger
than several tens of Kelvins within one timestep), leading to unrealistic condensation rates. In practice, at each model timestep,when the methane mass mixing ratio qCH4 is detected to exceed
saturation (or if methane ice is already present), one must simultaneously calculate the temperature at the end of the timestep, T , asinﬂuenced by the condensation/sublimation and the corresponding
saturation mixing ratio qsat CH4 (T ). For this purpose we numerically

Please cite this article as: F. Forget et al., A post-new horizons global climate model of Pluto including the N2 , CH4 and CO cycles, Icarus(2017), http://dx.doi.org/10.1016/j.icarus.2016.11.038

ARTICLE IN PRESS

JID: YICAR
6

F. Forget et al. / Icarus 000 (2017) 1–18

determine T by solving the following equation:

T  = T + [qCH4 − qsatCH4 (T  )]3. Model initialization and choice of key parameters

LCH4
cp

(15)

The change in CH4 gas and ice mass mixing ratios (kg/kg) arethen given by

δ qCH4 = −δ qice = (qsatCH4 (T  ) − qCH4 ),(16)

unless all the atmospheric CH4 ice is sublimed (and T is adjusted
accordingly).
Once the mass mixing ratio of CH4 ice qice is known, the ice isdistributed to form ice cloud particles around cloud condensation
nuclei (CCN). We assume that the number of cloud condensation
nuclei [CCN] per mass of atmosphere (kg−1 ) is constant throughout
the atmosphere. Assuming that the cloud particle size distribution
is monodispersed in each volume element, the cloud particle radius r is then given by:r=([m5G;January 2, 2017;11:46]

3qice
3+ r[CCN]) 1/3
4π ρice [CCN](17)with ρ ice the CH4 ice density (520 kg m−3 ), and r[CCN] the radiusof the CCN set to 0.2 μm.
Once r is known, the cloud particle sedimentation velocity is
calculated using Stokes law corrected for low pressure by the Cunningham slip-ﬂow correction (Rossow, 1978). The calculated particle radius, r, is also used to estimate the apparent opacity of theclouds. However, we neglected the radiative effect of the clouds inthis paper.
[CCN] is clearly a key parameter which directly controls theproperties of the clouds and their sedimentation. What is thepossible range of [CCN]? On the Earth, the number mixing ratioof activated cloud condensation nuclei in the troposphere ranges
between 106 kg−1 (for low saturation in clean polar air) and
1010 kg−1 (polluted air mass)[Hudson and Yun, 2002, Andreae,2009]. It is signiﬁcantly lower for icy cirrus clouds (< 104 kg−1 )[e.g. Demott et al. 2003] . On Pluto, it is likely that the organic
haze particles may serve as CCN. In Bertrand and Forget (2016a)
we discuss the possible range of the mass mixing ratio for theseparticles. However, the actual number mixing ratio strongly depends on the degree of aggregation of the monomers and on theiractivation, which is poorly known. In our baseline simulations, we
assumed [CCN] = 105 kg−1 .
2.9. CO cycle
The CO cycle is computed using the same parameterizations
than for methane, modiﬁed to use the CO properties: the CO latent heat is set to LCO = 2.74 × 105 J kg−1 and the saturation mass
mixing ratio q sat CO , is calculated as a function of temperature T (K)and pressure p (Pa) using the following expression derived fromFray and Schmitt (2009):

qsatCO = 0.1537 × 105 e

2.74×105R(1/68.1−1/T ) × MCO × 1

Mair

p

(18)

Here MCO /Mair is the ratio of molecular masses use to convert
volume mixing ratio into mass mixing ratio and R = 8.314/MCO =296.8 m2 s−2 K−1 the CO gas constant.
CO is almost as volatile as N2 and thus much more volatile thanCH4 . In practice,we found that CO only condenses when N2 ice ispresent at the surface, and never forms pure CO deposits. A key parameter controlling the CO is thus the CO mixing ratio in the surface N2 :CO ice solutions. This ratio has been estimated remotelyusing spectroscopic investigations of Pluto. Following the recent
analysis of Very Large Telescope observations by Merlin (2015), we
set this ratio to 0.3%.

Even if we had designed a perfect model of the processes atwork in the Pluto environment, simulating Pluto would remainchallenging. First, in spite of the New Horizons’ achievements, several key parameters remain too poorly known to be used “as observed” (e.g., the global topography). Second, unlike on Mars, theEarth or even Venus, the timescales involved in the evolution of
the climate system at Pluto are so long that it is diﬃcult to reach arealistic model state insensitive to the initial state, even after running the model for weeks of computer time. Here we describe how
we deal with these issues.
3.1. Topography
In our baseline simulations we assume a mostly ﬂat surface except that we placed a 3800 m-deep circular crater roughly at thelocation of Sputnik Planum (in agreement with Moore et al., 2016)as well as two smaller craters corresponding to the informallynamed Burney crater (10 0 0 m deep) and Guest crater (800 mdeep). See Fig. 1.
As discussed below, we also performed sensitivity runs witha perfectly ﬂat topography, and with two additional hypothetical
4 km-high, 800 km wide mountains that we put on the hemisphere opposite to the one better observed by New Horizons (inaddition to the three craters mentioned above).
3.2. Initial subsurface temperatures and ices distribution on thesurfaceOn Pluto, the distribution of surface ices and subsurface temperatures (which plays a key role in the Pluto environment) are
the outcome of thousand of years of evolution (Hansen and Paige,
1996; Young, 2013; Toigo et al., 2015). Running the GCM for such along duration is not feasible. However, initializing the model withprescribed subsurface temperatures and surface ice deposits unrelated to a natural surface evolution may be very unrealistic.
To deal with this issue, as described in Vangichith and Forget(2011) and like Toigo et al. (2015), we designed a reduced versionof the GCM in which the 3D atmopsheric transport and dynamics
are replaced by a simple global mixing function for N2 , CH4 and
CO. Such a model works well on Pluto because the surface energybalance is not signiﬁcantly sensitive to the atmospheric sensible
heat ﬂux and to the radiative transfer through the air. Without atmospheric dynamic and complex radiative transfer to deal with,
we can perform much faster numerical simulations spanning morethan 40,0 0 0 Earth years with the same horizontal grid, the same
subsurface model, and the same surface/atmosphere volatiles exchange parametrizations than with the full GCM.
The details of this reduced model, its validation and the results that we have obtained are described in a separate paper
Bertrand and Forget (2016b). The key ﬁnding is that when we assume a topography map as described above (Fig. 1) with a 3800 mdeep “Sputnik Planum”-like basin and a seasonal ground thermal
inertia larger than 800 J m−3 K−1 , after 40,0 0 0 Earth years theseasonal cycle repeats itself every year with all the nitrogen andCO ices trapped in the “Sputnik Planum”-like basin. This resultsfrom the fact that nitrogen preferentially condenses at lower altitude where the surface pressure is higher, inducing higher condensation temperature and thus enhanced thermal infrared cooling. In this model, methane still undergoes a seasonal cycle andmakes seasonal deposits in both hemispheres, except in an equatorial belt which remains frost-free. Using the set of parametersdescribed in Section 3.4 we establish a realistic, equilibrated initial state for the surface N2 , CH4 and CO deposits and subsurface
temperatures.

Please cite this article as: F. Forget et al., A post-new horizons global climate model of Pluto including the N2 , CH4 and CO cycles, Icarus(2017), http://dx.doi.org/10.1016/j.icarus.2016.11.038

JID: YICAR

ARTICLE IN PRESS[m5G;January 2, 2017;11:46]

F. Forget et al. / Icarus 000 (2017) 1–18

7

Fig. 1. Maps of surface ice distribution and topography at the beginning of the reference and alternative simulations presented in this paper. The black lines show theasssumed topography contours (km). a): Initial state of the reference simulation with no N2 condensation in the south polar region, in Earth year 1988. This state is theoutcome of a 40,0 0 0-Earth-years simulation performed with the reduced 2D model. b): Initial state of the alternative simulation with N2 condensation in the south polar
region, in Earth year 2005. This state is derived from the reference simulation state in 2005, with nitrogen added between 35°N and 48°N and subsurface temperature
poleward of 65°S decreased by 0.5 K.

3.3. Sensitivity to initial atmospheric temperatures and winds
Once the surface ices and subsurface temperatures have been
initialized with the reduced GCM, the full 3D GCM should be
run long enough to reach a realistic regime insensitive to the initial state assumed for the atmosphere. This is challenging becauseof the long radiative time-scale of the Pluto atmosphere (several
Earth years) and the time required to reach established methane
and CO cycles in equilibrium with the surface reservoir. Sensitivity
experiments performed with various initial temperatures, winds,
and atmospheric CH4 and CO contents showed that it takes about
20 years for two simulations initiated with two temperature proﬁles chosen at the end of the realistic possibilities (e.g. differing by
30 K) to differ by less than 2 K. On this basis, we start our simulations at the end of Earth year 1988 and analyse the results after
2010. The convergences of the CO and CH4 cycles are discussed in
Section 5.
3.4. Two kind of simulations
In this paper, we describe two kinds of simulations, with andwithout nitrogen condensation in the south polar region in 2015.
3.4.1. Reference simulation, without N2 condensation at the southpole
For the ﬁrst simulation, we directly use the initial state obtained for Earth date 1988 after 40,0 0 0 Earth years of simulated
climate history performed with the reduced model.
As described by Hansen and Paige (1996) and Young (2013),
the evolution of pressure is sensitive to the surface N2 ice radiative properties. Some tuning was performed to select a reference
value for the N2 ice albedo AN2 and emissivity ε N2 within therange of possible values. By choosing AN2 = 0.67 and εN2 = 0.85,
we obtained an evolution of pressure (shown in Fig. 2) in qualitative agreement with the available observations (Sicardy et al.,2016; Gladstone et al., 2016), reaching a mean surface pressure of
1.1 Pa in July 2015.
Fig. 1 shows the corresponding distribution of ice and subsurface temperature in 1988. In this simulation, the heat stored in
the southern hemisphere during the previous southern hemispheresummer keeps the surface temperature above the nitrogen frost
point, and nitrogen ice is only found in the “Sputnik Planum”-likebasin.
The albedo of the surface CH4 ice deposits was set to ACH4 = 0.5
and its emissivity to εCH4 = εN2 = 0.85. In 1988, Methane frost covers most of the planet except for an equatorial belt which remain

Fig. 2. Evolution of the pressure at r = 1215 km from the planet center (Top) andof the global mean surface pressure (Bottom) in the reference simulation with nosouth pole N2 condensation (red solid line) and in the alternative simulation withsouth pole N2 condensation (blue dashed line) starting at the end of 2005. The
black dots with error bars show the pressure data at r = 1215 km obtained by stellar occultations, as compiled by Sicardy et al. (2016). (For interpretation of the references to colour in this ﬁgure legend, the reader is referred to the web version ofthis article.)

frost free and dark (the albedo and emissivity of the ice-free surfaces were set to A = 0.15 and ε = 1) in agreement with the observations (Stern et al., 2015; Grundy et al., 2016).
3.4.2. Alternative simulation, with N2 condensation at the south pole
It is possible that nitrogen is condensing in the south polar
region in 2015. In that case, we show in this paper that Pluto’s
atmospheric circulation would be very different than without
winter condensation, because of the induced North-south condensation ﬂow. However, to be consistent with the evolution of
surface pressure inferred from the stellar occultations since 1988,

Please cite this article as: F. Forget et al., A post-new horizons global climate model of Pluto including the N2 , CH4 and CO cycles, Icarus(2017), http://dx.doi.org/10.1016/j.icarus.2016.11.038

JID: YICAR

ARTICLE IN PRESS

8[m5G;January 2, 2017;11:46]

F. Forget et al. / Icarus 000 (2017) 1–18

Fig. 3. Maps of surface temperature and winds at 20 m above the surface in July 2015 at different local times for 3 simulations: a) the reference simulations with no N2
condensation at the south pole, b) The same simulation with ﬂat topography (started from the reference run on Juanuary 1st, 2015, and analyzed on July 14, 2015) c) thealternative simulation with N2 condensation at the south pole. From top to bottom, the local time LT in the middle of the map (longitude 180°) is 3:0 0, 9:0 0, 15:0 0 and
21:00, with LT (hours) = [longitude (°) - subsolar point longitude (°)]/15 + 12.

this winter condensation must be balanced by sublimation of
nitrogen frost outside our modeled Sputnik Planum. In fact, New
Horizons observations suggest that mid-northern latitude nitrogen frost deposits were present on Pluto in 2015 (Grundy et al.,2016).
Within that context we designed an artiﬁcial, alternative simulation by taking a modeled state from the ﬁrst reference simulation
at the end of 2005, with two modiﬁcations. First, we added a layerof nitrogen ice in a latitudinal belt between 35°N and 48°N. Second, we decreased the subsurface temperature poleward of 65°S
by 0.5 K to induce nitrogen condensation. This value was chosen
in order to maintain an evolution of pressure similar to the ﬁrst
reference run, as shown in Fig. 2. All other modeled parameters
are the same as in the reference simulation.

4. Model results: temperatures and winds
4.1. Surface temperatures and low level winds
4.1.1. Surface temperatures
Fig. 3 shows maps of surface temperatures and winds at 20 m
above the surface at various times of the day for our different
simulations. The epoch corresponds to July 2015, the time of the
New Horizons encounter. In these simulations, surface temperatures range between 36.6 and 48 K. The lowest values correspond to the N2 frost point around 1 Pa. The highest temperatures are more model dependent, and vary with the assumed diurnal thermal inertia Iday . Daytime surface temperatures reach 57 K
in GCM runs, assuming Iday = 20 J s−1/2 m−2 K−1 (as reported byPlease cite this article as: F. Forget et al., A post-new horizons global climate model of Pluto including the N2 , CH4 and CO cycles, Icarus(2017), http://dx.doi.org/10.1016/j.icarus.2016.11.038

JID: YICAR

ARTICLE IN PRESS[m5G;January 2, 2017;11:46]

F. Forget et al. / Icarus 000 (2017) 1–18

9Fig. 4. Maps of surface temperature and winds at 20 m above the surface in July 2015 in the sub-charon hemisphere, where two artiﬁcial 40 0 0 m-high mountains has been
added to illustrate the formation of downward slope winds on Pluto. The topography is shown by white contours. The local time at longitude 0°E is 3:00 (nighttime) and
15:00 (daytime).

Fig. 5. Zonal-mean temperatures (K) in (a) the reference and (c) the alternative simulations in July 2015. The central plot (b) shows the global mean temperature proﬁles.

Lellouch et al., 2011b) instead of Iday = 50 J s−1/2 m−2 K−1 , as assumed in our baseline simulations.
4.1.2. Slope winds
On ﬂat surfaces and where nitrogen condensation-sublimation
ﬂows are negligible, wind velocities at 20 m remain well below
1 m s−1 . In particular surface temperature gradients do not induce signiﬁcant thermal circulations. As on Mars however, slopes
can create signiﬁcant downward katabatic winds resulting from thefact that the surface is much colder than the atmosphere. The air
close to the slopes is cooled and tends to ﬂow down because it isdenser than the air away from the slope at the same level. Fig. 4 illustrates the formation of such winds on two (hypothetical) 4-km
high, 800-km wide mountains. The wind at 20 m above the surface reaches 4 m s−1 . Because the atmosphere is always warmer
than the surface, and because of its long radiative timescale, thediurnal variations of surface temperature have a limited effect on
the katabatic winds which only increase by 20 % during the nightcompared to the day. Downward katabatic winds can also be observed on the modeled Burney and Guest craters at 45°N in Fig. 3,
left column.
4.1.3. Surface winds induced by condensation-sublimation ﬂows
Wind velocities larger than several meters per second can also
result from the condensation and sublimation of nitrogen. In ourreference circulation (with no condensation at the South pole), this
only occurs in the modeled “Sputnik Planum” area. If one assume aﬂat topography (Fig. 3, center column), intense inward ﬂows form

during the night when N2 condenses, and outward ﬂows are predicted when N2 sublimes during the afternoon. In a more realistic
simulation taking into account the topographic depression in Sputnik Planum (Fig. 3, left column), this effect is combined with theslope winds on the sides of the basin. During the night, when N2
condenses, both slope winds and condensation ﬂows contribute tocreate winds ﬂowing into the modeled Sputnik Planum. During theday, however, the outward sublimation ﬂow is damped by the opposite katabatic ﬂow.
In our alternative model (Fig. 3, right column), N2 condenses inthe south polar region and this sink is balanced by the sublimation of mid-northern latitude N2 deposits. This creates planetaryscale condensation ﬂows from the northern hemisphere toward thesouth pole, and from the dayside toward the nightside. The wind at
20 m reaches several meters per seconds over most of the planet.
In both hemisphere its direction is affected by the Coriolis force,
which prevents the atmosphere from ﬂowing directly southward.
4.2. Atmospheric temperatures
4.2.1. Zonal-mean temperatures
Fig. 5 presents the zonal-mean and global-mean atmospherictemperatures. As found by Toigo et al. (2015), the horizontal
gradients of temperature are very small because of the longradiative timescale. In particular, the meridional variations intemperatures are less than 1 K. In our reference simulation withno south pole N2 condensation, the atmospheric concentration of
methane is realistic (see Section 5.1), and the mean temperature

Please cite this article as: F. Forget et al., A post-new horizons global climate model of Pluto including the N2 , CH4 and CO cycles, Icarus(2017), http://dx.doi.org/10.1016/j.icarus.2016.11.038

JID: YICAR

ARTICLE IN PRESS

10[m5G;January 2, 2017;11:46]

F. Forget et al. / Icarus 000 (2017) 1–18

Fig. 6. Comparison of the two temperature proﬁles retrieved by the New Horizons REX experiment (Hinson et al., 2015b; Gladstone et al., 2016) at 193.5°E, 17.0°S and Local
time 16:31 (red) and 15.7°E, 15.1°N and Local time 04:42 (blue) with GCM results. The model data are taken at the same location and time, except for the proﬁle at latitude
17.0°S which is shifted to 7.5°N in order to locate it just within the modeled Sputnik Planum basin ﬁlled with N2 ice, as it is the case in reality (see text). (For interpretationof the references to colour in this ﬁgure legend, the reader is referred to the web version of this article.)

proﬁle is in acceptable agreement with available observations(Hinson et al., 2015b; Gladstone et al., 2016; Dias-Oliveira et al.,2015), except that above 160 km modeled temperatures are 10–
15 K higher than reported. The thermal structure produced in ouralternative simulation with South pole N2 condensation is evenwarmer, because of the excessive methane concentration in thissimulation, as explained in Section 12.
4.2.2. Comparison with the observed REX proﬁles
In Fig. 6, the simulated temperature proﬁles are compared in
more details with the New Horizons REX radio-occultation proﬁles
obtained at two locations on opposite sides of Pluto. The modeled
proﬁles are taken at the same location and time, except that the
ingress proﬁle is shifted from latitude 17.0°S to 7.5°N, in order tolocate it just inside the modeled Sputnik Planum basin. Indeed,on the real Pluto the ingress proﬁle corresponds to a location justabove the southern tip of the Sputnik Planum depression, above
a surface covered by nitrogen ice. At the same coordinates in ourmodel, we are outside the basin and the surface is frost free. However, we found that taking into account the low topography andN2 coverage is key to understand the differences between the two
REX proﬁles. We plot the modeled temperature proﬁles as a function of altitude above the surface. This creates an apparent shift
in temperatures (the proﬁles are much more similar when shownin pressure coordinates) which contributes to the apparent differences reported in the observations.
Of special interest are the lowest kilometers of the simulated
ingress proﬁles which exhibit a low temperature layer analogous
to the bottom of the observed ingress proﬁle. Which process creates this layer? To better understand this behaviour, and possiblyinterpret the observations, we show in Fig. 7 the diurnal evolutionof the atmospheric proﬁle in the lowest 4 km in different modeled
conﬁgurations. In the reference simulation, the atmospheric temperature in the Spunik Planum basin varies with local time, withcoldest temperatures in the afternoon. This results from the sublimation of nitrogen ice when the sun heats the area, as proposedby Hinson et al. (2015a). In fact, the volume of gas involved in thecondensation-sublimation cycle is considerable in our model. Fig. 8
shows the nitrogen ice budget in the modeled Sputnik Planumbasin at 7.5°N and 45°N. At this last position, about 230 g m−2 ofice sublimates every Pluto day in 2015. As shown on the right axisof Fig. 8, this corresponds to more than 2500 m3 of N2 gas per
square meter. At 7.5°N, the solar ﬂux is weaker in 2015 and thedaily N2 ice budget corresponds to a net gain in N2 ice (net condensation). Nevertheless, every afternoon the equivalent of 800 m3

per square meters is injected into the atmosphere. Moreover, in the
GCM the large amount of cold N2 gas produced at higher latitude(where the insolation is higher) is spread throughout the basin in
the lowest kilometers. In fact, in the alternative simulation thisprocess contributes to increasing the amount of cold air present
in the modeled Spunik Planum basin (Fig. 7b), adding the freshlysublimed cold N2 gas transported from the N2 ice belt at 35°N (asseen on Fig. 3, right column, local time 15:00 and 21:00).
Interestingly, as shown in Fig. 7c, a simulation performed without taking into account the topographic depression in the modeled Spunik-planum does not create a signiﬁcant cold layer. Twofacts explain that. First, the freshly-sublimed N2 gas is eﬃciently
transported away as discussed above (and as seen on Fig. 3, midcolumn). Second, in an atmosphere with radiative timescale aslong as Pluto, in a local topographic depression the temperature
lapse rate is not as steep as on average because temperatures tend
to be homogeneous at a given pressure level. This is illustrated on
Fig. 7d which shows the temperatures at the bottom of the basin
in a simulation with N2 condensation and sublimation completely
switched off. Without N2 sublimation, the air is not as cold as in
the reference simulation, but at a given altitude above the surface,temperatures in the basin remain 5 to 10 K below what they would
have been outside (compare Fig. 7c and d).
4.2.3. Thermal tides and waves
Stellar occultations have shown that vertical proﬁles of density ﬂuctuations in the atmosphere of Pluto often exhibit wave-likestructure (e.g. Elliot et al., 2003b; Person et al., 2008) with an amplitude of a few percent and vertical wavelengths of a few kilometers. On the basis of theoretical calculations, Toigo et al. (2010) suggested that such waves could correspond to the tidal response of
Plutos atmosphere to solar-induced sublimation breathing from N2
frost patches. Here we brieﬂy examine the type of wavelike structure present in the temperature proﬁles generated by our GCM.
Note, however, that the horizontal and vertical resolution used in
the GCM simulations is unlikely to capture waves with vertical
wavelengths smaller than ∼20 km.
Fig. 9a presents the 4-sols evolution of the difference between
instantaneous temperatures and 1-Pluto-day gliding averages at
0°E–0°N in our reference simulation. The observed temperature
excursions are lower than 0.2 K. Nevertheless, they are characteristic of upward atmospheric thermal tides, with, below 80 km,
diurnal, wavenumber = 1 thermal tides with a vertical wavelengtharound 20 km and a downward phase velocity. Above 150 km,
semi-diurnal wavenumber = 2 tide with much longer vertical

Please cite this article as: F. Forget et al., A post-new horizons global climate model of Pluto including the N2 , CH4 and CO cycles, Icarus(2017), http://dx.doi.org/10.1016/j.icarus.2016.11.038

JID: YICAR

ARTICLE IN PRESS
F. Forget et al. / Icarus 000 (2017) 1–18

[m5G;January 2, 2017;11:46]
11

Fig. 7. Diurnal variations of atmospheric temperature in the lower atmosphere at 193.5°E-7.5°N (at the bottom of the modeled Sputnik Planum basin) for (a) the reference
simulation (without South Pole N2 condensation), (b) the alternative simulation (with South Pole N2 condensation), (c) a version of the reference simulation with a ﬂat
topography, and (d) No N2 condensation/sublimation at all on the planet. The simulations with ﬂat topography and No N2 condensation/sublimation were started from thereference run initial state on January 1st, 2015, and analyzed on July 14, 2015.

Fig. 8. Diurnal variation of the surface N2 ice loading at two different latitudes in the modeled “Sputnik Planum” basin in July 2015. The right axis illustrates the corresponding volume of N2 gas, assuming a pressure of 1 Pa and a temperature of 40 K. The different line colours correspond to different kinds of simulations: reference (blue),
alternative with South pole N2 condensation (black, partly hidden by the blue line), and with a ﬂat topography (red). The curves do not loop (i.e. the values at 24:00 differ
from the values at 0:00) because every Pluto day the integrated surface budget corresponds to a net gain of N2 ice by condensation at 7.5°N and a net loss by sublimation
at 45°N, where the incident solar ﬂux is stronger than at 7.5°N. (For interpretation of the references to colour in this ﬁgure legend, the reader is referred to the web versionof this article.)

Fig. 9. Temperature anomaly (difference between instantaneous value and diurnal average) at 0°E–0°N in the reference and alternative simulations in July 2015. Thermaltides are clearly visible in the reference simulation, whereas the alternative simulation is characterized by atmospheric barotropic waves (see text).

Please cite this article as: F. Forget et al., A post-new horizons global climate model of Pluto including the N2 , CH4 and CO cycles, Icarus(2017), http://dx.doi.org/10.1016/j.icarus.2016.11.038

ARTICLE IN PRESS

JID: YICAR
12[m5G;January 2, 2017;11:46]

F. Forget et al. / Icarus 000 (2017) 1–18

Fig. 10. Zonal-mean zonal and meridional winds (m s−1 ) in the reference and alternative simulations in July 2015.

wavelengths start to dominate. As predicted by Toigo et al. (2010),
the source of the tides is the diurnal N2 condensation-sublimation
cycle of the N2 ice: Tidal amplitude are 4-times weaker if N2
condensation-sublimation processes are switched off.
Fig. 9 b presents the same anomaly plot in the alternative simulations with N2 condensation occuring at the south pole. The amplitude of the waves are signiﬁcantly larger, reaching more than ±1 K around 120 km. However, a careful examination of Fig. 9b reveals that the period of the stronger waves is not 1 nor 0.5 Plutoday. These are not thermal tides: the same waves are present insimulations forced by a diurnally-averaged insolation (no diurnal
cycle and no tides). These waves appear to be barotropic waves
produced by a southern polar jet, as described in Section 4.3.2.
4.3. Atmospheric circulation and waves
Fig. 10 shows cross-sections of the average zonal (west-east)and meridional (south-north) winds in our two baseline simulations.
4.3.1. Reference case without N2 condensation at the south pole,In the reference case with no condensation-ﬂow induced byN2 condensation at the south pole, the circulation is relatively
weak with slow retrograde zonal winds in the northern hemisphere and the equatorial regions (Fig. 10a). This circulation remains unchanged with a ﬂat topography, no diurnal cycle, or whenN2 condensation and sublimation processes are switched off. It canbe explained by the north-south latitudinal gradient of solar heating rates. It induces a very small temperature contrast between thespring and fall hemisphere and, in turn, forces weak zonal winds
corresponding to the thermal wind balance. Consistently, the weak
meridional circulation (Fig. 10c) is characterized by a cell centeredat the equator (where the Coriolis force is null) between 80 and
140 km, with the upper branch ﬂowing from the sunlit hemispheretoward the polar night hemisphere.

4.3.2. Alternative case with N2 condensation at the south pole
The circulation is profoundly inﬂuenced by the North-South
condensation ﬂow if N2 condenses in the South polar regions.
If N2 ice deposits were covering the entire northern polar regions (which is not observed) and the southern hemisphere condensation much more intense, the condensation ﬂow would bevery strong. As obtained in some of our past simulations (notshown) and as reported in some scenarios analysed by Toigo et al.
(2015) (see their Fig. 11 and 18), the meridional circulation would
be characterized by a global ﬂow from the northern hemisphere to
the southern hemisphere. In such conditions, the zonal circulationis characterized by a global “retro-superrotation” with retrograde
winds at most latitude. Such winds result from the conservation ofangular momentum of the air particles as they ﬂow from the sunlit
pole to the polar night above the equator, where they are farther
from the rotation axis than where they started from.
In our simulations however, the North-South condensation ﬂow
remains limited compared to this extreme case. We consider
that this is in better agreement with the observations because
1) outside Sputnik Planum the N2 ice frost deposits are limited to patches around 45–60°N (Grundy et al., 2016), and 2) because the south pole N2 condensation cannot be very intense in2015 since Pluto’s surface pressure has been increasing in recent
years.
With the realistic assumptions made in our “alternative” simulation, the meridional circulation remains weak (Fig. 10d) andstrongly modulated by waves (see below). The overall transport
pattern is southward, as revealed when analysing tracer transport(Bertrand and Forget, 2016a).
The zonal wind is charaterized by an intense prograde jetstream poleward of 40°S and a prograde superrotation at mostother latitudes (Fig. 10c). The high-latitude jet is a classical featureof terrestrial atmospheres, and likely result here from the poleward condensation ﬂow and the conservation of angular momentum. Superrotation is more surprising. It is observed on Venus andPlease cite this article as: F. Forget et al., A post-new horizons global climate model of Pluto including the N2 , CH4 and CO cycles, Icarus(2017), http://dx.doi.org/10.1016/j.icarus.2016.11.038

JID: YICAR

ARTICLE IN PRESS[m5G;January 2, 2017;11:46]

F. Forget et al. / Icarus 000 (2017) 1–18

13

Fig. 11. Characteristics of the barotropic waves present in the simulations with South Pole N2 condensation inducing a condensation ﬂow. a) Hövmoller diagram of thetemperature anomaly (difference between the local and the zonal-mean temperature) at 0°N. b) Same at 60°N. c) Zonal average of the root-mean-square standard deviationof the local meridional wind from the zonal-averaged meridional wind. d) Hövmoller diagram of the meridional wind anomaly (difference between the local and the zonalmean wind, in m s−1 ) at 0°N.

Titan and has been the subject of many studies (see, e.g. Lebonnois
et al., 2010, and references therein). In these cases, superrotation isconsidered to primarily result from the so-called Gierasch-RossowWilliams mechanism (from Gierasch, 1975; Rossow, 1979). In thismechanism, waves, possibly generated by barotropic instabilities
from the high-latitude jets, redistribute angular momentum equatorward. Preliminary analysis suggest that this is what is happening in our simulation. A study of the variations of the high-latitude
jet show that it is subject to instabilities that create a wavenumber 1 wave that propagates eastward with a 0.5–0.8 Pluto day
period. At 60°N, such waves are clearly visible at an altitude of
140 km in the temperature and meridional wind ﬁelds (Fig. 11b
and d). In Fig. 11c, the extension of this wave is mapped by plotting the meridional wind variability as a function of latitude andheight. One can see that it propagates to all latitudes, and notably
to the equator, where the signature in the thermal ﬁeld dominates
the temperature variability (Fig. 11a). Similar results are obtainedin model runs without a diurnal cycle or with a ﬂat topography.
In addition to the wind predictions published by Toigo et al.
(2015), already discussed, our results can be compared with the results from the other Pluto GCM proposed by Zalucha and Michaels(2013) and Zalucha (2016). The comparison with Zalucha andMichaels (2013) is diﬃcult to achieve because this version of their
GCM did not yet include nitrogen condensation and because theirmodeled thermal structure was very different than what was observed on Pluto by New Horizons. In fact the updated version presented by Zalucha (2016) yielded completely different results. Her“Case1”, in which a surface pressure of 0.8 Pa and 1% of CH4 is assumed, can be compared to our simulations. The zonal wind structure ressemble our reference simulation, suggesting that, for un-

known reasons, the condensation ﬂow is weak in this GCM in spiteof the fact that Pluto is assumed to be covered by nitrogen ice.
5. Model results: CH4 and CO cycles
5.1. Evolution and distribution of gaseous CH4
Fig. 12 shows the global-mean mixing ratio of methane (determined from the ratio of CH4 and N2 column densities) in ourbaseline simulations, and how this ratio varies over time. Fig. 12cshows the evolution of the global-mean mixing ratio of methane.
The three red curves correspond to reference simulations (without poleward condensation ﬂow) with methane volume mixing ratio inialized at 0.1%, 0.5% and 1% in 1988. One can see that in2015 the results are still sensitive to the initialization, although
the three simulations clearly converge toward a global mean value
near 0.5 %. Fig. 12a and b present the zonal-mean methane abundances as a function of latitude and altitude in 2010 (mid-point
between the 2008 and 2012 observations by Lellouch et al., 2015)
and 2015 (New Horizons). These ﬁgures show that methane is nothomogeneously distributed, notably because the high latitude deposits are increasingly heated and sublimed as the sub-solar point
moves northward with time. As a result methane tends to be enriched in the lower atmosphere at high northern latitudes compared to the rest of the planet, but is otherwise vertically well
mixed and near 0.5% at most altitudes. This is consistent with theobservation analysis of Lellouch et al. (2015) who concluded that
their data “imply a roughly uniform mixing ratio in at least theﬁrst 22 - 27 km of the atmosphere”, and that “high concentrationsof low-temperature methane near the surface can be ruled out”.

Please cite this article as: F. Forget et al., A post-new horizons global climate model of Pluto including the N2 , CH4 and CO cycles, Icarus(2017), http://dx.doi.org/10.1016/j.icarus.2016.11.038

JID: YICAR

ARTICLE IN PRESS

14[m5G;January 2, 2017;11:46]

F. Forget et al. / Icarus 000 (2017) 1–18

Fig. 12. a) - b): Zonal mean methane volume mixing ratio (%) in the reference simulation (without south pole N2 condensation and [CH4 ] initially at 0.5% in 1988) in 2010
and 2015. c): Evolution of the mean volume mixing ratio: globally averaged with different initialization (red), the apparent mixing ratio as seen from the Earth (green, seetext) and in the alternative simulation with south pole N2 condensation started in 2005 (dashed blue). d): Zonal mean methane volume mixing ratio (%) in the alternative
simulation (with south pole N2 condensation) in 2015. (For interpretation of the references to colour in this ﬁgure legend, the reader is referred to the web version of thisarticle.)To compare with Earth-based near-infrared observations, one mustnevertheless take into account the fact that such observations arebiased toward the methane column near the sub-Earth-subsolar
points for geometrical reason (Pluto is a sphere) and because thisis where the insolation is maximum. Taking into account that the
sub-Earth and subsolar points are always very close, we can estimate the apparent mixing ratio as seen from the Earth by performing an average of the local column mixing ratio weighted by
the square of the cosine of the solar zenith angle. The apparent
mixing ratio for the reference simulations started with 0.5 % CH4
is shown in green on Fig. 12c. The difference with the global-meanvalue remain small and has only become signiﬁcant recently withthe local increase of methane above the North pole.
On the same Fig. 12c, the blue dashed curve shows the evolution of the global-mean methane in the alternative scenario (withN2 condensing at the south pole) starting in 2005. Fig. 12d show
the corresponding methane abundances as a function of latitude
and altitude in 2015. One can see that the methane content islarger and still increasing in this simulation. This results from thestronger near-surface winds induced by the condensation ﬂow, andthe fact that the near-surface mixing is directly proportional to thehorizontal wind as formalized in Eq. (8) presented in Section 2.5.
5.2. Formation of CH4 ice clouds
Fig. 13 shows maps of methane ice clouds in our reference andalternative simulations at various local time in July 2015. In both
simulations, atmospheric condensation is induced by the subliming
nitrogen ice on the surface. On the dayside, freshly-sublimed nitrogen gas tends to cool the atmosphere nearby and trigger methane

condensation in the ﬁrst hundreds of meters above the surface, as
illustrated in Fig. 14. In the alternative simulations with surface N2
ice between 35°N and 48°N, the cold air and the clouds particles
are transported by the sublimation ﬂows (see Fig. 3, right column)
and can extend outside the N2 ice covered regions, reaching 20°N
and 75°N.
5.3. CO cycle
Fig. 15 shows the evolution of the carbon monoxide mixing ratio as a function of time since 1988. The red curves correspond to
the global-averaged mixing ratio for three different initial values(0%, 0.05%, 0.1%). Clearly, the three simulations have not converged
but one can estimate that the system evolves toward a mean mixing ratio near 0.03%. A mixing ratio of 0.03% is in acceptable agree.01
ment with the 0.05+0−0.025 % reported by Lellouch et al. (2011a) fromtelescopic observations performed in 2010, and of the same orderof magnitude as the 0.0515 ± 0.004% just retrieved by Lellouch
et al. (2016) using the ALMA interferometer on June 12–13, 2015.
In details, the CO cycle is dominated by a condensationsublimation cycle inside Sputnik Planum. For instance in 2015
there is a net ﬂux from the northern part and the center part of
Sputnik Planum to the southern part where nitrogen is condensing
along with CO. We do not show here the spatial distribution of CO
since we have found that CO is usually very well mixed with N2 .
As a consequence, the apparent CO mixing ratio as seen from the
Earth (green curve in Fig. 15) is very close to the global mean.
When the alternative simulation is started in 2005 with N2 condensing in the high southern latitudes (blue lines in Fig. 15), theCO mixing ratio rapidly decreases to reach values below 0.03%. This

Please cite this article as: F. Forget et al., A post-new horizons global climate model of Pluto including the N2 , CH4 and CO cycles, Icarus(2017), http://dx.doi.org/10.1016/j.icarus.2016.11.038

JID: YICAR

ARTICLE IN PRESS
F. Forget et al. / Icarus 000 (2017) 1–18

[m5G;January 2, 2017;11:46]
15

Fig. 13. Maps of methane ice clouds mass (10−6 kg per m2 ) in July 2015 for the reference and alternative simulations for different local times at center of the map (180°E).
The black contours show the atmospheric temperature 20 m above the surface.

Fig. 14. Methane clouds as a function of latitude and altitude above the surface, around July 14, 2015, in a) the reference simulation (N2 only in Sputnik Planum) at longitude
180°E and Local Time 15:00, and b) the alternative simulation (with surface N2 ice between 35°N and 48°N) at longitude 0°E and Local Time 15:00.

Please cite this article as: F. Forget et al., A post-new horizons global climate model of Pluto including the N2 , CH4 and CO cycles, Icarus(2017), http://dx.doi.org/10.1016/j.icarus.2016.11.038

JID: YICAR

ARTICLE IN PRESS

16[m5G;January 2, 2017;11:46]

F. Forget et al. / Icarus 000 (2017) 1–18

very grateful to two anonymous reviewers for their exceptionallydetailed comments which helped to improve the readibility of thearticle.
Appendix A. Computing mass, momentum and heat verticalﬂuxes induced by N2 condensation and sublimation in the GCM
vertical coordinates

Fig. 15. Evolution of the mean volume mixing ratio of gaseous carbon monoxide.
The red curves present the globally averaged values with different initialization.
The green curve shows the apparent mixing ratio as seen from the Earth. The blue
curves shows the global mean mixing ratio in the alternative simulation with southpole N2 condensation started in 2005 (see text). (For interpretation of the references to colour in this ﬁgure legend, the reader is referred to the web version ofthis article.)is even the case when we assume that all mid-northern latitude N2
frost deposits contains 0.3% of CO. In these conditions, the atmospheric CO appears to decrease below 0.03% because the ices thatcondense in the south polar cap tends to be enriched in CO, up to
0.05% at the pole.
In reality, the mid-latitude N2 frost deposits have been observed
by New Horizons to be strongly depleted of CO compared to Sputnik Planum (Grundy et al., 2016). If we take this into account andset the N2 :CO mixing ratio to zero in these deposits, we obtain
the evolution shown by the dashed blue line Fig. 15, with an additional decrease of atmospheric CO down to less than 0.01% in2015. One can guess that these values could be tuned up by increasing the assumed N2 :CO ice mixing ratio in Sputnik Planum.
This would still be consistent with the Merlin (2015)’s telescopic
measurements since they included both Sputnik Planum and the
mid-latitude deposits. Further work will be required to fully understand the long term CO equilibrium, its evolution, and the surface
N2 :CO mixing ratio.
5.4. Conclusions
The goal of this paper was to describe, for the ﬁrst time, our
new Global Climate Model of Pluto including the N2 , CH4 andCO cycles. We presented two baseline simulations which can shedlight on New Horizons observations, for instance to understand the
low atmosphere temperature proﬁles measured by REX and the
distribution of ices. However, this is just the beginning. One of our
key conclusions is that the Pluto climate system is extremely sensitive to the assumed model parameters, such as the ice properties,
the ground thermal inertia, or the topography. Many more studies
will have to be performed to better simulate the reality and understand the processes at work. It will also be very useful to perform
longer simulations, with higher model resolution, with a more realistic topography, etc... We hope that this GCM will be applied tomany speciﬁc studies regarding clouds, hazes, frost deposits, seasonal evolution, and paleoclimates.
Acknowledgments
The authors thank the NASA New Horizons team for their excellent work on a fantastic mission and their interest in this research. We also thank CNES for its support. Finally, the authors are

In the GCM, the changes in atmospheric mass due to the condensation and sublimation of nitrogen are taken into account bymodifying the surface pressure p0 at each timestep by: δ p0 =−g N
k=0 δ mk , with N the number of atmospheric model layers andδ mk the mass condensed (or sublimed if < 0) in layer k or at thesurface (k = 0), as described in Section 2.6. This ensures the conservation of the total mass of N2 (surface caps + atmosphere).
As described in Section 2.1, the vertical coordinate of eachmodel layer is deﬁned by its σl = pl /p0 coordinates. The changes
in p0 due to the N2 condensation-sublimation induce “artiﬁcial”
movements of the σ levels in the atmosphere. This must be reﬂected in the temperature and wind ﬁelds.
Consider a layer l delimited by the levels σl− 1 and σl+ 1 . At22each timestep, its mass Ml = g0 (σl− 1 − σl+ 1 ) (in kg m−2 ) variesp2

2

because of the global variation of p0 . Such a variation δ Ml is associated with transfers of mass between the layers (on which one
must add the sink corresponding to the local condensation −δ ml ).
The local mass balance may be written :δ Ml =

δ p0g(σl− 12 − σl+ 12 ) = Wl− 12 − Wl+ 12 − δ ml

(19)

where Wl− 1 is the air mass (kg m−2 ) “transfered” through the level2

σl− 1 (> 0 when up) during the timestep. Eqs. (19) can be rear2

anged to yield a recursive formula on W :Wl+ 1 = Wl− 1 − δ ml −
2

2

δ p0
g

(σl− 12 − σl+ 12 )(20)

with, in the ﬁrst layer:

W 1 = −δ m0

(21)

2

The knowledge of W can then be used to compute the exchangeof heat and momentum between the layers. For cp T (enthalpy), thelocal heat balance can be written :δ (Ml Tl ) = Wl− 12 T l− 12 − Wl+ 12 T l+ 12 − δ ml T cl

(22)

with T l− 1 the mean temperature of the gas transported through2

the σl− 1 interface. The calculation of T l− 1 is like in a classical22transport problem. We use the “Van-Leer I” ﬁnite volume transport scheme (Van Leer, 1977; Hourdin and Armengaud, 1999). Separately, one can also write :δ (Ml Tl ) = (Ml + δ Ml )δ Tl + Tl δ Ml

(23)

with δ Tl the correction to be applied at every timestep in eachlayer after the N2 condensation or sublimation. Eqs. (22) and(23) may be combined to obtain δ Tl

δ Tl =

1[W 1 (T 1 −T ) − Wl+ 1 (T l+ 1 −Tl ) − δ ml (T cl −Tl )]22
Ml + δ Ml l− 2 l− 2 l(24)

The ﬁrst two terms, with Wl− 1 and Wl+ 1 , correspond to the22

re-arrangement of the temperatures over the entire column due tothe pressure variations in σ coordinates. The last term δ ml (T cl −Tl ) is negligible when N2 condenses or partially sublimes since wethen have T cl = Tl . However, when the N2 totally sublimes , it becomes a cooling term accounting for the mixing of the newly sublimed mass −δ ml with the rest of the layer at Tl > Tcl .

Please cite this article as: F. Forget et al., A post-new horizons global climate model of Pluto including the N2 , CH4 and CO cycles, Icarus(2017), http://dx.doi.org/10.1016/j.icarus.2016.11.038

ARTICLE IN PRESS

JID: YICAR

F. Forget et al. / Icarus 000 (2017) 1–18

On the ground, if δ m0 > 0 (condensation), we set T 1 = T1 . As
2mentioned above, the near-surface cooling of the condensing N2
gas from T1 to T0 is then taken into account in the surface energy balance. If δ m0 < 0 (sublimation), we set T 1 = T0 . The term2

δ m0 (T0 − T1 ) then accounts for the cooling of the lowest level by
the freshly-sublimed nitrogen.
Similarly, the momentum distribution must be re-arranged. Fora wind component v, we shall simply write:

δvl =

1
[W 1 (v 1 − vl ) − Wl+ 1 (vl+ 1 − vl )]22
Ml l− 2 l− 2(25)

with, on the ground, v 1 = v1 if δ m0 > 0 and v 1 = 0 if δ m0 < 02

2(the velocity of the N2 gas that has just sublimed is zero).
References
Arakawa, A., Lamb, V.R., 1977. Computational design of the UCLA general circulation
model. Methods Comput. Phys. 17, 174–267.
Bertrand, T., Forget, F., 2016a. 3d modeling of organic haze in Pluto s atmosphere,
icarus, this issue.
Bertrand, T., Forget, F., 2016. Observed glacier and volatile distribution on Pluto fromatmosphere-topography processes. Nature 540, 86–89.
Boursier, C., Menard, J., Doyennette, L., Menard-Bourcin, F., 2003. Rovibrational relaxation of methane in CH4n2 mixtures: time-resolved IRIR doubleresonancemeasurements at 193 k and kinetic modeling. J. Phys. Chem. 107, 52805290.
Demott, P.J., et al., 2003. Measurements of the concentration and composition of
nuclei for cirrus formation. Proc. Natl. Acad. Sci. 100, 14655–14660.
Dias-Oliveira, A., Sicardy, B., Lellouch, E., Vieira-Martins, R., Assaﬁn, M., Camargo, J.I.B., Braga-Ribas, F., Gomes-Júnior, A.R., Benedetti-Rossi, G., Colas, F., Decock, A., Doressoundiram, A., Dumas, C., Emilio, M., Fabrega Polleri, J., Gil-Hutton, R., Gillon, M., Girard, J.H., Hau, G.K.T., Ivanov, V.D., Jehin, E., Lecacheux, J.,
Leiva, R., Lopez-Sisterna, C., Mancini, L., Manfroid, J., Maury, A., Meza, E.,
Morales, N., Nagy, L., Opitom, C., Ortiz, J.L., Pollock, J., Roques, F., Snodgrass, C.,
Soulier, J.F., Thirouin, A., Vanzi, L., Widemann, T., Reichart, D.E., LaCluyze, A.P.,
Haislip, J.B., Ivarsen, K.M., Dominik, M., Jørgensen, U., Skottfelt, J., 2015. Pluto s
atmosphere from stellar occultations in 2012 and 2013. Astrophys. J. 811, 53.
Douté, S., Schmitt, B., Quirico, E., Owen, T.C., Cruikshank, D.P., de Bergh, C.,
Geballe, T.R., Roush, T.L., 1999. Evidence for methane segregation at the surfaceof Pluto. Icarus 142, 421–444.
Elliot, J.L., Ates, A., Babcock, B.A., Bosh, A.S., Buie, M.W., Clancy, K.B., Dunham, E.W., Eikenberry, S.S., Hall, D.T., Kern, S.D., Leggett, S.K., Levine, S.E.,
Moon, D.-S., Olkin, C.B., Osip, D.J., Pasachoff, J.M., Penprase, B.E., Person, M.J.,
Qu, S., Rayner, J.T., Roberts, L.C., Salyk, C.V., Souza, S.P., Stone, R.C., Taylor, B.W.,
Tholen, D.J., Thomas-Osip, J.E., Ticehurst, D.R., Wasserman, L.H., 2003. The recent expansion of Pluto’s atmosphere. Nature 424, 165–168.
Elliot, J.L., Ates, A., Babcock, B.A., Bosh, A.S., Buie, M.W., Clancy, K.B., Dunham, E.W., Eikenberry, S.S., Hall, D.T., Kern, S.D., Leggett, S.K., Levine, S.E.,
Moon, D.-S., Olkin, C.B., Osip, D.J., Pasachoff, J.M., Penprase, B.E., Person, M.J.,
Qu, S., Rayner, J.T., Roberts, L.C., Salyk, C.V., Souza, S.P., Stone, R.C., Taylor, B.W.,
Tholen, D.J., Thomas-Osip, J.E., Ticehurst, D.R., Wasserman, L.H., 2003. The recent expansion of Pluto’s atmosphere. Nature 424, 165–168.
Elliot, J.L., Dunham, E.W., Bosh, A.S., Slivan, S.M., Young, L.A., Wasserman, L.H., Millis, R.L., 1989. Pluto’s atmosphere. Icarus 77, 148–170.
Elliot, J.L., Person, M.J., Gulbis, A.A.S., Souza, S.P., Adams, E.R., Babcock, B.A., Gangestad, J.W., Jaskot, A.E., Kramer, E.A., Pasachoff, J.M., Pike, R.E., Zuluaga, C.A.,
Bosh, A.S., Dieters, S.W., Francis, P.J., Giles, A.B., Greenhill, J.G., Lade, B., Lucas, R.,
Ramm, D.J., 2007. Changes in Pluto’s atmosphere: 1988–2006. Astrophys. J. 134,
1–13.
Forget, F., Haberle, R.M., Montmessin, F., Levrard, B., Head, J.W., 2006. Formationof glaciers on Mars by atmospheric precipitation at high obliquity. Science 311,
368–371.
Forget, F., Hourdin, F., Fournier, R., Hourdin, C., Talagrand, O., Collins, M., Lewis, S.R.,
Read, P.L., Huot., J.-P., 1999. Improved general circulation models of the martian atmosphere from the surface to above 80 km. J. Geophys. Res. 104 (24),155–176.
Forget, F., Hourdin, F., Talagrand, O., 1998. CO2 Snow fall on Mars: simulation witha general circulation model. Icarus 131, 302–316.
Fray, N., Schmitt, B., 2009. Sublimation of ices of astrophysical interest: a bibliographic review. Planet. Space Sci. 57, 2053–2080.
Gierasch, P.J., 1975. Meridional circulation and the maintenance of the Venus atmospheric rotation. J. Atmos. Sci. 32, 1038–1044.
Gladstone, G.R., Stern, S.A., Ennico, K., Olkin, C.B., Weaver, H.A., Young, L.A., Summers, M.E., Strobel, D.F., Hinson, D.A., Kammer, J.A., Parker, A.H., Steﬄ, A.J., Linscott, I.R., Parker, J.W., Cheng, A.F., Slater, D.C., Versteeg, M.H., Greathouse, T.K.,
Retherford, K.D., Throop, H., Cunningham, N.J., Woods, W.W., Singer, K.N.,
Tsang, C.C.C., Schindhelm, E., Lisse, C.M., Wong, M.L., Yung, Y.L., Zhu, X.,
Curdt, W., Lavvas, P., Young, E.F., Tyler, G.L., 2016. The atmosphere of Pluto as
observed by new horizons. Science 351, 1280.

[m5G;January 2, 2017;11:46]
17

González-Galindo, F., Forget, F., López-Valverde, M.A., Angelats i Coll, M., Millour, E.,
2009. A ground-to-exosphere martian general circulation model. 1. seasonal,
diurnal and solar cycle variation of thermospheric temperatures. J. Geophys.
Res.(Planets) 114 (E13), 4001.
Grundy, W.M., Binzel, R.P., Buratti, B.J., Cook, J.C., Cruikshank, D.P., DalleOre, C.M.,
Earle, M., Ennico, K., Howett, C.J.A., Lunsford, A.W., Olkin, C.B., Parker, A.H.,
Philippe, S., Protopapa, S., Quirico, E., Reuter, D.C., Schmitt, B., Singer, K.N.,
Verbiscer, A.J., Beyer, R.A., Buie, M.W., Cheng, A.F., Jennings, D.E., Linscott, I.R.,
Parker, J.W., Schenk, P.M., Spencer, J.R., Stansberry, J.A., Stern, S.A., Throop, H.B.,
Tsang, C.C.C., Weaver, H.A., Weigle II, G.E., Young, L.A., 2016. Surface compositions across Pluto and Charon. Science 351, 1283.
Grundy, W.M., Olkin, C.B., Young, L.A., Buie, M.W., Young, E.F., 2013. Near-infrared
spectral monitoring of Pluto’s ices: spatial distribution and secular evolution.
Icarus 223, 710–721.
Gurwell, M., Lellouch, E., Butler, B., Moullet, A., Moreno, R., Bockelée-Morvan, D.,
Biver, N., Fouchet, T., Lis, D., Stern, A., Young, L., Young, E., Weaver, H., Boissier,J., Stansberry, J., 2015. Detection of atmospheric CO on Pluto with ALMA.
AAS/division for planetary sciences meeting abstracts. vol. 47 of AAS/division
for planetary sciences meeting abstracts. p. 105.06.
Hansen, C.J., Paige, D.A., 1996. Seasonal nitrogen cycles on Pluto. Icarus 120,
247–265.
Hinson, D. P., Linscott, I., Tyler, G. L., Bird, M. K., Paetzold, M., Strobel, D. F., Summers, M. E., Woods, W. W., Stern, A., Weaver, Jr., H. A., Olkin, C., Young, L. A.,
Ennico Smith, K., Gladstone, R., Greathouse, T., Kammer, J., Parker, A. H., Parker,
J. W., Retherford, K. D., Schindhelm, E., Singer, K. N., Steﬄ, A., Tsang, C., Versteeg, M., 2015a. Radio occultation measurements of Pluto’s atmosphere withnew horizons. AGU fall meeting abstracts.
Hinson, D. P., Linscott, I., Tyler, L., Bird, M., Paetzold, M., Strobel, D., Summers, M.,
Woods, W., Stern, A., Weaver, H., Olkin, C., Young, L., Ennico, K., Gladstone, R.,
Greathouse, T., Kammer, J., Parker, A., Parker, J., Retherford, K., Schindhelm, E.,
Singer, K., Steﬄ, A., Tsang, C., Versteeg, M., 2015b. Radio occultation measurements of Pluto’s atmosphere with new horizons. AAS/division for planetary sciences meeting abstracts. vol. 47. p. 105.01.
Hourdin, F., Armengaud, A., 1999. Test of a hierarchy of ﬁnite-volume schemes fortransport of trace species in an atmospheric general circulation model. Mon.
Wea. Rev. 127, 822–837.
Hourdin, F., Forget, F., Talagrand, O., 1995. The sensitivity of the martian surfacepressure to various parameters: a comparison between numerical simulationsand viking observations. J. Geophys. Res. 100, 5501–5523.
Hubbard, W.B., Hunten, D.M., Dieters, S.W., Hill, K.M., Watson, R.D., 1988. Occultation evidence for an atmosphere on Pluto. Nature 336, 452–454.
Hubbard, W.B., Yelle, R.V., Lunine, J.I., 1990. Nonisothermal Pluto atmosphere models. Icarus 84, 1–11.
Lebonnois, S., Burgalat, J., Rannou, P., Charnay, B., 2012. Titan global climate model:
a new 3-dimensional version of the IPSL titan GCM. Icarus 218, 707–722.
Lebonnois, S., Hourdin, F., Eymet, V., Crespin, A., Fournier, R., Forget, F., 2010. Superrotation of Venus’ atmosphere analyzed with a full general circulation model. J.
Geophys. Res.(Planets) 115, 6006.
Lefèvre, F., Bertaux, J.-L., Clancy, R.T., Encrenaz, T., Fast, K., Forget, F., Lebonnois, S.,
Montmessin, F., Perrier, S., 2008. Heterogeneous chemistry in the atmosphere ofMars. Nature 454, 971–975.
Lellouch, E., de Bergh, C., Sicardy, B., Forget, F., Vangvichith, M., Käuﬂ, H.-U., 2015.
Exploring the spatial, temporal, and vertical distribution of methane in Pluto’s
atmosphere. Icarus 246, 268–278.
Lellouch, E., de Bergh, C., Sicardy, B., Käuﬂ, H.U., Smette, A., 2011. High resolution
spectroscopy of Pluto’s atmosphere: detection of the 2.3 μm CH4 bands andevidence for carbon monoxide. Astron. Astrophys. 530, L4.
Lellouch, E., Gurwell, M., Butler, B., Fouchet, T., Lavvas, P., Strobel, D. F., Sicardy, B.,
Moullet, A., Moreno, R., Bockelée-Morvan, D., Biver, N., Young, L., Lis, D., Stansberry, J., Stern, A., Weaver, H., Young, E., Zhu, X., Boissier, J., in press, D. o. C.,
Icarus, H. i. P. a. w. A., 2016.
Lellouch, E., Sicardy, B., de Bergh, C., Käuﬂ, H.-U., Kassi, S., Campargue, A., 2009.
Pluto’s lower atmosphere structure and methane abundance from high-resolution spectroscopy and stellar occultations. Astron. Astrophys. 495, L17–L21.
Lellouch, E., Stansberry, J., Emery, J., Grundy, W., Cruikshank, D.P., 2011. Thermal
properties of Pluto’s and Charon’s surfaces from spitzer observations. Icarus 214,
701–716.
Mellor, G.L., Yamada, T., 1982. Development of a turbulence closure model for geophysical ﬂuid problems. Rev. Geophys. 20 (4), 851–875.
Merlin, F., 2015. New constraints on the surface of Pluto. Astron. Astrophys. 582,
A39.
Moore, J.M., McKinnon, W.B., Spencer, J.R., Howard, A.D., Schenk, P.M., Beyer, R.A.,
Nimmo, F., Singer, K.N., Umurhan, O.M., White, O.L., Stern, S.A., Ennico, K.,
Olkin, C.B., Weaver, H.A., Young, L.A., Binzel, R.P., Buie, M.W., Buratti, B.J.,
Cheng, A.F., Cruikshank, D.P., Grundy, W.M., Linscott, I.R., Reitsema, H.J.,
Reuter, D.C., Showalter, M.P., Bray, V.J., Chavez, C.L., Howett, C.J.A., Lauer, T.R.,
Lisse, C.M., Parker, A.H., Porter, S.B., Robbins, S.J., Runyon, K., Stryk, T.,
Throop, H.B., Tsang, C.C.C., Verbiscer, A.J., Zangari, A.M., Chaikin, A.L., Wilhelms, D.E., 2016. The geology of Pluto and Charon through the eyes of newhorizons. Science 351, 1284.
Olkin, C.B., Young, L.A., Borncamp, D., Pickles, A., Sicardy, B., Assaﬁn, M., Bianco, F.B.,
Buie, M.W., de Oliveira, A.D., Gillon, M., French, R.G., Ramos Gomes, A., Jehin, E., Morales, N., Opitom, C., Ortiz, J.L., Maury, A., Norbury, M., Braga-Ribas, F.,
Smith, R., Wasserman, L.H., Young, E.F., Zacharias, M., Zacharias, N., 2015. Evidence that Pluto’s atmosphere does not collapse from occultations including
the 2013 May 04 event. Icarus 246, 220–225.

Please cite this article as: F. Forget et al., A post-new horizons global climate model of Pluto including the N2 , CH4 and CO cycles, Icarus(2017), http://dx.doi.org/10.1016/j.icarus.2016.11.038

JID: YICAR

ARTICLE IN PRESS

18[m5G;January 2, 2017;11:46]

F. Forget et al. / Icarus 000 (2017) 1–18

Owen, T.C., Roush, T.L., Cruikshank, D.P., Elliot, J.L., Young, L.A., de Bergh, C.,
Schmitt, B., Geballe, T.R., Brown, R.H., Bartholomew, M.J., 1993. Surface ices andthe atmospheric composition of Pluto. Science 261, 745–748.
Person, M.J., Elliot, J.L., Gulbis, A.A.S., Zuluaga, C.A., Babcock, B.A., McKay, A.J., Pasachoff, J.M., Souza, S.P., Hubbard, W.B., Kulesa, C.A., McCarthy, D.W., Benecchi, S.D., Levine, S.E., Bosh, A.S., Ryan, E.V., Ryan, W.H., Meyer, A., Wolf, J., Hill, J.,
2008. Waves in Pluto’s upper atmosphere. Astron. Jour. 136, 1510–1518.
Rossow, W.B., 1978. Cloud microphysics: analysis of the clouds of Earth, Venus,
Mars, and Jupiter. Icarus 36 (1), 50.
Rossow, W.B., 1979. Large–scale motion in the Venus stratosphere. J. Atmos. Sci. 36,
377–389.
Sicardy, B., Talbot, J., Meza, E., Camargo, J.I.B., Desmars, J., Gault, D., Herald, D.,
Kerr, S., Pavlov, H., Braga-Ribas, F., Assaﬁn, M., Benedetti-Rossi, G., Dias-Oliveira, A., Ramos-Gomes-Jr., A., Vieira-Martins, R., Berard, D., Kervella, P.,
Lecacheux, J., Lellouch, E., Beisker, W., Dunham, D., Jelinek, M., Duffard, R.,
Ortiz, J.L., Castro-Tirado, A.J., Cunniffe, R., Querel, R., Yock, P.A., Cole, A.A.,
Giles, A.B., Hill, K.M., Beaulieu, J.P., Harnisch, M., Jansen, R., Pennell, A., Todd, S.,
Allen, W.H., Graham, P.B., Loader, B., McKay, G., Milner, J., Parker, S., Barry, M.A.,
Bradshaw, J., Broughton, J., Davis, L., Devillepoix, H., Drummond, J., Field, L.,
Forbes, M., Giles, D., Glassey, R., Groom, R., Hooper, D., Horvat, R., Hudson, G.,
Idaczyk, R., Jenke, D., Lade, B., Newman, J., Nosworthy, P., Purcell, P., Skilton, P.F.,
Streamer, M., Unwin, M., Watanabe, H., White, G.L., Watson, D., 2016. Pluto’s atmosphere from the 29 June 2015 ground-based stellar occultation at the timeof the new horizons ﬂyby. Astrophys. J. Lett. 819, L38.
Stern, S.A., Bagenal, F., Ennico, K., Gladstone, G.R., Grundy, W.M., McKinnon, W.B.,
Moore, J.M., Olkin, C.B., Spencer, J.R., Weaver, H.A., Young, L.A., Andert, T., Andrews, J., Banks, M., Bauer, B., Bauman, J., Barnouin, O.S., Bedini, P., Beisser, K.,
Beyer, R.A., Bhaskaran, S., Binzel, R.P., Birath, E., Bird, M., Bogan, D.J., Bowman, A., Bray, V.J., Brozovic, M., Bryan, C., Buckley, M.R., Buie, M.W., Buratti, B.J., Bushman, S.S., Calloway, A., Carcich, B., Cheng, A.F., Conard, S., Conrad, C.A., Cook, J.C., Cruikshank, D.P., Custodio, O.S., Dalle Ore, C.M., Deboy, C.,
Dischner, Z.J.B., Dumont, P., Earle, A.M., Elliott, H.A., Ercol, J., Ernst, C.M., Finley, T., Flanigan, S.H., Fountain, G., Freeze, M.J., Greathouse, T., Green, J.L.,
Guo, Y., Hahn, M., Hamilton, D.P., Hamilton, S.A., Hanley, J., Harch, A.,
Hart, H.M., Hersman, C.B., Hill, A., Hill, M.E., Hinson, D.P., Holdridge, M.E.,
Horanyi, M., Howard, A.D., Howett, C.J.A., Jackman, C., Jacobson, R.A., Jennings, D.E., Kammer, J.A., Kang, H.K., Kaufmann, D.E., Kollmann, P., Krimigis, S.M., Kusnierkiewicz, D., Lauer, T.R., Lee, J.E., Lindstrom, K.L., Linscott, I.R.,
Lisse, C.M., Lunsford, A.W., Mallder, V.A., Martin, N., McComas, D.J., McNutt, R.L.,
Mehoke, D., Mehoke, T., Melin, E.D., Mutchler, M., Nelson, D., Nimmo, F.,
Nunez, J.I., Ocampo, A., Owen, W.M., Paetzold, M., Page, B., Parker, A.H.,
Parker, J.W., Pelletier, F., Peterson, J., Pinkine, N., Piquette, M., Porter, S.B., Protopapa, S., Redfern, J., Reitsema, H.J., Reuter, D.C., Roberts, J.H., Robbins, S.J.,
Rogers, G., Rose, D., Runyon, K., Retherford, K.D., Ryschkewitsch, M.G., Schenk, P.,
Schindhelm, E., Sepan, B., Showalter, M.R., Singer, K.N., Soluri, M., Stanbridge, D.,Steﬄ, A.J., Strobel, D.F., Stryk, T., Summers, M.E., Szalay, J.R., Tapley, M., Taylor, A., Taylor, H., Throop, H.B., Tsang, C.C.C., Tyler, G.L., Umurhan, O.M.,
Verbiscer, A.J., Versteeg, M.H., Vincent, M., Webbert, R., Weidner, S., Weigle, G.E., White, O.L., Whittenburg, K., Williams, B.G., Williams, K., Williams, S.,
Woods, W.W., Zangari, A.M., Zirnstein, E., 2015. The Pluto system: initial resultsfrom its exploration by new horizons. Science 350.

Strobel, D.F., Zhu, X., Summers, M.E., Stevens, M.H., 1996. On the vertical thermal
structure of Pluto’s atmosphere. Icarus 120, 266–289.
Toigo, A.D., French, R.G., Gierasch, P.J., Guzewich, S.D., Zhu, X., Richardson, M.I., 2015.
General circulation models of the dynamics of Pluto’s volatile transport on theeve of the new horizons encounter. Icarus 254, 306–323.
Toigo, A.D., Gierasch, P.J., Sicardy, B., Lellouch, E., 2010. Thermal tides on Pluto.
Icarus 208, 402–411.
Van Leer, B., 1977. Towards the ultimate conservative difference scheme : IV. a new
approach to numerical convection. J. Comput. Phys. 23, 276–299.
Vangichith, M., Forget, F., 2011. Climate Models of Triton and Pluto. New Horizons
Workshop on Icy Surface Processes. Lowell Observatory, Flagstaff Arizona. p.
1447
Yelle, R.V., Lunine, J.I., 1989. Evidence for a molecule heavier than methane in theatmosphere of Pluto. Nature 339, 288–290.
Young, E.F., French, R.G., Young, L.A., Ruhland, C.R., Buie, M.W., Olkin, C.B.,
Regester, J., Shoemaker, K., Blow, G., Broughton, J., Christie, G., Gault, D., Lade, B.,
Natusch, T., 2008. Vertical structure in Pluto’s atmosphere from the 2006 June12, stellar occultation. Astron. Jour. 136, 1757–1769.
Young, L.A., 2013. Pluto’s seasons: new predictions for new horizons. Astrophys. J.
766, L22.
Young, L.A., Elliot, J.L., Tokunaga, A., de Bergh, C., Owen, T., 1997. Detection of
gaseous methane on Pluto. Icarus 127, 258.
Zalucha, A.M., 2016. An atmospheric general circulation model for Pluto with predictions for new horizons temperature proﬁles. Mon. Not. R. Astron. Soc. 459,
902–923.
Zalucha, A.M., Gulbis, A.A.S., Zhu, X., Strobel, D.F., Elliot, J.L., 2011. An analysis of
Pluto occultation light curves using an atmospheric radiative-conductive model.
Icarus 211, 804–818.
Zalucha, A.M., Michaels, T.I., 2013. A 3d general circulation model for Pluto and Triton with ﬁxed volatile abundance and simpliﬁed surface forcing. Icarus 223,
819–831.

Please cite this article as: F. Forget et al., A post-new horizons global climate model of Pluto including the N2 , CH4 and CO cycles, Icarus(2017), http://dx.doi.org/10.1016/j.icarus.2016.11.038