# The New Titan Planetary Climate Model. II. Titan's Haze and Cloud Cycles

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

► To cite this version:

Bruno de Batz de Trenquelléon, Pascal Rannou, Jérémie Burgalat, Sébastien Lebonnois, Jan Vatant D'ollone. The New Titan Planetary Climate Model. II. Titan's Haze and Cloud Cycles. The Planetary Science Journal, 2025, 6 (4), pp.79. 10.3847/PSJ/adbb6c. hal-05041890

HAL Id: hal-05041890

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

Submitted on 22 Apr 2025

HAL is a multi-disciplinary open access
archive for the deposit and dissemination of sci-
entific research documents, whether they are pub-
lished 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>

# The New Titan Planetary Climate Model. II. Titan's Haze and Cloud Cycles

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

<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

<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

Received 2024 August 27; revised 2025 February 19; accepted 2025 February 26; published 2025 March 31

## Abstract

Observations of Titan have long revealed the presence of seasonal cycles, including haze, clouds, and gases. The lack of information on the different processes that govern these cycles prevents us from understanding all the phenomena taking place in Titan's atmosphere. In this work, we survey Titan's cycles exhaustively to understand their mechanisms, in particular the haze and cloud cycles that govern a large part of Titan's climate and its surface. The Titan Planetary Climate Model (Titan PCM)—an upgraded version of the IPSL Titan Global Climate Model— incorporates a microphysical model for haze and clouds, along with a methane cycle (evaporation, transportation, condensation, precipitation). The model also features improved radiative transfer and updated spectroscopic data on gases. The Titan PCM includes phenomena hitherto omitted from Titan 3D simulations, such as cloud nucleation and condensation. Furthermore, we have coupled microphysics with radiative transfer in order to take account of feedback from haze and clouds on the thermal structure. We find that the haze follows a seasonal cycle constrained by the pole-to-pole circulation in the stratosphere. Below, the cold regions between 30 and 90 km trigger the condensation of minor species, forming a condensate mist layer throughout the year. Methane clouds form near the surface (5–25 km), in the ascending branches of the circulation, while polar clouds form higher up (60–300 km) at the winter pole after enrichment in trace compounds. The Titan PCM now enables detailed exploration of the long-term variations in Titan's cycles, providing insights for future observations.

*Unified Astronomy Thesaurus concepts: Titan (2186); Planetary atmospheres (1244); Planetary climates (2184); Methane (1042); Atmospheric clouds (2180)*

## 1. Introduction

Observations of Titan by the Cassini-Huygens mission (13.6 yr mission) have mapped the winds (M. K. Bird et al. 2005; F. M. Flasar et al. 2005; R. K. Achterberg et al. 2008, 2011; E. Karkoschka 2016; J. Sharkey et al. 2021), thermal structure, and atmospheric composition (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) of this Saturnian moon over almost half a Titanian year, covering the northern winter to the northern summer solstice ($L_S = 293°–93°$). Analyses of these observations have revealed several seasonal variations in both the thermal structure and the atmospheric circulation. They confirmed the existence of a planetary-scale summer to winter circulation in the stratosphere due to intense temperature contrasts predicted by the model of F. Hourdin et al. (1995), with a rapid reversal of this circulation after the northern spring equinox in 2009 (R. K. Achterberg et al. 2008, 2011; P. J. Schinder et al. 2012; N. A. Teanby et al. 2012; S. Vinatier et al. 2015; R. A. West et al. 2018; B. Seignovert et al. 2021; J. Sharkey et al. 2021). In addition, Huygens measurements provided a unique wind vertical profile at 10°S during the northern winter period (M. K. Bird et al. 2005). A strong zonal wind shear was observed around 70 km, within a region 10 km high, with winds of a few meters per second

surrounded by zonal winds of tens of m s⁻¹ below and above. This signature is characteristic of two distinct circulation systems. Although it may be linked to the strong positive temperature gradient, its origin remains poorly understood (S. Lebonnois et al. 2012; J. Li et al. 2012).

These observations also highlight the existence of seasonal cycles in the atmospheric composition and in the haze layer. A pronounced latitudinal variation of trace molecules such as C₂H₂, C₂H₆, and HCN is observed, with an enrichment in the stratosphere at the winter pole (S. Vinatier et al. 2007, 2015, 2020; N. A. Teanby et al. 2009, 20010, 2012; A. Coustenis et al. 2016, 2018; C. Mathé et al. 2020). The haze layer is also modulated by the seasonal circulation and shows an interhemispheric asymmetry and a winter polar enrichment of the same nature as for minor species (N. A. Teanby et al. 2008a, 2017, 2019; R. A. West et al. 2018; J. Sharkey et al. 2020; B. Seignovert et al. 2021).

It has been shown that the polar enrichment in chemical species and haze is due to a thermally direct cell in the stratosphere, with an ascending circulation in the summer hemisphere and a descending circulation toward the winter pole, able to transport the air loaded in minor molecules in the mesosphere directly toward the winter pole (S. Lebonnois et al. 2001, 2003; P. Rannou et al. 2004; A. Crespin et al. 2008; S. Vinatier et al. 2015; N. A. Teanby et al. 2017; A. Coustenis et al. 2018; C. Mathé et al. 2020). The role of this enrichment in the radiative balance of the lower winter stratosphere, as well as its role in the strengthening of the polar vortex observed by the Composite Infrared Spectrometer (CIRS; N. A. Teanby et al. 2008a) and the Imaging Science Subsystem (ISS; R. A. West et al. 2016) onboard Cassini, was studied by S. Lebonnois et al.

<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. Any further distribution of this work must maintain attribution to the author(s) and the title of the work, journal citation and DOI.
(2003), P. Rannou et al. (2004), and N. A. Lombardo & J. M. Lora (2023).

A wide variety of clouds have been observed (C. A. Griffith et al. 2005, 2006, 2009; C. C. Porco et al. 2005; E. L. Schaller et al. 2006, 2009; S. Rodriguez et al. 2009, 2011; C. M. Anderson et al. 2010; M. E. Brown et al. 2010; C. M. Anderson & R. E. Samuelson 2011; E. P. Turtle et al. 2011, 2018; D. E. Jennings et al. 2012; S. Le Mouëlic et al. 2012, 2018; R. J. de Kok et al. 2014; R. A. West et al. 2016), with methane clouds associated with condensation in the ascending circulation in the troposphere and clouds of minor species related to descending branches in the stratosphere or more complex seasonal changes. Despite CIRS and ISS observations of clouds and studies of their composition (S. Vinatier et al. 2007, 2018; C. M. Anderson et al. 2010, 2018; R. J. de Kok et al. 2014; C. Mathé et al. 2020; L. E. Hanson et al. 2023), doubts remain about the exact nature (phase, composition, mixing) of the condensate clouds that form at the poles (D. E. Jennings et al. 2015; C. M. Anderson et al. 2018).

Observations of liquid methane and ethane in lakes and seas (E. R. Stofan et al. 2007; R. H. Brown et al. 2008; A. Hayes et al. 2008; S. P. D. Birch et al. 2017; M. Mastrogiuseppe et al. 2018, 2019) suggest the presence of a vigorous hydrological-type cycle of methane and ethane on Titan, similar to the water cycle on Earth, which regulates part of the atmospheric meteorology (J. I. Lunine & R. D. Lorenz 2009; H. G. Roe 2012; J. L. Mitchell & J. M. Lora 2016; S. M. Hörst 2017; A. G. Hayes et al. 2018; E. P. Turtle et al. 2018; J. M. Lora 2024). The fate of minor species once condensed is still unclear, as are the effects of rain and aerosol deposition on the surface (P. J. Schinder et al. 2012). Some elements of answers can be provided by studying their surface and tropospheric fluxes.

These observations display the complexity of atmospheric phenomena, revealing unexpected features such as the significant decrease of trace compounds near the stratospheric north pole in 2016 (A. Coustenis et al. 2016, 2018), well after the inversion of the meridional circulation in 2012 predicted by models (S. Lebonnois et al. 2012; J. M. Lora et al. 2015, 2019; N. A. Lombardo & J. M. Lora 2023). In addition, a similar delay, probably of the same origin, is also observed in the behavior of the seasonal cycle of the aerosol detached haze layer (DHL) (R. A. West et al. 2011, 2018; B. Seignovert et al. 2021).

Our understanding of Titan's atmosphere has improved significantly since the rise of dedicated circulation models (e.g., A. D. del Genio et al. 1993; F. Hourdin et al. 1995; T. Tokano et al. 1999) and the coupling of several radiatively active components (e.g., S. Lebonnois et al. 2001, 2003; P. Rannou et al. 2002, 2004, 2006). Today, four models continue to supply the scientific community with studies of Titan's atmosphere. The TitanWRF model (C. E. Newman et al. 2011, 2016), the Titan Atmospheric Model (TAM; J. M. Lora et al. 2015, 2022; N. A. Lombardo & J. M. Lora 2023), the Cologne Titan Global Climate Model (GCM; T. Tokano 2019, 2023), and the IPSL Titan GCM (F. Hourdin et al. 1995; S. Lebonnois et al. 2012). For instance, the Cologne Titan GCM, the TitanWRF, and the TAM incorporate a methane cycle, including detailed surface reservoir distributions (J. M. Lora et al. 2015; C. E. Newman et al. 2016), latent heat effects (T. Tokano 2009), and an idealized precipitation scheme.

The importance of the coupling between dynamics, microphysics, and chemistry in GCMs has been discussed in S. Lebonnois et al. (2009). The first couplings between different model components emerged in the early 2000s, with examples such as dynamics–microphysics coupling in P. Rannou et al. (2002, 2004) and photochemistry–transport coupling in S. Lebonnois et al. (2003). These coupling approaches contribute to a more realistic modeling of the interactions between the various atmospheric processes, but only a few models take them into account. To date, the modeling integration of microphysics for haze began initially for Titan in 1D (O. B. Toon et al. 1988; C. P. McKay et al. 1989; M. Cabane et al. 1992, 1993; J. Burgalat et al. 2014; J. Burgalat & P. Rannou 2017), then 2D models (P. Rannou et al. 2002, 2004, 2006), and is now extending to 3D models as the IPSL Titan GCM (S. Lebonnois et al. 2012). The IPSL Titan GCM (S. Lebonnois et al. 2012) is the only one of the four models mentioned above to include a microphysical model capable of studying the haze cycles. Additionally, none of the 3D GCMs have a microphysical model for clouds, unlike the 2D model of P. Rannou et al. (2006), which produces methane clouds in the ascending branches of circulation and ethane clouds in descending branches.

Despite the increasing development of models, some limitations remain. The Cologne Titan GCM is restricted mainly to the troposphere and does not reproduce significant superrotation, which hampers its ability to deal with seasonal variations in the stratospheric circulation, and probably alters predictions in lower layers. Both the TitanWRF model and the IPSL Titan GCM have difficulty reproducing temperature profiles consistent with observations in the stratosphere beyond 300 km, which has a significant impact on stratospheric circulation and mixing. The TAM reproduces the thermal structure of the stratosphere correctly, thanks to an advanced radiative transfer scheme (N. A. Lombardo & J. M. Lora 2023), which enables the stratopause to be modeled. Although the results are satisfying, it does not use coupling, but a thermal balance forced with the spatial and seasonal variations of haze opacity and minor species abundance interpolated from observations.

