F Forget
email: forget@lmd.jussieu.fr

T Bertrand

M Vangvichith

J Leconte

E Millour

E Lellouch

A post-New Horizons Global climate model of Pluto including the N 2 , CH 4 and CO cycles

We have built a new 3D Global Climate Model (GCM) to simulate Pluto as observed by New Horizons in 2015.

All key processes are parametrized on the basis of theoretical equations, including atmospheric dynamics and transport, turbulence, radiative transfer, molecular conduction, as well as phases changes for N 2 , CH 2 and CO. Pluto's climate and ice cycles are found to be very sensitive to model parameters and initial states. Nevertheless, a reference simulation is designed by running a fast, reduced version of the GCM with simplified atmospheric transport for 40,000 Earth years to initialize the surface ice distribution and sub-surface temperatures, from which a 28-Earth-year full GCM simulation is performed. Assuming a topographic depression in a Sputnik-planum (SP)-like crater on the anti-Charon hemisphere, a realistic Pluto is obtained, with most N 2 and CO ices accumulated in the crater, methane frost covering both hemispheres except for the equatorial regions, and a surface pressure near 1.1 Pa in 2015 with an increase between 1988 and 2015, as reported from stellar occultations. Temperature profiles are in qualitative agreement with the observations. In particular, a cold atmospheric layer is obtained in the lowest kilometers above Sputnik Planum, as observed by New Horizons's REX experiment. It is shown to result from the combined effect of the topographic depression and N 2 daytime sublimation. In the reference simulation with surface N 2 ice exclusively present in Sputnik Planum, the global circulation is only forced by radiative heating gradients and remains relatively weak. Surface winds are locally induced by topography slopes and by N 2 condensation and sublimation around Sputnik Planum. However, the circulation can be more intense depending on the exact distribution of surface N 2 frost. This is illustrated in an alternative simulation with N 2 condensing in the South Polar regions and N 2 frost covering latitudes between 35 • N and 48 • N. A global condensation flow is then created, inducing strong surface winds everywhere, a prograde jet in the southern high latitudes, and an equatorial superrotation likely forced by barotropic instabilities in the southern jet. Using realistic parameters, the GCM predict atmospheric concentrations of CO and CH 4 in good agreement with the observations. N 2 and CO do not condense in the atmosphere, but CH 4 ice clouds can form during daytime at low altitude near the regions covered by N 2 ice (assuming that nucleation is efficient enough). This global climate model can be used to study many aspects of the Pluto environment. For instance, organic hazes are included in the GCM and analysed in a companion paper (Bertrand and Forget, Icarus, this issue).

Introduction

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

Among them, Pluto presents a unique case, with an atmosphere significantly warmer than the underlying surface,

A C C E P T E D M

A N U S C R I P T long radiative timescales, and a circulation dominated by condensation/sublimation process of the main atmospheric component. Studying this exotic case can provide new insight into the physics of terrestrial atmosphere.

