# The New Titan Planetary Climate Model. I. Seasonal Variations of the Thermal Structure and Circulation in the Stratosphere

Bruno de Batz de Trenquelléon, Lucie Rosset, Jan Vatant D'ollone, Sébastien Lebonnois, Pascal Rannou, Jérémie Burgalat, Sandrine Vinatier

► To cite this version:

Bruno de Batz de Trenquelléon, Lucie Rosset, Jan Vatant D'ollone, Sébastien Lebonnois, Pascal Rannou, et al. The New Titan Planetary Climate Model. I. Seasonal Variations of the Thermal Structure and Circulation in the Stratosphere. The Planetary Science Journal, 2025, 6 (4), pp.78. 10.3847/PSJ/adbbe7. hal-05041885

HAL Id: hal-05041885

https://hal.sorbonne-universite.fr/hal-05041885v1

Submitted on 22 Apr 2025

**HAL** is a multi-disciplinary open access archive for the deposit and dissemination of scientific research documents, whether they are published or not. The documents may come from teaching and research institutions in France or abroad, or from public or private research centers.

L'archive ouverte pluridisciplinaire **HAL**, est destinée au dépôt et à la diffusion de documents scientifiques de niveau recherche, publiés ou non, émanant des établissements d'enseignement et de recherche français ou étrangers, des laboratoires publics ou privés.

<figure><img src="image_1.png" /></figure>

Distributed under a Creative Commons Attribution 4.0 International License
OPEN ACCESS

<figure><img src="image_2.png" /></figure>

CrossMark

# The New Titan Planetary Climate Model. I. Seasonal Variations of the Thermal Structure and Circulation in the Stratosphere

Bruno de Batz de Trenquelléon<sup>1,2 ID</sup>, Lucie Rosset<sup>2</sup>, Jan Vatant d'Ollone<sup>2</sup>, Sébastien Lebonnois<sup>2 ID</sup>, Pascal Rannou<sup>1 ID</sup>, Jérémie Burgalat<sup>1</sup>, and Sandrine Vinatier<sup>3 ID</sup>

## Abstract

Observations of Titan through Cassini's mission allowed Saturn's moon's stratospheric thermal structure and composition to be mapped over half a Titan year. Seasonal variations revealed various unexplained phenomena, such as mechanisms within the polar vortex, thermal structure evolution at high latitudes, and the impact of the enrichment in trace compounds during winter polar nights. We have developed the Titan Planetary Climate Model (Titan PCM)—an improved version of the IPSL Titan Global Climate Model (GCM)—including upgraded radiative transfer, now based on a flexible correlated-k method, updated spectroscopic data on gases, and integration of a new microphysics model for haze and clouds. Our photochemical solver extends computation of the composition above the top of the model up to 1300 km. The radiative transfer is now coupled with microphysics, consisting of the first full radiative coupling of microphysics within a Titan GCM. The model is presented along with its dynamics, thermal structure, and seasonal variations. Despite biases above 10 Pa due to upper boundary limitations, the Titan PCM leads to better modeling of the temperature profiles in the middle atmosphere. Consequently, it is now possible to address some scientific issues about thermal structure in polar regions. Investigations into fall and winter polar phenomena, particularly focusing on thermal structure control and equinoctial circulation reversal, were conducted. This study is linked to the radiative destabilization of the lower polar stratosphere, observed at the end of winter by Cassini radio-occultations. The results confirm a dynamic interplay between haze and gas distributions, which influence the thermal structure.

*Unified Astronomy Thesaurus concepts: Titan (2186); Planetary atmospheres (1244); Planetary climates (2184); Radiative transfer (1335); Planetary polar regions (1251)*

## 1. Introduction

The longevity of the Cassini mission has provided a gathering of data that enlightened our knowledge of Titan's atmosphere to a large degree. Among others, infrared observations by the Composite InfraRed Spectrometer (CIRS) on board Cassini (F. M. Flasar et al. 2004) enabled the seasonal evolution of the thermal structure and composition to be mapped over almost half a Titan year in the middle atmosphere (M. Fulchignoni et al. 2005; H. B. Niemann et al. 2005, 2010; S. Vinatier et al. 2007, 2015, 2020; N. A. Teanby et al. 2008b, 2009, 2010, 2012; P. J. Schinder et al. 2011, 2012; B. Bézard 2014; A. Coustenis et al. 2016, 2018; C. Mathé et al. 2020) and retrieve zonal winds (M. K. Bird et al. 2005; F. M. Flasar et al. 2005; R. K. Achterberg et al. 2008a, 2011; J. Sharkey et al. 2021). All these works allowed the global overturning circulation to be inferred, as well as the ongoing processes in the winter polar vortex (N. A. Teanby et al. 2008a, 2017, 2019; S. Vinatier et al. 2015, 2020; J. Sharkey et al. 2020) along with the detection of high-altitude condensates (D. E. Jennings et al. 2012, 2015; R. J. de Kok et al. 2014; C. M. Anderson et al. 2018; S. Le Moulicé et al. 2018; S. Vinatier et al. 2018).

Still, some features remained unexplained or partially understood (S. M. Hörst 2017). For example, S. Vinatier et al. (2015) and N. A. Teanby et al. (2017) suggested that radiative cooling by trace compounds advected by the descending branch of the circulation could compensate adiabatic heating at the south pole. This hypothesis could explain the observed evolution of the thermal structure following the 2009 equinox, particularly the unexpected temperature drop in early 2012. However, further investigation is needed to confirm this mechanism.

Deeper into the atmosphere, the lower stratosphere was probed by Cassini's radio-occultations for temperature (P. J. Schinder et al. 2011, 2020), and by CIRS nadir observations for both thermal structure and composition (M. Sylvestre et al. 2018, 2020). The destabilization of the low winter polar stratosphere remains a topic of discussion concerning its origin, whether it be haze, clouds (D. E. Jennings et al. 2015), or dynamics. In a broader range, this region is of particular interest regarding seasonal variations as it is the transition region below which the characteristic radiative relaxation timescales are long enough to mute the seasonal response of the atmosphere (B. Bézard et al. 2018). It is also the altitude level where the vertical zonal wind profile has a sharp minimum between 60 and 100 km, as observed at low latitudes by the Doppler Wind Experiment on board the Huygens probe (M. K. Bird et al. 2005). The proposed hypotheses to explain this minimum are still largely unclear.

<figure><img src="image_3.png" /></figure>