This study aims to improve our understanding of the atmospheric cycles in Titan's atmosphere using 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 microphysical model for haze and clouds was implemented in the model, along with a methane cycle. Furthermore, the chemistry and microphysics of the model are now coupled with radiative transfer and dynamics. The main objective is to acquire an exhaustive overview of all the fluxes that govern Titan's cycles of haze, gas, and clouds. This will provide a better understanding of their associated processes, such as their sources or sinks.

A general overview of the model with a detailed description of the haze and cloud microphysical model is given in Section 2. An analysis of the haze cycle carried out with the new model is done in Section 3. A detailed description of the cloud cycle, with emphasis on the methane cycle, is given in Section 4. The impact of coupling haze and clouds with radiative transfer on tropospheric temperature and surface temperature is discussed in Section 5.

In the companion paper (B. de Batz de Trenquelléon et al. 2025), hereafter referred to as “Paper I,” a detailed description
of the new radiative transfer scheme based on the correlated-k method and recent spectroscopic data is provided. This paper also includes a study of the mechanisms and feedback in the upper atmosphere, focusing on their control of the thermal structure and their consequences for atmospheric circulation, as investigated in N. A. Lombardo & J. M. Lora (2023). Particular attention is paid to the destabilization of the low winter polar stratosphere (M. Sylvestre et al. 2020), with a study of the mechanisms taking place in polar regions.

## 2. The Titan PCM

### 2.1. Overview of the Model

In recent years, significant improvements have been made to the IPSL Titan GCM (B. Charnay & S. Lebonnois 2012; S. Lebonnois et al. 2012; J. Vatant d'Ollone 2020), now referred to as the Titan PCM (see also Paper I). These include an improved photochemistry module based on recent spectroscopic data and an updated radiative transfer, both presented in Sections 2.2 and 2.3. In addition, the Titan PCM integrates an advanced microphysical model adapted for Titan's haze and clouds (J. Burgalat et al. 2014; J. Burgalat & P. Rannou 2017), coupled with radiative transfer, described in Section 2.4.

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 enstrophy for barotropic nondivergent flows and the total angular momentum for axisymmetric flows. The model solves the primitive equation of the meteorology with a constant gravity field g. It is used with the same resolution as in S. Lebonnois et al. (2012), with a horizontal resolution of 32 longitudes by 48 latitudes, each grid covering 11.25° × 3.75°. The vertical extension is also unchanged, with 55 levels extending from the surface to ~500 km. The dynamical time step is 3 minutes, which corresponds to 8000 calls in a Titan day, with a Matsuno-Leapfrog time integration scheme. In the top four layers of the GCM, a sponge layer is included, with Rayleigh friction damping horizontal winds to their mean value. Damping time constants are 10⁴ s for the last level, then 2 × 10⁴ s, 4 × 10⁴ s, and 8 × 10⁴ s, respectively, for the three levels below.

The physical parameterizations of the Titan PCM have been drawn from the previous IPSL Titan GCM (S. Lebonnois et al. 2001, 2012; A. Crespin et al. 2008; B. Charnay & S. Lebonnois 2012). The boundary layer scheme and surface parameters have been maintained in consistency with the approach outlined in S. Lebonnois et al. (2012). The only exception is the surface thermal inertia, which has been recalibrated to a value of $I = 335 \text{ J m}^{-2} \text{ s}^{-0.5} \text{ K}^{-1}$, taken from T. Tokano (2005). Tests with different values of thermal inertia showed that the near-surface environment of the model was not significantly affected by this parameter. Nonetheless, future studies could be conducted using the thermal inertia maps by S. M. MacKenzie et al. (2019). For the surface, we now include in the model the topography from P. Corlies et al. (2017), limited at the eighth order in the spherical harmonic expansion. The physical processes are called 200 times per Titan day, i.e., once every 2 hr.

This version of the model is subject to a loss of angular momentum in the troposphere, which causes momentum to leak from the stratosphere into the troposphere and a fall in

stratospheric superrotation. As detailed in Paper I, to overcome this problem, we are nudging the zonal wind over the first 23 layers of the atmosphere, from 0 to 65 km. We force the zonal wind on the profile averaged over the first year of our simulation, starting from the converged state of S. Lebonnois et al. (2012), for each latitude, with a relaxation time of one Titanian year ($2.8 \times 10^8$ s).

### 2.2. Photochemistry

The Titan PCM photochemistry solver is detailed in Paper I and will only be presented here briefly. The photochemistry has been updated by integrating the model of V. Vuitton et al. (2019), with only neutral chemistry. It is based on recent spectroscopic data from the HITRAN2012 (HIGH-resolution TRANsmission molecular absorption) database (L. S. Rothman et al. 2013), resulting in a more comprehensive representation of gas cycles. The revised reaction scheme now extends from the planetary surface up to an altitude of 1300 km on each 1D column, as described in S. Lebonnois et al. (2001). Above the dynamical core upper limit, we use a vertical grid discretization of 10 km built from the atmospheric structure derived by R. J. J. Vervack et al. (2004), with an eddy diffusion vertical profile taken from V. Vuitton et al. (2019). The model encompasses 43 distinct chemical species engaged in 323 chemical reactions. An altitude-dependent gravity field value has been incorporated into the chemistry framework, as this is imperative for accurately determining photodissociation rates, which are determined by the optical depth of gases at a given level (S. Lebonnois et al. 2001). Given the time it takes to calculate chemical processes, we only perform these calculations two times per Titan day. Furthermore, photochemical tendencies are computed using zonal means and then recast in the longitudinal dimension proportionally to the difference of each cell to the zonal mean.

The photochemical solver integrated into the Titan PCM yields mixing ratios that are approximately 10 times lower than those observed by CIRS in the lower stratosphere. To address this discrepancy, we enforce constant mixing ratios, close to the observed values, in the model's sponge layer for species involved in radiative transfer and microphysics (i.e., CH₄, C₂H₂, C₂H₆, HCN). This update results in more accurate mixing ratios in the stratosphere, significantly enhancing their radiative impact within the model and thereby improving its thermal and dynamic structures and cloud formation.

### 2.3. Radiative Transfer

The new radiative transfer of the Titan PCM is described in detail in Paper I, which is summarized here. A flexible radiative transfer scheme developed for the generic version of the PCM (R. D. Wordsworth et al. 2011; B. Charnay et al. 2013; F. Forget et al. 2013; S. Guerlet et al. 2014; M. Turbet et al. 2016) has been integrated into the Titan PCM. This update was deemed necessary, particularly because the methane band at 7.7 µm/1305 cm⁻¹ was not previously considered. This band has a significant impact on the radiative balance of the stratosphere, with a pronounced cooling effect. 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), similar to the radiative transfer of other models of Titan such as A. J. Friedson et al. (2009) and J. M. Lora et al. (2015). Because of the number of lines, this method has the strong advantage of requiring only the use of a
limited number of bands with enough accuracy for GCM applications. The Titan PCM narrow bands are 23 in the visible range and 23 in the infrared range, covering the full spectral range from 0.3 to 10,000 µm.

The absorption from CH₄, C₂H₆, C₂H₂, and HCN are included in the Titan PCM. 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 input. This database is particularly well suited to the study of planetary atmospheres and contains several million atomic transitions for around 50 molecules. An updated version of the database, HITRAN2020 (I. E. Gordon et al. 2022), is available and will be used for future studies. Spectra are generated on a 2D grid with 10 pressure bins (P = 10⁻¹⁻10⁶ Pa) and eight temperature bins (T = 60 K–200 K). To compensate for the lack of HITRAN2012 data in the near-infrared for methane, the band list is supplemented by the ab initio calculations of M. Rey et al. (2018) between 8000 and 12,500 cm⁻¹.

Collision-induced absorptions are transitory phenomena (10⁻¹³ s). According to Heisenberg's indeterminacy principle, this results in a very broad spectral signature, which is thus considered as a continuum in the model. Collision-induced absorptions are also tabulated in HITRAN2012 (C. Richard et al. 2012). The model also takes into account the Rayleigh scattering in the radiative transfer computation. Finally, the radiative transfer equations are solved in the plane-parallel approximation, with a two-stream source function method (C. P. McKay et al. 1989; O. B. Toon et al. 1989).

In addition, radiative transfer is now coupled with the chemistry (see Paper I) and microphysics (see Section 2.4.3) of the model. With this improvement, we now provide a better representation of the thermal structure of the upper stratosphere in the Titan PCM, including the modeling of the stratopause. To account for these improvements, we have increased the number of radiative transfer calls to 40 per Titan day, improving our ability to account for short-term opacity variability and provide a comprehensive description of the influence of clouds on Titan's thermal profile.

## 2.4. Microphysics

We present here the new microphysical model for haze and clouds implemented in the Titan PCM. This model simulates microphysical haze phenomena, including aerosol production, coagulation, and sedimentation, as detailed in Section 2.4.1. It also simulates the nucleation and condensation of hydrocarbon and nitrile clouds, phenomena that have never before been integrated into a 3D GCM of Titan, as detailed in Section 2.4.2.

### 2.4.1. Microphysical Model of Haze

In planetary environments, aerosols significantly influence climate dynamics, with Titan providing a prominent example. On Titan, aerosols are produced at high altitudes through the photodissociation of nitrogen and methane. The resulting particles sediment and produce spherical particles. These particles, called monomers, are defined as the primary unit of an aggregate or the largest available spherical particle (R. A. West & P. H. Smith 1991). Subsequent coagulation between monomers forms fractal aggregates (R. A. West & P. H. Smith 1991; M. Cabane et al. 1992), whose specific structure necessitates adaptations to conventional microphysical laws (M. Cabane et al. 1993).

<figure><img src="image_4.png" /><figcaption>**Figure 1.** Top: size distribution of aerosols used in the Titan PCM, derived from the 2D model of P. Rannou et al. (2004) and J. Burgalat & P. Rannou (2017) utilizing a bin representation. This distribution, characterized by a combination of power laws, is parameterized using only two variables: the normalization factor $C_0$ and the mean aerosol radius $r_c$. Bottom: intermoment function associated to the size distribution used in the Titan PCM. The function is presented for the interval $k \in [-2; 6]$ containing all the orders necessary for the development of the microphysical laws.</figcaption></figure>

J. Burgalat & P. Rannou (2017) developed a microphysical model specific to Titan's aerosols, integrating both spherical and fractal aggregates. This model, derived from the classical bin distribution approach, employs distribution moments to address computational challenges. It governs interactions between two aerosol populations with differing structures: spherical aerosols (referred to as S) smaller than monomers, and fractal aggregates (referred to as F) larger than monomers.

Aerosol mechanical properties are influenced by their fractal structure, which results in a size distribution on Titan that differs significantly from typical log-normal patterns observed in most natural environments (P. Rannou et al. 2004; T. Cours et al. 2011; E. J. L. Larson et al. 2014). In the Titan PCM, we use a unique type of aerosol size distribution, derived from the 2D Titan bin model of P. Rannou et al. (2004). The model of P. Rannou et al. (2004) describes the size distribution of Titan's aerosols using a radius grid, with 40 bins and a radius step of 2¹/³, covering a wide range of radii from 1.64 × 10⁻⁹ to 20.01 × 10⁻³ m (O. B. Toon et al. 1988). Figure 1 shows this aerosol size distribution controlled by two parameters linked to the aerosol total number and volume: a characteristic radius $r_c = 4.582 \times 10^{-7}$ m and a normalization factor $C_0 = 2.122 \times 10^{-3}$.