The observations made by the New Horizons spacecraft have revealed the nature of the surface of Pluto and have provided unprecedented constraints on the state of the atmosphere in 2015 (Stern et al., 2015;[START_REF] Gladstone | The Atmosphere of Pluto as observed by New Horizons[END_REF][START_REF] Grundy | Surface compositions across Pluto and Charon[END_REF][START_REF] Moore | The geology of Pluto and Charon through the eyes of New Horizons[END_REF]. Within that context, it is now interesting to test our ability to create a 3D numerical simulator of the Pluto climate system, analogous to the climate models already used on the Earth as well as on Mars, Venus and Titan. Conversely, the output of such a Global Climate Model is useful to interpret the available atmospheric measurements, and can even shed light on some geological observations.

Pluto's ices and atmosphere observations

The presence of a significant atmosphere on Pluto was demonstrated in 1988 by observing a stellar occultation by Pluto [START_REF] Hubbard | Occultation evidence for an atmosphere on Pluto[END_REF][START_REF] Elliot | Pluto's atmosphere[END_REF]. This atmosphere was predicted to be mainly composed of molecular nitrogen in vapour-pressure equilibrium with N 2 ice deposits observed on the surface. In 2015, New Horizons determined a pressure of about 10 µbar, i.e. 1 Pa [START_REF] Hinson | Radio Occultation Measurements of Pluto's Atmosphere with New Horizons[END_REF][START_REF] Gladstone | The Atmosphere of Pluto as observed by New Horizons[END_REF]. However, a series of stellar occultations conducted since 1988 have shown that the pressure at a specific reference level (e.g. 1275 km from Pluto's center) has increased by a factor of three during that period (Elliot et al., 2003a,b;[START_REF] Olkin | Evidence that Pluto's atmosphere does not collapse from occultations including the 2013 May 04 event[END_REF], suggesting a similar rise of the surface pressure.

Prior to the New Horizons flyby, spectroscopic observations had demonstrated that, in addition to N 2 ice, Pluto's surface is covered by patches of CH 4 and CO ices [START_REF] Owen | Surface ices and the atmospheric composition of Pluto[END_REF][START_REF] Douté | Evidence for Methane Segregation at the Surface of Pluto[END_REF][START_REF] Grundy | Near-infrared spectral monitoring of Pluto's ices: Spatial distribution and secular evolution[END_REF].

New Horizons was able to map these ices and revealed that the main reservoir was a thick ice cap informally named Sputnik Planum, with thinner N 2 frost covering the mid-northen latitudes and CH 4 frost possibly everywhere except in the equatorial dark regions [START_REF] Grundy | Surface compositions across Pluto and Charon[END_REF].

Accordingly, CH 4 and CO gas were observed from Earth to be present in present-day Pluto's atmosphere [START_REF] Young | Detection of Gaseous Methane on Pluto[END_REF][START_REF] Lellouch | Pluto's lower atmosphere structure and methane abundance from high-resolution spectroscopy and stellar occultations[END_REF]Lellouch et al., , 2011a[START_REF] Lellouch | Exploring the spatial, temporal, and vertical distribution of methane in Pluto's atmosphere[END_REF][START_REF] Lellouch | Detection of CO and HCN in Pluto's atmosphere with ALMA[END_REF], with, during the 2008-2012 period, volume mixing ratios near 0.05% for CO and 0.5% for CH 4 . (Lellouch et al., 2011a[START_REF] Lellouch | Exploring the spatial, temporal, and vertical distribution of methane in Pluto's atmosphere[END_REF]. CH 4 has also been oberved by New Horizons Alice spectrograph. in 2015 and estimated to range between 0.6 and 0.84 % in the lower atmosphere [START_REF] Gladstone | The Atmosphere of Pluto as observed by New Horizons[END_REF]. Using the hydrostatic equation, atmospheric temperature profiles have been derived from vertical density profiles retrieved from Earth-based stellar occultations [START_REF] Elliot | Pluto's atmosphere[END_REF](Elliot et al., , 2003b[START_REF] Elliot | Changes in Pluto's Atmosphere: 1988-2006[END_REF][START_REF] Person | Waves in Pluto's Upper Atmosphere[END_REF][START_REF] Young | Vertical Structure in Pluto's Atmosphere from the 2006 June 12 Stellar Occultation[END_REF][START_REF] Dias-Oliveira | Pluto s Atmosphere from Stellar Occultations in 2012 and 2013[END_REF][START_REF] Sicardy | Pluto's atmosphere from the 29 June 2015 ground-based stellar occultation at the time of the New Horizons flyby[END_REF], and from radio and solar occultations performed by New Horizons [START_REF] Hinson | Radio Occultation Measurements of Pluto's Atmosphere with New Horizons[END_REF][START_REF] Gladstone | The Atmosphere of Pluto as observed by New Horizons[END_REF]. The latest stellar occultations and New Horizon's data consistently show that the temperature profile is characterized by a steep temperature gradient in the lower atmosphere, with temperature increasing from surface values (38 to 55 K) at 0 km to about 110 K at 20 km. This

A C C E P T E D M

A N U S C R I P T has been interpreted as resulting from the absorption of near-infrared solar radiation by gaseous methane [START_REF] Yelle | Evidence for a molecule heavier than methane in the atmosphere of Pluto[END_REF][START_REF] Strobel | On the Vertical Thermal Structure of Pluto's Atmosphere[END_REF]. Above 30 km, the temperature appears to decrease with altitude to reach about 70-80 K around 200 km (Dias-Oliveira et al., 2015;[START_REF] Gurwell | Detection of Atmospheric CO on Pluto with ALMA[END_REF][START_REF] Gladstone | The Atmosphere of Pluto as observed by New Horizons[END_REF]. Such a structure requires infrared-cooling species acting only at a specific altitude range. C 2 H 2 and HCN -respectively detected by New Horizons [START_REF] Gladstone | The Atmosphere of Pluto as observed by New Horizons[END_REF] and from the ground [START_REF] Lellouch | Detection of CO and HCN in Pluto's atmosphere with ALMA[END_REF] -have been proposed, but the details of exactly how Plutos upper atmosphere is being cooled remains poorly understood [START_REF] Gladstone | The Atmosphere of Pluto as observed by New Horizons[END_REF].

Finally, stellar occultation observations suggest that the temperature profiles are affected by oscillation that can be related to gravity waves or thermal tides (Elliot et al., 2003b;[START_REF] Person | Waves in Pluto's Upper Atmosphere[END_REF][START_REF] Toigo | Thermal tides on Pluto[END_REF].

3D Modelling of the Pluto surface-atmosphere system

To improve our understanding of the complex Pluto surface-atmosphere system, we have built a new Global Climate Model (GCM) including a full simulation of the nitrogen, methane, and carbon monoxide cycles. This GCM computes the temporal evolution of the variables which control the meteorology and the climate of the planet in different points of a 3D grid covering the entire atmosphere. On the Earth, GCMs have been applied to weather forecasting and climate change projections. Because these models are almost entirely built on physical equations (rather than empirical parameters), several teams around the world have been able to succesfully adapt them to the other terrestrial planets or satellites that have a solid surface and a thick enough atmosphere. The Pluto GCM presented in this paper is derived from the LMD Global Climate Model of planet Mars [START_REF] Forget | Improved general circulation models of the Martian atmosphere from the surface to above 80 km[END_REF]) which has been used for numerous applications including simulating CO 2 ice caps analogous to Pluto's N 2 ice caps [START_REF] Forget | CO 2 snow fall on Mars: Simulation with a general circulation model[END_REF], the thermosphere (Gonzalez-Galindo et al., 2009), photochemistry (Lefevre et al. 2008) or paleoclimatology (e.g. [START_REF] Forget | Formation of Glaciers on Mars by Atmospheric Precipitation at High Obliquity[END_REF]). The LMD GCM has been adapted to Venus [START_REF] Lebonnois | Superrotation of Venus' atmosphere analyzed with a full general circulation model[END_REF] and Titan [START_REF] Hourdin | The sensitivity of the Martian surface pressure to various parameters: A comparison between numerical simulations and Viking observations[END_REF][START_REF] Lebonnois | Titan global climate model: A new 3-dimensional version of the IPSL Titan GCM[END_REF]. All these GCMs have been able to predict or accurately reproduce the observed thermal structure and circulation, giving us some confidence in its ability to predict the characteristics of the Pluto atmosphere in spite of the scarcity of observations. For Pluto, after the simplified General Circulation Model (without phase changes) presented by [START_REF] Zalucha | A 3D general circulation model for Pluto and Triton with fixed volatile abundance and simplified surface forcing[END_REF] for Pluto and Triton, a realistic model was developed by [START_REF] Toigo | General circulation models of the dynamics of Pluto's volatile transport on the eve of the New Horizons encounter[END_REF] a few months before the New Horizons encounter. This model includes a "robust treatment of nitrogen volatile transport", and initializes the full GCM using a two dimensional surface volatile exchange model and a one-dimensional radiative-conductiveconvective model. In this paper we present a new model with a different origin and which benefits from the New Horizons observations. We include an improved N 2 condensation-sublimation scheme, the full CO and CH 4 cycles, and explore the effect of topography. Nevertheless, we use an analogous strategy for the initialization.

A C C E P T E D M A N U S C R I P T

In Sections 2, we provide a detailed description of the different components of our Pluto Global Climate Model, and in Section 3 we discuss how the different model parameters were chosen and how the 3D GCM is initialized for our two baseline simulations. The model results for temperature and winds and for the CH 4 and CO cycles are then presented in Sections 4 and 5, before the conclusion.

Model description

Generalities

As mentioned above, our Pluto Global Climate Model is derived from the LMD Mars GCM [START_REF] Forget | Improved general circulation models of the Martian atmosphere from the surface to above 80 km[END_REF] , with several new parameterizations. Its core is a hydrodynamical code dedicated to the temporal and spatial integration of the equations of hydrodynamics, used to compute the large scale atmospheric motions and the transport. The equations are solved using a finite difference scheme on an "Arakawa C" grid [START_REF] Arakawa | Computational design of the UCLA general circulation model[END_REF]. Such a scheme is equally valid for the Earth, Mars or Pluto. Therefore the hydrodynamical core has not been modified for Pluto. While the estimated surface pressure on Pluto (10 µbar or 1 Pa) is much lower than on the Earth or even on Mars, the atmosphere is thick enough to be modeled with the primitive equations of meteorology used in the model.

In fact, it is generally found that GCMs dynamical cores are valid almost up to the exobase. For instance, on Mars our dynamical core has been used sucessfully up to the thermosphere at pressures lower than 10 -7 Pa [START_REF] González-Galindo | A Ground-to-Exosphere Martian General Circulation Model[END_REF].

In this paper, we present simulations with a horizontal grid of 32×24 points to cover the planet, that is a grid-point spacing of 7.5 • latitude by 11.25 • longitude. The corresponding spatial resolution is of about 150 km, which is equal or better to the typical resolution used in planetary GCMs, and which is sufficient to resolve possible planetary waves.

We also performed simulations with a doubled resolution (64×48) and even an experimental run with a 360×180 grid, and did not find any fundamental differences in the results that could change the conclusions of this paper. Their analysis is out of the scope of this paper and will be presented in a future article, in which we will take into account a more realistic topography. In the vertical, the model uses the terrain-following "sigma" coordinate system in finite difference form (i.e. each layer is defined by a constant value of the ratio pressure devided by surface pressure). 25 levels are typically used. In the baseline model, most of the levels are located in the first 15 km to obtain a good resolution close to the surface, in the boundary layer. The altitude of the first mid-layers are 7 m, 15 m, 25 m, 40 m, 80 m etc.. Above 10 km, the resolution is about one scale height, with the upper pressure level equal to 0.007 times surface pressure, i.e. up to 250 km. (Note that in a companion paper dedicated to the study of atmospheric hazes Bertrand and Forget (2016), the top of the model is extended to about 600 km to include the altitudes of methane photolysis).

A C C E P T E D M

A N U S C R I P T

Radiative transfer

The incident insolation upon each modeled atmospheric column is calculated at each timestep, taking into account the variation of the Pluto-Sun distance throughout its orbit, the seasonal inclination and the diurnal cycle.

While N 2 is the major constituent of the atmosphere of Pluto, its radiative effects are neglected in the lower atmosphere since N 2 is transparent at solar and infrared wavelengths. Nevertheless we account for (1) the radiative heating and cooling by CH 4 , which can vary in space and time depending of the results of the methane cycle model (see section 2.8) (2) cooling by the thermal infrared rotational lines of CO, which volume mixing ratio is prescribed at 0.05 % everywhere (Lellouch et al., 2011a[START_REF] Lellouch | Detection of CO and HCN in Pluto's atmosphere with ALMA[END_REF] and 3) the effect of other infrared emitting species in altitude.

Radiative transfer through CH 4 and CO

For CH 4 and CO, we use a correlated k-distribution radiative transfer model, with 17 spectral bands in the thermal infrared and 23 for solar wavelengths. The bands are designed to well represent the 1.6, 2.3 and 3.3 µm CH 4 vibrational bands in the near infrared as well as the 7.6 µm CH 4 emission band in the thermal infrared. To calculate the k absorption coefficients in each bands, high resolution line-by-line spectra combining CO and CH 4 were computed from the HITRAN 2012 database using the open-source "kspectrum" tool. Spectra and k coefficients were calculated to fill a look up table matrix (from which the k coefficients are interpolated by the GCM in each spectral band) comprising 8 temperatures × 7 log-pressure × 7 CH 4 volume mixing ratio grid, with T = {30, 40, 50, 70, 90, 110, 150, 200} K, p = {10 -4 , 10 -3 , 10 -2 , 10 -1 , 1, 10, 100} Pa, and [CH 4 ]= {10 -4 , 10 -3 , 5 × 10 -3 , 10 -2 , 5 × 10 -2 , 10 -1 , 5 × 10 -1 } kg/kg. We found that no less than 33 points were needed for the g-space integration to get accurate results throughout the matrix space (g is the cumulated distribution function of the absorption data for each band).

Non Local Thermal Equilibrium processes

In the low-pressure, low temperature Pluto environment, a major difficulty (and therefore uncertainty) in the radiative transfer calculations results from the fact that the methane lines can be far from Local Thermal Equilibrium (LTE). It is not the case of CO rotational lines which are assumed to remain in LTE for the pressure levels that we model in this paper.

To account for non-LTE effects for the 7.6 µm CH 4 band, we modify the LTE cooling rates obtained with the correlated k-distribution radiative transfer model as in [START_REF] Strobel | On the Vertical Thermal Structure of Pluto's Atmosphere[END_REF]. However, the total CH 4 cooling rates we obtain are found to be much lower than shown in [START_REF] Strobel | On the Vertical Thermal Structure of Pluto's Atmosphere[END_REF], and significantly smaller than the CO cooling rates. This is also found in recent models from the same authors (D. Strobel, personnal communication) The difference is thought to result from the updated spectroscopic database (HITRAN 2012vs HITRAN 1986) and the fact that the [START_REF] Strobel | On the Vertical Thermal Structure of Pluto's Atmosphere[END_REF] are larger than here. The uncertainties on the NLTE calculations for the 7.6 µm CH 4 band have thus a limited effect on our results.

A C C E P T E D M A N U S C R I P T temperatures used in
For the near-infrared solar bands, we first reproduced the calculations from [START_REF] Strobel | On the Vertical Thermal Structure of Pluto's Atmosphere[END_REF] updated by [START_REF] Zalucha | An analysis of Pluto occultation light curves using an atmospheric radiative-conductive model[END_REF] for each of the 2.3 and 3.3 µm bands. We had no information on the 1.6 µm band. Within that context, and given the overall uncertainty in the NLTE calculations [START_REF] Boursier | Rovibrational relaxation of methane in CH4N2 mixtures: Time-resolved IRIR doubleresonance measurements at 193 K and kinetic modeling[END_REF], we authorized some empirical modifications of the theoretical NLTE variations with atmospheric density (while keeping the theoretical shape) to adjust the heating rates in order to get temperatures closer to the thermal structure observed by New Horizons.

Therefore, the ability of our GCM to roughly reproduce the observed mean thermal structure should not be regarded as a success of our radiative transfer model. In practice, we multiply the total CH 4 heating rate provided by the LTE radiative transfer code by a vertically varying non-LTE efficiency coefficient ε NLTE .

ε NLTE = 0.1 + 0.9 1 + ρ .55 /ρ , (1)
with ρ the atmospheric density (kg m -3 ), and ρ .55 the reference density for which ε NLTE = 0.55. After tuning, we set ρ .55 = 2 × 10 -6 kg m -3 .

Additional radiative coolers

As mentioned in the introduction, the presence of radiatively cooling species at a specific altitude has been suggested to explain the decrease of temperature above 30 km [START_REF] Dias-Oliveira | Pluto s Atmosphere from Stellar Occultations in 2012 and 2013[END_REF][START_REF] Gladstone | The Atmosphere of Pluto as observed by New Horizons[END_REF]. Using the cooling-to-space approximation, we phenomenologically represent this effect with the following cooling rate for pressures below 0.12 Pa:

∂T ∂t = -5 × 10 -11 B(λ 0 , T ) (2)
with T the atmospheric temperature (K) and B(λ 0 , T ) the Planck function (in W m -2 µm -1 sr -1 ) at wavelength λ 0 .

We use λ 0 = 14 µm since the main emission bands of the most likely cooling species C 2 H 2 and HCN [START_REF] Gladstone | The Atmosphere of Pluto as observed by New Horizons[END_REF] are respectively centered at 13.7 and 14.05 µm (we here neglect the rotational bands of HCN at submillimeter wavelengths). The value -5 × 10 -11 was chosen to simulate a moderate cooling yielding temperatures below 90 K in our reference simulation.

Atmospheric molecular thermal conduction and viscosity

We account for the effect of molecular conduction on temperature and molecular viscosity on winds. Both processes are governed by similar equations. Assuming the plane-parallel approximation, for thermal conduction we

A C C E P T E D M A N U S C R I P T get: ∂T ∂t = 1 ρc p ∂ ∂z k ∂T ∂z ( 3
)
where T is the temperature (K), ρ the density (kg m -3 ) and k the thermal conduction coefficient (J m -1 s -1 K -1 ), expressed as k = k 0 T s , with k 0 = 5.63 × 10 -5 J m -1 s -1 K -(1+s) and s = 1.12 [START_REF] Hubbard | Nonisothermal Pluto atmosphere models[END_REF].

For molecular viscosity:

∂S ∂t = 1 ρ ∂ ∂z µ ∂S ∂z (4)
where S stands for the components of the horizontal wind (m s -1 ) and µ is the coefficient of molecular viscosity (kg m -1 s -1 ), that is related to the thermal conduction coefficient by k = 1 4 [9c p -5(c p -R)]µ. Given its similarity, both equations are discretized and solved using the same implicit numerical schemes.

Surface temperatures and thermal conduction in the subsurface

Surface temperature evolution T s is governed by the balance between solar insolation, thermal emisssion in the infrared, latent heat exchanges (see section 2.6), sensible heat flux from the atmosphere (usually negligible on Pluto, but taken into account in the model) and thermal conduction in the soil. On a weakly irradiated body like Pluto, the radiative fluxes are small compared to the internal heat stored in the ground. In particular, the subsurface heat stored during one season can play a major role in the control of the surface temperature at the opposite season.

The heat flux from and to the subsurface is computed using a classical model of the evolution of the subsurface temperatures T as a function of time t and depth z. It satisfies the following equation:

C ∂T ∂t = ∂ ∂z λ ∂T ∂z ( 5
)
where λ is the heat conductivity of the ground, (J s -1 m -1 K -1 ) and C the ground volumetric specific heat (J m -3

K -1
). This equation is solved using a finite differences approach and an implicit Euler scheme. The key parameter which controls the influence of the subsurface heat storage and conduction on the surface temperature is the thermal

Inertia I = √ λC.
In practice, we thus use I as the key model parameter, assuming a constant value for C=10 6 J m -3 K -1 and making λ vary accordingly.

On Pluto the discretization requires a special attention compared to the Earth or Mars because one need to simultaneously capture 1) the short period diurnal thermal waves in the near-surface, low thermal inertia terrain and 2) the much longer seasonal thermal waves which can penetrate deep in the high thermal inertia substrate. In this paper,

A C C E P T E D M A N U S C R I P T
we assumed a relatively low diurnal thermal inertia I day = 50 J s -1/2 m -2 K -1 , slightly higher than the 20 to 30 SI range reported by Lellouch et al. (2011b) from their Spitzer data analysis. For the seasonal thermal inertia, we set

I year = 800 J s -1/2 m -2 K -1
, which corresponds to a low porosity ice/rock-like substrate.

The skin depth of a thermal wave of period P (s) is:

δ P = I C P π (6)
The modeled diurnal and annual skin depths are thus 0.02 m and 40 m respectively. To represent this accurately, the subsurface is divided into N = 22 discrete layers, with a geometricaly stretched distribution of layers with higher resolution near the surface and a coarser grid at depth:

z k = z 1 2 k-1 (7)
where z 1 = 1.414 × 10 -4 m is the depth of the first layer. The deepest layer depth is thus near 300 m.

Mixing in the boundary layer

Turbulent mixing and convection are parameterized as in [START_REF] Forget | Improved general circulation models of the Martian atmosphere from the surface to above 80 km[END_REF]. In practice, the boundary layer dynamics is accounted for by a [START_REF] Mellor | Development of a turbulence closure model for geophysical fluid problems[END_REF] unstationary 2.5-level closure scheme, used to compute turbulent mixing coefficients induced by wind shears depending on the temperature profile stability and the evolution of turbulent kinetic energy. It is completed by a "convective adjustment" scheme which rapidly mixes the atmosphere in the case of unstable temperature profiles (rare on Pluto).

Turbulence and convection mix energy (potential temperature), momentum (wind), and tracers (gases and aerosols).

In the surface layer, the turbulent surface flux is given by

F = ρC d U 1 (q 1 -q 0 ), (8)
where q 1 and q 0 are the variable values in the first atmospheric layer and at the surface (q 0 = 0 for winds), U 1 is the horizontal wind velocity in the first layer, and C d is the drag coefficient. Because of the small depth of the first layer z 1 , we assume that the wind profile in the first meters above the surface is logarithmic and not influenced by stability, and simply use

C d =        κ ln z 1 z 0        2 (9)
where κ is the von Karman constant (κ = 0.4) and z 0 is the roughness coefficient, set to z 0 = 0.01 m everywhere like

A C C E P T E D M A N U S C R I P T
in the Mars GCM [START_REF] Forget | Improved general circulation models of the Martian atmosphere from the surface to above 80 km[END_REF].

Turbulent mixing is negligible outside the boundary layer (which is often shallow on Pluto because of the positive lapse rate above the surface). In our GCM there is no other vertical "eddy diffusion" process. In particular, species are only transported upwards by the large scale circulation.

N 2 Condensation and Sublimation

The condensation and sublimation of nitrogen ice must be carefully computed in the Pluto environment. The amount of energy and the relative mass of the atmosphere involved in phases changes at each timestep can be very significant. Locally, it not only changes the surface temperature and pressure, but it also modifies the structure of the boundary layer by "pumping" the air when condensation occurs on the surface, and by releasing large amount of cold, pure nitrogen (with no horizontal velocity) when N 2 sublimes. Our scheme is adapted from [START_REF] Forget | CO 2 snow fall on Mars: Simulation with a general circulation model[END_REF].

However, we found it necessary to make several changes in the equations to better represent the intense condensation and sublimation at the surface of Pluto.

The variation of the condensation temperature T c with nitrogen partial pressure P N2 is derived from the thermodynamic relations computed by Fray and Schmitt ( 2009 -1 (11) with L N2 = 2.5.10 5 J kg -1 the latent heat of condensation for nitrogen.

Surface Condensation and sublimation

The condensation and sublimation of nitrogen on the ground is primarily controlled by energy and mass conservation. At a given timestep, if the surface temperature predicted by radiative and conductive balance T * 0 falls below the condensation temperature at surface pressure T c 0, an amount δm 0 (kg m -2 ) of N 2 condenses, releasing the latent heat required to keep the solid-gas interface at the condensation temperature (T 0 = T c0 ):

δm 0 = c s (L N2 + c p (T 1 -T c0 )) (T c0 -T * 0 ) (12) A C C E P T E D M A N U S C R I P T c s is the surface heat capacity (in J m -2 K -1
), c p the air specific heat at constant pressure (set to 1040 J kg -1 K -1

for N 2 ) and L N2 the latent heat of N 2 (2.5 10 5 J kg -1 ).

The term c p (T 1 -T c0 ) (J kg -1 ) corresponds to the extra heat brought by the atmosphere (assumed to be at temperature T 1 in the first model layer) when cooled to the condensation temperature T c0 just above the surface. Because Pluto's lower atmosphere is a warm stratosphere lying just above a surface, we found that this term can be significant.

With T 1 typically 10 K above T c0 when N 2 condenses in the topics, it reaches 4% of the latent heat. Conversely, when surface N 2 ice predicted temperature T * 0 is above the frost point T c0 , N 2 sublimes and δm 0 is negative:

δm 0 = c s L N2 (T ca 0 -T * 0 ) (13)
We set T 0 = T c0 , unless all the local ground ice of mass m 0 (kg m -2 ) completely sublimes. We then set δm 0 = -m 0 and the new surface temperature is expressed as:

T 0 = T * 0 -L N2 m 0 /c s
The formation or disapearance of nitrogen ice on the substrate is taken into account in the calculations of the surface albedo and emissivity.

Atmospheric Condensation and sublimation

In the atmosphere, things are, in theory, more complex. The condensation of a gas involves various microphysical processes: supersaturation, nucleation, crystal growth, sedimentation, etc... In our model, we have kept the detailed Mars GCM CO 2 ice sheme described in the appendix of [START_REF] Forget | CO 2 snow fall on Mars: Simulation with a general circulation model[END_REF] and directly adapted it to N 2 ice. Supersaturation is neglected and atmospheric condensation and sublimation are computed using energy conservation principles as above. We do not simulate the growth and transport of nitrogen ice particles. Instead, after condensing at a given model level, we assume that N 2 ice falls through the atmospheric layers located below it (where it can sublimate), possibly down to the ground within a model timestep.

Because the atmosphere is warmer than the surface most of the time, we have found that atmospheric condensation is a processes of little importance on Pluto as we model it with a 150 km resolution. In reality, ascending motions induced by local slopes or gravity waves could trigger condensation in N 2 ice covered regions. We will explore that in future versions of the model.

Computing mass, momentum and heat vertical fluxes induced by N 2 condensation and sublimation

The condensation and sublimation of nitrogen induce significant transport of air (mass, heat, momentum, tracers) through the model layers as well as to and from the surface. These processes must be taken into account on Pluto where an atmospheric layer of several tens of meters thick can undergo a phase change at each timestep. The numerical resolution of these processes in the "σ" vertical coordinates used in the GCM (see Section 2.1) is given in the appendix.

A C C E P T E D M

A N U S C R I P T

Organic hazes

New Horizons revealed the presence of extensive hazes thought to be primarily composed of organic particles indirectly produced by methane photolysis. Our GCM includes a model of the formation and transport of these particles. This model and its outputs are described and analyzed in a companion paper by Bertrand and Forget (2016),

and not detailed here.

Methane cycle and CH 4 ice clouds

The 3D evolution of CH 4 on the surface and in gaseous and solid phase in the atmosphere is simulated taking into account 1) the condensation and sublimation at the surface and in the atmosphere (see below), 2) the transport by the general circulation using the "Van-Leer I" finite volume scheme from [START_REF] Hourdin | Test of a hierarchy of finite-volume schemes for transport of trace species in an atmospheric general circulation model[END_REF], 3) the mixing in the atmosphere by turbulent diffusion and possibly convection (see Section 2.5), and 4) the gravitational sedimentation of CH 4 ice particles (see below).

