Icarus 309 (2018) 277–296

Contents lists available at ScienceDirect

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

The nitrogen cycles on Pluto over seasonal andastronomical timescales
T. Bertrand a,b,∗, F. Forget a, O.M. Umurhan b, W.M. Grundy c, B. Schmitt d, S. Protopapa e,f,
A.M. Zangari f, O.L. White b, P.M. Schenk g, K.N. Singer f, A. Stern f, H.A. Weaver h, L.A. Young f,
K. Ennico b, C.B. Olkin faLaboratoire de Métórologie Dynamique, IPSL, Sorbonne Universits, UPMC Univ Paris 06, CNRS, 4 place Jussieu, Paris 75005, France
National Aeronautics and Space Administration (NASA), Ames Research Center, Space Science Division, Moffett Field, CA 94035, United States
Lowell Observatory, Flagstaff, AZ, United StatesdUniversit Grenoble Alpes, CNRS, Institut de Planétologie et Astrophysique de Grenoble, Grenoble F-38000, FranceeUniversity of Maryland, Department of Astronomy, College Park, MD 20742, United StatesfSouthwest Research Institute, Boulder, CO 80302, United StatesgLunar and Planetary Institute, 3600 Bay Area Blvd. Houston, TX 77058, United StateshJohns Hopkins University Applied Physics Laboratory, Laurel, MD 20723, United Statesbca r t i c l e

i n f oArticle history:
Received 2 October 2017
Revised 7 March 2018
Accepted 9 March 2018
Available online 15 March 2018
Keywords:
Pluto
Nitrogen
Paleo
Modeling
GCM
Sputnik Planitia

a b s t r a c t
Pluto’s landscape is shaped by the endless condensation and sublimation cycles of the volatile ices covering its surface. In particular, the Sputnik Planitia ice sheet, which is thought to be the main reservoirof nitrogen ice, displays a large diversity of terrains, with bright and dark plains, small pits and troughs,
topographic depressions and evidences of recent and past glacial ﬂows. Outside Sputnik Planitia, New
Horizons also revealed numerous nitrogen ice deposits, in the eastern side of Tombaugh Regio and at
mid-northern latitudes.
These observations suggest a complex history involving volatile and glacial processes occurring on
different timescales. We present numerical simulations of volatile transport on Pluto performed with amodel designed to simulate the nitrogen cycle over millions of years, taking into account the changesof obliquity, solar longitude of perihelion and eccentricity as experienced by Pluto. Using this model, we
ﬁrst explore how the volatile and glacial activity of nitrogen within Sputnik Planitia has been impacted
by the diurnal, seasonal and astronomical cycles of Pluto. Results show that the obliquity dominates theN2 cycle and that over one obliquity cycle, the latitudes of Sputnik Planitia between 25°S-30°N are dominated by N2 condensation, while the northern regions between 30°N and -50°N are dominated by N2
sublimation. We ﬁnd that a net amount of 1 km of ice has sublimed at the northern edge of Sputnik
Planitia during the last 2 millions of years. It must have been compensated by a viscous ﬂow of the thickice sheet. By comparing these results with the observed geology of Sputnik Planitia, we can relate theformation of the small pits and the brightness of the ice at the center of Sputnik Planitia to the sublimation and condensation of ice occurring at the annual timescale, while the glacial ﬂows at its easternedge and the erosion of the water ice mountains all around the ice sheet are instead related to the astronomical timescale. We also perform simulations including a glacial ﬂow scheme which shows that the
Sputnik Planitia ice sheet is currently at its minimum extent at the northern and southern edges. We
also explore the stability of N2 ice deposits outside the latitudes and longitudes of the Sputnik Planitia
basin. Results show that N2 ice is not stable at the poles but rather in the equatorial regions, in particular in depressions, where thick deposits may persist over tens of millions of years, before being trappedin Sputnik Planitia. Finally, another key result is that the minimum and maximum surface pressures obtained over the simulated millions of years remain in the range of milli-Pascals and Pascals, respectively.
This suggests that Pluto never encountered conditions allowing liquid nitrogen to ﬂow directly on its surface. Instead, we suggest that the numerous geomorphological evidences of past liquid ﬂow observed on∗
Corresponding author at: Laboratoire de Métórologie Dynamique, IPSL, Sorbonne Universits, UPMC Univ Paris 06, CNRS, 4 place Jussieu, Paris 75005, France.
E-mail address: tanguy.bertrand@lmd.jussieu.fr (T. Bertrand).

https://doi.org/10.1016/j.icarus.2018.03.012
0019-1035/© 2018 Elsevier Inc. All rights reserved.

278

T. Bertrand et al. / Icarus 309 (2018) 277–296

Pluto’s surface are the result of liquid nitrogen that ﬂowed at the base of thick ancient nitrogen glaciers,
which have since disappeared.
© 2018 Elsevier Inc. All rights reserved.