Once the distribution is set, the effects of microphysical processes on the concentration $C^a(r, z)$ of aerosols of radius r at altitude z can be modeled using a continuity equation:

$$ \frac{\partial}{\partial t} C^a(r, z) = Q(z) + \frac{\partial}{\partial t} C^a(r, z)|_{\text{trans.}} + \frac{\partial}{\partial t} C^a(r, z)|_{\text{collec.}} + \frac{\partial}{\partial t} C^a(r, z)|_{\text{nucl.}} \quad (1) $$
<table><caption>**Table 1**

General Properties of Aerosols in The Titan PCM</caption><thead><tr><th>Parameters</th><th>Values</th><th>References</th></tr></thead><tbody><tr><td>Fractal dimension</td><td>D<sub>f</sub> = 2.0</td><td>M. Cabane et al. (1993)</td></tr><tr><td>Monomer radius</td><td>r<sub>m</sub> = 50 nm</td><td>M. G. Tomasko et al. (2008)</td></tr><tr><td>Production radius</td><td>r<sub>p</sub> = 10 nm</td><td>P. Rannou et al. (2004)</td></tr><tr><td>Aerosol density</td><td>ρ<sub>p</sub> = 800 kg m<sup>-3</sup></td><td>M. G. Trainer et al. (2006)</td></tr><tr><td>Production rate</td><td>Q<sub>p</sub> = 2.1 × 10<sup>-13</sup> kg m<sup>-2</sup> s<sup>-1</sup></td><td>M. G. Tomasko et al. (2008)</td></tr><tr><td>Production altitude</td><td>z<sub>0</sub> = 415 km (1 Pa)</td><td>P. Rannou et al. (2004)</td></tr><tr><td>Gaussian standard deviation</td><td>Δz = 20 km</td><td>P. Rannou et al. (2004)</td></tr><tr><td>Electric charge</td><td>n<sub>e</sub><sup>-</sup> = 15 e<sup>-</sup> μm<sup>-1</sup></td><td>P. Lavvas et al. (2010)</td></tr></tbody></table>

where $Q(z)$ corresponds to the particle mass production function. The other three terms correspond to the evolution of the concentration under the effect of transportation, collection, and nucleation, respectively. In the Titan PCM, we address only the phenomena of aerosol production, sedimentation for transport processes, and Brownian coagulation for collection processes. Aerosol diffusion is treated into the Titan PCM along with other tracer quantities. The nucleation term will be discussed in Section 2.4.2 dedicated to clouds.

The model of J. Burgalat & P. Rannou (2017) describes these phenomena in terms of the distribution of moments. The kth moment of a distribution law of aerosols in radius is defined by

$$M_k(n) = \int_0^\infty r^k n(r) dr. \qquad (2)$$

To describe the evolution of particle distribution through microphysical processes in terms of moments, we integrate Equation (1) by applying Equation (2), as follows:

$$\frac{\partial M_k}{\partial t} = \frac{\partial}{\partial t} \int_0^\infty r^k \frac{C^a(r)}{\partial r} dr. \qquad (3)$$

In this model, we only deal with two selected moments for both modes (S and F): the number density of aerosols $M_0$ in $m^{-3}$, and the volume concentration of aerosols $M_3$ in $m^3 m^{-3}$. For some processes, formatting the calculation can produce moments of different orders in integrals. However, for a given type of size distribution, it is possible to relate the moments of different orders with a simple relation. For this purpose, we use the intermoment function $\alpha(k)$, which depends on the size distribution $n(r)$, as presented by J. Burgalat et al. (2014):

$$\frac{M_k}{M_0} = r_c^k \times \alpha(k), \qquad (4)$$

where $r_c$ is the characteristic radius of the distribution, and $k$ is the order of the moment. As the size distribution is fixed, the intermoment relationship function can be expressed once for all as a function of the order $k$ only, as shown in Figure 1.

The first term of the right side of Equation (1), $Q(z)$, represents the production of aerosols (number of particles per unit of volume and time) due to the production of photochemical aerosols at high altitude. Aerosols are only produced in the small spherical aerosol mode (S). The injection is performed at the top of the model with a simple production law using a Gaussian function (M. Cabane et al. 1992). We express this production in terms of the zero-order moments for

spherical aerosols:

$$\frac{\partial M_0^S}{\partial t} = \frac{Q_p}{\frac{4}{3}\pi \rho_p r_p^3} \times \frac{1}{\sqrt{2\pi} \Delta z} \exp\left[-\frac{(z-z_0)^2}{2\Delta z^2}\right], \qquad (5)$$

where $Q_p$ is the aerosol production rate, $\rho_p$ is the aerosol density, $r_p$ is the initial model radius that represents the smallest size of aerosols, $z_0$ is the production altitude, and $\Delta z$ is the standard deviation of the production zone (Table 1). The third-order moment evolution is obtained by directly applying Relation (4), using a characteristic radius $r_c = r_p$ to represent the radius of spherical aerosols produced by photochemistry:

$$\frac{\partial M_3^S}{\partial t} = \frac{\partial M_0^S}{\partial t} \times (r_p^3 \times \alpha(3)). \qquad (6)$$

The second term of the right side of Equation (1) describes the transport processes. 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. The equation describing the temporal evolution of the kth order moment is expressed in moments as

$$\frac{\partial M_k}{\partial t} = -\vec{\nabla} \cdot \vec{\Phi}_{\text{sed}} = -\frac{d}{dz} \left[ \frac{2\rho_p g}{9\eta E} \left( M_{\frac{D_f(k+3)}{D_f}} + \frac{A\lambda_a}{E} M_{\frac{D_f(k+3)-6}{D_f}} \right) \right], \qquad (7)$$

$\frac{D_f-3}{D_f}$

where $E = r_m^{\frac{D_f}{D_f-3}}$, $r_m$ is the monomer radius, and $D_f$ is the aerosol fractal dimension, which defines the evolution of the volume of matter $M$ contained within a sphere of increasing radius $R$ ($M \propto R^{D_f}$). $\rho_p$ is the aerosol density, $\eta$ is the air viscosity, and $\lambda_a$ is the atmospheric mean free path. $A$ is the coefficient used in the first-order slip-flow correction, with $A = 1.591$.

Finally, the main processes involved in aerosol microphysics are intramodal (SS and FF) and intermodal (SF) coagulations, included in the third term of the right side of Equation (1). The equation describing the time evolution of the moments through coagulation depends on the flow regime in which we are working. On Titan, the flow regime spans from the molecular regime in the upper stratosphere and the mesosphere to the continuous regime in the troposphere, so we also have to deal with the transitional regime. Details of the calculations are given in J. Burgalat & P. Rannou (2017); we will remind here only the final form of the equations for the couple ($M_0$, $M_3$) in
both modes:

$$ \frac{\partial M_0^S}{\partial t} = \beta_1 \times M_0^{S2} - \beta_2 \times M_0^S M_0^F, \quad (8) $$

$$ \frac{\partial M_0^F}{\partial t} = \gamma_1 \times M_0^{S2} - \gamma_2 \times M_0^F, \quad (9) $$

$$ \frac{\partial M_3^S}{\partial t} = \delta_1 \times M_3^{S2} - \delta_2 \times M_3^S M_3^F, \quad (10) $$

$$ \frac{\partial M_3^F}{\partial t} = -\delta_1 \times M_3^{S2} + \delta_2 \times M_3^S M_3^F, \quad (11) $$

where the factors $\beta_i$, $\gamma_i$, and $\delta_i$ define the coagulation kernels, which represent the collision frequency of particles. They are dependent on atmospheric properties, the function $\alpha(k)$ (Equation (4)), and the morphological parameters of aerosols. We observe that the sum of the two last equations is zero (Equations (10) + (11)), indicating the conservation of the system's total mass during the coagulation process.

When the aerosols are electrically charged with a number of electrons per length unit $n_e^-$, the coagulation kernels must be adjusted by a multiplicative factor to take into account the electrical interactions of the aerosols involved in the coagulation process. The aerosol charge is repulsive and can completely inhibit the coagulation. However, the parameter $n_e^-$ is not well constrained, with typical values ranging between 15 and $30 \ e^- \ \mu m^{-1}$ for Titan's aerosols (W. J. Borucki et al. 1987; C. P. McKay et al. 1989; P. Rannou et al. 2003). The general parameters for aerosols used in the Titan PCM are summarized in Table 1 with the references for each value.

### 2.4.2. Microphysical Model of Clouds

Knowledge of cloud modeling is crucial because of their impact on climate. In particular, due to their single-scattering albedo, which is higher than that of haze in the visible spectrum, clouds and condensed mist have a major influence on the distribution of solar energy in the lower atmosphere. Clouds are also associated with precipitation formation processes and interact with other atmospheric components, modulating their quantity and distribution. Finally, clouds are essential to Titan's atmospheric cycles, particularly the methane cycle.

Cloud microphysics involves two concurrent processes: nucleation and condensation. In addition, clouds are also subject to transport as aerosols. In the same way as aerosols, the general equation governing the evolution of condensates in the atmosphere can be written as follows:

$$ \frac{\partial}{\partial t} C^c(r, z) = \frac{\partial}{\partial t} C^c(r, z)|_{\text{nucl.}} + \frac{\partial}{\partial t} C^c(r, z)|_{\text{cond./evap.}} + \frac{\partial}{\partial t} C^c(r, z)|_{\text{trans.}} + \frac{\partial}{\partial t} C^c(r, z)|_{\text{collec.}} \quad (12) $$

Here, $C^c(r, z)$ is the concentration of condensates of radius $r$ at altitude $z$. The terms on the right side of Equation (12) correspond, respectively, to the evolution of condensates under the effect of nucleation, condensation/evaporation, transportation, and collection.

In the Titan PCM, we have adapted the microphysical cloud model of J. Burgalat et al. (2014), initially developed by F. Montmessin et al. (2002), to use with a bin description for

Mars and used for Titan by P. Rannou et al. (2006). The model of J. Burgalat et al. (2014) describes microphysical phenomena in terms of moments for nucleation, condensation/evaporation, and for sedimentation included in transport phenomena. We will therefore not be dealing with deep convection, nor with collection phenomena.

In environments where gaseous species become oversaturated, they can initiate the formation of small condensed material embryos on aerosols. If these embryos exceed a critical radius—which depends on the environment and the characteristics of the aerosol—they trigger spontaneous growth, nucleation. Subsequent to nucleation, the growth process due to phase change is referred to as condensation. To describe these processes, we require the size distribution of aerosols and cloud condensation nuclei (CCN), with two pairs of moments used for aerosols and CCN distributions: ($M_0^{\text{aer}}$, $M_3^{\text{aer}}$) and ($M_0^{\text{CCN}}$, $M_3^{\text{CCN}}$). For condensate materials, as the number of drops corresponds to the number of nuclei $M_0^{\text{CCN}}$, only volumes for each involved species are necessary, denoted as moments $M_3^{\text{cond}}$, where $\chi$ corresponds to the species. In the present model, we deal with methane ($CH_4$), ethane ($C_2H_6$), acetylene ($C_2H_2$), and hydrogen cyanide (HCN).