Surface condensation and sublimation.

The mass fluxes of methane to and from the atmosphere are computed using Eq. 8, with q 0 and q 1 the mass mixing ratios (kg/kg) just above the surface and in the middle of the atmospheric first layer, respectively. Note that an important consequence of Equation 8is that the sublimation rate of methane is proportional to the horizontal wind velocity in the lower atmosphere.

When pure methane is on the surface, q 0 is set equal to the saturation vapour pressure mass mixing ratio of methane q sat CH4 , calculated as a function ot temperature T (K) and pressure p using the following expression derived from Fray and Schmitt (2009):

q sat CH4 = 0.117 × 10 5 e 6.12×10 5 R (1/90.7-1/T ) × M CH4 M air × 1 p (14)
Here M CH4 /M air is the ratio of molecular masses use to convert volume mixing ratio into mass mixing ratio and R = 8.314/M CH4 = 519 m 2 s -2 K -1 the methane gas constant. When both methane and nitrogen ices are present on the surface and methane is subliming, we assume that methane is diluted in a solid solution N 2 :CH 4 with 0.3% of methane [START_REF] Merlin | New constraints on the surface of Pluto[END_REF]. Applying Raoult's law, we thus set q 0 = 0.005q sat CH4 If the total amount of methane on the surface is sublimed within a model timestep, the flux to the atmosphere is limited accordingly. If no methane ice is present on the surface, then q 0 = q 1 if q 1 < q sat CH4 (no condensation) and q 0 = q sat CH4 if q 1 > q sat CH4 (direct condensation onto the surface). The latent heat released by methane surface condensation and sublimation is taken into account in the surface energy budget assuming a latent heat L CH4 = 5.867 × 10 5 J kg -1 (Fray and Schmitt, 2009).