1. Introduction
1.1. Pluto’s ices observations
Among the most striking observations of Pluto made by New
Horizons in July 2015 is the prominent nitrogen ice sheet laying in
Sputnik Planitia1 (SP), which displays a highly diverse range of terrains, as described by White et al. (2017), Moore et al. (2017) and
McKinnon et al. (2016). First, bright nitrogen-rich plains (0°-30°N)
contrast with darker plains at higher latitudes (30°-40°N) having higher amounts of diluted methane, and even darker, moremethane and tholins rich plains (40°-50°N) at the northern edge
of SP (see composition maps Fig. 5.C in Protopapa et al. (2017) and
Figs. 13.2, 15 and 18 in Schmitt et al. (2017). Cellular patterns,
indicative of active solid-state convection (McKinnon et al., 2016;
Trowbridge et al., 2016), are observed in the northern part of SP(0°-40°N) but not in the southern part of SP (0°-25°S). The absence of convection cells coincides with the presence of hundred
meters deep pits on the surface of the ice sheet (Moore et al.,2017). Glacial ﬂow activity is observed through the valleys at
the eastern edge of SP (ﬂowing toward the center of SP fromthe uplands, 20°S-30°N) and at the northern edge of SP (ﬂowing outward the basin), as shown by Howard et al. (2017) andUmurhan et al. (2017). Finally, rugged water ice mountains surround the SP region (Al-Idrisi, Hillary and Tenzing Montes), suggesting that they have been eroded and shaped over time by theaction of glacial ﬂow (Stern et al., 2015; Howard et al., 2017), as
the terrestrial “Nunataks” in Greenland and Antarctica.
New Horizons also revealed bright deposits of nitrogen outside
SP. In the equatorial regions, the eastern side of Tombaugh Regio
is covered by numerous patches of nitrogen ice, which are also
observed further east in the depressions and valleys between thehigh altitude “bladed terrains” deposits (Moore et al., 2017) andfurther west at the bottom of deep craters (Schmitt et al., 2017). In
addition, bright patches of nitrogen ice are detected over a latitudinal band between 30°N and 60°N (Schmitt et al., 2017).
What is the history of Sputnik Planitia and the nitrogen deposits? Resurfacing by glacial ﬂow, solid-state convection, or nitrogen sublimation and condensation have been proposed to explain the formation and disappearance of the pits and polygonalcells within SP (White et al., 2017; Moore et al., 2017). In addition,
sublimation-condensation processes are thought to drive the difference in ice albedos, composition and distribution outside Sputnik
Planitia (White et al., 2017; Protopapa et al., 2017; Howard et al.,2017). However, the timescales and amounts of ice involved are notknown, which prevents us from distinguishing the roles of eachprocess and the nature of the reservoirs (perennial/seasonal).
1.2. Modeling of the present-day Pluto’s volatile cycle
The condensation and sublimation of Pluto’s volatiles and theirtransport over several Pluto seasons has been modeled under current orbital conditions by Bertrand and Forget (2016), that is withan obliquity of 119.6°, an orbit eccentricity of 0.2488 and a solar
longitude at perihelion of 3.8°. Solar longitude at perihelion (Ls peri )is the Pluto-Sun angle at perihelion, measured from the Northern
1
The place names mentioned in this paper include a mix of oﬃcially approved
names and informal names.

Hemisphere spring equinox (Ls peri =0° when the perihelion coincides with the northern spring equinox). It deﬁnes the link between Pluto season and the distance from the Sun (and thus also
the duration of the season). Results using a seasonal thermal inertia for the sub-surface within 50 0-150 0 J s−1/2 m−2 K−1 showedthat nitrogen and carbon monoxide (CO) are trapped inside theSputnik Planitia basin because its low altitude (its surface is located at ∼ 3 km below the surroundings terrains) induces a higher
condensation temperature of the ice (compared to the ice outside
the basin), leading to an enhanced thermal infrared cooling andcondensation rate at the bottom of the basin (Bertrand and Forget, 2016). Note that the reverse occurs on Earth, where the ice
caps form at the poles and at high altitude. In these simulations,methane also accumulates in SP but its low volatility allows it tocondense on warmer surfaces (where N2 and CO would instantly
sublimate) and form seasonal frosts of pure methane everywherein the fall-winter hemisphere except at the equator that tends toremain free of volatile ice. Results also showed that bright methane
frosts in the northern hemisphere could favour nitrogen condensation on it and thus lead to the formation of a seasonal nitrogen
polar cap. These polar deposits sublime in spring from the poleand in 2015, only a latitudinal band of nitrogen frost around 45°N
remained.

1.3. The changes of obliquity and orbital parameters overastronomical timescales
While Bertrand and Forget (2016) focused on the volatile cycles
in Pluto’s current orbital conditions, these latter are known to varyover timescales of 100 000 terrestrial years. Pluto’s high obliquity
varies between 104° and 127° over 2.8 Myrs. The solar longitudeof perihelion of Pluto’s orbit varies with a precession period of
3.7 Myrs, while its eccentricity oscillates between 0.222 and 0.266
with a 3.95 Myrs period (Earle et al., 2017; Binzel et al., 2017;
Hamilton et al., 2016). The variation of these parameters, known
as the Milankovitch mechanism, also occurs on the Earth, Mars
and Titan. The parameters combine to modulate the solar insolation and surface temperatures, forcing the volatile ices that form
glaciers, frost or lakes to migrate in different regions with time.
On Earth, obliquity changes are known to have played a criticalrole in pacing glacial and interglacial eras. Because of the absenceof a big moon and its proximity to Jupiter, Mars has an obliquity
that varies much more strongly than on Earth, and experienced periods when poles were warmer than equator, like on Pluto. Many
surface structures of Mars are thought to be the effect of orbital forcing or of the Milankovitch cycle on the climate of Mars.
For instance, climate modeling efforts showed that during highobliquity periods, ice can be deposited almost anywhere in the
mid-latitudes, explaining the evidences of glaciers and widespreadground-ice mantle in these regions, while during low obliquity ice
is transported back to the poles (Levrard et al., 2007; Mischna
et al., 2003; Forget et al., 2006; Madeleine et al., 2009, e.g.). On
Titan, the cooler summers in the north pole explain the lakes preference for the northern latitudes (Schneider et al., 2012). The differing solar insolation between both hemisphere would result fromthe eccentricity of Titans orbit and the obliquity of Saturn, coupled with Titans low inclination and obliquity (Aharonson et al.,2009). Climatic changes similar in scale to Earths climatic cycles

T. Bertrand et al. / Icarus 309 (2018) 277–296

are expected as the obliquity and orbital parameters of Titan varyon timescales of tens of thousands of years.
By comparison with the Earth, Mars and Titan, Pluto’s climate
is expected to be dictated by the universal Milankovitch mechanism as well. A few studies have explored the variation of insolation on Pluto caused by the changes of obliquity, solar longitudeof perihelion and eccentricity, and have shown that the obliquity
is the main driver of Pluto’s insolation (Earle et al., 2017; Binzel
et al., 2017; Stern et al., 2017). However these thermal models neglected the impact of seasonal thermal inertia (TI), which strongly
controls surface temperatures (see Section 3) and they did not address the transport of volatiles, necessary to fully understand how
Pluto’s ices evolved in the past.
1.4. Objectives of this paper
Our objective is to investigate the evolution and distribution of
nitrogen ice on Pluto over the past millions of Earth years (Myrs).
To do that, we extend Bertrand and Forget (2016) study and use
the Pluto’s volatile transport model, taking into account (1) thechanges of the astronomical cycles (obliquity, solar longitude ofperihelion and eccentricity) induced by the perturbation of the Sunon the Pluto-Charon binary system (Dobrovolskis et al., 1997), (2)realistic reservoirs of nitrogen ice, (3) the changes of ice thickness
induced by glacial ﬂow.
Pluto’s astronomical cycles are thought to have been stable
over at least the last 20 million years Binzel et al. (2017), andprobably even before. The ﬁrst reason is that Pluto’s orbit is ina relatively isolated region of the Solar System, never getting
within ∼ 11 AU of any major planet, and it is therefore subject
to very little perturbations (Dobrovolskis et al., 1997). In addition,
Binzel et al. (2017) state that the presence of ancient craters at theequator demonstrates a certain stability of the astronomical cycles,
which could extend back in time by hundreds of Myrs (otherwise
the craters would have been eroded away or completely buried). In
this paper, we assume that the astronomical cycles remained stable during the last 30 Myrs.
Here, we explore the impact of orbital and obliquity changes
on the nitrogen cycle only. Thus, all the simulations of this paper are performed without the cycles of methane and CO. This
choice is driven by the fact that the number of sensitivity parameters and initial states explored in this work already makes it startat a certain complexity level. In fact, we know from Bertrand andForget (2016) that CO “follows“ N2 and remains always trapped in
N2 ice. The case of CH4 ice is much more complex because it canform CH4 -rich deposits (with 3-5% N2 ) on Pluto’s surface and trigger N2 condensation if its albedo is high enough. Exploring thecycle of methane over astronomical timescales and its impact on
Pluto’s climate will be the topic of a separate paper.
In Section 2, we describe the Pluto volatile transport modeland its recent development allowing for the simulations of thenitrogen cycle over several astronomical cycles (in particular, themodel includes an ice redistribution algorithm, glacial ﬂow modeling, changes of topography, obliquity, and orbital parameters withtime).
In Section 3, we show the impact of the obliquity, the orbitalparameters and the thermal inertia on the surface temperatures
averaged over the past Myrs.
Then, in Section 4, we investigate how the past cycles of nitrogen sublimation and condensation (at diurnal, seasonal, astronomical timescales), as well as the glacial ﬂow of N2 ice may haveshaped Sputnik Planitia as it is observed today.
In Section 5, we explore possible steady state for nitrogen deposits outside Sputnik Planitia, by performing simulations over thelast 30 Myrs with different reservoirs of N2 ice initially placed at
the poles, at the equator or uniformly spread over the surface.

279

We also explore the maximum and minimum surface pressures
Pluto encountered during that time. We discuss these results in
Section 6.
2. Model description
We use the latest version of the Pluto volatile transport modelof the Laboratoire de Météorologie Dynamique (LMD) (Bertrand
and Forget, 2016, see Methods).
This model represents the physical processes that control thecondensation and sublimation of Pluto’s volatiles (insolation, suface thermal balance) and uses a simple global mixing function toparametrize the atmopsheric transport and dynamics. Note that inthis model, the atmosphere is considered transparent: there is no
atmospheric process taken into account aside from the condensation, sublimation and exchanges of latent heat with the surface.
Such a model works well on Pluto because the surface energy balance is not signiﬁcantly sensitive to the atmospheric sensible heatﬂux and to the radiative transfer through the air.
In this section, we describe the grid and surface properties usedin our simulations (Sections 2.1 and 2.2) and the recent developments performed in the code, allowing the simulation of the nitrogen cycles over astronomical timescales. These improvements concern a paleoclimate mode (Section 2.3), the implementation of thelatest topography of Pluto with a speciﬁc relief for the SP basin(Section 2.4) and a glacial ﬂow scheme (Section 2.5). Note that all
ﬁgures and maps are shown using the IAU convention, spin northsystem for deﬁnition of the north pole (Buie et al., 1997; Zangari,2015), that is with spring-summer in the northern hemisphere during the 21th Century.
2.1. Model grid
In this paper, the simulations investigating the stability of N2
ice deposits outside SP (Section 5) have been performed with ahorizontal grid of 32 × 24 points, that is a grid-point spacing of
7.5° latitude by 11.25° longitude and a spatial resolution of about
150 km. Simulations focusing on the N2 cycle within Sputnik Planitia (Section 4) have been performed with a twice higher spatial
resolution of 5.6° in longitude and 1.875° in latitude.
2.2. Model grid and surface properties
As in Bertrand and Forget (2016), the reference albedo of nitrogen ice is set to 0.7, and its emissivity to 0.8. The albedo andemissivity of the bare ground (free of N2 ice) are set to 0.1 and
1, respectively. In the soil model, the near-surface layers have a
low-TI to capture the short-period diurnal thermal waves, while
the deeper layers have a high-TI to capture the much longer seasonal thermal waves. The diurnal TI is set to 20 J s−1/2 m−2 K−1 (orSI), as inferred from Spitzer thermal observations (Lellouch et al.,2011). In this paper, simulations have been performed without adiurnal cycle (the insolation is averaged over the Pluto day) andtherefore the diurnal TI has no impact on the results. The reference seasonal thermal inertia TI is uniformly set to 800 SI, because
it corresponds to the best case simulation reproducing the threefold increase of surface pressure observed between 1988 and 2015
(Bertrand and Forget, 2016) and the 1-1.2 Pa value in 2015. Here
we assume that all terrains (water ice bedrock and nitrogen ice)have the same TI. In this paper, sensitivity simulations have also
been performed using 40 0 and 120 0 SI. The modeled diurnal andannual skin depths are 0.008 m and 20–60 m respectively.
To adequately resolve these scale lengths, the subsurface is divided into 24 discrete layers, with a geometrically stretched distribution of layers with higher resolution near the surface to capture
the short waves (the depth of the ﬁrst layer is z1 =1.4×10−4 m) and280

T. Bertrand et al. / Icarus 309 (2018) 277–296

a coarser grid for deeper layers and long waves (the deepest layer
depth is near 10 0 0 m). The depth of each layer is given by:

zk = z1 2k−1

(1)

Our simulations are performed assuming no internal heat ﬂux.
Adding an internal heat ﬂux of few mW m−2 , as suggested in
Robuchon and Nimmo (2011) for Pluto, does not change signiﬁcantly the results. Our tests show that the surface temperature increases by 0.2 K when taking into account an internal heat ﬂux of
3 mW m−2 (see the discussion in 2.5).
Note that the seasonal thermal inertia of N2 ice does not impact
the amount of condensed or sublimed N2 ice to ﬁrst order.
Indeed, as shown by equation 12 and 13 in Forget et al. (2017),
the variation of the exchanged mass between the surface and theatmosphere δ m0 is nearly proportional to a product involving thesurface heat capacity cs (in J m−2 K−1 ), which depends on the thermal inertia:

δ m0 ∝

cs
Ts
LN
2

with LN2(2)

the latent heat of N2 (2.5 105 J kg−1 ) and Ts the sur-face temperature. Yet Ts is nearly proportional to cFs to ﬁrst order,

with F the thermal ﬂux absorbed by the surface (W m−2 ). Hence:

δ m0 ∝F
LN

Pluto year, and the orbital parameters and the obliquity of Pluto
are changed according to the new date t+t.
The paleo-timestep must be small enough so that the changesof obliquity and orbital parameters allow the surface ice distribution to reach a steady state at each paleo timestep, but it mustbe high enough to reduce signiﬁcantly the computing time of thesimulation. In our simulations described here, we typically usedt = 50 0 0 0 Earth years, which corresponds to about 200 Pluto
orbits and a maximal change in its obliquity of 1° of latitude(Binzel et al., 2017).
The changes of obliquity, solar longitude of perihelion
and eccentricity are taken from Earle et al. (2017) andDobrovolskis et al. (1997) and are shown on Fig. 1. Assumingthat the obliquity and orbital parameters remain periodic withtime, we extrapolated the data back to 30 Myrs ago, which corresponds to the starting date of our simulations. Here the variationsof eccentricity (between emin =0.222 and emax =0.266) are taken
into account but not shown for the sake of simplicity. Note thatits impact on Pluto’s climate is negligible compared to the changeof the two other parameters (Hamilton et al., 2016).

(3)2Consequently, to ﬁrst order, δ m0 is independent of cs and ofthermal inertia.
2.3. The paleoclimate mode and ice equilibration algorithm
Because the solar ﬂux received by Pluto’s surface is very low(about 1 W m−2 ), and because Pluto makes a full orbit around theSun every 248 years, the modeled surface and subsurface temperatures and the surface ice distribution require simulations of several Pluto years to reach a steady state. Running the Pluto volatiletransport model with the N2 cycle only can take around 5 min ofcomputing time for one simulated Pluto year at the chosen resolution. Consequently, running climate evolution over 1 Myrs would
require 2 weeks of simulation and performing paleoclimate simulations over several astronomical cycles, e.g. over the last 30 Myrs,
would be prohibitively time-consuming. To resolve this problem,
we implemented in the model a paleoclimate mode containing
an ice iteration scheme. The algorithm is similar to the approachtaken for Mars simulations in Wordsworth et al. (2013): Startingfrom the initial surface ice distribution q− (kg m−2 ), the model isran normally for 5 Pluto years (Step 1), in order to reach repeatable annual condensation and sublimation rates as well as repeatable annual temperatures variations. In the last year, the annual
mean ice rate of change < dq/dt > is evaluated at each grid point
and then multiplied by a multiple-year timestep t to give the updated surface ice distribution q+ and reservoir (Step 2):

q+
= q−
+<
ice
ice

dqice
> t
dt

(4)

If the reservoir of ice at one gridpoint has entirely sublimed
during the paleo-timestep t, the amount of ice is set to 0 kg m−2 ,and after redistribution, the amount in each cell is normalized toconserve the total nitrogen mass (ice and vapour) in the system.
Then, the topography is updated according to the new thicknessof the deposits on the surface (Step 3). Finally, the orbital parameters and the obliquity are changed according to the new date of
the simulation (Step 4). Then the loop started again: the model isrun again for 5 Pluto years ; at each paleo-timestep t, the topography is updated according to the annual mean ice rate of the last

2.4. Topography
The simulations are performed with the latest topography dataof Pluto (Schenk et al., 2018a). In the south hemisphere, where
there is no data, we considered a ﬂat surface (at mean radius). If
we model a surface topography in the southern hemisphere varying of few kilometers, results remain unchanged and N2 ice willaccumulate at similar latitudes.
In the simulations performed with glacial ﬂow of N2 ice, we
modiﬁed the topography of Sputnik Planitia by placing the bedrockmuch deeper than the actual surface of SP, in order to represent SPwith realistic amounts of ice. We assume that the bedrock below
the center of SP is a 10 km deep elliptical basin, which is in therange of the estimates for the thickness of volatiles where polygonal cells are observed (Moore et al., 2016; McKinnon et al., 2016;
Trowbridge et al., 2016; Keane et al., 2016).
The elliptical basin covers the latitudes 10°S-50°N (Fig. 2), witha semi major axis of 1200 km and the foci F=(42°N,163°E) andF’=(1.75°N,177°W). The edges and the southern parts of the basin
are less deep, in accordance with typical impact basin shapes andthe absence of convective cells there. We have different cases
where we placed the bedrock in these areas at 4 km or 5 km
below the mean surface level, as shown by Fig. 2. In the simulations focusing on N2 ice inside SP (Section 4), this modeled basinis then ﬁlled with N2 ice up to 2.5 km below the mean surface
level, which corresponds to the observed altitude of the surface of
the SP ice sheet (in that case, the thickness of ice at the centre of
SP would be 10-2.5=7.5 km). Note that in the modeled topography we removed the water ice mountains from Al-Idrisi to Hillary
Montes for the sake of simplicity.
2.5. Ice viscous ﬂow modeling
Despite being completely frozen, Pluto’s surface is remarkably
active. Evidence of current and past ﬂowing nitrogen ice have been
observed by New Horizons all around Sputnik Planitia (rugged
eroded water ice mountains, glaciers with moraine-like ridges...)and at higher latitudes (Stern et al., 2017; Howard et al., 2017;
Umurhan et al., 2017).
In order to represent the nitrogen ice ﬂow in the model,
we use a laminar glacial ﬂow scheme, as presented in
Umurhan et al. (2017), which is based on the N2 ice rheology for
low surface temperatures described in Yamashita et al. (2010) anddepends on the thickness and temperature of the ice. The model
has several limitations. First, very little laboratory data is available

T. Bertrand et al. / Icarus 309 (2018) 277–296

281

Fig. 1. The astronomical cycles of Pluto during the past 14 Earth million years ( ∼ 5 obliquity cycles): cycles of obliquity (blue solid line, period: 2.8 Myrs) and Ls peri (redsolid line, period: 3.7 Myrs). High obliquity refers to values close to 90°, while ”low“ obliquity is used here to designate the periods with minimal obliquity, that is far
from 90°, although here it remains relatively high compared to other bodies in the Solar System. Present-day Pluto (obliquity=119.6°, that is 60.4° in retrograde rotation) is
subject to an intermediate obliquity (white background). A low obliquity period occurred 2.35 Myrs ago (blue background) and a high obliquity period 0.9 Myrs ago (redbackground). (For interpretation of the references to color in this ﬁgure legend, the reader is referred to the web version of this article.)

Fig. 2. Left: Topography of Sputnik Planitia as seen by New Horizons, ﬁlled by ice at 2.5 km below the mean surface (Schenk et al., 2016a; 2016b). Right: The modeled initial
topography of Sputnik Planitia with a 10 km deep bedrock (here not ﬁlled by nitrogen ice) assumed in the model.

under Plutonian surface conditions and the rheology of solid N2
has not been very well constrained to date (see Umurhan et al.,2017), not to mention the rheology of the possible ice mixtures
on Pluto and of the α crystalline structure. N2 ice behaves as
a viscous ﬂuid with a viscosity ranging between 108 and 1012
Pa s at stresses of 105 Pa (Durham et al., 2010), and ﬂows more
rapidly than water ice on Earth, despite the lower gravity of Pluto.
Secondly, the model strongly depends on the height of the columnof ice, its temperature and the N2 ice rheology properties, which
are poorly known (Umurhan et al., 2017). In addition, we are limited by the horizontal resolution, which prevents us to reproducewith precision small glacial ﬂows (e.g. in the narrow channels
east of SP or around Tenzing Montes in the southern part of SP).

Consequently, the model of glacial ﬂow used in this paper does
not intend to be quantitatively accurate given its simplicity andthe unknowns but instead intends to reproduce to ﬁrst order theactivity of the nitrogen glaciers (e.g. in Sputnik Planitia rapid ﬂows
in the centre of Sputnik Planitia and slow ﬂows on the edges).
We use the scheme described in Umurhan et al. (2017) under
the following assumptions. First, we consider that the ice within
SP ﬂows like pure nitrogen ice. Then, we consider the simple caseof a laminar ﬂow with an isothermal ice, without basal melting(thus we are in the case of a basally cold and dry glacier). Thisis typically valid for thin layers of ice (shallow ice-sheet approximation). In fact, a conductive temperature gradient of 15 K km−1
due to internal heat ﬂux on Pluto is suggested in (McKinnon et al.,

282

T. Bertrand et al. / Icarus 309 (2018) 277–296

2016), assuming no convection. In that case, basal melting would
occur below about 2 km of N2 ice, assuming a surface temperatureof 37 K. This 15 K km−1 value is probably an upper limit becauseof the convection within the ice. (McKinnon et al., 2016) suggestthat the ice layer is convecting in the so-called sluggish lid regime,
which involves the entire layer of ice in the overturn. Therefore
the ice temperature within the layer is likely colder than the temperature assumed without convection. (McKinnon et al., 2016) use
a Nusselt number of 3.2, which means that the effective temperature gradient (in the horizontal mean) is rather 5 K km−1 . In thatcase, basal melting would occur below about 5 km of N2 ice.
If the surface temperature of the ice has approached the melting temperature of 63 K at the triple-point of N2 in the past, assuggested in (Stern et al., 2017), then it is likely that the thin layers of ice at the edges of the ice sheet were “temperate” glaciersat this time. However, our results show that during an entire astronomical cycle (even during high obliquity periods), the surface temperatures of nitrogen ice remain low, below 40 K (seeSection 5.3), which reinforces our assumption of dry glacier at theedges of SP. Here we apply the case of a basally cold and dry
glacier to all encountered ice thicknesses (with no internal heatﬂux and no basal melting). This is acceptable to ﬁrst order because (1) our study focus on the edges of SP and on the glaciers
outside SP, whose thickness in the model remains thin ( < 1 km),(2) the large amount of ice in the centre of SP already ﬂows extremely rapidly; a basal melting here would allow for even more
rapid ﬂow which would not signiﬁcantly impact our results, (3) although the impact of internal heat ﬂux on soil and surface temperature is not negligible (+0.2 K at the surface), it has a small effecton the ﬂow compared to other parameters of the model, which
are not well constrained, such as the albedo and the thickness of
the ice (the glacial ﬂow modeling strongly depends on the depthof the bedrock). We tested the model assuming that the effective
temperature controlling the glacial ﬂow is the one that we would
obtain at the bottom of the glacier taking into account a conductive temperature gradient of 15 K km−1 (that is ∼ 55 K for 1 km
thick glacier), and it does not change signiﬁcantly the results ofthis paper.
Finally, we consider that the bedrock remains static and is notaltered by the glacial ﬂow. Consequently, in the model, the ice istransferred from one grid point to another one using the modiﬁed Arrhenius-Glen analytic function of the mass-ﬂux given in
Umurhan et al. (2017):

q0 = gQ exp[

HHa

1 + HHT

qglen = A(ρ g)n]qglenH n+2tann−1 (θ )nn + 2 (1 + tan2 (θ )) 2

(5)

(6)

With q0 in m2 s−1 , g the gravity at Pluto’s surface(0.6192 m s−1 ), H the ice thickness of the considered column ofice. ρ is the nitrogen ice density, set to 10 0 0 kg m−3 (Scott, 1976;
McKinnon et al., 2016; Umurhan et al., 2017). gQ is a corrective
factor given in Umurhan et al. (2017) and set to 0.5. Ha and HT
are parameters depending on the surface pressure and are given
by equation 14 in Umurhan et al. (2017). θ is the angle between
the two adjacent columns of ice (see Figure 7 in Umurhan et al.,2017), and is deﬁned by θ =arctan(H/Lref ), with H the differenceof altitude between both columns (computed from the bedrock
topography and the amount of ice) and Lref is the characteristic distance of the glacial ﬂow (distance between both columns,
that is both adjacent grid-points). The parameters A and n are
given in Umurhan et al. (2017) are only depend on the surface
temperature. By using this scheme in our model, we obtain thesame relaxation times for the ice than those shown on Figure 8-9
in Umurhan et al. (2017) (the corresponding relaxation time of a

50 km long channel, sloping at θ = 10° and initiated with 200 mof glacial ice is about 50 years). We adapted this scheme so that it
ﬁts on the model grid and so that each grid point can redistributethe correct amount of ice to the neighboring points.
3. Orbital, obliquity and TI changes as drivers of surface
temperatures
Obliquity is the main driver of insolation changes on Pluto. The
polar regions receive more ﬂux than the equator on annual average during high obliquity periods and about the same ﬂux thanthe equator during low obliquity periods. However, although one
could think the rule also applies for surface temperatures, it is not
systematically the case. The following paragraph explains why.
If one assumes that the surface temperature is only driven by
the absorbed ﬂux and the infrared cooling (neglecting the soil TI,
the latent and sensible heat ﬂux), then the surface temperature atequilibrium is given by:( 1 − A )F(7)
σwith A and  the surface albedo and emissivity respectively, F theabsorbed ﬂux and σ the Boltzmann constant. This equation shows
Teq =

4

that when the ﬂux F strongly increases, the surface temperature
only increases by a factor of F1/4 . In other words, the thermal infrared cooling limits the increase of surface temperatures.
As a result, the surface temperatures at the poles do not increase as much as the insolation during summer, and the poles canbe colder than the equator on annual average, even if the mean insolation is not. This is true for low thermal inertia ( < 800 SI), but
not for medium-to-high TI, which enables the poles to store theheat accumulated during summer and release it during winter, as
illustrated on Fig. 3. In the cases of TI between 800 and 1200 SI,
the poles and the equator have similar surface temperatures andthe coldest regions are around ± 30° latitude. The higher the TI,
the more equatorial are the coldest regions.
The variation of Ls peri is also signiﬁcant because of the higheccentricity of Pluto’s orbit. In the past Myrs, the Ls peri parameter has created a North-South asymmetry of annual mean insolation and surface temperatures, favouring a warmer southern hemisphere when Ls peri values were close to 90° and a warmer northernhemisphere when Ls peri values were close to 270° (Fig. 4.A). The
surface temperatures tend to be the same in both hemispheres for
Ls peri values close to 0° and 180°.
The surface temperatures averaged over the last 14 Myrs (whichcorresponds to the last 5 obliquity cycles) and over the last
2.8 Myrs are shown on Fig. 4.B. For the same reasons mentioned above, the equatorial regions are colder than the poles
for medium-to-high TI and warmer for low TI ( < 800 SI). For
medium-to-high TI, the lowest temperatures are obtained at 30°N45°N, which corresponds to the latitudes where a band of nitrogenice has been observed by New Horizons. It can also be noted thatthe northern hemisphere is in average over several Myrs sligthly
colder than the southern hemisphere in all TI cases. This is because the last high obliquity periods of Pluto’s past (during which
the poles receive the most of insolation) remained coupled witha solar longitude at perihelion close to 90°, thus favouring colder
northern latitudes during these periods and in average over several
Myrs.
Assuming that the evolution of obliquity and Ls peri remained
stable during the last billion of years, one can quantify the shift
between the obliquity and the Ls peri values. As shown by Fig. 5,
between 260 and 165 Myrs ago, the Ls peri during high obliquity
periods varied from 225° to 315°, which favoured colder southern latitudes in average over several Myrs. Between 165 and 70
Myrs ago, the Ls peri during high obliquity periods varied from −45°

T. Bertrand et al. / Icarus 309 (2018) 277–296

283

Fig. 3. Surface thermal model results. Left: Annual mean incident solar ﬂux for an obliquity of 104°, 116°, 119.6° and 127°. Generally speaking, the poles receive more ﬂux inaverage than the equatorial regions, except during the low obliquity periods (127°) where mid-latitudes receive less ﬂux in average. Right: Annual mean surface temperatures
obtained with the obliquity of 119.5° and the incident solar ﬂux shown on the left panel, a Ls peri of 0° and TI of 40 0 SI (red), 60 0 SI (purple), 80 0 SI (blue), 10 0 0 SI (orange),
1200 SI (green). The surface albedo is uniformly set to 0.1, and the emissivity to 1. The coldest points are the poles for the low TI case and the “low latitudes bands” at ±30° for the high TI case. (For interpretation of the references to color in this ﬁgure legend, the reader is referred to the web version of this article.)

Fig. 4. Surface thermal model results. Left: Annual mean surface temperatures obtained assuming uniform and constant surface conditions (surface albedo = 0.1, TI=800 SI)
for the maximum and minimum values of obliquities (red: 104°, green:127°) and Ls peri (solid line: 0°, dashed: 90°, dot-dashed: 270°). While the obliquity and the TI drive
the location of the coldest region (polar or equatorial regions), the Ls peri induces an asymmetry of temperatures with Ls=90° and Ls=270° leading to a colder and warmer
north pole respectively. Right: Surface temperatures averaged over the last 14 Myrs (last 5 obliquity cycles, solid lines) and 2.8 Myrs (last obliquity cycle, dashed lines) for
different cases of TI. The surface albedo is uniformly set to 0.1, and the emissivity to 1. The polar regions are warmer than the equatorial regions, except in the case of TI
lower than 800 SI. Our reference simulation (TI=800 SI) shows that the regions around ± 30° are colder in average.

to +45°, leading to symmetric surface temperatures between bothhemisphere in average over several Myrs. Finally, from 70 Myrs up
to now, the Ls peri during high obliquity periods varied from 45°
to 113°, which favoured colder northern latitudes in average over
several Myrs. The entire period of the cycle obliquity+Ls peri is 375
Myrs.

4. Exploring the changes of N2 ice thickness in Sputnik Planitia
In this section, we explore the past evolution of the N2
ice thickness within Sputnik Planitia using the volatile transportmodel in the conﬁguration as described above and with all of
the initial N2 ice reservoir sequestered in the deep Sputnik Planitia basin. We explore the changes of N2 ice thickness considering
its condensation and sublimation cycles, ﬁrst without glacial ﬂow(Section 4.1) and then with glacial ﬂow (Section 4.2). In this pa-

per we assume a compact N2 -rich ice so that 1 kg of ice per m2
corresponds to a thickness of 1 mm.
4.1. The cycles of condensation and sublimation
We ﬁrst start the simulation 30 Myrs ago with nitrogen ice sequestered in SP and let the amount of ice evolve at the surface. In
this section, we do not perform the simulations with glacial ﬂow.
Instead, we assume that the timescale for ice viscous ﬂow is very
short and that the ice sheet surface is effectively a level sheet, remaining ﬂat at all times, and we only evaluate the condensationand sublimation rates within SP. Note that in this simulation, no
nitrogen frost form outside SP.
Fig. 6shows the net change of N2 ice thickness obtained withthe model over four different timescales: A, one Pluto day in July
2015; B, one current Pluto year; C, during the last 50 0,0 0 0 Earthyears, which correspond to the estimated time of full resurfacing

284

T. Bertrand et al. / Icarus 309 (2018) 277–296

Fig. 5. Evolution of the solar longitude of perihelion (Ls peri ) at high obliquity (104°) during the last 10 0 0 Myrs (assuming astronomical cycles stable with time). The Ls peri athigh obliquity has been slowly shifted with time, e.g. from 110.8 to 113.4° during the last 6 Myrs, due to the sligth difference of periods between both Ls peri and obliquity
cycles. During the last 70 Myrs, the Ls peri value at high obliquity remained close to 90° and thus lead to an asymmetry of insolation and surface temperatures favouring asligthly warmer south hemisphere (see Fig. 4).

of SP by the action of the convection cells (McKinnon et al., 2016);
D, during the last 2.8 Myrs (last obliquity cycle). These results arecompared with geologic features observed by New Horizons within
SP (Fig. 6.E-F).
4.1.1. The current annual timescaleFig. 7shows the normalized diurnal variations of N2 ice thickness over a Pluto day in July 2015. Fig. 8 shows the evolution of
the diurnal mean condensation-sublimation rate (net change of N2
ice thickness) within Sputnik Planitia, over one current Pluto year,
since 1800, while Fig. 9 shows the normalized variations of N2 ice
thickness at same dates (gross change of N2 ice thickness).
Over one Pluto day, several tens of micrometres of nitrogen ice
move from the summer (North) to the winter (South) parts of SP(Fig. 6.A, Fig. 7), while in one current Pluto year, a net amount of
20–50 mm of ice accumulates around 30°N, 10–15 mm are lost in
the southern part ( < 10°N) and 20–50 mm are lost in the northernedge of SP ( > 40°N, Fig. 6.B).
In 2015, the regions above 15°N are in a sublimation-dominatedregime, while regions below 15°N are in a condensation-dominatedregime (Fig. 6.A, Fig. 8). The southern regions of SP entered thecondensation-dominated regime after the northern spring equinox
in 1988, where a fast transition of regime between the northernand southern regions occurred. Before 1988, the southern regions
had been in a sublimation-dominated regime since 1865. The net
variation of ice thickness after one Pluto year reaches tens of mm(Fig. 6.B, Fig. 8) but the sublimation-condensation during this period involves thicker layers of ice (by a factor 10-30, Fig. 9). Between 1865 and 1988 (7033 Pluto days), the southern regions lost
0.3-1 m of N2 ice. Between 1988 and 2015, the regions below 10°S
accumulated 0.15–0.25 m of N2 ice.
Association with the observed pits south of Sputnik Planitia Fig. 10 shows the annual mean condensation-sublimation rates
over the last 6 Myrs. The southern latitudes of SP (20°S-10°N) arein a sublimation-dominated regime since at least 100 000 years.
The region below 20°S is a sublimation-dominated regime since
1.3 Myrs and currently starts to enter a condensation-dominatedregime. The net loss of ice involved at the annual timescale inthese regions occurs in the model at the same latitudes where
the small pits are observed, explaining their formation there if
they formed by sublimation. Between 20°S and 10°N, if we assume relatively similar insolation conditions over the last 10 0 0 0 0
Earth years, with a mean net sublimation rate of 15 mm of N2 ice

per Pluto year (Fig. 6.B), then the total loss of ice in this region
could reach 6 m. Below 20°S, assuming a mean sublimation rate of
50 mm of N2 ice per Pluto year (Fig. 10) over 1.3 Myrs, the lossof ice reaches ∼ 260 m. These values are in accordance with theobserved depth of the pits (tens of meters, Moore et al. (2017)).
Other mechanisms not taken into account in this model, such asatmospheric winds, light reﬂection and deposition of dark materials at the bottom of the pits may further increase the sublimation
rate and favor deeper pits formation. This annual mean sublimation pattern could also explain the disappearance of the polygonalcells (if they are erased by sublimation), although this may also be
related to a lower ice thickness (too low to trigger convection), as
it is probably the case for all edges of SP.
Association with the ice albedo and composition Our results also show that the latitudes where N2 ice accumulates inaverage over one Pluto year (Fig. 6.B) correspond to the latitudes where bright N2 plains and a weaker amount of CH4 in
the N2 -CH4 mixture (both are correlated) are observed in SP(Protopapa et al., 2017; Schmitt et al., 2017; Buratti et al., 2017).
Protopapa et al. (2017) note that the abundance of N2 :CH4 ishigher at the center of Sputnik Planitia with respect to the northern area of the basin, contrary to the dilution content of CH4 in
the mixture (Fig. 11). They interpret this trend in the composition
maps as a possible north-south sublimation transport of nitrogenin Sputnik Planitia (indicated schematically by the arrow in panel
B of Fig. 11). This is now supported by our results showing a net
deposition of N2 ice in the middle of the basin over the seasonal
timescales (Fig. 6.B), and a recent deposition of few micrometersof N2 ice in the southern latitudes during the past 30 Earth years(Fig. 8).
Despite its net daily sublimation since about 30 Earth years(Fig. 8), the darker cellular plains are currently an area of net ice
accumulation on the annual timescale (about 4–8 m in the last
10 0 0 0 0 Earth years, Figs 6.B and 10), explaining why they also
remain relatively bright compared to the northern dark troughbounding plains (above 40°N), which are subjected to net annual
sublimation since almost 1.8 Myrs (Fig. 10).
4.1.2. The astronomical timescaleFig. 6.D shows that over the last obliquity cycle (2.8 Myrs ago
up to now), up to 300 m of N2 ice accumulated between 20°S and
30°N, while an intense loss of about 1 km of ice occurred at thenorthern edge of the ice sheet between 30°N and 50°N. In addition,

T. Bertrand et al. / Icarus 309 (2018) 277–296

285

Fig. 6. Sublimation-condensation rates of N2 (zoom at Sputnik Planitia) at different timescales (note the order of magnitude differences between the colorbars, from μm tom): (A) During one Pluto day in July 2015. (B) During one Pluto year in current orbital conditions. (C) During the last 0.5 Myrs. (D) During the last 2.8 Myrs (obliquity cycle).
(E) Geological map of the Sputnik Planitia region (a full resolution version can be found in White et al. (2017). (F) New Horizons mosaic of Sputnik Planitia, with recent
glacial activity indicated by the red area. The purple line indicates the extent of the N2 ice sourced for the glaciation, the cyan line indicates the current ice deposition limit,
and the red line indicates the inferred former ice deposition limit. The black arrows indicate the direction of the ﬂow. Originally shown as Fig. 6 in Howard et al. (2017).

at the southern edge of SP (below 20°S), a net loss of 150 m of ice
also occurred. As shown by Fig. 10, N2 sublimation at the northernedge of SP is the most intense during the periods of high obliquity(e.g 0.85 Myrs ago), and still occurs there during a large part ofthe obliquity cycle, for obliquities lower than 119° (higher than 61°
in retrograde rotation). As shown by Fig. 6.C, during the last 0.5
Myrs, the mean accumulation and loss of ice occurred at similar

latitudes than during the last 2.8 Myrs, except between 15°S-20°S
since these latitudes are currently transitioning to a sublimationdominated regime. During the last 0.5 Myrs, the center of SP accumulated up to 100 m of ice while the northern regions lost 200–
300 m of ice, that is one third of what they have lost in averageover the last obliquity cycle (2.8 Myrs). Note that a net loss of ice
continuously occurs at the northern edge of SP since the last 1.8

286

T. Bertrand et al. / Icarus 309 (2018) 277–296

Fig. 7. Variations of N2 ice thickness within Sputnik Planitia, during one Pluto day in July 2015, normalized to 0 at t=0. The data is taken at the longitude 180°, where theice covers the latitudes 30°S–50°N. As shown by Fig. 6, the ﬂux does not vary with longitude within Sputnik Planitia.

Fig. 8. Evolution of the diurnal mean condensation-sublimation rate within Sputnik Planitia (mm per Pluto day), in current orbital conditions, from 1800 to 2050 assumingthat the glacier remains ﬂat.

Fig. 9. Variations of N2 ice thickness within Sputnik Planitia (normalized to 0 at t=1800), in current orbital conditions, from 1800 to 2050 assuming that the glacier remains
ﬂat. The data is taken at the longitude 180°. Although the net budget of ice within one Pluto year varies by tens of mm (Fig. 6.B), the thickness of ice involved during thisyear reaches hundreds of mm.

T. Bertrand et al. / Icarus 309 (2018) 277–296

Fig. 10. Evolution of the annual mean condensation-sublimation rate of N2 ice withtime (mm per Pluto year), assuming that the glacier remains ﬂat (same as Fig. 8.
The vertical solid and dashed lines correspond to the periods of high (104°) andlow (127°) obliquity, respectively.

Myrs. During the same period, ice has been continuously condensing between 20°S and 25°N. We can associate several structures of
SP to the change of N2 ice thickness averaged over this astronomical timescale.
Depressions and outward glacial ﬂows at the northern andsouthern edge of Sputnik Planitia First, the latitudes where 1–2
km deep depressions are observed at the northern and southernboundaries of the ice sheet (see Figure 8 and 17 in Howard et al.,2872017) coincide with the latitudes where intense sublimation of ice
occurred in the last 2.8 Myrs. This loss of ice should tend to becompensated by glacial ﬂow, in line with the outward direction of
the ﬂow observed at these edges, and with the evidences of theparticularly strong erosion of the Al-Idrisi Montes at the northernedge of SP (Howard et al., 2017). Simulations with glacial ﬂow are
explored in Section 4.2.
Recent glacial activity at the eastern side of Sputnik Planitia Secondly, the recent glacial activity of ice ﬂowing westward
through the valleys of the eastern side of SP (pink color on Fig. 6.F,
the observations show that the ice at the eastern side of SP ﬂows
toward the center of SP) occurs at the same latitudes where nitrogen ice continuously accumulated during the last 1.8 Myrs (20°S30°N). We suggest that the accumulation of ice at these latitudes
created a topography gradient at the edge of SP as the thicklayer of ice far from the edge (closer to the center of SP) ﬂowed
more rapidly than the shallow one at the edge. These glacial ﬂows
should be reduced or disappear during the next hundreds of thousand years since these latitudes gradually enter a sublimationdominated regime (Fig. 10).
Why are such glacial ﬂows not observed on the western side ofSP, where the latitudes receive the same insolation? Glacial activity
on the western side of SP has occurred as evidenced by the manyerosional valleys, but these valleys are not ﬁlled with ﬂowing ice
like they are to the east. This could be due to the signiﬁcant difference of geology between the western and eastern side. As shownby Schenk et al. (2018a), a North-South giant fault system passes

Fig. 11. Modeling results from Protopapa et al. (2017) showing the abundance (A) and the path length (B) of the N2 -enriched. Panel C shows the dilution content of CH4 in
N2 .

288

T. Bertrand et al. / Icarus 309 (2018) 277–296

Fig. 12. Variations in elevation within SP at different locations, assuming that the bedrock below SP include a 6–9 km deep elliptical basin and edges at 4 km (solid lines)
and 5 km (dashed lines) below the mean surface level. The red dotted line corresponds to a case with edges at about 3 km below mean level. In that case, the ice has
been entirely sublimed during the last millions of years at Al-Idrisi, as the elevation shown during that time is the bedrock level at this location (ﬂat line at 2.55 km below
the mean surface). The vertical solid and dashed lines correspond to the periods of high (104°) and low (127°) obliquity respectively. (For interpretation of the references tocolor in this ﬁgure legend, the reader is referred to the web version of this article.)under the western edge of the ice sheet, which may explain the
fragmentation of water ice blocks and the presence of deep ridge,
faults, and cliffs observed on the western edge of SP. This topography may prevent the ice from ﬂowing easily through the west sideof SP and form large glacial ﬂows. In addition, the western side of
SP may correspond to a deeper bedrock than on the eastern side,
preventing strong gradients of nitrogen ice thickness and the inward ﬂow observed in the eastern raising valleys (the ice would
locally ﬂow faster on the western edge of SP).
The dark northern plains of Sputnik Planitia The dark andmethane rich aspect of the northern edge of Sputnik Planitia(40°N-50°N) is also consistent with our results, which show that
this area is a sublimation area at all timescales: the current diurnal, the current annual, the last 0.5 Myrs and the last 2.8 Myrs.
Why does this area not have small pits, like in the southern part ofSP? A suggestion is that the area lacks an intake of fresh nitrogenice necessary for the formation of pits. This may be because thearea is subject to a net loss of ice at all timescales and the layer of
ice became too shallow to undergo solid-state convection. In addition or alternatively, the methane rich composition and the size ofgrains of this dark ice may play a restricting role in the formationof pitted plains. Finally, the formation of pits is a process of erosion
by reﬂected light. If the ice albedo is too low, the direct absorptionof solar energy predominates the reﬂection and the ice sublimates
uniformly, inhibiting the pit formation (Moore et al., 2017).
4.2. The astronomical cycles with glacial ﬂow
In this section, we repeat the same simulations as in
Section 4.1 except that we turn on the glacial ﬂow scheme of N2
ice (described in Section 2.5), which enables the ice to ﬂow in themodeled Sputnik Planitia basin (Section 2.4). The basin is initially
ﬁlled with N2 ice up to 2.5 km below the mean surface level. We
assume that the edges of the basin are at about 3 km, 4 km or
5 km below the mean surface level (Fig. 2).
4.2.1. General overview
As a general tendency, our results show that over one obliquity cycle, only small variations of elevation up to 25 m are obtained in the centre of SP at 20°N (SP remains relatively ﬂat where
the ice is thick) while variations of elevations of 20 0–30 0 m are
obtained at the edges of SP (Figs. 12 and 13), in particular at thenorthern and southern edges. The variations are reduced if we assume the bedrock deeper below the ice sheet. These variations ofFig. 13. Maximal variation in elevation of N2 ice over the last obliquity cycle (last
2.8 Myrs), with a bedrock on the edges of SP at 4 km and an initial ﬁlling of SP at
2.5 km.

elevation obtained with our model are consistent with the depressions observed at the northern and southern boundaries of SP andwith the eroded mountain blocks observed west and east of Tenzing Montes, indicative of the presence of large amounts of ice there
in the past. Fig. 12 suggests that the ice sheet was at its maximal
North-South extension 1.5-2 Myrs ago, since the ice level in the
Al-Idrisi region and in the far south of SP was well above the levelof the centre of SP. Conversely, the last million of years coincides
with a period of minimal extension, which is consistent with theice ﬂowing outward from SP at its northern and southern edges(Howard et al., 2017) and the lower MVIC-derived topography observed north of the Al-Idrisi region (Schenk et al., 2018b).

T. Bertrand et al. / Icarus 309 (2018) 277–296

289

Fig. 14. Maps of N2 ice distribution on Pluto (m) for simulations starting 30 Myrs ago with a polar reservoir ( Polar1, Polar4, Polar8, Polar12). Results are shown after
1 Myrs (left panel), 5 Myrs (middle panel) and 30 Myrs (right panel).

4.2.2. Comparisons with observations
The Al-Idrisi Montes As shown by Fig. 13, the entire northern part of SP (above 40°N) is subject to variations of altitudesof 100-280 m, which is consistent with the intense sublimationand condensation of N2 ice occurring over one obliquity cycle atthese latitudes (Fig. 6.C-D). In particular, the latitudes of the AlIdrisi Montes displays one of the largest variations of elevations(up to 280 m), in agreement with the scenario of strong and endless erosion of the water ice blocks in this region (Figs. 12 and 13).
In the simulation with the less deep bedrock on the edges of SP,
the ice at the latitudes of the Al-Idrisi Montes sublimed during thelast 2 Myrs and revealed the bedrock (which is 2.55 km below themean surface at this location in that case, Fig. 12 red dotted line).
Currently, this region enters a regime of net accumulation over astronomical timescales.
The recent glacial activity at the eastern side of Sputnik
Planitia Fig. 12 also shows the evolution of the ice at the easternedge of SP, at 15°N. As predicted by Fig. 10, the area accumulated
up to 150 m of ice during the high obliquity periods, if the bedrockis at 4 km below mean surface (the variation is less if the bedrockis deeper). The elevation of this area decreases since 0.6 Myrs because the ﬂow of ice toward the center of SP overcomes the intake
of nitrogen from condensation. Note that this area is always higher
than the center of SP. Thus, the glacial ﬂows induced from the uplands to the center of SP, as observed by New Horizons, should
never stop.

The southern latitudes of Sputnik Planitia The south of SP
also displays strong variations of elevation, that are about 200 m
over one obliquity cycle, with a bedrock at 4 km below mean surface (Fig. 13). At the east of the Tenzing Montes (27°S-170°E), theelevation of the ice is higher than the center of SP during most of
the obliquity cycle.
5. Possible steady-state conditions for Pluto’s ices
In this section, we explored the stability of N2 ice deposits outside Sputnik Planitia. To do that, we performed several simulations
using different sensitivity parameters and initial states.
5.1. Simulation settings
We used the following settings for the simulations: (1) We performed the simulations over the last 30 Myrs, taking into accountthe obliquity and orbital changes over time described in Section 1.
(2) We used realistic reservoirs of N2 ice corresponding to a global
surface coverage of Pluto of 20 0 m, 50 0 m or 10 0 0 m of ice. The
case of 500 m corresponds to a basin of 120 0x10 0 0 km ﬁlled by∼ 7 km of ice, which is in the range of what is assumed for SP.
(3) At the beginning of these simulations 30 Myrs ago, the surface is not initialized with the entire N2 reservoir trapped inside
SP as in Section 4. Instead, the initial N2 reservoir is either globally
uniformly distributed (Simulations Glob), or placed at the equatorial regions between ± 30° latitude (Simulations Equa), or at the290

T. Bertrand et al. / Icarus 309 (2018) 277–296

Fig. 15. Maps of N2 ice distribution on Pluto (m) for simulations starting 30 Myrs ago with a global reservoir ( Glob1, Glob4, Glob8, Glob12). Results are shown after
1 Myrs (left panel), 5 Myrs (middle panel) and 30 Myrs (right panel).

poles above 50° latitude (Simulations Polar). As an example, an
initial global reservoir of 500 m redistributed over the equatorial
regions between ± 30° latitude corresponds to an initial equatorial reservoir of ∼ 1 km of ice. (4) We used the latest topography
data from New Horizons coupled with a deep bedrock for SP (upto 10 km deep), as described in Section 2.4, and the glacial ﬂow
scheme described in Section 2.5.
The sensitivity parameters of the simulations are the following: (1) Seasonal thermal inertia of 40 0, 80 0 and 1200 SI are used,similar than those used in Bertrand and Forget (2016). The diurnal
thermal inertia remains ﬁxed at 20 SI (Lellouch et al., 2011), as in
Bertrand and Forget (2016). (2) The reference N2 albedo and emissivity used are set to 0.7 and 0.8 respectively, while those for bareground are set to 0.1 and 1 respectively, which is in the range ofwhat has been used in Bertrand and Forget (2016), assuming thewater ice bedrock is covered by dark tholins. We also explored thecase of an albedo of 0.4 for N2 ice.
5.2. Simulation results
The results are summarized in Table 1 and illustrated by
Figs. 14 and 15.
5.2.1. Overall outcome
As a general rule, N2 ice quickly accumulates in the Sputnik
Planitia basin and in the equatorial regions (preferentially at latitudes around ± 30°), with stronger condensation rates inside SP

and inside other depressions because of the stronger infrared cooling effect, as detailed in Section 1 and Bertrand and Forget (2016).
However, in many of the simulations, large deposits also remainin the equatorial regions outside SP after 30 Myrs, and even beyond as they seem to remain relatively stable with time. The simulation Polar8 described in Table 1 and shown on Fig. 14 provides
a typical example. Starting with an initial global reservoir of 500 m
conﬁned at the poles and a thermal inertia of 1200 SI, the ice migrates toward the equatorial regions by forming latitudinal bands
which get closer to the equator with time. The basin SP is progressively ﬁlled by N2 ice, with a decreasing rate with time, because
the ice outside SP migrates toward the more stable equatorial regions, leading to lower condensation rates inside SP, and also because as the surface of the basin becomes less deep, the infrared
SP indicatescooling effect becomes less eﬃcient. The parameter τ95%
the time needed to ﬁll SP at 95% of its ﬁnal state (in Myrs). It depends on the TI, the reservoir and the initial state. The lower thistime, the more stable are the deposits outside SP. As an example,
in the simulation Polar8, the basin is already in a relatively stablestate after 11.90 Myrs. After 30 Myrs, it is ﬁlled by ice up to 2350
m below the mean surface (Table 1). Outside SP, N2 ice remainedat the equator forming 60 0-80 0 m deposits. After these 30 Myrs,
the ice still migrates in the basin because of the infrared cooling
effect but at a very slow rate. Typically, in 1 Myrs, the equatorial
deposits lost 5–10 m of ice. Consequently, the 60 0–80 0 m deposits
outside SP should end trapped inside SP after at least 60 additional
Myrs.

T. Bertrand et al. / Icarus 309 (2018) 277–296

291

Table 1
Settings and results of the simulations performed from 30 Myrs ago to present-day. From left to right, settings are: name of the run (the runs marked by ∗ are illustrated by
Figs. 14 and 15), thermal inertia, N2 ice reservoir (globally averaged), N2 ice albedo. Results are: latitudes between which N2 ice deposits are obtained outside SP, latitudes
where the N2 ice deposits are obtained outside SP peak, altitude of the N2 ice deposits inside SP, altitude of N2 ice outside SP, time needed to ﬁll SP with N2 ice at 95% of
its ﬁnal state, minimum and maximum surface pressures obtained during the last 15 Myrs
RunPolar1∗
Polar2
Polar3
Polar4∗
Polar5
Polar6
Polar7
Polar8∗
Polar9
Polar10
Polar11
Polar12∗
Equa1Equa2Equa3Equa4
Equa5
Equa6
Equa7
Equa8Equa9Equa10Equa11Equa12Glob1∗
Glob2
Glob3
Glob4∗
Glob5Glob6
Glob7
Glob8∗Glob9Glob10
Glob11Glob12∗

Thermal inertia

Reservoir

AlbedoLatitude N2

Peak N2

Hice

Hice

(J s−1/2 m−2 K−1 )(kg m−2 )AN2outside SPoutside SP

in SP

outside SP

(m)

(m)

(Myrs)

Pmin

Pmax

400
400
400
800
800
800
1200
1200
1200
400
800
1200
400
400
400
800
800
800
1200
1200
1200
400
800
1200
400
400
400
800
800
800
1200
1200
1200
400
800
1200

200
500
10 0 0200
500
10 0 0200
500
10 0 0
500
500
500
200
500
10 0 0200
500
10 0 0200
500
10 0 0
500
500
500
200
500
10 0 0200
500
10 0 0200
500
10 0 0
500
500
500

0.7
0.7
0.7
0.7
0.7
0.7
0.7
0.7
0.7
0.4
0.4
0.4
0.7
0.7
0.7
0.7
0.7
0.7
0.7
0.7
0.7
0.4
0.4
0.4
0.7
0.7
0.7
0.7
0.7
0.7
0.7
0.7
0.7
0.4
0.4
0.4

15°N-22.5°N
30°S-30°N
50°S-50°N
No ice30°S-30°N
50°S-50°N
No ice30°S-30°N
50°S-50°N
30°S-37.5°N
30°S-37.5°N
30°S-37.5°N
No ice30°S-30°N
50°S-50°N
20°S-20°N
37.5°S-37.5°N
50°S-50°N
20°S-20°N
30°S-30°N
50°S-50°N
30°S-37.5°N
37.5°S-45°N
30°S-37.5°N
No ice30°S-30°N
50°S-50°N
20°S-20°N
30°S-30°N
50°S-50°N
20°S-20°N
30°S-30°N
50°S-50°N
30°S-37.5°N
30°S-37.5°N
30°S-30°N

18°N
15°S
10°S
/
EquatorEquator/EquatorEquator
5°N or locally in the depressions
Equator or locally in the depressions
Equator or locally in the depressions
/
10°S
18°S
10°S 10°N, less ice at the equator
EquatorEquator
10°S 10°N, less ice at the equator
EquatorEquator
5°N or locally in the depressions
Equator or locally in the depressions
Equator or locally in the depressions
/
10°S
18°S
10°S 10°N, less ice at the equator
EquatorEquator
10°S 10°N, less ice at the equator
EquatorEquator
5°N or locally in the depressions
Equator or locally in the depressions
Equator or locally in the depressions

−3960
−1540
440
−3290
−2110
440
−3290
−2350
430
−1740
−1820
−1850
−3290
−1460
450
−5290
−2030
460
−5340
−2160
440
−1740
−1820
−1850
−3290
−1500
470
−5200
−2040
450
−5010
−2250
450
−1740
−1830
−1860

330
40 0–60 0
10 0 0–120 00
500–750
10 0 0–120 00
50 0–80 0
10 0 0–120 0
60 0–80 0
60 0–80 0
70 0–80 00
50 0–70 0
10 0 0–120 0
30 0–40 0
60 0–70 0
10 0 0–120 0
30 0–40 0
70 0–80 0
10 0 0–120 0
60 0–80 0
60 0–80 0
70 0–80 00
500–650
10 0 0–120 0
250–400
650–7501200
200–350
70 0–80 0
1200
60 0–80 0
60 0–80 0
70 0–80 0

22.45
10.15
6.55
26.55
27.05
12.30
35.00
11.90
12.50
7.45
7.50
7.55
19.15
9.75
7.40
7.95
27.85
12.25
7.70
30.35
12.20
9.65
9.40
8.10
17.15
11.45
6.65
6.10
25.05
12.00
7.40
26.10
12.30
8.70
7.90
7.95

0.00560.01660.02340.02290.0371
0.05230.14980.05130.06790.8433
1.82242.59420.00730.01600.02320.02410.03200.05270.03750.04850.07040.84301.82202.59150.00730.01600.02220.02580.03200.05330.03720.0461
0.06690.84291.82192.5917

4.945
3.977
2.371
2.355
2.101
1.527
2.862
1.393
1.162
88.48570.00562.690
4.919
4.192
2.397
2.218
1.975
1.536
1.381
1.448
1.184
88.48370.00662.693
4.754
4.147
2.391
2.190
1.932
1.532
1.373
1.373
1.170
88.48770.01162.722

We note that the ice in the equatorial regions outside SP is
sligthly less stable for Ls peri values close to 90° and 270° (values favouring an asymmetry of surface temperatures between bothhemisphere). Finally, N2 ice is never stable at the poles and any initial polar deposit up to 1 km thick is entirely sublimed after lessthan 2 Myrs.
5.2.2. Sensitivity to the ice reservoir
The larger the ice reservoir, the faster the glacial ﬂow and themore easily the ice reaches the equatorial regions (outside SP),
where it is able to form relatively stable deposits several hundred
meters deep. In addition, larger reservoirs of ice lead to larger N2
ice deposits outside SP, spread from the equator toward higher latitudes.
In our simulations, equatorial deposits outside SP inevitably
form as soon as the initial reservoir is equal to or larger than
500 m, independently of any initial distribution as long as it is assumed that the initial reservoir is outside SP. If the initial reservoir
is lower, the ice does not ﬂow easily and the presence of equatorial deposits outside SP depends on the TI and the initial state (seeSection 5.2.3).
As an example, in the simulations Polar4 and Polar7 performed with a reservoir of 200 m of ice, all the ice ends in the
SP basin after 30 Myrs (see Fig. 14) and ﬁlls it up to 3.29 km below the mean level. If the 200 m of ice are initially present at theequator (Simulation Equa and Glob), and if the TI is higher thanSurface Pressure
SP
τ95%

(Pa)

80 0 SI then 20 0-40 0 m thick ice deposits can persist in the equatorial regions after 30 Myrs (e.g. simulations Glob4, Glob7, Equa4,
Equa7, see Table 1).
While in the simulations performed using a reservoir of 200 m
the SP basin is ﬁlled by ice up to 3.29 km below the mean level
after 30 Myrs, in all simulations using a reservoir of 500 m of ice,
the SP basin is ﬁlled up to 150 0–230 0 m below the mean level. In
all simulations using a reservoir of 1 km of ice, the basin is entirely ﬁlled with ice. In these cases, the basin ﬁlls up very rapidly
because a large amount of ice is able to ﬂow directly inside thebasin.
5.2.3. Sensitivity to thermal inertia and initial state
Our results are sensitive to the assumed thermal inertia. As
shown by Figs. 3 and 4, the lower the TI, the less equatorial are
the cold points on Pluto’s surface (in average over the last Myrs). In
particular, if TI is lower than 800 SI, the equator becomes warmer
in average than the other latitudes. In the simulations using areservoir of 200 m of ice, different ﬁnal states are obtained after
30 Myrs depending on the TI and the initial state.
(1) If the initial reservoir is distributed at the poles ( Polar1,
Polar4, Polar7), we found that the ice subsists outside SP after
30 Myrs only if the TI is equal to or lower than 400 SI, where arelatively stable latitudinal band of 330 m N2 ice can form at 20°N
( Polar1, Fig. 14). If the TI is larger ( Polar4, Polar7), the entire
reservoir is quickly trapped in SP (Fig. 14).

292

T. Bertrand et al. / Icarus 309 (2018) 277–296

Fig. 16. Evolution of maximum (solid lines) and minimum (dashed lines) annual surface pressure over the last 15 Myrs for simulations starting with a polar reservoir:
Polar1 (blue), Polar4 (red), Polar8 (black), Polar12 (green). The present-day maximum surface pressure is ∼ 1.1 Pa (Stern et al., 2015; Gladstone et al., 2016; Hinson et al.,2017). (For interpretation of the references to color in this ﬁgure legend, the reader is referred to the web version of this article.)(2) If some of the initial reservoir is distributed at the equator ( Equa, Glob), the ice can remain outside SP after 30 Myrs
only if the TI is equal to or larger than 800 SI. The ice forms latitudinal bands of 30 0-80 0 m N2 ice at the latitudes ± 10°, withhigher amounts in the local depressions (see Fig. 15, Glob4 andGlob7). If the TI is lower ( Glob1, Equa1), the ice is less stable atthe equator and ends trapped in the SP basin (Fig. 15). Generally
speaking, results obtained using an initial ice reservoir distributedat the equator and over the globe are similar and therefore we
only show the latter results on Fig. 15. Note that when larger TI
are used, larger amounts of ice remain in the equatorial regions.
In addition, the deposits tend to peak at the equator for large TI,
while they tend to peak at higher latitudes ( ± 10°, or ± 20°) for
low TI (see e.g. Polar1, Glob4).
5.2.4. Sensitivity to the albedo
The simulations using a reservoir of 500 m of ice have also
been performed using an albedo of 0.4 for N2 ice ( Polar10-12,
Equa10-12, Glob10-12), instead of 0.7 (reference value). In allthese low albedo simulations, the ice sublimates very rapidly fromthe poles and accumulates in the equatorial regions between 37.5°S
and 37.5°N (the ice is sligthly more spread than in the cases withan albedo of 0.7), forming stable deposits about 60 0–80 0 m deep.
These results are found to be relatively independent of the initialstate and of the thermal inertia. This is because the lower albedo
enables the ice to be warmer and to gain greater mobility (by bothglacial ﬂow and condensation-sublimation ﬂux) to reach the coldest point on Pluto’s surface.
In these simulations, SP is rapidly ﬁlled by ice and reaches a
SP =7-9 Myrs, the elevation of SP surrelatively stable level. After τ95%
face then increases by only a few meter every Myrs due to N2 ice
condensation.
5.3. Minimum and maximum surface pressures
Fig. 16shows the evolution of the maximum and minimum annual surface pressures during the last 15 Myrs obtained in simulations Polar1, Polar4, Polar8 and Polar12. Generally speaking, thesurface pressures (and surface temperatures) remain within 10−2 10 Pa (31-40 K) in all simulations using an albedo for N2 ice of0.7, and within the range 1–100 Pa (39-45 K) in all simulations

using an albedo for N2 ice of 0.4 (in the model we are always in
the global-atmosphere regime). Higher maximum surface pressures
could be obtained, if we lower the albedo below 0.4 or the TI below 400 SI, but such values seem quite distant from reality. We
note that (1) there are two peaks of maximum surface pressure per
obliquity cycle (2) Maximum pressures are lower during the highobliquity periods (104°). This is because the main reservoirs of ice
are located in the equatorial regions, which receive less ﬂux on annual average during these periods (Fig. 3) (3) The surface pressuresin the simulations with N2 ice deposits outside SP are sligthly lessthan in the simulation without deposits outside SP.
6. Discussion
In this paper, we do not seek to reproduce precisely how the SPbasin ﬁlled with N2 ice, since many parameters are unknown (e.g.
the origin of the basin, or the obliquity and the orbital conditions
at the time it formed). In addition we do not take into account thereorientation of the rotation axis (Keane et al., 2016) and, last but
not least, the methane and CO cycles and the presence of methane
and CO ices which can strongly affect the surface albedo, emissivity, temperatures and the rheology of the N2 ice (and its sublimation if CH4 -rich ice layers form on the N2 ice). Instead, we seek toevaluate if N2 deposits outside SP could have remained for a longtime on Pluto and form perennial deposits, and if yes, at which
latitudes.
In our previous paper Bertrand and Forget (2016), we showedthat any condensed N2 ice on Pluto’s surface tends to end inside
the Sputnik Planitia basin. Here, we reproduced similar simulations by taking into account large reservoirs of N2 ice able to sublimate, condense and ﬂow over several Myrs through the changes of
obliquity and orbital parameters of Pluto. We found again that any
large N2 ice deposits outside SP would accumulate in SP and ﬁll
the basin with several kilometres of ice. However, this would takeseveral tens of Myrs during which transient states exist for the deposits. Indeed, assuming that the basin formed initially without N2
ice inside, our results show that large deposits of several hundredsof meters of N2 ice, placed at the poles, are not stable there, andwould inevitably accumulate ﬁrst at mid-latitudes over an entire
latitudinal band after few Myrs, and then, in some cases, in moreequatorial regions after tens of Myrs. We estimate that the basin

T. Bertrand et al. / Icarus 309 (2018) 277–296

would be ﬁlled up by several kilometers of ice in few Myrs. In the
mid-latitude and equatorial regions, the deposits are relatively stable and may remain there during 10-100s of Myrs before to end in
Sputnik Planitia, depending on the thermal inertia, the albedo of
the ice, the local topography, etc. These results raise discussions
about the impact of such glaciers outside SP on the geology of
Pluto and on the surface pressures encountered in Pluto’s past.
First, parts of the equatorial regions of Pluto (and in particular in Cthulhu region) are covered by numerous geologically old
craters which do not seem particularly eroded by ancient deposition of N2 ice. Is it possible that nitrogen ice accumulated in thisregion and did not eroded the bedrock? We believe that cold/dry
based glaciation has a good erosive ability on Pluto and therefore
that it unlikely that hundreds of meters of ice accumulated in thisregion in the past. Although the erosive properties of nitrogen ice
at these temperatures are unknown (nothing has been publishedon this issue yet), we are guided (1) by the erosive mechanismsthat exist on the Earth, where dry/cold based glaciation has beenshown to be possible (Atkins et al., 2002), although it is diﬃcult
to show if it is eﬃcient or not under Pluto’s conditions, and (2)by the fact that the water ice bedrock has been strongly eroded
around Sputnik Planitia, possibly involving dry glaciation (since on
the edges of the ice cap, the ice layer is thin and dry basal ﬂow
should dominate).
In addition, if large N2 deposits existed outside SP, they may nothave been large enough to ﬂow toward the equatorial regions (likein our simulations started with a global reservoir of ice less than
200 m). Or, the equatorial regions may have been already warmer
than the higher latitudes, due to a low thermal inertia (less than800 SI) or due to albedo gradients (dark tholins at the equator andbright methane ice at higher latitudes, see also Earle et al., 2018).
In some of our simulations, relatively stable deposits of N2 ice
are obtained outside SP at higher latitudes, without any ice atthe equator. This is for example the case for simulations with areservoir of 200 m, such as Polar1, Polar4, where the ice does
not ﬂow toward the equator but forms latitudinal bands between
25°S-45°S and 25°N-45°N. Interestingly, several surface features on
Pluto have been interpreted as evidence for past liquid ﬂow, andthey are all observed around the latitudes ± 30-60° (Stern et al.,2017). These latitudes correspond to the regions where ice accumulates in our model (outside SP), ﬁrst as a transient state (thickglaciers over latitudinal bands) and then as a ﬁnal state at midlatitudes around 10–30°N (swallower glaciers because most of theice is trapped inside SP), in particular if a low TI and N2 reservoir
are considered. It has been suggested that epochs with higher atmospheric pressure occurred in Pluto’s geologic past and enabled
the nitrogen ice to be much warmer, perhaps even to turn to liquid, and to ﬂow on the surface leading to the formation of thesefeatures (Stern et al., 2017). However, here our simulations demonstrate that surface pressures higher than 100 Pa are unlikely tohave occurred in Pluto’s past, because the large reservoirs of nitrogen ice are located in the cold and stable equatorial regions andbecause the relatively high thermal inertia and albedo of the ice
limit the sublimation and condensation ﬂuxes. We even found thatpressures are the lowest during the high obliquity periods, because
N2 ice is never stable at the poles at the scale of the astronomical cycles and therefore not available for intense polar sublimation during these periods. Note that in 2015, which corresponds tonorthern spring on Pluto, the observed surface of the north polar
regions of Pluto was free of N2 ice (Schmitt et al., 2017; Grundy
et al., 2016). There will be no N2 ice available at the pole to sublimate during summer and increase the pressure. In addition, in our
simulations, the maximum surface pressures raise up to ∼ 100 Pa
if the ice albedo is set to 0.4, which is a very low value for an ice
as mobile as N2 on Pluto. Therefore we propose that the paleoliquids - and other terrains thought to have been shaped and altered293

by liquid ﬂows - are the results of past liquid ﬂows which occurredat the base of massive nitrogen glaciers (basal ﬂow), which accumulated at the mid-latitudes because they are the coldest points
on Pluto in average (an effect depending on thermal inertia, as
shown by Fig. 4). These glaciers may have remained at these latitudes for millions of years before they disappear, the ice ending
inside SP.
What would trigger the formation of perennial N2 ice deposits
on Pluto outside SP ? Since the astronomical cycles of Pluto are
relatively stable, we can make the hypothesis that the perennial
ice deposits on present-day Pluto reached a steady-state. In thatcase the entire reservoir of N2 ice should be trapped in SP, as suggested by our model results showing that N2 ice deposits outside
SP still accumulates in SP after 30 Myrs, losing about 10 m per
Myrs. However, other processes could help to maintain perennial
N2 ice deposits and feed Pluto’s surface with N2 ice outside SP,
such as cryovolcanism or bright methane deposits enabling N2 ice
to condense on it (see below).
What is the nature (seasonal or perennial) of the different
reservoirs of N2 ice observed in 2015 by New Horizons? How did
they form? Observationally, it is diﬃcult to know because we do
not know the thickness of these reservoirs, although they do notlook like several hundreds of meters deep. The amounts of diluted
CH4 and CO vary in these deposits, which is indicative of volatile
evolution processes (Protopapa et al., 2017; Schmitt et al., 2017).
In Bertrand and Forget (2016), we show that regions covered bydark tholins do not favour N2 condensation on it, while surfaces
covered by bright methane frost do. In fact, the latitudinal band of
nitrogen observed by New Horizons between 30°N-60°N has been
reproduced by the volatile transport model when high methane
albedo ( > 0.65) were considered (see Figure 3 in Bertrand andForget, 2016). In this scenario, the latitudinal band of N2 ice is seasonal since it forms on the cold methane polar frost in winter andsublimates during spring from the pole. However, if the thermal
inertia is lower than the 800 SI assumed in this scenario, then our
results suggest that the ice may be more stable at these latitudes
and the latitudinal band of N2 ice may be perennial, continuously
fed by seasonal frosts. In other words, bright methane frosts may
have helped to maintain the latitudinal bands of massive N2 deposits as a perennial reservoir (e.g. as the one obtained in the case
Polar1 or Polar4 in Fig. 14). Similar arguments apply for the region East of Tombaugh region: bright methane deposits coupledwith relatively low-altitude terrains may favour the accumulationof N2 ice there, which can remain relatively stable over time, especially if the TI is high and thus favouring more stable deposits close
to the equator. In this paper, thin seasonal polar nitrogen frosts
have been obtained in most of the simulations. Although we notedthat a lower thermal inertia favour thicker deposits at the poles,
simulations taking into account bright methane deposits are necessary to fully investigate the evolution of polar frosts, and will be
the topic of future studies.
Finally, as predicted by the model, N2 ice is more stable in the
depressions than in higher terrains. In fact, a limited number ofspots of N2 -rich ice have been observed in the dark equatorial region of Cthulhu, in particular in the Oort and Edgeworth craters(Schmitt et al., 2017). Note that preferential deposition of N2 ice atthe latitudes ± 10° or ± 20° (with the equator free of ice) would
be consistent with our results obtained with TI between 40 0-80 0
SI showing latitudinal bands of stable deposits at these latitudes( Polar1, Equa4, Equa7, Glob4, Glob7). The lack of data makes
it diﬃcult to assess, and low-resolution data in the sub-Charon
hemisphere is currently under processing and analysis. Groundbased telescopic observations rule out the presence of large expanses of N2 ice on the sub-Charon hemisphere, but not the presence of small patches, which are impossible to see from the ground(Grundy et al., 2013).

294

T. Bertrand et al. / Icarus 309 (2018) 277–296

7. Conclusions
The Pluto volatile transport model has been used to investigate
the cycles of nitrogen on Pluto over diurnal, seasonal and astronomical timescales, taking into account the changes of obliquity,
longitude of perihelion and eccentricity and the ﬂow of N2 ice andthe changes of topography induced (following the rheology andglacial ﬂow equations as described in Umurhan et al., 2017).
Our ﬁrst conclusion is that Pluto’s climate is impacted by theuniversal Milankovitch mechanism, as the Earth, Mars and Titan.
The changes of obliquity and orbital parameters lead to differencesof surface temperatures between poles and equator, and asymmetries in the season. We described in this paper how the mostvolatile ice of Pluto, N2 ice, is impacted by these changes overtime.
We ﬁrst focused on the nitrogen cycles within the Sputnik
Planitia basin, considering that it is the only known perennial
reservoir of nitrogen ice on Pluto. The results suggest that Sputnik
Planitia has a complex history, related to sublimation, condensation, and glacial ﬂow involved at different timescales. High obliquity periods induce intense polar summers and thus intense sublimation rates in the northern part of the ice sheet. During the last
2 million years, this part would have lost up to 1 km of ice bysublimation. On the other hand, low obliquity periods favour sublimation in the center of Sputnik Planitia and condensation at thenorth and south extremities, of up to 300 m of ice in 1 million
years. The glacial ﬂow activity (ice ﬂowing toward the center of
Sputnik Planitia) observed at the eastern edge of the ice sheet canthus be related to the intense condensation of nitrogen ice which
occurred at these latitudes during the past 2 million years, while
the methane-enriched N2 ice dark plains are linked to the intense
sublimation which occurred north of Sputnik Planitia during thesame period. The deep pits observed in the south of Sputnik Planitia may have started to form 10 0,0 0 0 years ago, when the southern latitudes of the ice sheet entered a net sublimation-dominatedregime. The bright plains in the center of Sputnik Planitia can be
explained by the current seasonal accumulation of ice there. Finally, the depressions observed north and south of the ice sheet,as well as the strong erosion of the Al-Idrisi Montes, are consistent with the simulated glacial activity of Sputnik Planitia, with
continuous variation of elevations at the edges of the ice sheet
up to 300 m every obliquity cycle. The results also show that incurrent epoch, the ice sheet is close to its minimal extension (inthe model in current epoch the center of SP has a higher elevation
than the northern and southern edges of SP), which is consistent
with the observations showing evidences of strong erosion further
north (Al-Idrisi) and south (West and East of Tenzing Montes) of
the ice sheet.
We also explored the stability of N2 ice deposits outside Sputnik Planitia. Our simulations show that nitrogen ice tends to end
inside Sputnik Planitia but if large deposits are formed outside SP,
they should accumulate and persist in the mid-latitude and equatorial regions for several tens of million years. In particular, N2 ice
accumulates in the depressions. For instance, in most of the simulations involving N2 ice in the equatorial regions, no N2 ice has
been obtained in the Tartarus Dorsa region, featuring the high altitude bladed terrains. The latitudes where N2 ice accumulates depends on the seasonal thermal inertia (the higher it is, the moreequatorial are the deposits), the ice albedo, the initial distribution
and probably other parameters not taken into account in this paper
such as the methane ice distribution. Our simulations support thecase of low to medium thermal inertia (40 0-80 0 SI) for several reasons. It enables to reproduce the evolution of pressure since 1988
(Bertrand and Forget, 2016). In some cases, it enables formation of
perennial deposits at mid-latitudes but not at the equator, which
remains free of volatile ice.

Geomorphological evidences of past liquid ﬂows have been observed at Pluto’s surface at the same mid-latitude. Therefore wesuggest that they formed by liquid nitrogen ﬂows at the base ofancient thick nitrogen glaciers instead of formed by liquid nitrogen
ﬂows directly at Pluto’s surface during higher pressure epochs in
Pluto’s geologic past, as suggested by Stern et al. (2017). This is reinforced by our results showing that the minimum and maximum
surface pressures obtained in our simulations always remain in
the range of milli-Pascals and tens of Pascals, respectively. Therefore surface temperatures never reach the triple point of nitrogen.
It is not possible to reach higher pressures in Pluto’s past with
our model because the sublimation-condensation ﬂux are limitedby the medium to high thermal inertia and the relatively brightalbedo assumed for the N2 ice ( > 0.4).
Finally, the cycle of nitrogen ice on Pluto can be impacted byother processes, not taken into account in the simulations of thispaper. In particular, methane ice is known to play a complex rolesince it can cold trap nitrogen ice if its albedo is high enough(Bertrand and Forget, 2016; Earle et al., 2018), which could explain
why many patches of nitrogen ice are observed outside Sputnik
Planitia. A study taking into account both cycles of methane andnitrogen, over all timescales, is in preparation and should help tobetter understand how these ices evolve on Pluto.
Acknowledgment
We acknowledge the Centre National d’Etudes Spatiales for its
ﬁnancial support through its “Système Solaire” program. This work
was supported by the CNES. It is based on observations with LEISA
and LORRI, embarked on New Horizons. The authors thank thewhole NASA New Horizons instrument and scientiﬁc team for their
excellent work on a fantastic mission and their interest in this research.
References
Aharonson, O., Hayes, A.G., Lunine, J.I., Lorenz, R.D., Allison, M.D., Elachi, C., 2009. An
asymmetric distribution of lakes on Titan as a possible consequence of orbital
forcing. Nat. Geosci. 2, 851–854. doi:10.1038/ngeo698.
Atkins, C.B., Barrett, P.J., Hicock, S.R., 2002. Cold glaciers erode and deposit: evidence
from Allan Hills, Antarctica. Geology 30, 659. doi:10.1130/0091-7613(2002)
0300659:CGEADE2.0.CO;2.
Bertrand, T., Forget, F., 2016. Observed glacier and volatile distribution on Pluto fromatmospheretopography processes. Nature 987. doi:10.1038/nature19337.
Binzel, R.P., Earle, A.M., Buie, M.W., Young, L.A., Stern, S.A., Olkin, C.B., Ennico, K.,
Moore, J.M., Grundy, W., Weaver, H.A., Lisse, C.M., Lauer, T.R., 2017. Climate
zones on Pluto and Charon. Icarus 287, 30–36. doi:10.1016/j.icarus.2016.07.023.
Buie, M.W., Tholen, D.J., Wasserman, L.H., 1997. Separate lightcurves of Pluto andCharon. Icarus 125, 233–244. doi:10.1006/icar.1996.5624.
Buratti, B.J., Hofgartner, J.D., Hicks, M.D., Weaver, H.A., Stern, S.A., Momary, T.,
Mosher, J.A., Beyer, R.A., Verbiscer, A.J., Zangari, A.M., Young, L.A., Lisse, C.M.,
Singer, K., Cheng, A., Grundy, W., Ennico, K., Olkin, C.B., 2017. Global albedos of
Pluto and Charon from LORRI new horizons observations. Icarus 287, 207–217.
doi:10.1016/j.icarus.2016.11.012.
Dobrovolskis, A.R., Peale, S.J., Harris, A.W., 1997. Dynamics of the Pluto-Charon Binary. In: Stern, S.A., Tholen, D.J. (Eds.), Pluto and Charon, University of Arizona
Press, Tucson, 159-190.
Durham, W.B., Prieto-Ballesteros, O., Goldsby, D.L., Kargel, J.S., 2010. Rheological andthermal properties of icy materials. Space Sci. Rev. 153, 273–298. doi:10.1007/
s11214- 009- 9619- 1.
Earle, A.M., Binzel, R.P., Young, L.A., Stern, S.A., Ennico, K., Grundy, W., Olkin, C.B.,
Weaver, H.A., 2017. Long-term surface temperature modeling of Pluto. Icarus
287, 37–46. doi:10.1016/j.icarus.2016.09.036.
Earle, A.M., Binzel, R.P., Young, L.A., Stern, S.A., Ennico, K., Grundy, W., Olkin, C.B.,
Weaver, H.A.New Horizons Surface Composition Theme, 2018. Albedo matters:
understanding runaway albedo variations on Pluto. Icarus 303, 1–9. doi:10.1016/
j.icarus.2017.12.015.
Forget, F., Bertrand, T., Vangvichith, M., Leconte, J., Millour, E., Lellouch, E., 2017. A
post-new horizons global climate model of Pluto including the N2 , CH4 and CO
cycles. Icarus 287, 54–71. doi:10.1016/j.icarus.2016.11.038.
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. doi:10.1126/science.1120335.
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.P., 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.,

T. Bertrand et al. / Icarus 309 (2018) 277–296
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., Bagenal, F., Grundy, W.M., McKinnon, W.B., Moore, J.M., Spencer, J.R., 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., Conard, S., Conrad, C.A., Cook, J.C., Cruikshank, D.P., Custodio, O.S.,
Ore, C.M.D., 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., 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., Holdridge, M.E., Horanyi, M., Howard, A.D.,
Howett, C.J.A., Jackman, C., Jacobson, R.A., Jennings, D.E., Kang, H.K., Kaufmann, D.E., Kollmann, P., Krimigis, S.M., Kusnierkiewicz, D., Lauer, T.R., Lee, J.E.,
Lindstrom, K.L., 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., 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., Ryschkewitsch, M.G., Schenk, P., Sepan, B., Showalter, M.R.,
Soluri, M., Stanbridge, D., Stryk, T., Szalay, J.R., Tapley, M., Taylor, A., Taylor, H.,
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., Zangari, A.M., Zirnstein, E., 2016. The atmosphere of Pluto as observed by New Horizons. Science 351, aad8866. doi:10.1126/science.aad8866.
Grundy, W.M., Binzel, R.P., Buratti, B.J., Cook, J.C., Cruikshank, D.P., Dalle Ore, C.M.,
Earle, A.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, G.E., Young, L.A., 2016. Surface compositions
across Pluto and Charon. Science 351, aad9189. doi:10.1126/science.aad9189.
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. doi:10.1016/j.icarus.2013.01.019.
Hamilton, D.P., Stern, S.A., Moore, J.M., Young, L.A., Binzel, R.P., Buie, M.W., Buratti, B.J., Cheng, A.F., Ennico, K., Grundy, W.M., Linscott, I.R., McKinnon, W.B.,
Olkin, C.B., Reitsema, H.J., Reuter, D.C., Schenk, P., Showalter, M.R., Spencer, J.R.,
Tyler, G.L., Weaver, H.A., 2016. The rapid formation of Sputnik Planitia early in
Pluto’s history. Nature 540, 97–99. doi:10.1038/nature20586.
Hinson, D.P., Linscott, I.R., Young, L.A., Tyler, G.L., Stern, S.A., Beyer, R.A., Bird, M.K.,
Ennico, K., Gladstone, G.R., Olkin, C.B., Pätzold, M., Schenk, P.M., Strobel, D.F.,
Summers, M.E., Weaver, H.A., Woods, W.W., 2017. Radio occultation measurements of Pluto’s neutral atmosphere with New Horizons. Icarus 290, 96–111.
doi:10.1016/j.icarus.2017.02.031.
Howard, A.D., Moore, J.M., Umurhan, O.M., White, O.L., Anderson, R.S., McKinnon, W.B., Spencer, J.R., Schenk, P.M., Beyer, R.A., Stern, S.A., Ennico, K.,
Olkin, C.B., Weaver, H.A., Young, L.A., 2017. Present and past glaciation on Pluto.
Icarus 287, 287–300. doi:10.1016/j.icarus.2016.07.006.
Keane, J.T., Matsuyama, I., Kamata, S., Steckloff, J.K., 2016. Reorientation and faulting of Pluto due to volatile loading within Sputnik Planitia. Nature 540, 90–93.
doi:10.1038/nature20120.
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. doi:10.1016/j.icarus.2011.05.035.
Levrard, B., Forget, F., Montmessin, F., Laskar, J., 2007. Recent formation andevolution of northern Martian polar layered deposits as inferred from a
global climate model. J. Geophys. Res. (Planets) 112 (E11), 6012. doi:10.1029/
20 06JE0 02772.
Madeleine, J.-B., Forget, F., Head, J.W., Levrard, B., Montmessin, F., Millour, E., 2009.
Amazonian northern mid-latitude glaciation on Mars: A proposed climate scenario. Icarus 203, 390–405. doi:10.1016/j.icarus.2009.04.037.
McKinnon, W.B., Nimmo, F., Wong, T., Schenk, P.M., White, O.L., Roberts, J.H.,
Moore, J.M., Spencer, J.R., Howard, A.D., Umurhan, O.M., Stern, S.A., Weaver, H.A.,
Olkin, C.B., Young, L.A., Smith, K.E., Beyer, R., Buie, M., Buratti, B., Cheng, A.,
Cruikshank, D., Dalle Ore, C., Gladstone, R., Grundy, W., Lauer, T., Linscott, I.,
Parker, J., Porter, S., Reitsema, H., Reuter, D., Robbins, S., Showalter, M., Singer, K.,
Strobel, D., Summers, M., Tyler, L., Banks, M., Barnouin, O., Bray, V., Carcich, B.,
Chaikin, A., Chavez, C., Conrad, C., Hamilton, D., Howett, C., Hofgartner, J., Kammer, J., Lisse, C., Marcotte, A., Parker, A., Retherford, K., Saina, M., Runyon, K.,
Schindhelm, E., Stansberry, J., Steﬄ, A., Stryk, T., Throop, H., Tsang, C., Verbiscer, A., Winters, H., Zangari, A.G. a. I. T. T. New Horizons Geology, 2016. Convection in a volatile nitrogen-ice-rich layer drives Pluto’s geological vigour. Nature
534, 82–85. doi:10.1038/nature18289.
Mischna, M.A., Richardson, M.I., Wilson, R.J., McCleese, D.J., 2003. On the orbital
forcing of Martian water and CO2 cycles: a general circulation model study withsimpliﬁed volatile schemes. J. Geophys. Res. (Planets) 108. 16–1
Moore, J.M., Howard, A.D., Umurhan, O.M., White, O.L., Schenk, P.M., Beyer, R.A.,
McKinnon, W.B., Spencer, J.R., Grundy, W.M., Lauer, T.R., Nimmo, F., Young, L.A.,
Stern, S.A., Weaver, H.A., Olkin, C.B., Ennico, K., 2017. Sublimation as a landformshaping process on Pluto. Icarus 287, 320–333. doi:10.1016/j.icarus.2016.08.025.
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.R., Bray, V.J., Chavez, C.L., Howett, C.J.A., Lauer, T.R.,

295

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., Bagenal, F., Gladstone, G.R., Andert, T., Andrews, J., Banks, M.,
Bauer, B., Bauman, J., Barnouin, O.S., Bedini, P., Beisser, K., Bhaskaran, S., Birath, E., Bird, M., Bogan, D.J., Bowman, A., Brozovic, M., Bryan, C., Buckley, M.R.,
Bushman, S.S., Calloway, A., Carcich, B., Conard, S., Conrad, C.A., Cook, J.C., Custodio, O.S., Ore, C.M.D., 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., Jackman, C., Jacobson, R.A., Jennings, D.E.,
Kammer, J.A., Kang, H.K., Kaufmann, D.E., Kollmann, P., Krimigis, S.M., Kusnierkiewicz, D., Lee, J.E., Lindstrom, K.L., 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., Nunez, J.I., Ocampo, A., Owen, W.M., Paetzold, M., Page, B.,
Parker, J.W., Pelletier, F., Peterson, J., Pinkine, N., Piquette, M., Protopapa, S., Redfern, J., Roberts, J.H., Rogers, G., Rose, D., Retherford, K.D., Ryschkewitsch, M.G.,
Schindhelm, E., Sepan, B., Soluri, M., Stanbridge, D., Steﬄ, A.J., Strobel, D.F.,
Summers, M.E., Szalay, J.R., Tapley, M., Taylor, A., Taylor, H., Tyler, G.L., Versteeg, M.H., Vincent, M., Webbert, R., Weidner, S., Weigle, G.E., Whittenburg, K.,
Williams, B.G., Williams, K., Williams, S., Woods, W.W., Zirnstein, E., 2016. The
geology of Pluto and Charon through the eyes of New Horizons. Science 351,
1284–1293. doi:10.1126/science.aad7055.
Protopapa, S., Grundy, W.M., Reuter, D.C., Hamilton, D.P., Dalle Ore, C.M., Cook, J.C.,
Cruikshank, D.P., Schmitt, B., Philippe, S., Quirico, E., Binzel, R.P., Earle, A.M.,
Ennico, K., Howett, C.J.A., Lunsford, A.W., Olkin, C.B., Parker, A., Singer, K.N.,
Stern, A., Verbiscer, A.J., Weaver, H.A., Young, L.A., 2017. Pluto’s global surface composition through pixel-by-pixel Hapke modeling of New Horizons
Ralph/LEISA data. Icarus 287, 218–228. doi:10.1016/j.icarus.2016.11.028.
Robuchon, G., Nimmo, F., 2011. Thermal evolution of Pluto and implications for surface tectonics and a subsurface ocean. Icarus 216, 426–439. doi:10.1016/j.icarus.
2011.08.015.
Schenk, P., Beyer, R. A., McKinnon, W. B., Moore, J. M., Spencer, J. R., White, O. M.,
Singer, K., Robbins, S., Umurhan, O. M., Nimmo, F., Lauer, T. D., Grundy, W., Stern,A., Weaver, H. A. Jr., Young, L. A., Olkin, C., the New Horizons Geology, Geophysics Investigation Team 2018a. Basins, fractures and volcanoes: global cartography and topography of Pluto from New Horizons. Submitted to Icarus.
Schenk, P., Beyer, R.A., Moore, J.M., Spencer, J.R., McKinnon, W.B., Weaver, H.A.,
2016a. Regional topographic properties of Pluto and Charon from New Horizons.
AGU Fall Meeting Abstracts.
Schenk, P., Beyer, R.A., Moore, J.M., Spencer, J.R., McKinnon, W.B., Weaver Jr., H.A.,
Young, L.A., Olkin, C., Ennico Smith, K., Stern, A., 2018b. Topography of Sputnik Planitia Basin on Pluto: What We Know and Don’t Know. AGU Fall Meeting
Abstracts http://adsabs.harvard.edu/abs/2017AGUFM.P24C..08S.
Schenk, P., Singer, K., Robbins, S., Bray, V., Beyer, R., Moore, J., McKinnon, W.B.,
Spencer, J., Runyon, K., Stern, S.A., Young, L.A., Olkin, C., Ennico, K., Weaver, H.A.,
2016b. Topography of Pluto and Charon: Impact Cratering. In: Lunar and Planetary Science Conference. In: Lunar and Planetary Science Conference, 47,
p. 2795.
Schmitt, B., Philippe, S., Grundy, W.M., Reuter, D.C., Cote, R., Quirico, E., Protopapa, S., Young, L.A., Binzel, R.P., Cook, J.C., Cruikshank, D.P., Dalle Ore, C.M.,
Earle, A.M., Ennico, K., Howett, C.J.A., Jennings, D.E., Linscott, I.R., Lunsford, A.W.,
Olkin, C.B., Parker, A.H., Parker, J.W., Singer, K.N., Spencer, J.R., Stansberry, J.A.,
Stern, S.A., Tsang, C.C.C., Verbiscer, A.J., Weaver, H.A., 2017. Physical state anddistribution of materials at the surface of Pluto from New Horizons LEISA imaging spectrometer. Icarus 287, 229–260. doi:10.1016/j.icarus.2016.12.025.
Schneider, T., Graves, S.D.B., Schaller, E.L., Brown, M.E., 2012. Polar methane accumulation and rainstorms on Titan from simulations of the methane cycle. Nature
481, 58–61. doi:10.1038/nature10666.
Scott, T.A., 1976. Solid and liquid nitrogen. Phys. Rep. 27, 89–157. doi:10.1016/
0370- 1573(76)90032- 6.
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.,

296

T. Bertrand et al. / Icarus 309 (2018) 277–296

Woods, W.W., Zangari, A.M., Zirnstein, E., 2015. The Pluto system: initial results from its exploration by New Horizons. Science 350, aad1815. doi:10.1126/
science.aad1815.
Stern, S.A., Binzel, R.P., Earle, A.M., Singer, K.N., Young, L.A., Weaver, H.A., Olkin, C.B.,
Ennico, K., Moore, J.M., McKinnon, W.B., Spencer, J.R., 2017. Past epochs of signiﬁcantly higher pressure atmospheres on Pluto. Icarus 287, 47–53. doi:10.1016/
j.icarus.2016.11.022.
Trowbridge, A.J., Melosh, H.J., Steckloff, J.K., Freed, A.M., 2016. Vigorous convection
as the explanation for Pluto’s polygonal terrain. Nature 534, 79–81. doi:10.1038/
nature18016.
Umurhan, O.M., Howard, A.D., Moore, J.M., Earle, A.M., White, O.L., Schenk, P.M.,
Binzel, R.P., Stern, S.A., Beyer, R.A., Nimmo, F., McKinnon, W.B., Ennico, K.,
Olkin, C.B., Weaver, H.A., Young, L.A., 2017. Modeling glacial ﬂow on and onto
Pluto’s Sputnik Planitia. Icarus 287, 301–319. doi:10.1016/j.icarus.2017.01.017.

White, O.L., 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., Stern, S.A., Ennico, K.,
Olkin, C.B., Weaver, H.A., Young, L.A., Cheng, A.F., Bertrand, T., Binzel, R.P.,
Earle, A.M., Grundy, W.M., Lauer, T.R., Protopapa, S., Robbins, S.J., Schmitt, B.,
2017. Geological mapping of Sputnik Planitia on Pluto. Icarus 287, 261–286.
doi:10.1016/j.icarus.2017.01.011.
Wordsworth, R., Forget, F., Millour, E., Head, J.W., Madeleine, J.-B., Charnay, B., 2013.
Global modelling of the early martian climate under a denser CO2 atmosphere:
water cycle and ice evolution. Icarus 222, 1–19. doi:10.1016/j.icarus.2012.09.036.
Yamashita, Y., Kato, M., Arakawa, M., 2010. Experimental study on the rheologicalproperties of polycrystalline solid nitrogen and methane: Implications for tectonic processes on Triton. Icarus 207, 972–977. doi:10.1016/j.icarus.2009.11.032.
Zangari, A., 2015. A meta-analysis of coordinate systems and bibliography of theiruse on Pluto from Charon’s discovery to the present day. Icarus 246, 93–145.
doi:10.1016/j.icarus.2014.10.040.