Cloud droplets (CCN + condensates) are subject to a similar sedimentation process to aerosols and are treated as shown in Equation (7), considering the droplets as spherical particles ($D_f = 3.0$). The model of J. Burgalat et al. (2014) gives the nucleation equation expressed in moments as

$$ \frac{\partial}{\partial t} M_k^{\text{CCN}} = 4\pi J_{\text{het}} r_m^{-1} \times M_{k+3}^{\text{aer}}, \quad (13) $$

where $J_{\text{het}}$ is the heterogeneous nucleation rate given by R. G. Keese (1989) and H. R. Pruppacher et al. (1998). In the Titan PCM, nucleation is effective only in the low stratosphere and the troposphere, where only fractal aerosols remain. In Equation (13), we therefore consider that all the moments associated with aerosols are those of the fractal mode. Moreover, since nucleation involves an exchange between aerosols and CCN, the k-order moment evolution of aerosols is the opposite of that for CCN:

$$ \left. \frac{\partial M_k^{\text{aer}}}{\partial t} \right|_{\text{nucl.}} = - \left. \frac{\partial M_k^{\text{CCN}}}{\partial t} \right|_{\text{nucl.}} \quad (14) $$

Condensation, evaporation, and sublimation control the evolution of droplet size. They affect the volume of condensate making up the droplet, while the number of nuclei remains constant. In the current model, each of the condensable species is considered independent, and no mixing between condensed species is assumed. However, in order to take into account nitrogen dissolution in methane in the model, we assume a constant dissolution in the troposphere, equal to 80%, and apply a factor of 0.8 in the calculation of the saturated vapor pressure of methane. The evolution of the volume of condensate $M_3^{\chi, \text{cond}}$ associated with the species $\chi$ corresponds to the product of the evolution of the volume of an effective droplet and the total number of available nuclei $M_0^{\text{CCN}}$:

$$ \frac{dM_3^{\chi, \text{cond}}}{dt} = 4\pi r_g \times \frac{S^{\chi} - S^{\chi*}}{R_{\text{diff}} + R_{\text{des}}} M_0^{\text{CCN}}, \quad (15) $$

where $r_g$ is the droplet radius, $S^{\chi}$ corresponds to the saturation ratio of the species $\chi$, i.e., the ratio between its partial pressure
and its saturation vapor pressure, and $S^{\lambda*}$ is the critical saturation, i.e., its saturation at equilibrium (H. Pruppacher & J. Klett 2004). $R_{\text{diff}}$ and $R_{\text{des}}$ are resistance terms that slow down condensation, resulting from molecular diffusion and the release of latent heat during phase change, respectively. We use constant values in the model: $R_{\text{diff}} = 1.519 \times 10^{-21} \text{ s m}^{-2}$ and $R_{\text{des}} = 1.519 \times 10^{-20} \text{ s m}^{-2}$ (P. Rannou et al. 2019).

### 2.4.3. Coupling Microphysics and Radiative Transfer

Radiative-convective and climate models (C. P. McKay et al. 1989; P. Rannou et al. 2004) have demonstrated that Titan's thermal structure is primarily influenced by the absorption of sunlight by the stratospheric haze, the opacity of gases, and radiative cooling in polar regions. For instance, during the polar night, the accumulation of haze, in the absence of incident solar radiation, acts as a powerful infrared emitter, effectively cooling the polar night atmosphere. In addition, clouds can influence tropospheric and surface temperatures. This feedback enhances the latitudinal contrasts in temperature, especially between tropics and the winter polar region. This strengthens the atmospheric circulation compared to simulations with a uniform haze layer (P. Rannou et al. 2004). Therefore, it is essential to include the effects of both haze and clouds in radiative transfer models for a more accurate depiction of the thermal structure of the atmosphere.

The radiative coupling of haze and clouds has been implemented using a lookup table, allowing for the derivation of optical properties directly from the $M_0$ and $M_3$ of haze and clouds. The tables were created for three types of particles: spherical aerosols, fractal aerosols, and cloud droplets. For each type of particle, they provide the extinction cross section, single-scattering albedo, and asymmetry parameter as functions of their radius. These properties are calculated within the spectral range of the radiative transfer of the PCM (comprising 46 channels; see Paper I) for 33 modal radii, ranging with geometric steps from $10^{-9}$ to $10^{-5}$ m for aerosols in both modes, and from $10^{-7}$ to $10^{-3}$ m for cloud droplets.

The characteristic radius of each type of particle is calculated in each cell of the model using the following equation:

$$r^i = \left( \frac{1}{\alpha(3)} \times \frac{M_3^i}{M_0^i} \right)^{1/3}, \quad (16)$$

where $i$ corresponds to the type of particle (spherical aerosol, fractal aerosol, or cloud droplet). In this way, we can directly link the modeled $M_0$ and $M_3$ to the total optical thickness $\Delta\tau$ of aerosols and droplets in every cell, using

$$\Delta\tau(\lambda) = \sum_{i} \sigma^{i}(r^{i}, \lambda) \times M_{0}^{i}, \quad (17)$$

where $\sigma^i$ is the extinction cross section of the $i$-type particle averaged in the spectral channel $\lambda$ and in the radius bin centered on $r^i$ given by Equation (16). We can then proceed in the same way to obtain the total single-scattering albedo $\bar{\omega}$ and the total asymmetry parameter $\bar{g}$ for aerosols and droplets:

$$\begin{aligned} \bar{\omega}(\lambda) &= \frac{\sum_i \bar{\omega}^i(\lambda) \cdot \Delta\tau^i(\lambda)}{\sum_i \Delta\tau^i(\lambda)} \\ \bar{g}(\lambda) &= \frac{\sum_i \bar{g}^i(\lambda) \cdot \bar{\omega}^i(\lambda) \cdot \Delta\tau^i(\lambda)}{\sum_i \bar{\omega}^i(\lambda) \cdot \Delta\tau^i(\lambda)}, \end{aligned} \quad (18)$$

where $\bar{\omega}^i$ and $\bar{g}^i$ are, respectively, the single-scattering albedo and asymmetry parameter of the $i$-type particle averaged in the narrow band centered on $\lambda$ and in the radius band centered on $r^i$.

We finally obtain properties that depend on the type and size of the particles, which, in turn, lead to variations in the optical properties of the atmosphere (such as albedo, refractive index, etc.) according to altitude and latitude, following the distribution of gases and haze. We now have a comprehensive extinction spectrum for Titan, encompassing haze and clouds, including the troposphere, where aerosol haze merges with condensate mist and clouds.

With the coupling of microphysics to radiation and chemistry to radiation, the model considers now the interaction between aerosols, gases, clouds, and radiative transfer for a thorough investigation of feedback mechanisms on the thermal equilibrium field. The validation and importance of these couplings are discussed in Paper I for gases and in the next sections for haze and clouds.

### 2.5. Surface Process

The methane cycle operates on a seasonal scale. Within a single season, methane evaporates from lakes and seas at the surface, forms clouds, precipitates, accumulates at the surface and in underground reservoirs, and then evaporates back into the atmosphere (F. M. Flasar 1998; T. Tokano et al. 2001; T. Tokano 2005; S. K. Atreya et al. 2006; P. Rannou et al. 2006; E. R. Stofan et al. 2007; A. Hayes et al. 2008; J. L. Mitchell 2008; J. M. Lora & J. L. Mitchell 2015; J. L. Mitchell & J. M. Lora 2016; S. P. D. Birch et al. 2017; E. P. Turtle et al. 2018; S. P. Faulk et al. 2020; J. M. Lora 2024). We have incorporated modules into the Titan PCM to simulate surface-atmosphere exchanges (evaporation and surface precipitation).

We have incorporated a standard bucket model (J. L. Mitchell 2008; J. M. Lora & J. L. Mitchell 2015) into the Titan PCM. This allows for the consideration of surface-atmosphere exchanges (evaporation and precipitation) at every point on the surface. Consequently, methane evaporation in the model originates from this surface reservoir and is directly injected into the first atmospheric layer through surface-atmosphere equilibrium flux. To replicate the distribution of lakes and seas localized at Titan's poles, we initially filled this reservoir only at high latitudes ($\ge 70^{\circ}$) (J. L. Mitchell 2008; J. M. Lora & J. L. Mitchell 2015; S. P. Faulk et al. 2020). The evaporation of methane is influenced by several factors: the methane saturation level in the first layer of the atmosphere $Q_{\text{sat}}$, a term related to horizontal surface winds $\omega_s$, and a term associated with the source of turbulent kinetic energy at the surface $C_d$. The surface-atmosphere methane flux $F_{\text{CH}_4}$ in the model is governed by a modified version of the standard bulk aerodynamic formula adapted from T. Tokano et al. (2001):

$$F_{\text{CH}_4} = f_{\text{lake}} \cdot \rho_{\text{air}} \cdot C_d \cdot \omega_s \times (h \cdot Q_{\text{sat}} - Q_{\text{CH}_4}), \quad (19)$$

where $\rho_{\text{air}}$ and $Q_{\text{CH}_4}$ are, respectively, the density of the air and the mole fraction of methane in the first layer of the atmosphere. $f_{\text{lake}}$ and $h$ are two free parameters related to the evaporation of the lakes. Lakes are represented in the model using a subgrid parameterization, with $f_{\text{lake}}$ representing the fraction of the lake surface covered within each grid cell. We set $f_{\text{lake}}$ to 10% in the northern polar region and 1% in the southern polar region (A. Hayes et al. 2008; A. G. Hayes
2016). This simplification allows us to represent the total lake fraction in polar regions at the local grid scale. However, it may introduce some uncertainty in local evaporation fluxes, which will be investigated in further studies. *h* is the saturation correction account for the dissolution of nitrogen in liquid methane (fixed to 80%; W. R. Thompson et al. 1992; M. J. Malaska et al. 2017; J. Hartwig et al. 2018).

## 2.6. Initial Conditions and Simulation Setting

The initial state of the simulation presented below is derived from the converged results of the model developed by S. Lebonnois et al. (2012). The reference simulation for this work was conducted, taking into account both haze (J. Burgalat & P. Rannou 2017) and clouds (J. Burgalat et al. 2014), coupled with radiative transfer, making the temperature profile and the winds highly sensitive to the input parameters of aerosols (see Table 1) and clouds. The simulation is also coupled with the new version of photochemistry (V. Vuitton et al. 2019; J. Vatant d'Ollone 2020) and the new radiative transfer scheme (J. Vatant d'Ollone 2020).

We ran the model over a period exceeding 25 Titan years to achieve a state of convergence. Clouds were activated from the first year, starting from an initial state without clouds. Their formation began in the first year, and it took approximately 15 yr for them to reach a stable state. The following 10 yr led to complete convergence of the model. When convergence is obtained, the temperature profile remains very stable from one year to the next. It oscillates around an average profile without drifting, with a maximum relative variation between the last 2 yr being less than 0.75%. The zonal wind and extinction profiles exhibit higher sensitivity to model parameters compared to temperature. In the stratosphere, they remain relatively stable but can fluctuate around an average profile up to 10% in the troposphere, where clouds form. However, when these variations are averaged over the last 5 yr, the maximum relative variation for the zonal wind drops below a few percent, while it remains higher for total extinction in the troposphere around 6%. In the troposphere, the extinction due to clouds varies considerably and does not repeat identically year after year.