A C C E P T E D M

A N U S C R I P T

Atmospheric condensation and CH 4 cloud formation

Methane can also condense (and then sublimate) in the atmosphere when the CH 4 mixing ratio exceeds the saturation mixing ratio provided by Equation 14. We do not know if CH 4 can easily nucleate or if large super-saturation is required. Organic particles resulting from the photochemistry in the upper atmosphere probably offer condensation nuclei suitable for heterogeneous condensation. In the GCM we assume that all atmospheric methane in excess of saturation condenses to form ice cloud particles.

The amount of latent heat released by methane condensation or sublimation is far from being negligible. We find that it can locally change the atmospheric temperature by more than 10 K. Moreover, latent heating actually limits the amount of methane that condenses when the atmosphere is supersaturated. If CH 4 condensation is calculated without simultaneously taking into account latent heat release, or using an explicit numerical scheme, the model predicts very unrealistic temperatures (e.g. changes larger than several tens of Kelvins within one timestep), leading to unrealistic condensation rates. In practice, at each model timestep, when the methane mass mixing ratio q CH4 is detected to exceed saturation (or if methane ice is already present), one must simultaneously calculate the temperature at the end of the timestep, T , as influenced by the condensation/sublimation and the corresponding saturation mixing ratio q sat CH4 (T ). For this purpose we numerically determine T by solving the following equation:

T = T + [q CH4 -q sat CH4 (T )] L CH4 c p (15)
The change in CH 4 gas and ice mass mixing ratios (kg/kg) are then given by δq CH4 = -δq ice = (q sat CH4 (T ) -q CH4 ), ( 16) unless all the atmospheric CH 4 ice is sublimed (and T is adjusted accordingly).

Once the mass mixing ratio of CH 4 ice q ice is known, the ice is distributed to form ice cloud particles around cloud condensation nuclei (CCN). We assume that the number of cloud condensation nuclei [CCN] per mass of atmosphere (kg -1 ) is constant throughout the atmosphere. Assuming that the cloud particle size distribution is monodispersed in each volume element, the cloud particle radius r is then given by:

r = ( 3q ice 4πρ ice [CCN] + r 3 [CCN] ) 1/3 (17)
with ρ ice the CH 4 ice density (520 kg m -3 ), and r [CCN] the radius of the CCN set to 0.2 µm.

Once r is known, the cloud particle sedimentation velocity is calculated using Stokes law corrected for low pres-

A C C E P T E D M A N U S C R I P T
sure by the Cunningham slip-flow correction [START_REF] Rossow | Cloud microphysics: Analysis of the clouds of Earth, Venus, Mars, and Jupiter[END_REF]. The calculated particle radius, r, is also used to estimate the apparent opacity of the clouds. However, we neglected the radiative effect of the clouds in this paper.

[CCN] is clearly a key parameter which directly controls the properties of the clouds and their sedimentation. What is the possible range of [CCN]? On the Earth, the number mixing ratio of activated cloud condensation nuclei in the troposphere ranges between 10 6 kg -1 (for low saturation in clean polar air) and 10 10 kg -1 (polluted air mass) [START_REF] Sicardy | Pluto's atmosphere from the 29 June 2015 ground-based stellar occultation at the time of the New Horizons flyby[END_REF]Yun, 2002, Andreae, 2009]. It is significantly lower for icy cirrus clouds (<10 4 kg -1 ) [e.g. Demott et al. 2003] . On Pluto, it is likely that the organic haze particles may serve as CCN. In Bertrand and Forget (2016) we discuss the possible range of the mass mixing ratio for these particles. However, the actual number mixing ratio strongly depends on the degree of aggregation of the monomers and on their activation, which is poorly known. In our baseline simulations, we assumed [CCN]=10 5 kg -1 .

CO cycle

The CO cycle is computed using the same parameterizations than for methane, modified to use the CO properties:

the CO latent heat is set to L CO = 2.74 × 10 5 J kg -1 and the saturation mass mixing ratio q sat CO , is calculated as a function of temperature T (K) and pressure p (Pa) using the following expression derived from Fray and Schmitt (2009):

q sat CO = 0.1537 × 10 5 e 2.74×10 5 R (1/68.1-1/T ) × M CO M air × 1 p ( 18
)
Here M CO /M air is the ratio of molecular masses use to convert volume mixing ratio into mass mixing ratio and R = 8.314/M CO = 296.8 m 2 s -2 K -1 the CO gas constant.

CO is almost as volatile as N 2 and thus much more volatile than CH 4 . In practice,we found that CO only condenses when N 2 ice is present at the surface, and never forms pure CO deposits. A key parameter controlling the CO is thus the CO mixing ratio in the surface N 2 :CO ice solutions. This ratio has been estimated remotely using spectroscopic investigations of Pluto. Following the recent analysis of Very Large Telescope observations by [START_REF] Merlin | New constraints on the surface of Pluto[END_REF], we set this ratio to 0.3%.

Model initialization and choice of key parameters

Even if we had designed a perfect model of the processes at work in the Pluto environment, simulating Pluto would remain challenging. First, in spite of the New Horizons' achievements, several key parameters remain too poorly known to be used "as observed" (e.g., the global topography). Second, unlike on Mars, the Earth or even Venus, the timescales involved in the evolution of the climate system at Pluto are so long that it is difficult to reach a

A C C E P T E D M A N U S C R I P T
realistic model state insensitive to the initial state, even after running the model for weeks of computer time. Here we describe how we deal with these issues.

Topography