Original content from this work may be used under the terms of the [Creative Commons Attribution 4.0 licence](https://creativecommons.org/licenses/by/4.0/). Any further distribution of this work must maintain attribution to the author(s) and the title of the work, journal citation and DOI.

<footnote n="1"><sup>1</sup> Groupe de Spectrométrie Moléculaire et Atmosphérique (GSMA), Université de Reims Champagne-Ardenne, CNRS, 51687 Reims cedex, France;
bruno.de-batz-de-trenquelleon@univ-reims.fr</footnote>

<footnote n="2"><sup>2</sup> Laboratoire de Météorologie Dynamique (LMD/IPSL), Sorbonne Université, ENS, PSL Research University, Ecole Polytechnique, IP Paris, CNRS, 4 Place Jussieu, 75252 Paris Cedex 05, France</footnote>

<footnote n="3"><sup>3</sup> LESIA, Observatoire de Paris, Université PSL, CNRS, Sorbonne Université, Univ. Paris Diderot, Sorbonne Paris Cité, 5 place Jules Janssen, 92195 Meudon, France
Received 2024 August 27; revised 2025 January 14; accepted 2025 January 28; published 2025 March 31</footnote>
Another topic open to question is the delay between the predicted reversal of the mean meridional circulation at the equinox and the observed enrichment at the winter pole. This also includes the asymmetry of this delay between the decrease of trace compounds at high northern latitudes, observed only since 2015, and the enrichment at high southern latitudes, which has increased considerably since 2012 (S. Vinatier et al. 2015, 2020; N. A. Teanby et al. 2017, 2019; A. Coustenis et al. 2018; C. Mathé et al. 2020; J. Sharkey et al. 2020).

Global climate models (GCMs) are very powerful tools to address these problems. Throughout the Cassini mission, a buildup of the observations provided a lot of constraints for these models and enabled their accuracy to be progressively strengthened. Subsequently, these models facilitate a deeper understanding of the observations, offering interpretations of the observed mechanisms. Four GCMs were recently used to bring concrete new insights into the behavior of Titan's atmosphere: the TitanWRF model (C. E. Newman et al. 2011, 2016), the IPSL Titan GCM (B. Charnay & S. Lebonnois 2012; S. Lebonnois et al. 2012), the Titan Atmospheric Model (TAM; J. M. Lora et al. 2015, 2022; N. A. Lombardo & J. M. Lora 2023), and the Cologne Titan GCM (T. Tokano 2019, 2023).

Most of these models use a two-stream radiative transfer model to determine heating rates at visible and infrared wavelengths. Gas and haze optical properties are determined from the scheme, described by C. P. McKay et al. (1989), as in the case of the TitanWRF model (C. E. Newman et al. 2011), the IPSL Titan GCM (S. Lebonnois et al. 2012), and the Cologne Titan GCM (T. Tokano et al. 1999). Among the four models mentioned, only the TAM does away with the radiative transfer of C. P. McKay et al. (1989) by opting for an alternative radiative solver and using recent spectroscopic data, with the help of imposed fields (N. A. Lombardo & J. M. Lora 2023). It distinguishes itself by obtaining a satisfactory thermal structure in the stratosphere, including a stratopause. The model by A. J. Friedson et al. (2009), based on the Titan NCAR Community Atmosphere Model, also does not use the C. P. McKay et al. (1989) radiative heating model, but it is unable to reproduce the superrotation. With the exception of the TAM (N. A. Lombardo & J. M. Lora 2023), no other GCM has a radiative transfer scheme allowing feedback from trace compound variations that could permit induced seasonal variations to be dealt with. A glimpse of these feedbacks can be found in S. Lebonnois et al. (2003), but this study was done with idealized tracers in the model. On these aspects, a review of the importance of the coupling between dynamics, microphysics, and chemistry has been performed in S. Lebonnois et al. (2009).

Within the scope of this work, we present the new Titan Planetary Climate Model (hereafter the Titan PCM), an improved version of the IPSL Titan GCM (S. Lebonnois et al. 2012). A new flexible radiative transfer scheme based on the correlated-k method and recent spectroscopic data was implemented. Furthermore, the Titan PCM is now fully coupled, encompassing chemistry, microphysics, radiative transfer, and dynamics. This approach leads to more accurate temperature profiles simulated in the middle atmosphere, and enables an overview study of the mechanisms and feedback processes involved in the upper atmosphere. This is crucial for understanding the control of the thermal structure and its implications for atmospheric circulation.

A general presentation of the model is given in Section 2, with particular attention on a description of radiative transfer in Section 2.4. Three-dimensional simulations that were performed with the Titan PCM are described and analyzed in terms of thermal and dynamical structure in Section 3. A study of the mechanisms taking place in polar regions is conducted in Section 4.

In the companion paper (B. de Batz de Trenquelléon et al. 2025), hereafter referred to as “Paper II,” the new microphysical model for haze and clouds implemented in the Titan PCM is presented, along with its coupling with radiative transfer. The study of haze and cloud cycles draws parallels with the thermal and dynamic structure presented in this paper.

## 2. The Titan PCM

### 2.1. Overview

The Titan PCM used in the present work is based on the reference version of the IPSL Titan GCM (B. Charnay & S. Lebonnois 2012; S. Lebonnois et al. 2012) with an updated dynamical core (F. Hourdin et al. 2013). Major changes are included in the physics part, with an updated radiative transfer, featuring a correlated-k solver (see Section 2.4), coupled to new photochemistry (see Section 2.2) and microphysics (see Section 2.3) modules.

We are using the dynamical core LMDZB, described in its most recent version in F. Hourdin et al. (2013). This dynamical core is based on a finite-difference discretization scheme that conserves among other things both the potential entrophy for barotropic nondivergent flows and the total angular momentum for axisymmetric flows. A horizontal spatial resolution of 32 longitudes by 48 latitudes (11.25° × 3.75°) is used in the presented simulations, and for a standard configuration, the dynamical time step is roughly 3 minutes, computed 8000 times per Titan day. The vertical grid is left unchanged compared to the previous version of the model, with 55 levels extending from the surface to the $10^{-1}$ Pa level (~500 km). In the top four layers of the Titan PCM, a sponge layer is included, with Rayleigh friction damping horizontal winds to their mean value. Damping time constants are $10^4$ s for the last level, then $2 \times 10^4$, $4 \times 10^4$, and $8 \times 10^4$ s, respectively, for the three levels below.

Boundary layer scheme and surface parameters are the same as in S. Lebonnois et al. (2012), apart from surface thermal inertia, which is adjusted to $I = 335 \text{ J m}^{-2} \text{ s}^{-0.5} \text{ K}^{-1}$. This value is taken from that estimated by the Cologne Titan GCM (T. Tokano 2005) for a soil consisting of porous icy regoliths. Tests with different thermal inertia values showed that only the near-surface environment of the model was sensitive to this parameter. The results presented here are not impacted by thermal inertia variations, but future studies focusing on the troposphere could be conducted using the thermal inertia maps by S. M. MacKenzie et al. (2019). Topography from P. Corlies et al. (2017) is now included in the model. We point out that the large-scale topography limited at the eighth order in the spherical harmonic expansion is used in the present study. Indeed, the higher resolution maps, whether they be in R. D. Lorenz et al. (2013) or P. Corlies et al. (2017), reveal stitches due to radar swathes that lead to irrelevant steep large-scale slopes between regions of extrapolated data that generate noisy stationary patterns in the large-scale circulation in all of the troposphere (J. M. Lora et al. 2019). In any case, the
inclusion of topography does not affect our results much in the stratosphere.

The physical parameterizations are computed200 times per Titan day, i.e., once every 2 hr. We have increased the number of radiative transfer calls to 40 per Titan day, enhancing our ability to account for the new scheme with a better description of the thermal profile (see Section 3), and to provide a comprehensive description of the influence of clouds on Titan's tropospheric thermal profile (see Paper II). Nevertheless, given the time it takes to calculate chemical processes, we only perform these calculations two times per Titan day.

## 2.2. Photochemistry

The photochemical solver implemented is the same as in the former version of the IPSL Titan 2D GCM (S. Lebonnois et al. 2001; A. Crespin et al. 2008) and used in the 3D reference version (S. Lebonnois et al. 2012). In this work, the reaction scheme has been updated based on recent modeling efforts by V. Vuitton et al. (2019), along with technical improvements, as described hereafter.

The setup used for the present work includes 43 species and 323 reactions, with only neutral chemistry, even though this is known to be a serious caveat for species such as benzene (V. Vuitton et al. 2008) or hydrogen cyanide (P. Lavvas et al. 2008; V. Vuitton et al. 2019). For the present study, which will focus on the radiative impact of ethane ($C_2H_6$), acetylene ($C_2H_2$), and hydrogen cyanide (HCN), the fact that ionospheric photochemistry is not included should mostly be a source of uncertainty for the production rate of HCN, which includes substantial production pathways through ionic processes (P. Lavvas et al. 2008; V. Vuitton et al. 2019).

In former studies with photochemistry coupled to advection, the upper boundary condition was external input fluxes of species at 500 km (0.1 Pa), previously generated offline from 1D chemistry. In the present work, the photochemistry is now solved from the surface up to 1300 km on each 1D column, as described in S. Lebonnois et al. (2001). In the upper atmosphere—above the dynamical core upper limit at 0.1 Pa—we use a vertical grid discretization of 10 km built from the atmospheric structure derived by R. J. J. Vervack et al. (2004), assuming hydrostatic equilibrium. An altitude-dependent gravity field is included, as it is compulsory to have correct photodissociation rates, which are determined by the optical depth at a given level (S. Lebonnois et al. 2001). Vertical transport by eddy diffusion is parametrized by a vertical profile—taken from V. Vuitton et al. (2019) for the upper atmosphere and then stopped at the upper limits of the GCM, as below the mixing is performed by advection of the tracers within the dynamical core. At the interface with the rest of the physics, at a given altitude, the obtained production rates are then recast on the grid of the dynamical core below 500 km. As in the previous version of the model (S. Lebonnois et al. 2012), due to the high computational costs, photochemical processes are conducted using the 3D fields zonally averaged. The resulting tendencies are then converted back to 3D fields proportionally to the corresponding tracer field. This implies that we have no sensitivity to the photochemical diurnal cycle.

The photodissociation rate branching ratios have been entirely updated from V. Vuitton et al. (2019), as well as the chemical scheme, but keeping the same limited number of species. These updates significantly improve the production rate of trace compounds in the upper atmosphere, as the former

scheme led to depleted profiles compared to observations in this region. However, the photochemical solver integrated into Titan's PCM consistently produces mixing ratios that are approximately 10 times lower than those observed by CIRS in the lower stratosphere. This results in an underestimation of mixing ratios through all lower layers of the model, thereby limiting stratospheric cooling during the polar night. To address this issue, we impose a constant mixing ratio in the model sponge layer for species involved in radiative transfer. Specifically, we apply a forcing equivalent to 3 times the photochemical production of acetylene, 5 times that of ethane, and 3 times that of hydrogen cyanide in the sponge layer. The application of these factors yields more satisfactory mixing ratios in the stratosphere, significantly enhancing their radiative impact within the model and consequently improving its thermal and dynamical structures.

Among other improved features, the saturation profiles are now computed from N. Fray & B. Schmitt (2009) on each grid point with local atmospheric conditions and no more than once at the beginning of the simulations from the atmospheric profile from E. Lellouch et al. (1989). In former versions of the model, this was done to avoid problems due to the nonrealistic temperature profiles above poles (F. Hourdin et al. 2004), but as the thermal structure above the poles is now more satisfactory—as presented in Section 3—this constraint is released, which enables us to study, for instance, the impact of seasonal variation of temperature in the polar low stratosphere on the rates and altitude of gas condensation, as well as the feedback on temperature in radiatively active variable species configurations.

## 2.3. Microphysics

Considering the computational cost imposed by the bin formalism in the previous microphysical scheme based on P. Rannou et al. (2004) and the need to pursue long-term simulations (tens of Titan years) led to a change of the microphysical model in the present study compared to previous versions of the model (S. Lebonnois et al. 2012). We have now integrated in the Titan PCM an advanced microphysical model adapted to Titan for haze (J. Burgalat & P. Rannou 2017) and clouds (J. Burgalat et al. 2014). Diverging from traditional bin distribution models, this model has been developed using moments to address computational resource challenges associated with bin models. It also simulates previously unrepresented phenomena in 3D simulations, including the nucleation and condensation of clouds. As the microphysical model is described in detail in Paper II, we will only summarize the main points here.

The fractal structure of aerosols on Titan leads to size distributions that differ significantly from typical log-normal patterns observed in natural environments (P. Rannou et al. 2004; M. G. Tomasko et al. 2008; T. Cours et al. 2011; E. J. L. Larson et al. 2014). As the size distribution of aerosols is complicated to predict, we use in the Titan PCM the size distribution from the bin model of P. Rannou et al. (2004) composed of 40 sections for an accurate description. In the models of J. Burgalat et al. (2014) and J. Burgalat & P. Rannou (2017), each population (spherical aerosols, fractal aggregates, and cloud condensation nuclei) is described by only two moments, $M_0$ and $M_3$, respectively the total number and total volume of particles, to a constant factor. This enables us to adequately describe the haze with only four tracers, and two
additional tracers for cloud nuclei and one tracer for each condensed (ice or liquid) species to account for. That makes $6 + N_{\text{species}}$ tracers to account for, compared to about 100 tracers for a classical 50-bin model.

In our model, the injection of spherical aerosols due to the production of photochemical aerosols at a high altitude is performed at the top of model with a simple production law using a Gaussian function (M. Cabane et al. 1992). The model for haze by J. Burgalat & P. Rannou (2017) solves the intramodal and intermodal Brownian coagulation equations for spherical and fractal aerosols. The sedimentation process is implemented in the model using a flux method, as presented by O. B. Toon et al. (1988), with the first-order approximation of the slip-flow correction. With the properties of aerosols being loosely constrained, we adjusted the values in S. Lebonnois et al. (2012) for the aerosol density and the aerosol production rate to $800 \text{ kg m}^{-3}$ and $2.1 \times 10^{-13} \text{ kg m}^{-2} \text{ s}^{-1}$, respectively, in order to obtain lighter aerosols and a more pronounced detached haze layer. The altitude of production is set at 1 Pa, that is around 400 km, as in previous versions of the microphysical model. This altitude is limited by the vertical extent of the model and can lead to biases in the distribution of haze and its detached layer.

We have also implemented in the Titan PCM the cloud model of J. Burgalat et al. (2014). It takes into account the nucleation phenomenon of cloud condensation nuclei and the cloud growth process called condensation. Clouds are subject to the same sedimentation as aerosols, with the difference that drops here are considered to be perfect spheres. In the current model, we deal with the condensation of CH$_4$, C$_2$H$_6$, C$_2$H$_2$, and HCN. Each of these condensable species is considered independent, and there is no mixing between condensed species.

These improvements enable us to fully model the haze, its main layer as well as its detached layer, in Titan's stratosphere. They have also made it possible to represent clouds in Titan's troposphere, paving the way for in-depth exploration of the hydrocarbon cycle. Furthermore, the microphysical model has been coupled with the Titan PCM radiative transfer, as detailed in Section 2.4.4. This new microphysical model is computationally more efficient than the bin model and therefore no longer necessitates us to have microphysical tendencies computed assuming zonal means, as was the case before (S. Lebonnois et al. 2012).

## 2.4. Radiative Transfer

### 2.4.1. Motivations and Overall Description

One of the main issues in the preexisting model (S. Lebonnois et al. 2012) was that the radiative transfer did not include all the opacities needed for stratospheric studies, such as that of HCN. In addition, the abundances used were based neither on recent Cassini observations nor on up-to-date spectroscopic data, with no radiative transfer calculations at wavelengths between 2.3 and $10.6 \mu\text{m}$. Indeed, based on C. P. McKay et al. (1989), which was motivated by tropospheric studies, it needed to be updated and especially to take into account the methane $\nu_4$ band at $7.7 \mu\text{m}/1305 \text{ cm}^{-1}$, which is prominent in the radiative balance of the stratosphere with a strong cooling effect. This band exchanges a small amount of flux at tropospheric temperatures but shall not be neglected at warmer stratospheric temperatures (see Figure 1).

<figure><img src="image_4.png" /><figcaption>**Figure 1.** Infrared absorption spectra of Titan's atmosphere at 1 Pa and 180 K computed with *kspectrum* (V. Eymet et al. 2016). Planck blackbody functions at 80 K (blue) and 180 K (red) are presented for comparison. This enlightens the necessity to take into account the CH$_4\nu_4$ band (7.7 µm/1305 cm$^{-1}$) and, to a lesser extent, C$_2$H$_2$, C$_2$H$_6$, and HCN between 600 and 900 cm$^{-1}$ at stratospheric temperatures.</figcaption></figure>

Apart from the TAM model (J. M. Lora et al. 2015), radiative transfer schemes were commonly based on C. P. McKay et al. (1989), which can explain commonly encountered difficulties to reproduce a correct thermal structure in the stratosphere. Motivated by this statement, along with the availability of recent spectroscopic values for methane as well as for trace compounds, we chose the approach to implement a new radiative transfer.

In this work, we use the flexible radiative transfer scheme developed for the Generic PCM, which has been used, among others, for early Mars (F. Forget et al. 2013), giants planets (S. Guerlet et al. 2014), Archean Earth (B. Charnay et al. 2013), and exoplanets (R. D. Wordsworth et al. 2011; M. Turbet et al. 2016). This radiative transfer scheme is based on a *k*-distribution method within the correlated-*k* approximation (R. Goody et al. 1989; Q. Fu & K. N. Liou 1992). When dealing with a high number of lines, this method has the strong advantage of requiring only the use of a limited number of narrow bands with enough accuracy for GCM applications. Here, we use 23 narrow bands in the visible range and 23 for infrared, covering the full spectral range [0.3–10000] µm. We use the cumulated distribution function *g* of the absorption within each narrow band. Integration is done on *k*(*g*) instead of *k*(*ν*) within a two-stream approximation radiative solver with multiple scattering (O. B. Toon et al. 1989). The detailed implementation of this method in the Titan PCM can be found in J. Vatant d'Ollone (2020).

### 2.4.2. Gaseous Absorption

Absorption from CH$_4$, C$_2$H$_2$, C$_2$H$_6$, and HCN are included. As the methane cycle is not included in these simulations, a vertical profile of the methane-mixing ratio is taken from H. B. Niemann et al. (2010) within the radiative transfer. For trace gases, mixing ratios are given by the chemical tracers computed by the photochemical solver and transported by the model in case they are treated as variable species (see Section 2.4.3). There is also the possibility to maintain each species as forced within the radiative transfer, to a reference profile derived from recent observational values. This may be used for simulations to study sensitivity to the variations of
<table><caption>**Table 1**
Vertical Profiles of the Radiatively Active Species Mixing Ratios in the Fixed-composition Radiative Transfer Configuration</caption><thead><tr><th>Pressure<br>(Pa)</th><th>10<sup>-1</sup></th><th>10<sup>0</sup></th><th>10<sup>1</sup></th><th>10<sup>2</sup></th><th>10<sup>3</sup></th><th>10<sup>3.3</sup></th><th>10<sup>3.7</sup></th><th>10<sup>4</sup></th><th>10<sup>5</sup></th><th>10<sup>6</sup></th></tr></thead><tbody><tr><td>CH<sub>4</sub></td><td></td><td></td><td></td><td></td><td>0.0148</td><td></td><td></td><td></td><td></td><td>0.0565</td></tr><tr><td>C<sub>2</sub>H<sub>6</sub></td><td></td><td>2.0 × 10<sup>-5</sup></td><td></td><td></td><td></td><td>2.0 × 10<sup>-6</sup></td><td>2.0 × 10<sup>-7</sup></td><td></td><td></td><td>10<sup>-8</sup></td></tr><tr><td>C<sub>2</sub>H<sub>2</sub></td><td></td><td>3.0 × 10<sup>-6</sup></td><td></td><td></td><td></td><td>3.0 × 10<sup>-7</sup></td><td>3.0 × 10<sup>-8</sup></td><td></td><td></td><td>10<sup>-11</sup></td></tr><tr><td>HCN</td><td></td><td>5.0 × 10<sup>-7</sup></td><td></td><td></td><td></td><td>5.0 × 10<sup>-9</sup></td><td></td><td></td><td>10<sup>-20</sup></td><td></td></tr></tbody></table>

each species independently. In this simulation, reference profiles are based on S. Vinatier et al.'s (2015) equinoctial values at the equator retrieved from CIRS observations for the stratosphere, and completed below the condensation level by values from previous model simulations. Profiles are smoothed vertically so that they do not show temporary local oscillations present in single observed profiles. The vertical distributions used in the fixed-composition radiative transfer configuration are summarized in Table 1.

First, high-resolution absorption spectra are computed using the *kspectrum* tool (V. Eymet et al. 2016) with spectroscopic data from the HITRAN2012 database (L. S. Rothman et al. 2013) as the input. These spectra are generated on a 2D grid with 10 bins in pressure ($P = [10^{-1}-10^{6}]$ Pa) times eight bins in temperature ($T = [60-200]$ K). To deal with the insufficient data from HITRAN in the near-infrared for methane, the line list is completed by ab initio calculations from M. Rey et al. (2018). An air-broadened Voigt line shape is used, corrected for methane in the troposphere by sub-Lorentzian profiles proposed by A. Campargue et al. (2012) for critical transparency windows, which is necessary to have a correct radiative balance in this region, as pointed out by M. Hirtzig et al. (2013). An updated version of the database, HITRAN2020 (I. E. Gordon et al. 2022), is available and will be used for future studies.

We use a modified 16-point Gauss-Legendre quadrature with eight Gauss points between 0 and 0.995 and eight Gauss points between 0.995 and 1 in order to have strong accuracy on high values of $g(k)$ (A. A. Lacis & V. Oinas 1991), which correspond to the intensity peak in line cores, especially crucial at low pressures in the upper layers of the model.

Collision-induced absorptions (CIAs)—also from the HITRAN2012 database (C. Richard et al. 2012)—are included for $N_2-N_2$, $N_2-CH_4$, $N_2-H_2$, $CH_4-CH_4$, and $H_2-H_2$, and are added as continua. The HITRAN2016 database (I. E. Gordon et al. 2017) brought some improvements for the methane line list but less than the M. Rey et al. (2018) one, and did not add any modifications to the CIAs used here (T. Karman et al. 2019).

### 2.4.3. Treatment of Variable Species

Since one of the main goals of this study is to investigate the sensitivity of the model to the variations of opacities induced by seasonal changes in the distribution of opacity sources, we chose to implement a method enabling an effective coupling: a method of recombination of the correlated-*k*. This method, first proposed by A. A. Lacis & V. Oinas (1991), is already applied in GCMs for exoplanets (P. Mollière et al. 2015; D. S. Amundsen et al. 2017). Using different values of several minor components in premixed spectra would induce a huge

increase in the size of the database to be computed upstream. The chosen method is much easier to produce and store, and also allows for more flexibility in terms of exploration of the sensitivity to variations of composition (the thermal structure in the low stratosphere being sensitive to the altitude of condensation, for instance).

Considering the mean transmittance $T_{\bar{\nu}}(u_x)$ in a given narrow band $\bar{\nu}$, for a single absorbing species *x* of surface density $u_x$ (density integrated over the considered layer), the strong hypothesis used in the recombination method is to neglect the potential cross correlation between the absorption of different species within a spectral interval, leading to

$$T_{\bar{\nu}}(u_x, u_y) = T_{\bar{\nu}}(u_x) \times T_{\bar{\nu}}(u_y). \quad (1)$$

This hypothesis leads to satisfying results in the case of random overlap of the different lines over a given narrow band (R. T. Pierrehumbert 2010). The narrowband discretization takes this into account. The random overlap approximation is statistically valid for low-pressure situations, while at high pressure, the large broadening of the lines, almost constant in the spectral interval, is also a situation satisfying this approximation. Consequences of correlations are also negligible in the case of very different optical depth between the different species, for example, when one minor species dominates its absorption band (J. Li & H. W. Barker 2005).

In the fully coupled simulations, the mixing ratio of $C_2H_2$, $C_2H_6$, and HCN are provided online by the advected tracers in the radiative transfer. The mixing ratio of $CH_4$ is set to H. B. Niemann et al.'s (2010) reference values. On every cell, the opacities of each compound are given by the local abundance of tracers, and a total gaseous opacity is given by successive recombination on top of the methane and rebinning down to a smaller number of Gauss points, where we end up having an equivalent local mixing. To allow computational cost saving, the operation is performed only for spectral intervals where the opacity of the added species is nonnegligible, which given the spectral discretization is the case for only a few narrow bands.

The limitations of this approach are illustrated in Figure 2. Heating and cooling rates are plotted for a given column of the Titan PCM (180°W, equator, at noon, northern spring equinox) at the first time step of a test run (so that temperature is the same for all cases), in four different configurations. In the “uncoupled” and “haze-coupled” cases, the premixed correlated-*k* tables for a fixed composition are used. In the “full-coupled” and “gas-coupled” cases, the recombination method is used, but the composition is set to the same profiles as for the premixed tables (i.e., according to Table 1). The differences between the two methods used for the gas opacity computation is mostly visible in the infrared range, with significant
<figure><img src="image_5.png" /><figcaption>**Figure 2.** Comparison of heating (left) and cooling (right) rates (in $10^{-5}$ K s$^{-1}$) computed in four different configurations, at the first step of a test run ($L_s = 0^{\circ}$, $0^{\circ}$N, $180^{\circ}$W). Differences are only due to the method used for computing the gas opacity with identical gas composition. Premixed method: uncoupled (green) and haze coupled (red). Recombination method: gas coupled (blue) and full coupled (black). The altitude scale on the right axis is approximate.</figcaption></figure>

differences ranging up to 10%–15%. However, this error is small compared to the uncertainty on the modeled abundances in the coupled simulations, where the abundances can be off by up to 1 order of magnitude when compared to observations. Thorough investigations to improve the application of this method in the case of Titan's stratosphere will be done in future studies.

### 2.4.4. Haze and Clouds

The importance of the coupling between microphysics and radiative transfer has been detailed in P. Rannou et al. (2004) and S. Lebonnois et al. (2003, 2009). As the details of the coupling between microphysics and radiative transfer in the Titan PCM are described in Paper II, we will only give a summary here.

For radiative coupling, optical properties of hazes and clouds (extinction cross section, asymmetry factor, and single-scattering albedo) are computed offline and read in a lookup table. The tables were created for spherical aerosols, fractal aerosols, and cloud droplets, calculating their properties over 33 radii, ranging from $10^{-9}$ to $10^{-5}$ m for aerosols in both modes, and from $10^{-7}$ to $10^{-3}$ m for cloud droplets. This implies a variation in properties depending on the type and size of the particles. Consequently, this leads to a variation in the optical properties of the atmosphere (albedo, refractive index, etc.) as a function of altitude and latitude, depending on the distribution of haze and clouds. Adding the gas contribution, we derive the total opacity, asymmetry factor, and single-scattering albedo for the whole cells.

The simulations uncoupled to haze were performed with a mean vertical profile of haze properties (opacity, asymmetry factor, and single-scattering albedo) used uniformly without horizontal or temporal variations from P. Lavvas et al. (2010).

### 2.5. Simulations Settings

Initial conditions were taken from simulations performed with the former version of the model (S. Lebonnois et al. 2012). An initial simulation coupled with the new version of the photochemistry (V. Vuitton et al. 2019), radiative transfer (A. A. Lacis & V. Oinas 1991; C. Richard et al. 2012; L. S. Rothman et al. 2013; M. Rey et al. 2018), and microphysical model of haze (J. Burgalat & P. Rannou 2017)

<figure><img src="image_6.png" /><figcaption>**Figure 3.** Vertical profile of the zonal wind in the Titan PCM at $L_s = 305^{\circ}$ and $10^{\circ}$S in black. The nudging profile is shown in red for $10^{\circ}$S. Huygens-observed profiles and zonal wind profiles from S. Lebonnois et al. (2012) are indicated, respectively, in dashed black and dotted black for comparison.</figcaption></figure>

and clouds (J. Burgalat et al. 2014) was run over 25 Titan years at standard resolution ($32 \times 48 \times 55$) to equilibrate the thermal structure.

We then ran four additional 5 yr simulations with different parameters: uncoupled simulation, labeled U (does not include any radiative feedback of trace compound variations nor haze variations), coupled gas simulation, labeled G (including radiative feedback of trace compound variations, but fixed haze opacities), coupled haze simulation, labeled H (including radiative feedback of haze variations, but with fixed gas opacities), and fully coupled gas and haze simulation, labeled F (haze and gas opacity variations are taken into account, as in the reference simulation of the first 25 yr). The results presented in this paper are the outputs of the last year of each simulation, in diurnal and zonal average. For the initial 25 yr run, the simulation has stabilized and shows little interannual variability. For instance, the temperature fluctuates around an average value with a maximum variation of 0.75% in the equatorial temperature profile during the same period over the past 2 yr. For the 5 yr sensitivity study runs, interannual stability is also generally satisfactory in the last couple of years.

A study of angular momentum within the model revealed that the integral of the angular momentum trends (which must be equal to 0) is negative in the lower troposphere, reducing the whole superrotation. Therefore, this version of the model is subject to a loss of angular momentum in the troposphere, which could lead momentum to leak from the stratosphere into the troposphere. To overcome this problem, we imposed a Newtonian relaxation of the zonal wind profiles over the first 23 layers of the atmosphere, from 0 to 65 km. The relaxation profiles are calculated for each latitude, based on the annual average tropospheric zonal wind of the Titan PCM with a factor of 2 to achieve an equatorial profile similar to that observed by Huygens during its descent. The relaxation time used is one Titanian year ($2.8 \times 10^8$ s). This forcing compensates for the numerical imbalance of angular momentum in the troposphere, allowing us to maintain a correct balance of momentum between the stratosphere and troposphere.

Figure 3 shows the tropospheric zonal wind profile at $L_s = 305^{\circ}$ and $10^{\circ}$S with the corresponding nudging profile. This approach is in line with expectations, given the very low seasonal variation of the zonal wind in the troposphere, resulting from relatively insignificant seasonal temperature changes. This approach also enables us to obtain a latitude-
<figure><img src="image_7.png" /><figcaption>**Figure 4.** Temperature profiles, in zonal and diurnal average, near the season and latitude of Huygens' landing ($L_s = 305°$; $10°S$) for all configurations (U: orange; F: red; H: green; G: blue), compared with the HASI profile of M. Fulchignoni et al. (2005). The altitude scale on the right axis is approximate.</figcaption></figure>

dependent wind profile, since zonal wind in the troposphere is more intense at the equator and midlatitudes than at the poles. In this way, nudging maintains a constant zonal wind in the troposphere and preserves superrotation in the stratosphere for long-term simulations.

### 3. Thermal and Dynamical Structure

The results presented in this work focus on the stratosphere. The thermal structure of the troposphere and the phenomena taking place there are not dealt with here, but are detailed in Paper II. A detailed study of the near-surface atmospheric region of a previous version of the IPSL Titan GCM has also been performed within a model intercomparison exercise by J. M. Lora et al. (2019).

#### 3.1. Temperature Structure and Seasonal Variations

To demonstrate the impact of the different couplings, the temperature profiles modeled for all the configurations are illustrated in Figure 4, around the season and latitude of the Huygens lander's arrival ($L_s = 305°$; $10°S$). In the fully coupled configuration (simulation **F**), recent improvements described in Section 2, particularly concerning radiative transfer, have significantly improved the temperature profile in the stratosphere. The stratopause is now simulated by the model at 30 Pa, at the same altitude as observed by the Huygens probe. However, the temperature predicted by the model is about 15 K colder than the temperature observed by Huygens (M. Fulchignoni et al. 2005) and 10 K below that found by CIRS at $15°$S (S. Vinatier et al. 2007). In the simulations coupled with haze and clouds (simulations **H** and **F**), the mesopause is observed around 10 Pa, which is significantly lower than expected (M. Fulchignoni et al. 2005). At $L_s = 305°$, the modeled mesosphere is highly compressed. This compression is caused by the detached haze layer of the model, which is too low compared to the observed one (see Paper II). It is the detached layer that produces the temperature inversion observed at 500 km in the observed thermal profile, and thus creates the mesopause. We will see below that, for other solar longitudes, the modeled mesosphere can be more extended or, conversely, can disappear completely where the stratopause merges with the warmer region of the detached layer.

The impact of coupling the gas composition is mainly visible between $2 \times 10^3$ and 10 Pa. The two simulations with decoupled gases (simulations **U** and **H**) are notably colder in the lower stratosphere, between $2 \times 10^3$ and $2 \times 10^2$ Pa, with respect to the coupled ones (**G** and **F**). The main differences between the coupling and decoupling of the gases are the mixing ratios used to compute their radiative feedback. The coupled cases use the mixing ratios computed by the Titan PCM, and the decoupled cases use fixed ratios based on observations (Table 1). The consequences of the use of computed composition are discussed in Section 3.3.

Below $2 \times 10^3$ Pa, the main difference is observed in the simulation coupled with gases (simulation **G**). This simulation shows a troposphere that is 5 K cooler than the other simulations, highlighting the importance of incorporating a cloud microphysics model in radiative transfer processes (see Paper II).

The seasonal evolution of the latitudinal temperature characteristics for simulation **F** is presented in Figure 5 and shows noticeable similarities with observations (R. K. Achterberg et al. 2008b; N. A. Teanby et al. 2019; S. Vinatier et al. 2020). Temperature variations are particularly pronounced at high latitudes ($>60°$). In midfall, the subsidence due to the descending branch of the global circulation causes adiabatic heating at high altitudes. The modeled temperature maps reproduce this polar heating: from the southern midfall to the end of southern winter, the high-latitude southern mesosphere and stratopause are warmer with respect to mid- and low latitudes. Lower in the atmosphere (below $10^2$ Pa), the meridional temperature gradient is reversed, with the polar latitudes colder than the equator.

At the equinoxes, the inversion of the circulation occurs. At $L_s = 180°$, the vertical gradient in temperature at high southern latitudes has substantially decreased. Before the following solstice, as the northern hemisphere is in midfall, the thermal structure has reversed, with the subsidence now in the northern hemisphere. Although these latitudinal features are similar to observations, it should be noted that the polar contrasts and structures in the winter stratosphere appear to be located at lower altitudes, specifically in the $10^0–10^2$ Pa range, compared to the observed data. This was already the case in S. Lebonnois et al. (2012), for a similar reason: the inaccurate representation of the haze production region, due to the too low top of the model, and the associated atmospheric structure above 10 Pa (roughly 300 km).

Despite these limitations, the altitude of the stratopause and its seasonal variations can be identified, as illustrated in Figure 6 over the last third of the southern fall, the entire southern winter, and almost all of the southern spring. During this period, the altitude of the stratopause changes with time and latitude. It is higher at high southern latitudes compared to the equator, until the southern spring equinox. After the reversal of the circulation (from $L_s = 200°$ and onward), the temperature of the stratopause increases at northern high latitudes, more and more as fall goes on. During the southern winter, the stratopause is consistently around 20 K warmer at the poles than at the equator, which compares positively with the observations collected by S. M. Hörst (2017).

In terms of latitudinal temperature contrasts, the thermal structure does not compare to observations as well as in the older 2D versions of the model of P. Rannou et al. (2004). Temperature contrasts are less marked than in the Voyager
<figure><img src="image_8.png" /><figcaption>**Figure 5.** Temperature maps, in zonal and diurnal average, as a function of latitude and pressure in simulation F. The altitude scale on the right axis is approximate.</figcaption></figure>

data, particularly around the equinox, where $\Delta T \approx 10$ K in the model against 15 K in the retrieved values by A. Coustenis & B. Bézard (1995). A more detailed comparison with the CIRS data is provided in Section 4.1.

## 3.2. Stratospheric Circulation

The seasonal variations of the zonally and diurnally averaged zonal wind predicted by the Titan PCM are shown in Figure 7, together with the meridional stream function. These seasonal variations are similar to the seasonal cycle described in S. Lebonnois et al. (2012), with some differences that are highlighted below.

First, the maximum velocity of the zonal wind jet ($L_s \approx 35^\circ$ south and $L_s \approx 215^\circ$ north) is stronger in the Titan PCM simulations compared to the previous version of the model (S. Lebonnois et al. 2012). The maximum predicted speed is approximately $240 \text{ m s}^{-1}$ ($L_s = 213^\circ$), corroborating observations that indicate a jet speed of around 190–250 m s$^{-1}$ (W. B. Hubbard et al. 1993; T. Kostiuk et al. 2005; B. Sicardy et al. 2006; R. K. Achterberg et al. 2008a; M. A. Cordiner et al. 2020; J. Sharkey et al. 2021), with peaks reaching up to $280 \text{ m s}^{-1}$ toward the end of the Cassini mission (S. Vinatier

et al. 2020). The intensity of the superrotation is highly sensitive to the model's parameters and configuration. Consequently, as we will discuss further, the new model couplings now enable the production of a strong jet ($\ge 200 \text{ m s}^{-1}$), consistent with observations. With the question of the maintenance of superrotation, which led us to impose tropospheric nudging, further investigations of the angular momentum budget for the Titan PCM will be necessary.

The zonal wind minimum observed in the lower stratosphere (M. K. Bird et al. 2005) was already predicted by the model (S. Lebonnois et al. 2012). However, this minimum is more pronounced compared to previous versions (see Figure 3). As this phenomenon occurs on spatial scales similar to the vertical resolution of our model (~10 km), we are still struggling to reproduce it. In addition, it should be noted that this result must be interpreted with caution, as the zonal wind forcing below 65 km may contribute to this outcome. This improvement may also be explained by a more accurate description of the temperature profile of this area provided by the microphysical cloud model (see Paper II). This region is characterized by very long radiative and dynamical timescales, and will be investigated in detail in a specific study focusing on this question.
<figure><img src="image_9.png" /><figcaption>**Figure 6.** Qualitative position of the stratopause (white dotted line) from before the summer solstice to the winter solstice ($L_s = 60°$ to $L_s = 260°$) in simulation **F**, plotted over a zonally and diurnally averaged temperature field, between 1 and 200 Pa.</figcaption></figure>

Above 10 Pa, the anomaly in the thermal structure due to the presence of the detached haze layer (DHL) induces a shear between the stratospheric circulation cell and a mesospheric one, here again more clearly marked around the 4–10 Pa layer than in our previous study (S. Lebonnois et al. 2012).

As demonstrated in previous works (P. Rannou et al. 2004; S. Lebonnois et al. 2009), the zonal circulation is sensitive to coupling. Figure 8 presents the seasonal variation of the zonal wind at two different altitudes, namely 1 and $10^2$ Pa, for the four configurations. These correspond roughly to the center of the circulation cells.

The impacts of coupling on the circulation are quite similar to those on the thermal structure. The more contrasted the thermal structure, the higher the intensity of the wind. Previous studies, like P. Rannou et al. (2004), S. Lebonnois et al. (2009), and J. Vatant d'Ollone (2020), have pointed out that successive couplings strengthen the zonal wind, and especially the coupling of the haze. Here, in line with the results for the thermal structure, the full coupling is indeed the one with the strongest impact on the reinforcement of the circulation, especially when compared to the coupling of the gases. There

is a difference of $25 \text{ m s}^{-1}$ between the peak velocity displayed at 1 Pa for simulation **F** (Figure 8(d)) versus simulation **G** (Figure 8(b)), and $80 \text{ m s}^{-1}$ when compared to simulation **U** (Figure 8(a)). When comparing simulations **H** (Figure 8(c)) and **F** (Figure 8(d)), we observe that the coupling with haze and gases intensifies temperature contrasts. This results in an increase of several $\text{m s}^{-1}$ in the peak zonal wind speed compared to coupling with haze only.

As detailed in studies based on different GCMs (C. E. Newman et al. 2011; S. Lebonnois et al. 2012; N. T. Lewis et al. 2023), wave activity is an important aspect of the stratospheric superrotation, especially present at high latitudes during fall and winter. The intensity of the circulation and wave activity are linked, and as such coupling also reinforces the wave activity. This impact is clearly visible in Figure 8.

### 3.3. *Seasonal Variations of Composition*

Figure 9 illustrates the seasonal evolution of the main trace gases included in the Titan PCM ($C_2H_2$, $C_2H_6$, HCN). These
<figure><img src="image_10.png" /><figcaption>**Figure 7.** Seasonal evolution of the zonally and diurnally averaged zonal wind in simulation F. The stream function of the meridional circulation is represented with dark contours, with solid lines indicating a clockwise circulation and dotted lines a counterclockwise circulation. The altitude scale on the right axis is approximate.</figcaption></figure>

chemical compounds are produced photochemically in the upper layers of the atmosphere, and their vertical distribution is controlled by vertical diffusion (>500 km). They are then transported down into the stratosphere by dynamical transport. In the sponge layer (~500 km), these species are kept constant at values close to those observed, as described in Section 2.2. In the vertical domain covered by the Titan PCM (0–500 km), they are also transported by the meridional circulation, leading to enrichment over the winter pole, as detailed in previous works (S. Lebonnois et al. 2001; F. Hourdin et al. 2004).

The temperature below 100 km triggers the condensation of species. Consequently, condensation occurs in the lower stratosphere, close to the tropopause, explaining the marked decrease observed around $10^3$ Pa in the mixing ratio of $C_2H_2$, $C_2H_6$, and HCN. Ethylene is the only species considered in the model that does not condense (E. L. Barth 2017), but it is not a species of interest in this study.

Figure 10 shows the mixing ratio profiles of the radiatively active species in the Titan PCM ($C_2H_2$, $C_2H_6$, HCN) for the four configurations, compared to corresponding CIRS data (S. Vinatier et al. 2015). $CH_4$ is considered in the calculation of radiative transfer with the profile from H. B. Niemann et al. (2010). Without the forcing of photochemistry in the sponge layer, all species would be underestimated, except for HCN between 0.1 and 1 Pa, which aligns well with observations in all configurations. Now, species are overestimated in the sponge layer due to forcing (above 1 Pa), but they match observations in the lower layers. The main discrepancy in the profile shapes occurs between 1 and 20 Pa: for all compounds and simulations, a steep vertical gradient appears in the modeled profiles but is absent from the observations. This discrepancy could be due to meridional transport taking effect at the top of the PCM. The vertical extension of the meridional transport cell is limited by the top of the model at 500 km. As a result, close to the equator, the maximum meridional flux of
<figure><img src="image_11.png" /><figcaption>(a) Simulation **U**, uncoupled</figcaption></figure>

<figure><img src="image_12.png" /><figcaption>(b) Simulation **G**, gas coupled</figcaption></figure>

<figure><img src="image_13.png" /></figure>

<figure><img src="image_14.png" /><figcaption>(c) Simulation **H**, haze coupled</figcaption></figure>

<figure><img src="image_15.png" /></figure>

<figure><img src="image_16.png" /><figcaption>(d) Simulation **F**, fully coupled</figcaption></figure>

<figure><img src="image_17.png" /></figure>

<figure><img src="image_18.png" /><figcaption>**Figure 8.** Evolution of the zonally and diurnally averaged zonal wind as a function of latitude and solar longitude at 110 Pa (left) and 1 Pa (right), for the four coupling configurations.</figcaption></figure>

species occurs between 1 and 10 Pa. If the model were extended vertically, the cell would stretch beyond 500 km, and the maximum meridional flux would be observed at a higher altitude. The lack of observations at these altitudes prevents us from confirming or refuting this hypothesis. However, the observational constraints on the latitudinal composition

contrasts confirm that such a flux exists and is well reproduced by the model (Figure 9). However, this bias could also be related to the interface between the photochemistry and the model dynamics. The vertical mixing could not be strong enough to compensate for the decrease in winds at the top of the model due to the sponge layer. Consequently, the dynamics
<figure><img src="image_19.png" /><figcaption>**Figure 9.** Seasonal evolution of zonally and diurnally averaged mixing ratios for minor species radiatively active in the Titan PCM, in simulation **F**. From top to bottom: C₂H₂, C₂H₆, HCN. The stream function of the meridional circulation is represented with dark contours, with solid lines indicating a clockwise circulation and dotted lines a counterclockwise circulation. The black dashed horizontal line signals the maximum altitude for the dynamical transport; above, there is only vertical diffusion.</figcaption></figure>

do not sufficiently mix the species, explaining why they do not reach the lower layers and thus leading to a decrease in mixing ratios.

These compounds are infrared radiators that act as cooling agents (B. Bézard et al. 1995, 2018). This explains the shift observed in Figure 4 between the profiles with coupled gases (G and F) and those decoupled (U and H). When the gases are decoupled, the abundances used for radiative feedback calculations are based on those observed by CIRS (S. Vinatier et al. 2015). These observed abundances are higher than those predicted by the model below 10 Pa (300 km), resulting in more efficient cooling. This bias is compensated through the adjustment of the optical properties of the haze. This should be corrected when the simulations will be able to better reproduce the species' observed profiles. Still, taking into account the seasonal variations of the composition is important because of the enrichments seen in the winter polar stratosphere, as demonstrated by N. A. Lombardo & J. M. Lora (2023).

## 4. Peculiar Features in the Polar Regions

This section focuses specifically on the polar regions, where the seasonal variations are the most pronounced. In particular, the moment of the reversal of the circulation is studied in more detail. This period is characterized by a complex evolution of the thermal structure. Thanks to the Cassini mission, a number

of observations are available to confront the results of the model. Here, the evolution of the thermal structure, from high to low stratosphere, as well as the presence of a polar shoulder in the observed temperature profiles of the winter stratosphere, are investigated.

### 4.1. *Comparison to Observations*

The passage of the vernal equinox has been extensively documented during the Cassini mission, notably through CIRS limb observations. Figures 11 and 12 compare these observations to the results obtained by the Titan PCM (simulation **F**) for the thermal structure and the C₂H₂ distribution, respectively, as retrieved by S. Vinatier et al. (2015, 2020). The maps selected here span from the end of the northern winter to the summer solstice.

During the northern winter (Figure 11—Ls = 301°–344°), the altitude and temperature of the stratopause in the observations match what has been detailed in the previous section: a warmer stratopause positioned a pressure level an order of magnitude smaller at the winter pole with respect to the other latitudes. Then, from the equinox, the reversal of the circulation occurs. The temperature and the altitude of the northern stratopause decrease while the subsidence initiates in the southern hemisphere, with an adiabatic heating starting in 2015 (Figure 11—Ls = 66°–67°). By the time of the solstice,
<figure><img src="image_20.png" /><figcaption>**Figure 10.** Modeled mixing ratio profiles of C₂H₂ (left), C₂H₆ (center), and HCN (right), close to the vernal equinox, at the equator (top row, 4°N) and low latitude (second row, 25°N), for four simulation settings: **U** (orange), **H** (green), **G** (blue), and **F** (red). The profiles are compared to CIRS observations at the same season and latitudes from S. Vinatier et al. (2015). The altitude scale on the right axis is approximate.</figcaption></figure>

the circulation has completely reversed. This is corroborated by the C₂H₂ distribution maps, with a clear enrichment visible at the north pole up to 2012 February (Figure 12—Ls = 29°–31°), and traces of enrichment apparent in the southern hemisphere as soon as 2011 May (Figure 12—Ls = 21°–23°). Following the conclusion of S. Vinatier et al. (2015), the circulation reverses, from one cell in the northern winter, to a two-cell temporary configuration that persists up to 2010 June, and finally a completely reversed one-cell circulation is established in 2011 June.

On the model side, in spite of the biases detailed in Section 3.1, the model reproduces quite well the passage of the equinox in terms of thermal structure. We observe a decrease in temperature above the north pole and a significant cooling at the south pole at a period similar to that of the observations (Ls = 21°–31°). Starting from Ls = 66°–67°, the stratopause becomes visible. Although too cold in the northern hemisphere, its behavior is similar to observations, with a higher altitude at the south pole (0.9–8 Pa) compared to the rest of the atmosphere (~30 Pa). The model also reproduces the intense cooling in the lower southern stratosphere at high latitudes evident in 2015 (Ls = 66°–67°). The modeled temperatures align here with those from CIRS data. However, the vertical development of this cold region is limited in the model, constrained by the pressure level of the stratopause, which is situated lower than observed.

Regarding the composition, the polar enrichments in the model remain slightly weaker compared to observations, despite the forcing in the sponge layer. A notable difference is the appearance of enrichment at the south pole in the model as early as Ls = 4°–5°, whereas in the observations, this occurs only in 2011 May (Ls = 21°–23°). Nevertheless, the persistence of an enriched cell at the north pole in spring, below the 10 Pa level, is well replicated in Titan's PCM. This sustained enrichment is explained by the presence of a wind cell that confines the species to the pole even after the circulation reversal. However, this cell appears lower in the model than in the observations.

The enrichment in the model also follows the same latitudinal expansion as in the observations: limited above 80°S at the beginning of southern fall, it expands to reach 60°S by the time of the solstice, at a rate closely following the data. This expansion of the enriched area is also consistent with the development of the polar vortex described by N. A. Teanby et al. (2017), responsible for trapping the compounds at high latitudes.

## 4.2. Temperature Minimum in the High Stratosphere

The model is now capable of replicating the observed cooling in the upper atmosphere of the south pole, above 3 Pa, between 2011 May and 2012 February (see Figure 11). This temperature decrease coincides with the increase in the
<figure><img src="image_21.png" /><figcaption>**Figure 11.** Evolution of the thermal structure during northern spring in CIRS data and in the Titan PCM. The temperature fields obtained by CIRS limb observations (left) are compared to the modeled zonally and diurnally averaged temperature maps (right) from Titan PCM simulation F. CIRS observations from northern winter, before 2009 August to 2012 February, come from S. Vinatier et al. (2015), and the remaining figures are adapted from S. Vinatier et al. (2020). The altitude scale on the right axis is approximate.</figcaption></figure>

concentration of photochemical compounds in the same region and period (see Figure 12). S. Vinatier et al. (2015) and N. A. Teanby et al. (2017) suggest that the dramatic increase in concentrations is indeed responsible for this cooling and that the radiative effect of these compounds is strong enough to surpass the adiabatic warming caused by the subsidence. Thus, the model results support this hypothesis. We propose an additional cause for this marked cooling, which could reinforce

the effect of compound enrichment and intensify this phenomenon. The collapse of the DHL and its small absorbing aerosols following the equinox reduces absorption and thus heating. The disappearance of this temperature minimum from 2015 March onward ($L_s = 66°$) is attributed to the continuous strengthening of subsidence through the southern fall, and thus adiabatic heating, which outweighs the cooling effect due to the infrared emission from chemical species.
<figure><img src="image_22.png" /><figcaption>**Figure 12.** Evolution of the C₂H₂ distribution during northern spring in CIRS data and in the Titan PCM. The distributions obtained by CIRS limb observations (left) are compared to the modeled zonally and diurnally averaged mole fraction maps (right) from Titan PCM simulation F. CIRS observations from northern winter, before 2009 August to 2012 February, come from S. Vinatier et al. (2015), and the remaining figures are adapted from S. Vinatier et al. (2020). The stream function is indicated in black contours on the PCM maps, along with the upper limit of the dynamical transport in line marks. The altitude scale on the right axis is approximate.</figcaption></figure>

Furthermore, the cooling of the stratosphere below 300 km at high latitudes is reproduced. The detail and driving force behind this cooling can thus be investigated and compared to what occurs in the high stratosphere. Figure 13 shows, from top to bottom, the evolution of temperature, heating rates, and extinction of the haze (much more dominant in the visible spectrum than gases) during southern fall. The focus is on the period $L_s = 0°–90°$ and area 20–80 Pa, $> 60°S$ where the

observed minimum temperature is most pronounced. By comparing these quantities, it is possible to qualitatively determine which process is driving the temperature variations observed. A distinction is made between the dynamic heating rate and the radiative heating rate, itself decomposed into its short- and long-wave components. Physically, the dynamic heating rate is related to the subsidence: the stronger the subsidence, the stronger the adiabatic heating, and therefore the
<figure><img src="image_23.png" /><figcaption>**Figure 13.** Seasonal evolution between $L_S = 0°$ and $L_S = 90°$ (southern fall) of the temperature (top), heating rates (middle), and haze extinction in the visible wavelength (bottom), averaged zonally and diurnally and between 20 and 80 Pa for latitudes higher than 60°S. The heating rates are expressed in K/Titan day. The radiative heating rate (red solid line) and its decomposition in the visible domain (short wave, red dotted line) and infrared (long wave, red dashed line), is compared to the dynamic heating rate (blue). The net radiative and long-wave heating rates are negative and are shown here in absolute values.</figcaption></figure>

stronger the dynamic heating. On the other hand, the short-wave component of the radiative heating rate corresponds to solar heating, while the long-wave component represents the cooling contribution of infrared-radiating elements such as haze and gases. As such, the net radiative and long-wave heating rates are negative, but as the interest here lies in their intensity, they are presented in absolute values.

The short-wave component is of little interest: as fall progresses, the southern high latitudes enter the polar night, and solar heating in the pressure domain considered decreases to reach zero around midseason. The global decrease in temperature during the first half of the season is surely in part due to this, but this is not the only reason. Up to roughly $L_S = 38°$, there is a stark decrease in the temperature (around 30 K), which coincides with a steady increase in the extinction and an increasing long-wave (and radiative) heating rate. The latter is the sign that haze and gases are accumulating in the high-latitude region, transported by the subsidence. The radiative heating rate is slightly superior to the dynamic one, so the radiative cooling effect due to the accumulation of haze and gas is stronger than the subsidence (that has also just begun). The temperature thus decreases.

Between $L_S = 38°$ and $L_S = 58°$, the extinction stabilizes and exhibits long-term oscillations. The radiative heating rate and temperature also tend to stabilize during this period. The radiative heating rate continues to slightly exceed the dynamic heating rate, except for intermittent and sporadic peaks. These peaks are caused by planetary wave activity that redistributes angular momentum, coinciding with temperature spikes. During these episodes, subsidence temporarily surpasses the radiative effect, leading to an increase in temperature.

Starting from $L_S = 58°$, the subsidence intensifies, and overall extinction decreases as the amount of haze follows the subsidence and moves deeper into the stratosphere. The dynamic heating rate surpasses the radiative rate most of the time, and thus begins leading the evolution of the temperature. The subsidence has become stronger than the radiative effect,

causing temperature increases due to adiabatic heating. The dynamic heating rate continues to rise until the end of the season, indicating that subsidence is continually strengthening. The period during which the radiative effect, particularly from long wavelengths, is stronger than the adiabatic heating from subsidence corresponds to the period when the minimum temperature is observed at high altitudes (Figure 11). It is likely that the same phenomenon is at play here, supporting the hypothesis that this temperature minimum is due to the radiative effect of compounds transported by subsidence.

### 4.3. Polar Shoulder in the Low Stratosphere

In the low polar stratosphere, strong cooling is observed during winter. At the end of winter, P. J. Schinder et al. (2012) studied radio-occultation profiles, showing a region of destabilization, also named “polar shoulder,” in the vertical temperature profile, present at high latitudes only. Figure 14 shows these profiles, compared to the results of the Titan PCM simulations at a similar time period and latitude. The radio-occultations show a strong inversion in the temperature profile around $10^3$ Pa at these polar latitudes, slightly more pronounced at $80°$N than at $75°$N. However, at the season and latitude of the observations, the model does not reproduce this destabilization, no matter the couplings used. A temperature inversion is present in the simulation coupled with haze (simulation **H**). However, it is not sufficiently pronounced and occurs too low in altitude (~2 × $10^3$ Pa) to allow for a thorough analysis. Cassini/CIRS data in the far-infrared show broad emission features at $220 \text{ cm}^{-1}$, as yet unidentified, which could have a strong effect on the radiative cooling of the winter polar lower stratosphere (R. de Kok et al. 2007; C. M. Anderson et al. 2018). C. M. Anderson et al. (2018) suggest that this emission originates from a stratospheric ice cloud, potentially composed of pure HCN or complex nitriles. Our model reproduces such a stratospheric ice cloud composed of hydrocarbons and HCN (see Paper II), but not the $220 \text{ cm}^{-1}$ emission, which could partly explain why it does not reproduce the polar shoulder.

The radio-occultation measurements available at high latitudes only cover the very end of the northern winter. CIRS data, however, comprise a wide range of solar longitudes during the southern fall and winter, even though they do not probe the atmosphere as deep as the radio-occultations. Based on S. Vinatier et al. (2010) and C. Mathé et al. (2020), we compiled an array of temperature profiles spanning from the northern vernal equinox to the following solstice, corresponding to southern fall. The focus was put on latitudes similar to those where radio-occultation data were available, i.e., $75°$S and $80°$S. Those observational profiles are shown in Figure 15. The radio-occultation profiles are also plotted as a point of reference for the end of the winter, even though they are from the opposite hemisphere.

Few profiles are available at $75°$S corresponding to the end of fall. They are discussed below together with the radio-occultation profiles at $80°$S. At $L_S = 5°$, above the stratopause visible at 20 Pa, the mesosphere is subject to heating located at 1 Pa, certainly associated with the beginning of subsidence. The two following profiles ($L_S = 25°$ and $L_S = 30°$) present only small shoulders visible at $6 \times 10^1$ and $4 \times 10^2$ Pa. It is unclear if they are related to the polar shoulder. The strong cooling discussed previously is visible above 10 Pa at $L_S = 30°$.
<figure><img src="image_24.png" /><figcaption>**Figure 14.** Polar shoulder at high latitudes at the end of northern winter from Cassini radio-occultations in black (P. J. Schinder et al. 2012), compared to modeled thermal profiles for the four configurations: **U** in orange, **H** in green, **G** in blue, and **F** in red. Only the two latitudes in which the polar shoulder is visible in the observations are represented: 80°N (solid lines) and 75°N (dashed-dotted lines). The altitude scale on the right axis is approximate.</figcaption></figure>

<figure><img src="image_25.png" /><figcaption>**Figure 15.** Temperature profiles extracted from radio-occultation and spectroscopic data from Cassini around 75° (top) and 80° (bottom). Radio-occultation profiles (black dashed-dotted lines) from P. J. Schinder et al. (2012) correspond to the end of northern winter. CIRS profiles (colored solid lines) from S. Vinatier et al. (2010) and C. Mathé et al. (2020) correspond to the southern fall. The altitude scale on the right axis is approximate.</figcaption></figure>

After that, there is another gap of a third of a season in the observations, as the next batch of profiles focuses between $L_s = 66°$ and the solstice at $L_s = 90°$. These profiles display significant cooling of the stratosphere, and a distinct polar shoulder. A drop in altitude of this polar shoulder can be clearly followed, starting out at $10^2$ Pa for $L_s = 66°$, and reaching an altitude close to that of the radio-occultation ($\sim 10^3$ Pa) at $L_s = 90°$. The profile at $L_s = 72°$ should be considered carefully, however, as it is closer to $85°$S, and a $5°$ gap in latitude may lead to nonnegligible variations. Also, because of the very low temperature in the deep stratosphere (deeper than $10^2$ Pa), the information on temperature retrieved from CIRS limb spectra is very limited.

From these profiles, it may be inferred that the temperature minimum, referred to as the polar shoulder, forms after the opacity enrichment starts to cool the mesosphere and stratopause area, followed by the cooling of the stratosphere.

<figure><img src="image_26.png" /><figcaption>**Figure 16.** Evolution of the zonally and diurnally averaged temperature in the low stratosphere ($10^1-10^3$ Pa), in the northern hemisphere (a), and in the southern hemisphere (b), during fall and winter (simulation **F**). The altitude scale on the right axis is approximate.</figcaption></figure>

As the season progresses, its altitude decreases. However, the destabilization is not as pronounced in CIRS observations at the solstice as in the radio-occultations at the end of winter. The vertical resolution of CIRS observations is approximately 30 km, with vertical profile smoothing applied during inversions. It is not possible to obtain as pronounced a polar shoulder as seen in radio-occultation data.
<figure><img src="image_27.png" /><figcaption>**Figure 17.** Seasonal evolution of the temperature in the stratosphere at 80°N, in zonal and diurnal average. Each panel corresponds to a coupling configuration: upper left = **U**, upper right = **H**, lower left = **G**, lower right = **F**. The altitude scale on the right axis is approximate.</figcaption></figure>

<figure><img src="image_28.png" /></figure>

4.4. Proposed Interpretation

The modeled seasonal evolution of the temperature in the low stratosphere region during both fall and winter, for both hemispheres, is presented in Figure 16 for simulation **F**. Regions showing a temperature minimum are initiated in the upper stratosphere, just before the middle of the season (Ls = 210°–220° in the north and Ls = 30°–40° in the south). This temperature minimum decreases rapidly in altitude during fall, which corroborates the rapid descent of this inversion in the stratosphere suggested from the observations.

A second, less intense cooling occurs around 40 Pa, taking place at Ls = 240° in the northern hemisphere and at Ls = 60° in the southern hemisphere. This secondary cooling follows the trend of the first and descends into the cold layers of the lower stratosphere as fall progresses. A slightly warmer layer becomes trapped between these two cold layers, creating a temperature inversion that extends into the lower stratosphere. This phenomenon could explain the thermal destabilization that produces the polar shoulder. However, the model does not maintain this minimum until the end of winter, as it disappears around the solstice, explaining the absence of an inversion in the modeled profile during the radio-occultation season at the end of winter (Figure 14). These maps also show the descent in altitude of the stratopause following the cooling of the upper layers.

The influence of the various couplings is illustrated in
Figure 17, which shows the stratospheric temperature evolution

with season at 80°N during a full year. It appears clearly that
the transition occurring after the fall equinox, as previously
described, is strongly reinforced by couplings, especially
coupling with the haze distribution. More pronounced cooling
and the presence of a temperature minimum are observed in
simulations **H** and **F**.

The scenario deduced from the model behavior, together
with observational data sets, is that after the fall equinox and
the reversal of the circulation, the descent of opacity sources
(haze and minor compounds) in the high polar stratosphere
initiates cooling of the temperature in the stratopause region.
This creates a temperature minimum that descends quickly as
fall progresses, as the opacity enrichment descends deeper in
the stratosphere. A second cooling phase associated with a
subsequent wave of enrichment appears in the latter half of fall,
contributing to a further decrease in temperature. This situation
tends to stabilize after winter solstice, with a temperature
inversion that should be maintained in the lower polar
stratosphere during all winter. Unfortunately, due to its vertical
extent biases, the Titan PCM does not maintain it past winter
solstice.

Test simulations using a higher horizontal resolution (96
longitudes × 96 latitudes) confirm this behavior, with a clearer
temperature minimum that tends to be maintained past the
winter solstice.
## 5. Conclusions

The Titan PCM, the new version of the IPSL Titan GCM, was presented in this work. Details on the modifications to physical processes were given, with a particular emphasis on the new radiative transfer scheme. Based on the *k*-distribution method within the correlated-*k* approximation, it allows a better representation of the full spectrum of gas opacities, and allows us to take into account composition variability. The reaction scheme of the photochemical solver has been updated, along with technical improvements, as the photochemistry is now being solved from the surface up to 1300 km. The companion paper (Paper II) presents the new microphysics scheme for haze and clouds included in the Titan PCM, along with a study of Titan's haze and cloud cycles.

Improvements have been made regarding the thermal structure and circulation of the stratosphere. Based on different coupling configurations (including or not coupling with minor gases and/or coupling with haze), this study focused on the middle atmosphere, from $10^3$ to $10^{-1}$ Pa, and on the processes occurring in the polar regions during fall and winter.

1. The temperature structure has been enhanced compared to previous model versions (S. Lebonnois et al. 2012; B. Charnay et al. 2013), with the stratopause now visible for much of the year in the model up to approximately 10 Pa. However, the detached haze layer is too low compared to observations (see Paper II), inducing heating above the stratopause and an artificially compressed atmosphere, a bias already present in previous versions of the IPSL Titan GCM. Nonetheless, the seasonal variations of the stratopause appear to be well represented.

2. The use of a coupled model with radiative feedback leads to a representation of the equinoctial circulation inversion at high altitudes that is in better agreement with observations. The feedback from the haze and photochemical species amplifies temperature contrasts, especially during the polar night when they act as an infrared radiator. At the same time, these temperature contrasts enhance atmospheric dynamics. The Titan PCM super-rotation jet now reaches values close to observations, exceeding $200 \text{ m s}^{-1}$. The improved description of the circulation allows us to obtain a distribution of atmospheric components that better matches observations, leading to a more precise study of the haze cycle and its detached layer (see Paper II).

3. At high altitudes, the heating induced by aerosols in the detached layer maintains the atmosphere in a configuration resembling that of the late season. Hence, the radiative feedback from the haze could explain the delay highlighted by CIRS observations (S. Vinatier et al. 2015; A. Coustenis et al. 2018) toward the equinox.

4. The cooling of the mesosphere and stratopause region observed in early 2012 ($L_s = 30^\circ$) is for the first time reproduced in a GCM. The mechanism that may explain it is demonstrated in the upper stratosphere: the increase of opacity and infrared emission that occur at that time cool the atmosphere, despite an increase of adiabatic heating due to subsidence.

5. Using some CIRS profiles obtained during fall in the southern polar regions and the behavior of the Titan PCM during this season, it has been demonstrated that the polar shoulder, observed in the lower stratosphere at the end of

winter with radio-occultations, may be present during the entire winter and may start in the upper stratosphere in the first half of fall. Our study suggests that the significant cooling of the stratopause at the beginning of fall is linked to the polar shoulder observed at the end of winter. This connection is due to the descent of opacity sources, causing cooling from the mesosphere and stratopause into the lower stratosphere during fall, creating an inversion that lasts through winter.

6. As in previous works (S. Lebonnois et al. 2003, 2009; P. Rannou et al. 2004), the role of couplings in modeling the atmospheric system is once more highlighted.

This study represents a significant step forward in understanding Titan's complex stratospheric thermal structure. Our investigation into seasonal variations, particularly in polar regions, has revealed valuable insights into previously unexplained phenomena, such as the destabilization of the lower polar winter stratosphere. While our model reproduces these phenomena to a certain extent, it also highlights the need for further refinement, particularly in increasing the model's upper limit. Removing the thin layer approximation appears to be a key point on which modeling efforts should be focused to keep progressing in the apprehension of the seasonal effects in the stratosphere and above. Addressing this limitation would require extensive work rewriting the primitive meteorological equations integrated into the dynamical core.

Including ionic chemistry is also a possible way to refine the present photochemical model toward vertical extent purposes. However, even if the previous challenges were to be tackled, vertical extent would still come up against the strongly limiting parallel-plan radiative transfer, leading to polar night not being able to be lit up at high altitudes in the model, which will likely last without a consequent leap forward in modeling. A project to develop a full 3D radiative transfer model, surpassing the limitations of plane-parallel calculations, is currently underway and will be the subject of future studies.

In conclusion, this study has exploited the new Titan PCM model to analyze previously poorly understood processes, such as the overall control of the thermal structure, the impact of the coupling on equinoctial circulation reversal, and the radiative destabilization of the lower polar stratosphere. By carrying out 3D simulations with the Titan PCM, we also present an analysis of the thermal structure, while comparing our results with the data collected by the Cassini mission. These investigations provide essential constraints for future observations of Titan. A study combining the PCM with a mesoscale model is underway to constrain the atmospheric and surface properties of the Dragonfly mission site.

## Acknowledgments

The authors thank Juan Lora, Bruno Bézard, Sandrine Guerlet, and Martin Turbet for their very useful discussions about radiative transfer modeling, as well as Ehouarn Millour for technical support on the GCM. This work has benefited from computation facilities provided by the Grand Équipement National de Calcul Intensif (GENCI) where Titan PCM simulations have been performed on High-Performance Computing (HPC) facilities under allocation A0120110391. The mesocenters ESPRI from the Institut Pierre-Simon Laplace and ROMEO from the Université de Reims Champagne-Ardenne were also used for Titan PCM simulations.
ORCID iDs

Bruno de Batz de Trenquelléon https://orcid.org/0000-0001-8988-8613

Sébastien Lebonnois https://orcid.org/0000-0002-2390-8164

Pascal Rannou https://orcid.org/0000-0003-0836-723X

Sandrine Vinatier https://orcid.org/0000-0001-5541-2502

References

Achterberg, R. K., Conrath, B. J., Gierasch, P. J., Flasar, F. M., & Nixon, C. A. 2008a, *Icar*, 194, 263

Achterberg, R. K., Conrath, B. J., Gierasch, P. J., Flasar, F. M., & Nixon, C. A. 2008b, *Icar*, 197, 549

Achterberg, R. K., Gierasch, P. J., Conrath, B. J., Flasar, F. M., & Nixon, C. A. 2011, *Icar*, 211, 686

Amundsen, D. S., Tremblin, P., Manners, J., Baraffe, I., & Mayne, N. J. 2017, *A&A*, 598, A97

Anderson, C. M., Samuelson, R. E., & Nna-Mvondo, D. 2018, *SSRv*, 214, 125

Barth, E. L. 2017, *P&SS*, 137, 20

Bézard, B. 2014, *Icar*, 242, 64

Bézard, B., Coustenis, A., & McKay, C. P. 1995, *Icar*, 113, 267

Bézard, B., Vinatier, S., & Achterberg, R. K. 2018, *Icar*, 302, 437

Bird, M. K., Allison, M., Asmar, S. W., et al. 2005, *Natur*, 438, 800

Burgalat, J., & Rannou, P. 2017, *JAerS*, 105, 151

Burgalat, J., Rannou, P., Cours, T., & Rivière, E. D. 2014, *Icar*, 231, 310

Cabane, M., Chassefière, E., & Israel, G. 1992, *Icar*, 96, 176

Campargue, A., Wang, L., Mondelain, D., et al. 2012, *Icar*, 219, 110

Charnay, B., Forget, F., Wordsworth, R., et al. 2013, *JGRE*, 118, 10414

Charnay, B., & Lebonnois, S. 2012, *NatGe*, 5, 106

Cordiner, M. A., Garcia-Berrios, E., Cosentino, R. G., et al. 2020, *ApJL*, 904, L12

Corlies, P., Hayes, A. G., Birch, S. P. D., et al. 2017, *GeoRL*, 44, 11754

Cours, T., Burgalat, J., Rannou, P., et al. 2011, *ApJL*, 741, L32

Coustenis, A., & Bézard, B. 1995, *Icar*, 115, 126

Coustenis, A., Jennings, D. E., Achterberg, R. K., et al. 2016, *Icar*, 270, 409

Coustenis, A., Jennings, D. E., Achterberg, R. K., et al. 2018, *ApJL*, 854, L30

Crespin, A., Lebonnois, S., Vinatier, S., et al. 2008, *Icar*, 197, 556

de Batz de Trenquelléon, B., Rannou, P., Burgalat, J., Lebonnois, S., & Vatant d'Ollone, J. 2025, *PSI*, 6, 79

de Kok, R., Irwin, P. G. J., Teanby, N. A., et al. 2007, *Icar*, 191, 223

de Kok, R. J., Teanby, N. A., Maltagliati, L., Irwin, P. G. J., & Vinatier, S. 2014, *Natur*, 514, 65

Eymet, V., Coustet, C., & Piaud, B. 2016, *JPhCS*, 676, 012005

Flasar, F. M., Achterberg, R. K., Conrath, B. J., et al. 2005, *Sci*, 308, 975

Flasar, F. M., Kunde, V. G., Abbas, M. M., et al. 2004, *SSRv*, 115, 169

Forget, F., Wordsworth, R. D., Millour, E., et al. 2013, *Icar*, 222, 81

Fray, N., & Schmitt, B. 2009, *P&SS*, 57, 2053

Friedson, A. J., West, R. A., Wilson, E. H., Oyafuso, F., & Orton, G. S. 2009, *P&SS*, 57, 1931

Fu, Q., & Liou, K. N. 1992, *JAtS*, 49, 2139

Fulchignoni, M., Ferri, F., Angrilli, F., et al. 2005, *Natur*, 438, 785

Goody, R., West, R., Chen, L., & Crisp, D. 1989, *JQSRT*, 42, 539

Gordon, I. E., Rothman, L. S., Hargreaves, R. J., et al. 2022, *JQSRT*, 277, 107949

Gordon, I. E., Rothman, L. S., Hill, C., et al. 2017, *JQSRT*, 203, 3

Guerlet, S., Spiga, A., Sylvestre, M., et al. 2014, *Icar*, 238, 110

Hirtzig, M., Bézard, B., Lellouch, E., et al. 2013, *Icar*, 226, 470

Hörst, S. M. 2017, *JGRE*, 122, 432

Hourdin, F., Grandpeix, J.-Y., Rio, C., et al. 2013, *CIDy*, 40, 2193

Hourdin, F., Lebonnois, S., Luz, D., & Rannou, P. 2004, *JGRE*, 109, E12005

Hubbard, W. B., Sicardy, B., Miles, R., et al. 1993, *A&A*, 269, 541

Jennings, D. E., Achterberg, R. K., Cottini, V., et al. 2015, *ApJL*, 804, L34

Jennings, D. E., Anderson, C. M., Samuelson, R. E., et al. 2012, *ApJL*, 761, L15

Karman, T., Gordon, I. E., van der Avoird, A., et al. 2019, *Icar*, 328, 160

Kostiuk, T., Livengood, T. A., Hewagama, T., et al. 2005, *GeoRL*, 32, L22205

Lacis, A. A., & Oinas, V. 1991, *JGR*, 96, 9027

Larson, E. J. L., Toon, O. B., & Friedson, A. J. 2014, *Icar*, 243, 400

Lavvas, P., Coustenis, A., & Varvadas, I. M. 2008, *P&SS*, 56, 27

Lavvas, P., Yelle, R. V., & Griffith, C. A. 2010, *Icar*, 210, 832

Lebonnois, S., Burgalat, J., Rannou, P., & Charnay, B. 2012, *Icar*, 218, 707

Lebonnois, S., Hourdin, F., Rannou, P., Luz, D., & Toublanc, D. 2003, *Icar*, 163, 164

Lebonnois, S., Rannou, P., & Hourdin, F. 2009, *RSPTA*, 367, 665

Lebonnois, S., Toublanc, D., Hourdin, F., & Rannou, P. 2001, *Icar*, 152, 384

Lellouch, E., Coustenis, A., Gautier, D., et al. 1989, *Icar*, 79, 328

Le Mouélic, S., Rodriguez, S., Robidel, R., et al. 2018, *Icar*, 311, 371

Lewis, N. T., Lombardo, N. A., Read, P. L., & Lora, J. M. 2023, *PSJ*, 4, 149

Li, J., & Barker, H. W. 2005, *JAtS*, 62, 286

Lombardo, N. A., & Lora, J. M. 2023, *Icar*, 390, 115291

Lora, J. M., Battalio, J. M., Yap, M., & Baciocco, C. 2022, *Icar*, 384, 115095

Lora, J. M., Lunine, J. I., & Russell, J. L. 2015, *Icar*, 250, 516

Lora, J. M., Tokano, T., Vatant d'Ollone, J., Lebonnois, S., & Lorenz, R. D. 2019, *Icar*, 333, 113

Lorenz, R. D., Stiles, B. W., Aharonson, O., et al. 2013, *Icar*, 225, 367

MacKenzie, S. M., Lora, J. M., & Lorenz, R. D. 2019, *JGRE*, 124, 1728

Mathé, C., Vinatier, S., Bézard, B., et al. 2020, *Icar*, 344, 113547

McKay, C. P., Pollack, J. B., & Courtin, R. 1989, *Icar*, 80, 23

Mollière, P., van Boekel, R., Dullemond, C., Henning, T., & Mordasini, C. 2015, *ApJ*, 813, 47

Newman, C. E., Lee, C., Lian, Y., Richardson, M. I., & Toigo, A. D. 2011, *Icar*, 213, 636

Newman, C. E., Richardson, M. I., Lian, Y., & Lee, C. 2016, *Icar*, 267, 106

Niemann, H. B., Atreya, S. K., Bauer, S. J., et al. 2005, *Natur*, 438, 779

Niemann, H. B., Atreya, S. K., Demick, J. E., et al. 2010, *JGRE*, 115, E12006

Pierrehumbert, R. T. 2010, Principles of Planetary Climate (Cambridge: Cambridge Univ. Press),

Rannou, P., Hourdin, F., McKay, C. P., & Luz, D. 2004, *Icar*, 170, 443

Rey, M., Nikitin, A. V., Bézard, B., et al. 2018, *Icar*, 303, 114

Richard, C., Gordon, I. E., Rothman, L. S., et al. 2012, *JQSRT*, 113, 1276

Rothman, L. S., Gordon, I. E., Babikov, Y., et al. 2013, *JQSRT*, 130, 4

Schinder, P. J., Flasar, F. M., Marouf, E. A., et al. 2011, *Icar*, 215, 460

Schinder, P. J., Flasar, F. M., Marouf, E. A., et al. 2012, *Icar*, 221, 1020

Schinder, P. J., Flasar, F. M., Marouf, E. A., et al. 2020, *Icar*, 345, 113720

Sharkey, J., Teanby, N. A., Sylvestre, M., et al. 2020, *Icar*, 337, 113441

Sharkey, J., Teanby, N. A., Sylvestre, M., et al. 2021, *Icar*, 354, 114030

Sicardy, B., Colas, F., Widemann, T., et al. 2006, *JGRE*, 111, E11S91

Sylvestre, M., Teanby, N. A., Vatant d'Ollone, J., et al. 2020, *Icar*, 344, 113188

Sylvestre, M., Teanby, N. A., Vinatier, S., Lebonnois, S., & Irwin, P. G. J. 2018, *A&A*, 609, A64

Teanby, N. A., Bézard, B., Vinatier, S., et al. 2017, *NatCo*, 8, 1586

Teanby, N. A., de Kok, R., Irwin, P. G. J., et al. 2008a, *JGRE*, 113, E12003

Teanby, N. A., Irwin, P. G. J., de Kok, R., et al. 2008b, *Icar*, 193, 595

Teanby, N. A., Irwin, P. G. J., de Kok, R., & Nixon, C. A. 2009, *RSPTA*, 367, 697

Teanby, N. A., Irwin, P. G. J., de Kok, R., & Nixon, C. A. 2010, *ApJL*, 724, L84

Teanby, N. A., Irwin, P. G. J., Nixon, C. A., et al. 2012, *Natur*, 491, 732

Teanby, N. A., Sylvestre, M., Sharkey, J., et al. 2019, *GeoRL*, 46, 3079

Tokano, T. 2005, *Icar*, 173, 222

Tokano, T. 2019, *Icar*, 317, 337

Tokano, T. 2023, *Icar*, 389, 115253

Tokano, T., Neubauer, F. M., Laube, M., & McKay, C. P. 1999, *P&SS*, 47, 493

Tomasko, M. G., Doose, L., Engel, S., et al. 2008, *P&SS*, 56, 669

Toon, O. B., McKay, C. P., Ackerman, T. P., & Santhanam, K. 1989, *JGR*, 94, 16287

Toon, O. B., Turco, R. P., Westphal, D., Malone, R., & Liu, M. S. 1988, *JAtS*, 45, 2123

Turbet, M., Leconte, J., Selsis, F., et al. 2016, *A&A*, 596, A112

Vatant d'Ollone, J. 2020, Thesis, Sorbonne Université, https://theses.hal.sciences/tel-03561601

Vervack, R. J. J., Sandel, B. R., & Strobel, D. F. 2004, *Icar*, 170, 91

Vinatier, S., Bézard, B., de Kok, R., et al. 2010, *Icar*, 210, 852

Vinatier, S., Bézard, B., Fouchet, T., et al. 2007, *Icar*, 188, 120

Vinatier, S., Bézard, B., Lebonnois, S., et al. 2015, *Icar*, 250, 95

Vinatier, S., Mathé, C., Bézard, B., et al. 2020, *A&A*, 641, A116

Vinatier, S., Schmitt, B., Bézard, B., et al. 2018, *Icar*, 310, 89

Vuitton, V., Yelle, R. V., & Cui, J. 2008, *JGRE*, 113, E05007

Vuitton, V., Yelle, R. V., Klippenstein, S. J., Hörst, S. M., & Lavvas, P. 2019, *Icar*, 324, 120

Wordsworth, R. D., Forget, F., Selsis, F., et al. 2011, *ApJL*, 733, L48