# 3. The Seasonal Cycle of Haze

## 3.1. Distribution of Haze

Figure 2 shows the meridional distribution of haze predicted by the Titan PCM in the stratosphere for the four seasons and the equinoxes. It highlights the seasonal variation of the main haze layer, located between 3 × 10³ and 5 × 10² Pa (80–130 km), as well as its detached layer, located at around 3 Pa (325 km). The stream function, representing the current lines of the meridional circulation, is plotted on the haze distribution map. It is characterized in the stratosphere by a large pole-to-pole cell, extending from the summer pole to the winter pole. Ascending winds predominate in the summer hemisphere ($L_s = 40°–150°$ at the north pole), while descending winds occur at the winter pole ($L_s = 220°–330°$ at the north pole).

Haze follows the stratospheric circulation, being carried from the summer pole to the winter pole. This leads to an enrichment of the winter pole in aerosols between 3 × 10³ and 10² Pa, to the detriment of the summer pole, as illustrated in

Figure 2(c) during the northern summer, and Figure 2(f) during the northern winter. The accumulation of haze during the polar night leads to significant cooling of the stratosphere by emission of aerosols in the infrared (see Paper I). On the other hand, at the summer pole, the solar flux warms the atmospheric column. This increases the latitudinal temperature gradient (equator to winter) that finally strengthens all the circulation (P. Rannou et al. 2004).

Below the main layer of haze ($>3 \times 10^3$ Pa), the temperature allows the minor species to reach condensation, and aerosols serve as a support for nucleation. The lower stratosphere and the troposphere are thus depleted of “dry” aerosols, but are made of wet aerosols onto which supersaturated species condense, forming a mist of condensates (see Section 4.2).

The Titan PCM predicts a reversal of the atmospheric circulation on a timescale of three Earth years, around the equinox (Paper I). This inversion induces a reversal in the haze distribution, which is transported toward the opposite pole. It enhances the presence of the haze in the fall hemisphere and depletes it in the spring hemisphere. This reversal of the haze layer is completed after an entire season following the equinox, i.e., over nearly eight Earth years (Figures 2(a)–(c) and (d)–(f)).

## 3.2. The Detached Haze Layer

The DHL is also a tracer of seasonal variations in atmospheric dynamics. Observations throughout the Cassini mission (R. A. West et al. 2011, 2018; B. Seignovert et al. 2021), as well as models (P. Rannou et al. 2002; E. J. L. Larson et al. 2015), suggest that the existence and location of the DHL is not the result of a static process. Its collapse after the equinox and its merging with the main layer, followed by its gradual reappearance later in the season, are the result of complex and dynamic seasonal processes (R. A. West et al. 2018; B. Seignovert et al. 2021) and were readily predicted by R. A. West et al. (2011).

Figure 3 illustrates the seasonal evolution of the DHL extinction at the equator between $L_s = 300°$ and $L_s = 90°$, as predicted by the Titan PCM. This map can be directly compared to Figure 15 of B. Seignovert et al. (2021), revealing a very good correlation between model predictions and Cassini observations. In the model, the detached layer dominates mainly at an altitude around 325 km (2–4 Pa), in line with the lower latitudinal temperature contrasts and weaker winds observed in this zone (Paper I). It then collapses shortly after the equinox and completely merges with the main layer below 200 km around $L_s = 50°$. This collapse explains the disparity in altitude of the DHL between the Cassini observations (500 km, $L_s = 300°–360°$; R. A. West et al. 2011) and the Voyager 2 observations (350 km, $L_s = 9°$; K. Rages & J. B. Pollack 1983). The Titan PCM subsequently predicts the reappearance of the DHL at a slightly higher altitude than before its collapse, around 375 km, before returning to its original altitude as the solstice approaches. Despite behavior similar to that of observations, the Titan PCM’s DHL is reproduced at too low an altitude compared to Cassini observations. This discrepancy is explained by the top of the model and the injection of small aerosols that are too low in altitude. Aerosols are injected at 1 Pa (~400 km) to avoid the sponge layer of the model where winds can be altered (see Section 2.4.1).

In the Titan PCM, the timing of the DHL descent coincides with the inversion of the atmospheric circulation around equinoxes. The model predicts that the ascending branch of
<figure><img src="image_5.png" /><figcaption>(a) Vernal equinox North < L$_{S}$ = 0° ></figcaption></figure>

<figure><img src="image_6.png" /><figcaption>(b) Spring North < L$_{S}$ = 60° ></figcaption></figure>

<figure><img src="image_7.png" /><figcaption>(c) Summer North < L$_{S}$ = 120° ></figcaption></figure>

<figure><img src="image_8.png" /><figcaption>(d) Fall equinox North < L$_{S}$ = 180° ></figcaption></figure>

<figure><img src="image_9.png" /><figcaption>(e) Autumn North < L$_{S}$ = 240° ></figcaption></figure>

<figure><img src="image_10.png" /><figcaption>(f) Winter North < L$_{S}$ = 300° ></figcaption></figure>

**Figure 2.** Haze distribution for (a) the vernal equinox ($L_s$ = 0°), (b) the northern spring ($L_s$ = 60°), (c) the northern summer ($L_s$ = 120°), (d) the fall equinox ($L_s$ = 180°), (e) the northern fall ($L_s$ = 240°), and (f) the northern winter ($L_s$ = 300°). The field presented is the extinction coefficient of aerosols at 0.671 µm, in zonal and diurnal average. The meridional stream function is shown as black contours (the solid line is clockwise rotation and the dashed line is counterclockwise).

wind, connecting the two Hadley cells, is initially located
above the summer polar region. It then gradually moves toward
the opposite polar region, crossing the equator at the equinoxes,
and reaches the opposite pole after Ls = 25°. As it crosses the
equator, the ascending branch lifts aerosols and maintains them
at high altitudes, forming a particularly pronounced detached
layer around Ls = 0° and Ls = 180°, reaching 3 Pa (Figure 3)
and extending from the north pole to the south pole
(Figures 2(a) and (d)). Just after the equinoxes, the passage
of the ascending branch from one pole to the other corresponds
to a weakening of the overall meridional circulation, and, in
particular, the weakening of the vertical wind in the ascending
region causes the collapse of this layer in the model. The
reappearance of the DHL occurs as the solstice approaches,
coinciding with the reestablishment of the pole-to-pole cell.

This is consistent with Cassini's observations, which noted the
descent of the detached layer after the 2010 equinox, and its
reappearance around 450 km in 2016, prior to the solstice
(R. A. West et al. 2018).