In our baseline simulations we assume a mostly flat surface except that we placed a 3800 m-deep circular crater roughly at the location of Sputnik Planum (in agreement with [START_REF] Moore | The geology of Pluto and Charon through the eyes of New Horizons[END_REF] as well as two smaller craters corresponding to the informally-named Burney crater (1000 m deep) and Guest crater (800 m deep). See Fig. 1.

As discussed below, we also performed sensitivity runs with a perfectly flat topography, and with two additional hypothetical 4 km-high, 800 km wide mountains that we put on the hemisphere opposite to the one better observed by New Horizons (in addition to the three craters mentioned above).

Initial Subsurface temperatures and ices distribution on the surface

On Pluto, the distribution of surface ices and subsurface temperatures (which plays a key role in the Pluto environment) are the outcome of thousand of years of evolution [START_REF] Hansen | Seasonal Nitrogen Cycles on Pluto[END_REF][START_REF] Young | Pluto's Seasons: New Predictions for New Horizons[END_REF][START_REF] Toigo | General circulation models of the dynamics of Pluto's volatile transport on the eve of the New Horizons encounter[END_REF].

Running the GCM for such a long duration is not feasible. However, initializing the model with prescribed subsurface temperatures and surface ice deposits unrelated to a natural surface evolution may be very unrealistic.

To deal with this issue, as described in [START_REF] Vangichith | Climate models of Triton and Pluto[END_REF] and like [START_REF] Toigo | General circulation models of the dynamics of Pluto's volatile transport on the eve of the New Horizons encounter[END_REF], we designed a reduced version of the GCM in which the 3D atmopsheric transport and dynamics are replaced by a simple global mixing function for N 2 , CH 4 and CO. Such a model works well on Pluto because the surface energy balance is not significantly sensitive to the atmospheric sensible heat flux and to the radiative transfer through the air. Without atmospheric dynamic and complex radiative transfer to deal with, we can perform much faster numerical simulations spanning more than 40,000 Earth years with the same horizontal grid, the same subsurface model, and the same surface/atmosphere volatiles exchange parametrizations than with the full GCM.

The details of this reduced model, its validation and the results that we have obtained are described in a separate paper Bertrand and Forget (2016). The key finding is that when we assume a topography map as described above (Fig. 1) with a 3800 m-deep "Sputnik Planum"-like basin and a seasonal ground thermal inertia larger than 800 J m -3 K -1 , after 40,000 Earth years the seasonal cycle repeats itself every year with all the nitrogen and CO ices trapped in the "Sputnik Planum"-like basin. This results from the fact that nitrogen preferentially condenses at lower altitude where the surface pressure is higher, inducing higher condensation temperature and thus enhanced thermal infrared cooling. In this model, methane still undergoes a seasonal cycle and makes seasonal deposits in both hemispheres, except in an equatorial belt which remains frost-free. Using the set of parameters described in Section 3.4 we establish a realistic, equilibrated initial state for the surface N 2 , CH 4 and CO deposits and subsurface temperatures.

A C C E P T E D M

A N U S C R I P T

Sensitivity to initial atmospheric temperatures and winds

Once the surface ices and subsurface temperatures have been initialized with the reduced GCM, the full 3D GCM should be run long enough to reach a realistic regime insensitive to the initial state assumed for the atmosphere.

This is challenging because of the long radiative time-scale of the Pluto atmosphere (several Earth years) and the time required to reach established methane and CO cycles in equilibrium with the surface reservoir. Sensitivity experiments performed with various initial temperatures, winds, and atmospheric CH 4 and CO contents showed that it takes about 20 years for two simulations initiated with two temperature profiles chosen at the end of the realistic possibilities (e.g.

differing by 30 K) to differ by less than 2 K. On this basis, we start our simulations at the end of Earth year 1988 and analyse the results after 2010. The convergences of the CO and CH 4 cycles are discussed in Section 5.

Two kind of simulations

In this paper, we describe two kinds of simulations, with and without nitrogen condensation in the south polar region in 2015.

Reference simulation, without N 2 condensation at the south pole

For the first simulation, we directly use the initial state obtained for Earth date 1988 after 40,000 Earth years of simulated climate history performed with the reduced model.

As described by [START_REF] Hansen | Seasonal Nitrogen Cycles on Pluto[END_REF] and [START_REF] Young | Pluto's Seasons: New Predictions for New Horizons[END_REF], the evolution of pressure is sensitive to the surface N 2 ice radiative properties. Some tuning was performed to select a reference value for the N 2 ice albedo A N2 and emissivity ε N2 within the range of possible values. By choosing A N2 = 0.67 and ε N2 = 0.85, we obtained an evolution of pressure (shown in Fig. 2) in qualitative agreement with the available observations [START_REF] Sicardy | Pluto's atmosphere from the 29 June 2015 ground-based stellar occultation at the time of the New Horizons flyby[END_REF][START_REF] Gladstone | The Atmosphere of Pluto as observed by New Horizons[END_REF], reaching a mean surface pressure of 1.1 Pa in July 2015.

Fig 1 shows the corresponding distribution of ice and subsurface temperature in 1988. In this simulation, the heat stored in the southern hemisphere during the previous southern hemisphere summer keeps the surface temperature above the nitrogen frost point, and nitrogen ice is only found in the "Sputnik Planum"-like basin.

The albedo of the surface CH 4 ice deposits was set to A CH4 = 0.5 and its emissivity to ε CH4 = ε N2 = 0.85.

In 1988, Methane frost covers most of the planet except for an equatorial belt which remain frost free and dark (the albedo and emissivity of the ice-free surfaces were set to A = 0.15 and ε = 1) in agreement with the observations (Stern et al., 2015;[START_REF] Grundy | Surface compositions across Pluto and Charon[END_REF].

Alternative simulation, with N 2 condensation at the south pole

It is possible that nitrogen is condensing in the south polar region in 2015. In that case, we show in this paper that Pluto's atmospheric circulation would be very different than without winter condensation, because of the induced North-south condensation flow. However, to be consistent with the evolution of surface pressure inferred from the stellar occultations since 1988, this winter condensation must be balanced by sublimation of nitrogen frost outside our modeled Sputnik Planum. In fact, New Horizons observations suggest that mid-northern latitude nitrogen frost deposits were present on Pluto in 2015 [START_REF] Grundy | Surface compositions across Pluto and Charon[END_REF].

Within that context we designed an artificial, alternative simulation by taking a modeled state from the first reference simulation at the end of 2005, with two modifications. First, we added a layer of nitrogen ice in a latitudinal belt between 35 • N and 48 • N. Second, we decreased the subsurface temperature poleward of 65 • S by 0.5 K to induce nitrogen condensation. This value was chosen in order to maintain an evolution of pressure similar to the first reference run, as shown in Fig. 2. All other modeled parameters are the same as in the reference simulation. by the sublimation of mid-northern latitude N 2 deposits. This creates planetary scale condensation flows from the northern hemisphere toward the south pole, and from the dayside toward the nightside. The wind at 20 m reaches several meters per seconds over most of the planet. In both hemisphere its direction is affected by the Coriolis force, which prevents the atmosphere from flowing directly southward.

Model results: Temperatures and winds

Surface temperatures and low level winds

Atmospheric temperatures 4.2.1. Zonal-mean temperatures

Fig. 5 presents the zonal-mean and global-mean atmospheric temperatures. As found by [START_REF] Toigo | General circulation models of the dynamics of Pluto's volatile transport on the eve of the New Horizons encounter[END_REF], the horizontal gradients of temperature are very small because of the long radiative timescale. In particular, the meridional variations in temperatures are less than 1 K. In our reference simulation with no south pole N 2 condensation, the atmospheric concentration of methane is realistic (see Section 5.1), and the mean temperature profile is in acceptable agreement with available observations [START_REF] Hinson | Radio Occultation Measurements of Pluto's Atmosphere with New Horizons[END_REF][START_REF] Gladstone | The Atmosphere of Pluto as observed by New Horizons[END_REF][START_REF] Dias-Oliveira | Pluto s Atmosphere from Stellar Occultations in 2012 and 2013[END_REF], except that above 160 km modeled temperatures are 10 to 15 K higher than reported. The thermal structure produced in our alternative simulation with South pole N 2 condensation is even warmer, because of the excessive methane concentration in this simulation, as explained in Section 12.

Comparison with the observed REX profiles

In Fig. 6, the simulated temperature profiles are compared in more details with the New Horizons REX radiooccultation profiles obtained at two locations on opposite sides of Pluto. The modeled profiles are taken at the same location and time, except that the ingress profile is shifted from latitude 17.0 • S to 7.5 • N, in order to locate it just inside the modeled Sputnik Planum basin. Indeed, on the real Pluto the ingress profile corresponds to a location just above the southern tip of the Sputnik Planum depression, above a surface covered by nitrogen ice. At the same coordinates in our model, we are outside the basin and the surface is frost free. However, we found that taking into account the low topography and N 2 coverage is key to understand the differences between the two REX profiles. We plot the modeled temperature profiles as a function of altitude above the surface. This creates an apparent shift in temperatures (the profiles are much more similar when shown in pressure coordinates) which contributes to the apparent differences reported in the observations. the Spunik Planum basin varies with local time, with coldest temperatures in the afternoon. This results from the sublimation of nitrogen ice when the sun heats the area, as proposed by Hinson et al. (2015a). In fact, the volume of gas involved in the condensation-sublimation cycle is considerable in our model. Fig. 8 shows the nitrogen ice budget in the modeled Sputnik Planum basin at 7.5 • N and 45 • N. At this last position, about 230 g m -2 of ice sublimates every Pluto day in 2015. As shown on the right axis of Fig. 8, this corresponds to more than 2500 m 3 of N 2 gas per square meter. At 7.5 • N, the solar flux is weaker in 2015 and the daily N 2 ice budget corresponds to a net gain in N 2 ice (net condensation). Nevertheless, every afternoon the equivalent of 800 m 3 per square meters is injected into the atmosphere. Moreover, in the GCM the large amount of cold N 2 gas produced at higher latitude (where the insolation is higher) is spread throughout the basin in the lowest kilometers. In fact, in the alternative simulation this process contributes to increasing the amount of cold air present in the modeled Spunik Planum basin (Fig. 7b), adding the freshly-sublimed cold N 2 gas transported from the N 2 ice belt at 35 • N (as seen on Fig. 3, right column, local time 15:00 and 21:00).

Interestingly, as shown in Fig. 7c, a simulation performed without taking into account the topographic depression in the modeled Spunik-planum does not create a significant cold layer. Two facts explain that. First, the freshlysublimed N 2 gas is efficiently transported away as discussed above (and as seen on Fig. 3, mid-column). Second, in an atmosphere with radiative timescale as long as Pluto, in a local topographic depression the temperature lapse rate is not as steep as on average because temperatures tend to be homogeneous at a given pressure level. This is illustrated on Fig. 7d which shows the temperatures at the bottom of the basin in a simulation with N 2 condensation and sublimation completely switched off. Without N 2 sublimation, the air is not as cold as in the reference simulation, but at a given altitude above the surface, temperatures in the basin remain 5 to 10 K below what they would have been outside (compare Fig. 7c and Fig 7d).

Thermal tides and waves

Stellar occultations have shown that vertical profiles of density fluctuations in the atmosphere of Pluto often exhibit wave-like structure (e.g. Elliot et al., 2003b;[START_REF] Person | Waves in Pluto's Upper Atmosphere[END_REF] with an amplitude of a few percent and vertical wavelengths of a few kilometers. On the basis of theoretical calculations, [START_REF] Toigo | Thermal tides on Pluto[END_REF] suggested that such waves could correspond to the tidal response of Plutos atmosphere to solar-induced sublimation breathing from N 2 frost patches. Here we briefly examine the type of wavelike structure present in the temperature profiles generated by our GCM. Note, however, that the horizontal and vertical resolution used in the GCM simulations is unlikely to capture waves with vertical wavelengths smaller than ∼ 20 km. Fig. 9b presents the same anomaly plot in the alternative simulations with N 2 condensation occuring at the south pole. The amplitude of the waves are significantly larger, reaching more than ±1 K around 120 km. However, a careful examination of Fig. 9b reveals that the period of the stronger waves is not 1 nor 0.5 Pluto day. These are not thermal tides: the same waves are present in simulations forced by a diurnally-averaged insolation (no diurnal cycle and no tides). These waves appear to be barotropic waves produced by a southern polar jet, as described in Section 4.3.2.

Atmospheric circulation and waves

Fig. 10 shows cross-sections of the average zonal (west-east) and meridional (south-north) winds in our two baseline simulations.

Reference Case without N 2 condensation at the south pole,

In the reference case with no condensation-flow induced by N 2 condensation at the south pole, the circulation is relatively weak with slow retrograde zonal winds in the northern hemisphere and the equatorial regions (Fig. 10a).

This circulation remains unchanged with a flat topography, no diurnal cycle, or when N 2 condensation and sublimation processes are switched off. It can be explained by the north-south latitudinal gradient of solar heating rates. It induces a very small temperature contrast between the spring and fall hemisphere and, in turn, forces weak zonal winds

A C C E P T E D M A N U S C R I P T
corresponding to the thermal wind balance. Consistently, the weak meridional circulation (Fig. 10c) is characterized by a cell centered at the equator (where the Coriolis force is null) between 80 and 140 km, with the upper branch flowing from the sunlit hemisphere toward the polar night hemisphere.

Alternative case with N 2 condensation at the south pole,

The circulation is profoundly influenced by the North-South condensation flow if N 2 condenses in the South polar regions.

If N 2 ice deposits were covering the entire northern polar regions (which is not observed) and the southern hemisphere condensation much more intense, the condensation flow would be very strong. As obtained in some of our past simulations (not shown) and as reported in some scenarios analysed by [START_REF] Toigo | General circulation models of the dynamics of Pluto's volatile transport on the eve of the New Horizons encounter[END_REF] (see their Fig. 11 and 18), the meridional circulation would be characterized by a global flow from the northern hemisphere to the southern hemisphere. In such conditions, the zonal circulation is characterized by a global "retro-superrotation" with retrograde winds at most latitude. Such winds result from the conservation of angular momentum of the air particles as they flow from the sunlit pole to the polar night above the equator, where they are farther from the rotation axis than where they started from.

In our simulations however, the North-South condensation flow remains limited compared to this extreme case.

We consider that this is in better agreement with the observations because 1) outside Sputnik Planum the N 2 ice frost deposits are limited to patches around 45-60 • N [START_REF] Grundy | Surface compositions across Pluto and Charon[END_REF], and 2) because the south pole N 2 condensation cannot be very intense in 2015 since Pluto's surface pressure has been increasing in recent years.

With the realistic assumptions made in our "alternative" simulation, the meridional circulation remains weak (Fig. 10d) and strongly modulated by waves (see below). The overall transport pattern is southward, as revealed when analysing tracer transport (Bertrand and Forget, 2016).

The zonal wind is charaterized by an intense prograde jet-stream poleward of 40 • S and a prograde superrotation at most other latitudes (Fig. 10c). The high-latitude jet is a classical feature of terrestrial atmospheres, and likely result here from the poleward condensation flow and the conservation of angular momentum. Superrotation is more surprising. It is observed on Venus and Titan and has been the subject of many studies (see, e.g. Lebonnois et al., 2010, and references therein). In these cases, superrotation is considered to primarily result from the so-called Gierasch-Rossow-Williams mechanism (from Gierasch, 1975;[START_REF] Rossow | Large-scale motion in the Venus stratosphere[END_REF]. In this mechanism, waves, possibly generated by barotropic instabilities from the high-latitude jets, redistribute angular momentum equatorward. Preliminary analysis suggest that this is what is happening in our simulation. A study of the variations of the high-latitude jet show that it is subject to instabilities that create a wavenumber 1 wave that propagates eastward with a 0.5-0.8 Pluto day period. At 60 • N, such waves are clearly visible at an altitude of 140 km in the temperature and meridional wind fields (Fig. ures show that methane is not homogeneously distributed, notably because the high latitude deposits are increasingly heated and sublimed as the sub-solar point moves northward with time. As a result methane tends to be enriched in the lower atmosphere at high northern latitudes compared to the rest of the planet, but is otherwise vertically well mixed and near 0.5% at most altitudes. This is consistent with the observation analysis of [START_REF] Lellouch | Exploring the spatial, temporal, and vertical distribution of methane in Pluto's atmosphere[END_REF] who concluded that their data "imply a roughly uniform mixing ratio in at least the first 22-27 km of the atmosphere", and that "high concentrations of low-temperature methane near the surface can be ruled out". To compare with Earthbased near-infrared observations, one must nevertheless take into account the fact that such observations are biased toward the methane column near the sub-Earth-subsolar points for geometrical reason (Pluto is a sphere) and because this is where the insolation is maximum. Taking into account that the sub-Earth and subsolar points are always very close, we can estimate the apparent mixing ratio as seen from the Earth by performing an average of the local column mixing ratio weighted by the square of the cosine of the solar zenith angle. The apparent mixing ratio for the reference simulations started with 0.5 % CH 4 is shown in green on Fig. 12c. The difference with the global-mean value remain small and has only become significant recently with the local increase of methane above the North pole.

On the same Figure 12c, the blue dashed curve shows the evolution of the global-mean methane in the alternative scenario (with N 2 condensing at the south pole) starting in 2005. Fig. 12d show the corresponding methane abundances as a function of latitude and altitude in 2015. One can see that the methane content is larger and still increasing in this simulation. This results from the stronger near-surface winds induced by the condensation flow, and the fact that the near-surface mixing is directly proportional to the horizontal wind as formalized in Equation 8 presented in Section 2.5.

A C C E P T E

CO cycle

Fig. 15 shows the evolution of the carbon monoxide mixing ratio as a function of time since 1988. The red curves correspond to the global-averaged mixing ratio for three different initial values (0%, 0.05%, 0.1%). Clearly, the three simulations have not converged but one can estimate that the system evolves toward a mean mixing ratio near 0.03%.

A mixing ratio of 0.03% is in acceptable agreement with the 0.05 +0.01 0.025 % reported by Lellouch et al. (2011a) from telescopic observations performed in 2010, and of the same order of magnitude as the 0.0515 ± 0.004% just retrieved by [START_REF] Lellouch | Detection of CO and HCN in Pluto's atmosphere with ALMA[END_REF] using the ALMA interferometer on June 12-13, 2015.

In details, the CO cycle is dominated by a condensation-sublimation cycle inside Sputnik Planum. For instance in 2015 there is a net flux from the northern part and the center part of Sputnik Planum to the southern part where nitrogen is condensing along with CO. We do not show here the spatial distribution of CO since we have found that CO is usually very well mixed with N 2 . As a consequence, the apparent CO mixing ratio as seen from the Earth (green curve in Figure 15) is very close to the global mean.

When the alternative simulation is started in 2005 with N 2 condensing in the high southern latitudes (blue lines in Figure 15), the CO mixing ratio rapidly decreases to reach values below 0.03%. This is even the case when we assume that all mid-northern latitude N 2 frost deposits contains 0.3% of CO. In these conditions, the atmospheric CO appears to decrease below 0.03% because the ices that condense in the south polar cap tends to be enriched in CO, up to 0.05% at the pole.

In reality, the mid-latitude N 2 frost deposits have been observed by New Horizons to be strongly depleted of CO compared to Sputnik Planum [START_REF] Grundy | Surface compositions across Pluto and Charon[END_REF]. If we take this into account and set the N 2 :CO mixing ratio to zero in these deposits, we obtain the evolution shown by the dashed blue line Figure 15, with an additional decrease of atmospheric CO down to less than 0.01% in 2015. One can guess that these values could be tuned up by increasing the assumed N 2 :CO ice mixing ratio in Sputnik Planum. This would still be consistent with the Merlin (2015)'s telescopic measurements since they included both Sputnik Planum and the mid-latitude deposits. Further work will be required to fully understand the long term CO equilibrium, its evolution, and the surface N 2 :CO mixing ratio.

A C C E P T E D M

A N U S C R I P T

Conclusions

The goal of this paper was to describe, for the first time, our new Global Climate Model of Pluto including the N 2 , CH 4 and CO cycles. We presented two baseline simulations which can shed light on New Horizons observations, for instance to understand the low atmosphere temperature profiles measured by REX and the distribution of ices.

However, this is just the beginning. One of our key conclusions is that the Pluto climate system is extremely sensitive to the assumed model parameters, such as the ice properties, the ground thermal inertia, or the topography. Many more studies will have to be performed to better simulate the reality and understand the processes at work. It will also be very useful to perform longer simulations, with higher model resolution, with a more realistic topography, etc... We hope that this GCM will be applied to many specific studies regarding clouds, hazes, frost deposits, seasonal evolution, and paleoclimates.

APPENDIX: Computing mass, momentum and heat vertical fluxes induced by N 2 condensation and sublimation in the GCM vertical coordinates

In the GCM, the changes in atmospheric mass due to the condensation and sublimation of nitrogen are taken into account by modifying the surface pressure p 0 at each timestep by: δp 0 = -g N k=0 δm k , with N the number of atmospheric model layers and δm k the mass condensed (or sublimed if < 0) in layer k or at the surface (k = 0), as described in Section 2.6. This ensures the conservation of the total mass of N 2 (surface caps + atmosphere).

As described in Section 2.1, the vertical coordinate of each model layer is defined by its σ l = p l /p 0 coordinates. The changes in p 0 due to the N 2 condensation-sublimation induce "artificial" movements of the σ levels in the atmosphere. This must be reflected in the temperature and wind fields.

Consider a layer l delimited by the levels σ l-1 2 and σ l+ 1 2 . At each timestep, its mass M l = p 0 g (σ l-1 2σ l+ 1 2 ) (in kg m -2 ) varies because of the global variation of p 0 . Such a variation δM l is associated with transfers of mass between the layers (on which one must add the sink corresponding to the local condensation -δm l ). The local mass balance may be written :

δM l = δp 0 g (σ l-1 2 -σ l+ 1 2 ) = W l-1 2 -W l+ 1 2 -δm l ( 19
)
where W l-1 2 is the air mass (kg m -2 ) "transfered" through the level σ l-1 2 (> 0 when up) during the timestep. Equations 19 can be rearanged to yield a recursive formula on W :

W l+ 1 2 = W l-1 2 -δm l - δp 0 g (σ l-1 2 -σ l+ 1 2 ) (20)

  ), taking into account the transition from the α to the β crystalline form near 35.61 K (corresponding to P N2 = 0.53 Pa): if P N2 < 0.53 Pa :

Figure 1 :

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

4. 1 . 1 .Figure 2 :

 112 Fig 3 shows maps of surface temperatures and winds at 20 m above the surface at various times of the day for our different simulations. The epoch corresponds to July 2015, the time of the New Horizons encounter. In these simulations, surface temperatures range between 36.6 and 48 K. The lowest values correspond to the N 2 frost point around 1 Pa. The highest temperatures are more model dependent, and vary with the assumed diurnal thermal inertiaI day . Daytime surface temperatures reach 57 K in GCM runs, assuming I day = 20 J s -1/2 m -2 K -1 (as reported byLellouch et al., 2011b) instead of I day = 50 J s -1/2 m -2 K -1 , as assumed in our baseline simulations.

1 Figure 3 :Figure 4 :

 134 Figure 3: Maps of surface temperature and winds at 20 meters above the surface in July 2015 at different local times for 3 simulations: a) the reference simulations with no N 2 condensation at the south pole, b) The same simulation with flat topography (started from the reference run on Juanuary 1st, 2015, and analyzed on July 14, 2015) c) the alternative simulation with N 2 condensation at the south pole. From top to bottom, the local time LT in the middle of the map (longitude 180 • ) is 3:00, 9:00, 15:00 and 21:00, with LT (hours) = [longitude ( • ) -subsolar point

  Of special interest are the lowest kilometers of the simulated ingress profiles which exhibit a low temperature layer analogous to the bottom of the observed ingress profile. Which process creates this layer? To better understand this behaviour, and possibly interpret the observations, we show in Fig 7 the diurnal evolution of the atmospheric profile in the lowest 4 km in different modeled configurations. In the reference simulation, the atmospheric temperature in

Fig. 9a presents 1 Figure 6 :Figure 7 :Figure 8 :Figure 9 :

 16789 Fig.9apresents the 4-sols evolution of the difference between instantaneous temperatures and 1-Pluto-day gliding averages at 0 • E -0 • N in our reference simulation. The observed temperature excursions are lower than 0.2 K. Nev-

  11b

Figure 11 : 1 .

 111 Figure 11: Characteristics of the barotropic waves present in the simulations with South Pole N 2 condensation inducing a condensation flow. a) Hövmoller diagram of the temperature anomaly (difference between the local and the zonal-mean temperature) at 0 • N. b) Same at 60 • N. c) Zonal average of the root-mean-square standard deviation of the local meridional wind from the zonal-averaged meridional wind. d) Hövmoller diagram of the meridional wind anomaly (difference between the local and the zonal-mean wind, in m s -1 ) at 0 • N.

  Figure 12: a) -b): Zonal mean methane volume mixing ratio (%) in the reference simulation (without south pole N 2 condensation and [CH 4 ] initially at 0.5% in 1988) in 2010 and 2015. c): Evolution of the mean volume mixing ratio: globally averaged with different initialization (red), the apparent mixing ratio as seen from the Earth (green, see text) and in the alternative simulation with south pole N 2 condensation started in 2005 (dashed blue). d): Zonal mean methane volume mixing ratio (%) in the alternative simulation (with south pole N 2 condensation) in 2015.

Figure 14 :

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

Figure 15 :

 15 Figure 15: Evolution of the mean volume mixing ratio of gaseous carbon monoxide. The red curves present the globally averaged values with different initialization. The green curve shows the apparent mixing ratio as seen from the Earth. The blue curves shows the global mean mixing ratio in the alternative simulation with south pole N 2 condensation started in 2005 (see text).

Acknowledgments

The authors thank the NASA New Horizons team for their excellent work on a fantastic mission and their interest in this research. We also thank CNES for its support. Finally, the authors are very grateful to two anonymous reviewers for their exceptionally detailed comments which helped to improve the readibility of the article.

and 11d). In Fig. 11c, the extension of this wave is mapped by plotting the meridional wind variability as a function of latitude and height. One can see that it propagates to all latitudes, and notably to the equator, where the signature in the thermal field dominates the temperature variability (Fig. 11a). Similar results are obtained in model runs without a diurnal cycle or with a flat topography.

In addition to the wind predictions published by [START_REF] Toigo | General circulation models of the dynamics of Pluto's volatile transport on the eve of the New Horizons encounter[END_REF], already discussed, our results can be compared with the results from the other Pluto GCM proposed by [START_REF] Zalucha | A 3D general circulation model for Pluto and Triton with fixed volatile abundance and simplified surface forcing[END_REF] and [START_REF] Zalucha | An atmospheric general circulation model for Pluto with predictions for New Horizons temperature profiles[END_REF]. The comparison with Zalucha and Michaels ( 2013) is difficult to achieve because this version of their GCM did not yet include nitrogen condensation and because their modeled thermal structure was very different than what was observed on Pluto by New Horizons. In fact the updated version presented by [START_REF] Zalucha | An atmospheric general circulation model for Pluto with predictions for New Horizons temperature profiles[END_REF] yielded completely different results. Her "Case 1", in which a surface pressure of 0.8 Pa and 1% of CH 4 is assumed, can be compared to our simulations. The zonal wind structure ressemble our reference simulation, suggesting that, for unknown reasons, the condensation flow is weak in this GCM in spite of the fact that Pluto is assumed to be covered by nitrogen ice.