Significant improvements have been made over previous
simulations (S. Lebonnois et al. 2012; J. Vatant
d'Ollone 2020), particularly with regard to the more
pronounced meridional contrasts of the main layer, and the
higher position and intensity of the DHL. The coupling of
microphysics with radiative transfer has led to improved
temperature contrasts, reinforcing the dynamics and favoring
more clearly defined haze cycles (see Paper I). In addition, the
reduction in aerosol density to 800 kg m⁻³, compared to
1000 kg m⁻³ in the study of S. Lebonnois et al. (2012),
generates lighter aerosols, facilitating their pole-to-pole
<figure><img src="image_11.png" /><figcaption>**Figure 3.** Map of haze extinction as a function of time and altitude at the equator as predicted in the Titan PCM. The field shows the extinction coefficient of aerosols at 0.671 µm. The white stars represent the position of the DHL, identified as the local maximum of the vertical extinction profile. The white and black vertical dotted lines represent, respectively, the period of the equinox and the reappearance of the DHL. This figure can be directly compared to Figure 15 of B. Seignovert et al. (2021).</figcaption></figure>

transport and favoring their retention in the detached layer. The DHL still remains less pronounced than in observations (R. A. West et al. 2018; B. Seignovert et al. 2021), with a contrast less than 1 order of magnitude in the extinction profile. The difference between the modeled and observed positions of the DHL has repercussions on the position and amplitude of the zonal wind jet (see Paper I). A vertical extension of the model could improve this aspect of the simulations. However, its implementation is not easy due to the atmospheric thickness, which is already well beyond the acceptable limit for thin-layer approximation.

# 4. Seasonal Distribution of Hydrocarbons

## 4.1. The Methane Cycle

The seasonal variation in methane evaporation on Titan PCM's surface is depicted in Figure 4. The evaporation peaks during the latter half of spring ($L_s = 50°–100°$ in the northern hemisphere and $L_s = 230°–270°$ in the southern hemisphere), while the evaporation reaches its minimum in the fall in the northern hemisphere ($L_s = 180°–250°$) and in the winter in the southern hemisphere ($L_s = 120°–180°$). Near the poles, the surface winds predicted by Titan's PCM are notably strong during spring, which coincides with the periods of maximum methane evaporation. This suggests that the surface winds play a significant role in the model's methane evaporation process.

Once evaporated, polar methane is transported in the lower troposphere by intense circulation, below 5 km, where eddy cells dominate (B. Charnay & S. Lebonnois 2012; J. M. Lora & J. L. Mitchell 2015; J. M. Lora 2024). Unlike the results of J. M. Lora & J. L. Mitchell (2015), where eddy cells extend to low latitudes (~30°), here the surface topography of the model confines the cells to midlatitudes (~60°; see Figure 4). This has the effect of reducing methane mixing in the lower troposphere, and generates a region at ±30° where the methane-mixing ratio is almost uniform, reaching 5.5% close to the surface at the equator (Figure 4). The uniform methane-mixing ratio simulated by the Titan PCM is discrepant compared to observations of the zonal mean meridional variation of tropospheric methane humidity (M. Ádámkovics et al. 2016; J. M. Lora & M. Ádámkovics 2017), although those observations from 2014 are likely only representative of one

<figure><img src="image_12.png" /><figcaption>**Figure 4.** Left: seasonal evolution of modeled methane evaporation at the surface, in zonal and diurnal average. Right: map of modeled amount of methane, in zonal and annual average, in the atmosphere below 25 km, as a function of latitude. The meridional stream function is shown as black contours.</figcaption></figure>

<figure><img src="image_13.png" /><figcaption>**Figure 5.** Top: modeled vertical profile of the methane-mixing ratio, in zonal and diurnal average, at $L_s = 305°$ and $10°S$. The vertical profiles from Huygens' observations are shown in color for comparison (H. B. Niemann et al. 2010; T. Gautier et al. 2024). Bottom: methane cloud extinction at 0.671 µm as a function of latitude and altitude, averaged zonally, for $L_s = 27.56°$. The meridional stream function is displayed as black contours, indicating the ascending currents of circulation.</figcaption></figure>

season. Tropospheric mixing in the model could be enhanced by incorporating moist convection, which would influence it in the mid and upper troposphere, where clouds form, as well as above lakes in methane-rich regions that may potentially become oversaturated, especially in the northern hemisphere.

The top panel of Figure 5 shows the vertical profile of methane in the troposphere at the time and location of Huygens' landing ($L_s = 305°$; $10°S$). The amount of methane predicted at Titan's surface by the Titan PCM is 5.6%. This concentration slightly increases over the first few kilometers, reaching 5.75% at 3 km. Two nonexclusive possibilities could explain this rise: on the one hand, atmospheric transport is more intense around 3 km, leading to greater mixing of methane in this region compared to the surface. On the other hand, the area between 3 and 5 km is located below the cloud base, which could result in methane accumulation when the
<figure><img src="image_14.png" /><figcaption>**Figure 6.** Seasonal evolution of methane cloud opacity at 0.671 µm predicted by the Titan PCM, in zonal mean. The field corresponds to the average over the last five simulated years. Opacity peaks correspond to the presence of large methane clouds in the atmospheric column. The black contours indicate the maximum methane surface precipitation (>0.1 mm T<sub>day</sub><sup>-1</sup>).</figcaption></figure>

clouds reevaporate. Then, the methane-mixing ratio curve drops above 5 km to 1.4% around 25 km where it crosses its saturation vapor pressure curve. The trend of the methane profile predicted by the model aligns well with the observed one (H. B. Niemann et al. 2010; T. Gautier et al. 2024). However, the model predicts a lower methane amount above 11 km than observed, due to a too cold modeled tropospheric temperature, which reduces the saturation threshold.

When methane evaporates and enriches the layers near the surface, it is then transported by the ascending flow of tropospheric circulation to the upper regions of the troposphere. This phenomenon is particularly pronounced during the spring and summer, when an upward current forms around 60°. The enrichment of methane in the colder layers of the troposphere creates a region oversaturated with methane. Methane then begins to condense, forming clouds at altitudes between 5 and 25 km. Figure 5 illustrates the formation of a spring methane cloud between 6 and 20 km in the ascending branch of the circulation (40°–70°S). Tropospheric cloud systems at a similar season have been observed, though at slightly lower latitudes than those predicted by the model (C. C. Porco et al. 2005; E. L. Schaller et al. 2006; S. Rodriguez et al. 2009). This phenomenon persisted until the end of Titan's southern summer upon the arrival of Cassini in 2004.

Figure 6 shows the seasonal evolution of methane cloud opacity predicted by the Titan PCM during a Titan year. Methane clouds appear very opaque and sporadic, localized mostly at low and midlatitudes. Significant cloud activity occurs at midlatitude (40°–75°) between the second half of spring and the first half of summer, both in the north ($L_s = 40°–110°$) and the south ($L_s = 230°–280°$). Short episodes of cloud formation are also observed at low latitudes (±30°) around the solstices and equinoxes. The supersaturation region is fed by vertical mixing, and the clouds that form at its lower boundary precipitate out of this region. These conditions are favorable to form frequent sporadic clouds at summer midlatitudes (Figure 6), with their lifespan limited by their sedimentation timescale of a few Earth days. The presence of less opaque and more persistent clouds is visible in the winter midlatitudes (north: $L_s = 70°–180°$; south: $L_s = 250°–360°$), due to colder regions, which is similar to the formation of clouds at the winter tropical region as found by P. Rannou et al. (2006).

The distribution and temporal evolution of methane clouds are valuable indicators of the tropospheric circulation and thermal structure. The fluctuations visible in Figure 6 are due to wave activity in the troposphere, which influences atmospheric conditions and consequently cloud formation (J. L. Mitchell et al. 2011; J. M. Battalio & J. M. Lora 2021). The Cassini mission detected clouds at southern midlatitudes between 2004 and 2010, during the southern summer, and then at northern midlatitudes after 2016, at the onset of northern summer (E. P. Turtle et al. 2018). This is consistent with the model predicting clouds at summer midlatitude. However, only ethane clouds were observed in the winter hemisphere by Cassini (C. A. Griffith et al. 2006; M. E. Brown et al. 2010; S. Le Mouélic et al. 2012). This could be explained by the fact that winter hemisphere clouds are less opaque and form beneath a denser haze layer present during winter (see Section 3), thus preventing their detection by Cassini.

Figure 6 illustrates the methane surface precipitation as predicted by the model over a Titan year. With radii of several hundred micrometers, the clouds precipitate and in most cases evaporate below 5 km in undersaturated layers. Thus, the formation of clouds does not necessarily imply precipitation, which is consistent with observations, which have detected only two precipitation events at the surface (E. P. Turtle et al. 2009, 2011; J. W. Barnes et al. 2013). However, when clouds form low enough, with droplets large enough, precipitation could reach the surface. This is the case at summer midlatitudes, shortly after the formation of large clouds, forming at lower altitudes (5–10 km). Another event is visible around the equator and the southern spring equinox, at the time of the reversal of circulation, where vertical winds are less strong. The Titan PCM precipitation appears to align with that of other models (S. P. Faulk et al. 2017; J. M. Lora et al. 2019, 2022), although it is less intense (0.1–1 mm T<sub>day</sub><sup>-1</sup>).

The north pole region is a special case: the model predicts precipitation over the northern polar lakes for most of the year. The evaporation of lakes produces methane supersaturation in the model's lower layers. As a result, the methane in this region condenses below 1 km into a mist, which precipitates to the surface. Since the mist forms close to the surface and is less opaque than clouds, it does not appear on the opacity maps (Figure 6).

Without moist convection into the Titan PCM, methane condensation limits the exchange between the troposphere and stratosphere. Consequently, methane concentration remains fairly stable over time and latitude in the stratosphere, with no annual or seasonal fluctuations, as shown in Figure 7. In contrast, there is a seasonal variation in methane quantity within the troposphere, with higher concentrations in the summer hemisphere compared to the winter hemisphere.

Figure 7 also illustrates the meridional methane fluxes within different atmospheric layers. These fluxes are predominantly below 10 km, where eddy cells dominate (B. Charnay & S. Lebonnois 2012; J. M. Lora & J. L. Mitchell 2015; J. M. Lora 2024), and where most exchanges occur. Methane evaporates from the lakes, enriching the summer region, and is then redistributed by meridional fluxes toward the winter pole. In the region where clouds form (10–40 km), the fluxes are also directed toward the winter pole, except in an area at 10°N where the fluxes are directed toward the summer pole. This region promotes the enrichment of the midlatitudes and the formation of large summer clouds. Above 40 km, these fluxes
<figure><img src="image_15.png" /></figure>

**Figure 7.** Top: modeled latitudinal mixing ratios of methane in the troposphere (solid line) and the stratosphere (dashed line). Bottom: meridional mass fluxes of methane in the Titan PCM, below 10 km (black), between 10 and 40 km (red), and above 40 km (blue). Positive flux values indicate a northward flux, while negative values indicate a southward flux. The left panels represent the average during the northern summer and the right panels during the southern summer. The figures show no latitudinal variation in methane quantity in the stratosphere and a higher quantity of methane in the troposphere during the summer.

become negligible. These results corroborate the existence of a
hydrological-type cycle of methane located in Titan's tropo-
sphere, where lakes act as a source for atmospheric methane,
and precipitation as a sink.

## 4.2. The Case of Minor Species

The cycle of minor species on Titan differs from that of
methane. Trace compounds such as hydrocarbons (C₂H₂,
C₂H₆) and nitriles (HCN) are produced in the upper layers of
the atmosphere through photodissociation. These compounds
then settle and are transported by atmospheric circulation to the
stratosphere (see Paper I).

In the Titan PCM, cold regions beneath the main haze layer ($<3 \times 10^3$ Pa) create a supersaturation zone for minor species. The mixing ratio of these species intersects with their saturation profile, initiating condensation onto aerosols at $\sim$70 km for hydrocarbons (C₂H₂, C₂H₆) and at $\sim$90 km for HCN. Consequently, below 90 km, the abundance of trace compounds decreases significantly, as illustrated for C₂H₆ in Figure 8.

This process depletes the atmosphere of gaseous chemical
species and dry aerosols below the main haze layer, replacing
them with a mist of condensed minor species. This layer is
most pronounced between 30 and 60 km, where saturation is
highest, and spans from one pole to the other, as shown in
Figure 9. Below 30 km, temperatures rise enough to exit the
supersaturation zone. The chemical species evaporate below
this altitude, enriching regions closer to the surface with trace
compounds (see Figure 8).

Similar to haze, trace compounds enrich the winter pole,
causing significant stratospheric cooling during the polar night
through infrared emission (see Paper I). The model predicts a
stratosphere enriched with C₂H₆ up to 5 times more in winter
than in summer, attributed to pole-to-pole circulation transport
(Figure 8). With the combination of higher enrichment and
colder temperatures, C₂H₆ condenses at higher altitudes in
winter and reaches higher concentration rates, triggering the

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

**Figure 8.** Ethane mole fraction predicted by the Titan PCM at the north pole (black line) and south pole (red line) for (a) southern winter and (b) northern winter. Mole fractions computed using CIRS data have been plotted with their error bar at the same periods for comparison (C. Mathé et al. 2020; S. Vinatier et al. 2020). The altitude scale on the right axis is approximate.

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

**Figure 9.** Annual and zonal average of the atmospheric extinction coefficient at 0.671 µm. The extinction above 80 km corresponds to the extinction of the main haze layer composed of dry aerosols (see Figure 2). The extinction between 20 and 60 km is associated with the mist of condensate. Peaks of extinction below 25 km correspond to methane clouds, as described in Section 4.1.

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

**Figure 10.** Temporal evolution of the mass fraction of trace compound condensates (C₂H₂, C₂H₆, and HCN) modeled by the Titan PCM between 60 and 350 km.

formation of polar clouds. C₂H₂ and HCN exhibit behavior
similar to C₂H₆, with the difference that HCN condenses much
higher in winter, up to 300 km. These polar clouds have been
identified at high latitudes up to 300 ± 70 km by VIMS during
fall and winter, with spectral signatures of C₂H₆, C₆H₆, HCN,
or HC₃N (C. A. Griffith et al. 2006; N. A. Teanby et al. 2008b;
S. Rodriguez et al. 2011; S. Le Mouélic et al. 2012, 2018; R. J.
de Kok et al. 2014; R. A. West et al. 2016; S. Vinatier et al.
2018).

Figure 10 illustrates the mass fraction of condensed C₂H₂, C₂H₆, and HCN as predicted by the Titan PCM between 60 and 350 km. The visible structures confined to the winter pole are
<figure><img src="image_19.png" /><figcaption>**Figure 11.** (a) The third-order moment profiles for haze (black) and mist of condensed species (red). (b) Radius profiles of aerosols (black) and cloud droplets (red). (c) Opacity profiles at 0.671 µm for haze (black) and total (red). The Huygens/DISR vertical profile from L. R. Doose et al. (2016) is plotted for comparison (dashed black line). (d) Temperature profile predicted by the Titan PCM (solid line) compared to the Huygens/HASI vertical profile from M. Fulchignoni et al. (2005; dashed line). All profiles correspond to the time and latitude of Huygens' landing ($L_s = 305^\circ$; 10°S). The two horizontal lines represent the altitude of the modeled DHL (4 Pa) and the main haze layer ($1.5 \times 10^3$ Pa). The altitude scale on the right axis is approximate.</figcaption></figure>

similar to the polar clouds observed by Cassini (C. A. Griffith et al. 2006; N. A. Teanby et al. 2008b; S. Rodriguez et al. 2011; S. Le Mouélic et al. 2012, 2018; R. J. de Kok et al. 2014; R. A. West et al. 2016; S. Vinatier et al. 2018). In the model, polar cloud condensation begins in midfall ($L_s = 40^\circ$ in the south, $L_s = 240^\circ$ in the north) at high altitudes (~300 km) and is confined to the poles by the polar vortex (>80°). As fall progresses, the polar vortex weakens, and the enrichment of trace compounds descends deeper into the atmosphere. Consequently, the clouds develop at lower altitudes, between 60 and 100 km, and extend to midlatitudes (60°) during fall and winter. The tracking of these massive clouds by VIMS observations (S. Le Mouélic et al. 2018) shows a very good match with the model predictions (Figure 10).

During winter, the droplets that emerge in this environment measure a few tens of micrometers, resulting in sedimentation timescales spanning several hundred Earth days. This allows the cloud to be maintained throughout the season. Following the equinoxes, a secondary wind cell forms near the pole exiting winter, extending to 60° and rotating opposite to the main cell (see Paper I). This cell traps the gases and clouds accumulated during winter, which persist in the polar region into spring. Observations by Cassini/ISS and Cassini/CIRS confirm the existence of this persistent polar cloud (N. A. Teanby et al. 2008b; R. A. West et al. 2016). In the second half of spring, the polar cloud dissipates, reintroducing most minor species into the atmosphere through reevaporation and eliminating the rest by precipitation at the surface.

Finally, an asymmetry between the north and south polar clouds is visible in the model. The northern polar cloud is less pronounced due to higher winter temperatures in the north compared to the south. It results in lower supersaturation, faster reevaporation of chemical species, and a shorter cloud lifespan.

# 5. Thermal Structure of the Lower Atmosphere

In this section, we examine the effects of coupling haze and cloud microphysics with radiative transfer on tropospheric and surface temperatures. A study of stratospheric temperature was carried out in Paper I.

## 5.1. Atmospheric Vertical Structure

Figure 11 shows the structure of the atmospheric components around the region and period of Huygens' landing as predicted by the Titan PCM ($L_s = 305^\circ$; 10°S). The vertical

profile of haze and mist distribution is depicted in Figure 11(a). The represented field is the third-order moment—their volume ratio—of aerosols (haze) and CCN with condensed species (mist). Their corresponding radius is shown in Figure 11(b) and their corresponding opacity in Figure 11(c). Four distinct layers can be identified, corresponding to local maxima of $M_3$: the DHL at approximately 4 Pa, the main haze layer at around $1.5 \times 10^3$ Pa, the mist of minor condensed species between $3 \times 10^3$ and $4 \times 10^4$ Pa, and the layer where methane clouds form at $10^5$ Pa.

The extinction profile is compared with that of L. R. Doose et al. (2016) using DISR data, showing disparities between the model and observations. The excessive extinction above 5 Pa is attributed to the presence of the DHL in the model at 325 km, a lower altitude than that observed at 500 km (R. A. West et al. 2018; B. Seignovert et al. 2021). On the other hand, the extinction of the main haze layer above $2 \times 10^3$ Pa is lower compared to Doose's profile. Below, a slope change in the opacity profile is observed where the mist layer forms (Figure 11(c)), as also noted by L. R. Doose et al. (2016). Mist particles are slightly larger than those comprising haze, measuring between 4 and 5 µm, compared to 0.5 µm for aerosols (Figure 11(b)). Although the opacity profile remains underestimated compared to observations, incorporating the optical properties of droplets improves the model profile and brings it closer to observations. Around $10^5$ Pa, a methane cloud is visible in the volume and opacity profiles, made up of large particles.

This new extinction profile plays an important role in the temperature profile (Figure 11(d)). The improved radiative transfer (see Paper I), along with haze and cloud coupling, has significantly enhanced the atmospheric temperature profile. The presence of a stratopause (local temperature maximum) around 30 Pa and a mesopause (local temperature minimum) at approximately 10 Pa is noticeable. The mesopause is attributed to the DHL made up of small, highly absorbent aerosols that heat the upper atmosphere (Figures 11(b) and (d)). This haze layer, causing the model's temperature inversion at 10 Pa, forms at a lower altitude than the one observed at 500 km (R. A. West et al. 2018; B. Seignovert et al. 2021). Around the main haze layer, the temperature is overestimated, warmer by nearly 15 K at 100 km. We observe a strong correlation between extinction and temperature profiles, highlighting the model's temperature sensitivity to these couplings. It is noteworthy that, due to their short lifespan, methane clouds have a minimal impact on the tropospheric temperature profile in the model.

## 5.2. Tropospheric Temperature

The Titan PCM's tropospheric temperature profile is now constrained by a comprehensive extinction profile extending to the surface (Figure 11), accounting for the opacity of gases, haze, and clouds. As a result, the tropospheric temperature profile now closely replicates the profile observed by HASI (M. Fulchignoni et al. 2005), as illustrated in Figure 12. We identify a tropopause of 68.75 K at $1.4 \times 10^4$ Pa near the equator, and of 66.75 K at $8.5 \times 10^3$ Pa and 67.5 K at $1.4 \times 10^4$ Pa at the north and south poles, respectively.

The lower stratosphere and the troposphere exhibit radiative relaxation timescales that are sufficiently long to dampen the atmospheric seasonal response (B. Bézard et al. 2018). Consequently, tropospheric temperatures display minimal
<figure><img src="image_20.png" /></figure>

**Figure 12.** Tropospheric temperature profiles in the Titan PCM, averaged zonally and diurnally, around the season of the Huygens' landing. Three latitudes are shown: the equator in black, the north pole (85°N) in red, and the south pole (85°S) in blue. The vertical profile from HASI (M. Fulchignoni et al. 2005) is plotted for comparison in black dotted lines. The altitude scale on the right axis is approximate.

seasonal variation (Figure 12). The poles exhibit a cooler troposphere than the equatorial region, with the winter pole being approximately 0.75 K colder than the summer pole and 2 K colder than the equator. These colder regions at the poles promote the development of mist, as observed in Figure 9.

### 5.3. Surface Temperature

Figure 13(a) illustrates the annual variation in surface opacity, showing the solar flux absorbed by the atmospheric column. Surface opacity in the visible spectrum is primarily governed by the haze column, which follows atmospheric circulation and accumulates at the winter pole. Consequently, higher opacity is observed around the winter poles ($L_s = 120°–240°$ at the south pole and $L_s = 300°–60°$ at the north pole). There is a half-season delay in the model between the reversal of circulation after equinoxes ($L_s = 0°$ and $L_s = 180°$) and the decrease in opacity caused by the haze shifting from one pole to the other ($L_s = 60°$ in the north and $L_s = 240°$ in the south). This delay was observed after the northern spring equinox in 2009, but only for a quarter of the season.

Methane clouds also contribute to the overall opacity due to their size, especially at high latitudes at the beginning of summer, as observed at 60°–90°S between $L_s = 240°$ and $L_s = 300°$, or at midlatitudes (~30°) in the north for $L_s = 300°–360°$, and in the south $L_s = 90°–120°$. Additionally, polar clouds appear at the beginning of fall (south at $L_s = 50°$; north at $L_s = 230°$).

The surface temperature follows a seasonal cycle, influenced by surface opacity, cooling in areas where opacity is highest, and warming where the solar flux is greatest (Figure 13). There is a slight asymmetry between the northern and southern winter poles. The southern pole experiences a more pronounced winter, which can be attributed to Titan's eccentricity (with the southern winter coinciding with aphelion) and the increased presence of clouds that take longer to dissipate compared to the northern pole. Additionally, the radiative effect of methane clouds at midlatitudes is negligible and has no impact on surface temperature.

The surface temperature predicted by the PCM, approximately 0.5 K lower than the recorded value of 93.65 K by Huygens

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

**Figure 13.** (a) Annual variation in surface opacity ($\lambda = 0.671 \mu m$) predicted by the Titan PCM. (b) Corresponding variation in surface temperature. The black crosses represent the maximum insolation received by the surface.

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

**Figure 14.** Latitudinal surface temperature predicted by the Titan PCM in black (a) in the middle of the northern winter, (b) at the early northern spring, (c) and at the northern summer solstice. For comparison, the CIRS profiles from D. E. Jennings et al. (2019) are plotted in red.

(M. Fulchignoni et al. 2005), shows slightly underestimated contrasts in both the polar winter regions and equatorial regions compared to observations (D. E. Jennings et al. 2009, 2019). The model contrasts are around 3 K (see Figure 14), whereas observations show contrasts up to 4 K. However, the contrasts obtained at the northern summer solstice ($L_s = 90°$) align with observations from Voyager 1/InfraRed Interferometer Spectrometer (Figure 14(c)).

The updated value of surface thermal inertia results in slight improvements in the temporal variations of surface temperature compared to previous versions (S. Lebonnois et al. 2012). This change reveals atmospheric planetary waves on the maps (Figure 13), and the temporal evolution of the maximum surface temperature is corrected, with the temperature peak at $L_s = 310°$ moving toward the midsouthern latitudes. Ultimately, the coupling with the haze microphysics significantly enhances the surface temperature description compared to the previous model. The contrasts are more pronounced, and the surface temperature is, on average, 1.5 K warmer, aligning more closely with observed temperatures.

## 6. Conclusions

Based on the IPSL Titan GCM, we have presented the new Titan PCM, including a new microphysical model for haze (J. Burgalat & P. Rannou 2017) and clouds (J. Burgalat et al. 2014), and updates of the radiative transfer scheme and reaction scheme of the photochemical solver (J. Vatant d'Ollone 2020 and Paper I). The microphysical model offers a comprehensive analysis of the distribution of haze and clouds within Titan's atmosphere, and now establishes an extensive coupling with radiative transfer. This integration also delivers a detailed depiction of the thermal and dynamical structure of the lower atmosphere.

Several aspects related to haze and cloud cycles and their impact on the thermal profile of the atmosphere were highlighted through this simulation:
1. As in previous works (e.g., P. Rannou et al. 2006; J. L. Mitchell 2008; T. Tokano 2009; J. M. Lora & J. L. Mitchell 2015; S. P. Faulk et al. 2020), the Titan PCM has shown that Titan's troposphere is dominated by the complex interaction between methane, photochemical species, clouds, condensate mist, and liquid reservoirs at the surface.

2. A better description of the haze's seasonal variations is obtained. The DHL is present in the Titan PCM at 325 km for most of the year, collapsing and disappearing after the equinox because of the shift of the ascending region from one pole to the other, crossing the equator during the equinox.

3. We now have access to a better description of the methane cycle (flux balances both in Titan's atmosphere and on its surface). The methane evaporates from the lakes and condenses in the ascending branches of circulation between 5 and 25 km, at spring and summer midlatitude. Subsequently, it returns to the surface through precipitation.

4. Between 20 and 60 km, the atmosphere is characterized by a persistent mist consisting of minor condensed species. This condensation depletes the troposphere of dry aerosols and alters the optical properties of the atmosphere within this region.

5. The minor species reach the cold polar region following the descending branches of the circulation during winter. This process leads to the formation of hydrocarbon and nitrile clouds between 60 and 300 km. The persistence of a circulation cell after the spring equinox traps trace compounds at the level of the tropopause, thereby sustaining the polar cloud through part of the spring season.

6. The coupling of microphysics with radiative transfer considerably improves the thermal description of the atmosphere. The role of the DHL in the heating of the upper atmosphere is highlighted. Taking into account the new optical properties of clouds provides a tropospheric temperature profile that reproduces observations very well.

7. In the stratosphere, the climate is mainly dominated by the cycles of haze and photochemical species and by strong feedback with atmospheric circulation, whereas the methane cycle takes place in the troposphere and regulates the climate.

The model provides information on the sources and sinks of the various elements in the atmosphere (haze, gas, and clouds), as well as on the fluxes that govern these cycles. Nevertheless several improvement avenues can be considered for the model. The observations show convective methane clouds until 2010 October, in agreement with the precipitation distributions predicted by models incorporating moist convection (e.g., J. L. Mitchell et al. 2006, 2011; T. Schneider et al. 2012; J. M. Lora et al. 2015; J. M. Battalio et al. 2022), but not taken into account by the Titan PCM. Moist convection is expected to play a crucial role in the global methane cycle. This phenomenon could modify the structure and distribution of methane clouds and might explain the variation in methane observed in the stratosphere (e.g., E. Lellouch et al. 2014; P. Rannou et al. 2021, 2022). On the other hand, although models predict that ethane is the main product of methane photolysis, production rates in the upper atmosphere limit its abundance, and it saturates in the troposphere at much lower vapor pressures than methane, making it a minor species.

However, its abundance in Lake Ontario may be comparable to that of methane (R. H. Brown et al. 2008; M. Mastrogiuseppe et al. 2018). Thus, an improvement in microphysics would consist of taking into account the interactions between methane and ethane through condensation processes favoring the precipitation of ethane at the surface on long-term simulations. One final consideration would be to solve the deep-atmosphere equations to study the upper atmosphere with a nonconstant gravity field. This would make it possible to increase the upper limit of the model and improve the injection of aerosols at the top of the model, which is currently too low, leading to a layer of a DHL that is too deep.

Finally, the new Titan PCM enables an in-depth analysis of long-term variations in haze, methane cycles, and trace species cycles. In Paper I, we also present an analysis of the thermal structure, comparing our results with data collected during the Cassini mission. These investigations open the way for a deeper understanding of Titan's atmospheric cycles and establish solid foundations for future studies in this field, while providing essential constraints for future observations of Titan.

## Acknowledgments

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

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

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

## References

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

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

Ádámkovics, M., Mitchell, J. L., Hayes, A. G., et al. 2016, Icar, 270, 376

Anderson, C. M., & Samuelson, R. E. 2011, Icar, 212, 762

Anderson, C. M., Samuelson, R. E., Bjoraker, G. L., & Achterberg, R. K. 2010, Icar, 207, 914

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

Atreya, S. K., Adams, E. Y., Niemann, H. B., et al. 2006, P&SS, 54, 1177

Barnes, J. W., Buratti, B. J., Turtle, E. P., et al. 2013, PISci, 2, 1

Battalio, J. M., & Lora, J. M. 2021, NatAs, 5, 1139

Battalio, J. M., Lora, J. M., Rafkin, S., & Soto, A. 2022, Icar, 373, 114623

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

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

Birch, S. P. D., Hayes, A. G., Dietrich, W. E., et al. 2017, Icar, 282, 214

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

Borucki, W. J., Levin, Z., Whitten, R. C., et al. 1987, Icar, 72, 604

Brown, M. E., Roberts, J. E., & Schaller, E. L. 2010, Icar, 205, 571

Brown, R. H., Soderblom, L. A., Soderblom, J. M., et al. 2008, Natur, 454, 607

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

Cabane, M., Rannou, P., Chassefière, E., & Israel, G. 1993, P&SS, 41, 257

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

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

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., 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., Rosset, L., Vatant d'Ollone, J., et al. 2025, *PSI*, 6, 78

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

del Genio, A. D., Zhou, W., & Eichler, T. P. 1993, *Icar*, 101, 1

Doose, L. R., Karkoschka, E., Tomasko, M. G., & Anderson, C. M. 2016, *Icar*, 270, 355

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

Faulk, S. P., Lora, J. M., Mitchell, J. L., & Milly, P. C. D. 2020, *NatAs*, 4, 390

Faulk, S. P., Mitchell, J. L., Moon, S., & Lora, J. M. 2017, *NatGe*, 10, 827

Flasar, F. M. 1998, *Natur*, 395, 541

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

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

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

Gautier, T., Serigano, J., Das, K., et al. 2024, *A&A*, 690, A165

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

Griffith, C. A., Penteado, P., Baines, K., et al. 2005, *Sci*, 310, 474

Griffith, C. A., Penteado, P., Rannou, P., et al. 2006, *Sci*, 313, 1620

Griffith, C. A., Penteado, P., Rodriguez, S., et al. 2009, *ApJL*, 702, L105

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

Hanson, L. E., Waugh, D., Barth, E., & Anderson, C. M. 2023, *PSJ*, 4, 237

Hartwig, J., Meyerhofer, P., Lorenz, R., & Lemmon, E. 2018, *Icar*, 299, 175

Hayes, A. G., Aharonson, O., Callahan, P., et al. 2008, *GeoRL*, 35, L09204

Hayes, A. G. 2016, *AREPS*, 44, 57

Hayes, A. G., Lorenz, R. D., & Lunine, J. I. 2018, *NatGe*, 11, 306

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

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

Hourdin, F., Talagrand, O., Sadourny, R., et al. 1995, *Icar*, 117, 358

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

Jennings, D. E., Flasar, F. M., Kunde, V. G., et al. 2009, *ApJL*, 691, L103

Jennings, D. E., Tokano, T., Cottini, V., et al. 2019, *ApJL*, 877, L8

Karkoschka, E. 2016, *Icar*, 270, 326

Keesee, R. G. 1989, *JGR*, 94, 14683

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

Larson, E. J. L., Toon, O. B., West, R. A., & Friedson, A. J. 2015, *Icar*, 254, 122

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

Le Mouélic, S., Rannou, P., Rodriguez, S., et al. 2012, *P&SS*, 60, 86

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

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., Bézard, B., Flasar, F. M., et al. 2014, *Icar*, 231, 323

Li, J., Liu, D., Coustenis, A., & Liu, X. 2012, *P&SS*, 63, 150

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

Lora, J. M. 2024, *Icar*, 422, 116241

Lora, J. M., & Ádámkovics, M. 2017, *Icar*, 286, 270

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., & Mitchell, J. L. 2015, *GeoRL*, 42, 6213

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

Lunine, J. I., & Lorenz, R. D. 2009, *AREPS*, 37, 299

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

Malaska, M. J., Hodyss, R., Lunine, J. I., et al. 2017, *Icar*, 289, 94

Mastrogiuseppe, M., Hayes, A. G., Poggiali, V., et al. 2018, *Icar*, 300, 203

Mastrogiuseppe, M., Poggiali, V., Hayes, A. G., et al. 2019, *NatAs*, 3, 535

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

Mitchell, J. L. 2008, *JGRE*, 113, E08015

Mitchell, J. L., Ádámkovics, M., Caballero, R., & Turtle, E. P. 2011, *NatGe*, 4, 589

Mitchell, J. L., & Lora, J. M. 2016, *AREPS*, 44, 353

Mitchell, J. L., Pierrehumbert, R. T., Frierson, D. M. W., & Caballero, R. 2006, *PNAS*, 103, 18421

Montmessin, F., Rannou, P., & Cabane, M. 2002, *JGRE*, 107, 5037

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

Porco, C. C., Baker, E., Barbara, J., et al. 2005, *Natur*, 434, 159

Pruppacher, H., & Klett, J. 2004, *Microphysics of Clouds and Precipitation* (New York: Springer)

Pruppacher, H. R., Klett, J. D., & Wang, P. K. 1998, *AerST*, 28, 381

Rages, K., & Pollack, J. B. 1983, *Icar*, 55, 50

Rannou, P., Coutelier, M., Rey, M., & Vinatier, S. 2022, *A&A*, 666, A140

Rannou, P., Coutelier, M., Rivière, E., et al. 2021, *ApJ*, 922, 239

Rannou, P., Curtis, D., & Tolbert, M. A. 2019, *A&A*, 631, A151

Rannou, P., Hourdin, F., & McKay, C. P. 2002, *Natur*, 418, 853

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

Rannou, P., McKay, C. P., & Lorenz, R. D. 2003, *P&SS*, 51, 963

Rannou, P., Montmessin, F., Hourdin, F., & Lebonnois, S. 2006, *Sci*, 311, 201

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

Rodriguez, S., Le Mouélic, S., Rannou, P., et al. 2009, *Natur*, 459, 678

Rodriguez, S., Le Mouélic, S., Rannou, P., et al. 2011, *Icar*, 216, 89

Roe, H. G. 2012, *AREPS*, 40, 355

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

Schaller, E. L., Brown, M. E., Roe, H. G., & Bouchez, A. H. 2006, *Icar*, 182, 224

Schaller, E. L., Roe, H. G., Schneider, T., & Brown, M. E. 2009, *Natur*, 460, 873

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, *JQSRT*, 113, 1020

Schneider, T., Graves, S. D. B., Schaller, E. L., & Brown, M. E. 2012, *Natur*, 481, 58

Seignovert, B., Rannou, P., West, R. A., & Vinatier, S. 2021, *ApJ*, 907, 36

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

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

Stofan, E. R., Elachi, C., Lunine, J. I., et al. 2007, *Natur*, 445, 61

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

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. 2018, *GeoRL*, 46, 3079

Thompson, W. R., Zollweg, J. A., & Gabis, D. H. 1992, *Icar*, 97, 187

Tokano, T. 2005, *AdSpR*, 36, 286

Tokano, T. 2009, *Icar*, 204, 619

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

Tokano, T., Neubauer, F. M., Laube, M., & McKay, C. P. 2001, *Icar*, 153, 130

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, 16287

Trainer, M. G., Pavlov, A. A., Dewitt, H. L., et al. 2006, *PNAS*, 103, 18035

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

Turtle, E. P., Perry, J. E., Barbara, J. M., et al. 2018, *GeoRL*, 45, 5320

Turtle, E. P., Perry, J. E., Hayes, A. G., et al. 2011, *Sci*, 331, 1414

Turtle, E. P., Perry, J. E., McEwen, A. S., et al. 2009, *GeoRL*, 36, L02204

Vatant d'Ollone, J. 2020, Thesis, Sorbonne Univ. https://theses.hal.science/tel-03561601

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

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., Klippenstein, S. J., Hörst, S. M., & Lavvas, P. 2019, *Icar*, 324, 120

West, R. A., Balloch, J., Dumont, P., et al. 2011, *GeoRL*, 38, L06204

West, R. A., Del Genio, A. D., Barbara, J. M., et al. 2016, *Icar*, 270, 399

West, R. A., Seignovert, B., Rannou, P., et al. 2018, *NatAs*, 2, 495

West, R. A., & Smith, P. H. 1991, *Icar*, 90, 330

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