A C C E P T E D M A N U S C R I P T

with, in the first layer:

The knowledge of W can then be used to compute the exchange of heat and momentum between the layers. For c p T (enthalpy), the local heat balance can be written :

with T l-1 2 the mean temperature of the gas transported through the σ l-1 2 interface. The calculation of T l-1 2 is like in a classical transport problem. We use the "Van-Leer I" finite volume transport scheme [START_REF] Van Leer | Towards the ultimate conservative difference scheme : IV. a new approach to numerical convection[END_REF][START_REF] Hourdin | Test of a hierarchy of finite-volume schemes for transport of trace species in an atmospheric general circulation model[END_REF]. Separately, one can also write :

with δT l the correction to be applied at every timestep in each layer after the N 2 condensation or sublimation.

Eqs 22 and 23 may be combined to obtain δT l

The first two terms, with W l-1 2 and W l+ 1 2 , correspond to the re-arrangement of the temperatures over the entire column due to the pressure variations in σ coordinates. The last term δm l (T c l -T l ) is negligible when N 2 condenses or partially sublimes since we then have T c l = T l . However, when the N 2 totally sublimes , it becomes a cooling term accounting for the mixing of the newly sublimed mass -δm l with the rest of the layer at T l > T c l .

On the ground, if δm 0 > 0 (condensation), we set T 1 2 = T 1 . As mentioned above, the near-surface cooling of the condensing N 2 gas from T 1 to T 0 is then taken into account in the surface energy balance. If δm 0 < 0 (sublimation), we set T 1 2 = T 0 . The term δm 0 (T 0 -T 1 ) then accounts for the cooling of the lowest level by the freshly-sublimed nitrogen.

Similarly, the momentum distribution must be re-arranged. For a wind component v, we shall simply write:

with, on the ground, v 1 2 = v 1 if δm 0 > 0 and v 1 2 = 0 if δm 0 < 0 (the velocity of the N 2 gas that has just sublimed is zero).