Icarus 218 (2012) 707–722

Contents lists available at SciVerse ScienceDirect

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

Titan global climate model: A new 3-dimensional version of the IPSL Titan GCM
Sébastien Lebonnois a,⇑, Jérémie Burgalat b, Pascal Rannou b,c, Benjamin Charnay aaLaboratoire de Météorologie Dynamique, IPSL, CNRS/UPMC, Box 99, F-75252 Paris Cedex 05, France
GSMA, UMR CNRS 6089, Université de Reims Champagne-Ardenne, Francec
LATMOS, IPSL, CNRS/UVSQ, Guyancourt, Franceb

a r t i c l e

i n f oArticle history:
Received 14 April 2011
Revised 17 November 2011
Accepted 25 November 2011
Available online 8 December 2011
Keywords:
Titan
Atmospheres, Dynamics

a b s t r a c t
We have developed a new 3-dimensional climate model for Titan’s atmosphere, using the physics of the
IPSL Titan 2-dimensional climate model with the current version of the LMDZ General Circulation Model
dynamical core. Microphysics and photochemistry are still computed as zonal averages. This GCM covers
altitudes from surface to 500 km altitude, with barotropic waves now being resolved and the diurnal
cycle included. The boundary layer scheme has been changed, yielding a strong improvement in the tropospheric zonal wind proﬁle modeled at Huygens descent position and season. The potential temperature
proﬁle is fairly consistent with Huygens observations in the lowest 10 km. The latitudinal proﬁle of thenear-surface temperature is close to observed values. The minimum of zonal wind observed by the Huygens probe just above the tropopause is also present in these simulations, and its origin is discussed bycomparing solar heating and dynamical transport of energy. The stratospheric temperature and windﬁelds are consistent with our previous works. Compared to observations, the zonal wind peak is too weak(around 120 m/s) and too low (around 200 km). The temperature structures appear to be compressed inaltitude, and depart strongly from observations in the upper stratosphere. These discrepancies are correlated, and most probably related to the altitude of the haze production. The model produces a detachedhaze layer located more than 150 km lower than observed by the Cassini instruments. This low production altitude is due to the current position of the GCM upper boundary. However, the temporal behaviourof the detached haze layer in the model may explain the seasonal differences observed between Cassini
and Voyager 1. The waves present in the GCM are analyzed, together with their respective roles in the
angular momentum budget. Though the role of the mean meridional circulation in momentum transport
is similar to previous work, and the transport by barotropic waves is clearly seen in the stratosphere, asigniﬁcant part of the transport at high latitudes is done all year long through low-frequency troposphericwaves that may be baroclinic waves.
Ó 2011 Elsevier Inc. All rights reserved.

1. Introduction
The Cassini mission to Saturn and Titan continues to send many
data, that greatly improves our knowledge of Titan’s atmosphereand surface. As the years go by, the season in the Saturnian system
has been moving from northern winter at Cassini’s arrival andHuygens mission (Ls  300°) into northern spring at present. This
allows the retrieval of seasonal effects in the distribution oftemperatures, tropospheric clouds, stratospheric haze andchemical compounds.
The dominant features of the temperature and dynamics ofTitan’s atmosphere observed recently with the CIRS instrument
on board the Cassini mission have been reviewed in Flasar andAchterberg (2009). In the temperature structure at this season, a
latitudinal contrast is observed in the stratosphere, with the equa⇑ Corresponding author.
E-mail address: Sebastien.Lebonnois@lmd.jussieu.fr (S. Lebonnois).
0019-1035/$ - see front matter Ó 2011 Elsevier Inc. All rights reserved.
doi:10.1016/j.icarus.2011.11.032

tor slightly warmer than the southern (summer) polar region, andthe northern (winter) polar region colder by 20 K. The warmest
region of the stratopause is located over the winter polar region(where the stratopause occurs at higher altitudes than at lowlatitudes, 0.01 hPa versus 0.1 hPa), indicative of a subsidence over
the winter pole. From the temperature structure, the cyclostrophic
winds are derived and show a strong jet in the winter hemisphere,with stratospheric zonal winds up to 190 m/s (Achterberg et al.,2008b). Another feature deduced from the CIRS dataset is the factthat the polar vortices appear to be tilted by a few degrees relative
to the spin axis of Titan (Achterberg et al., 2008a).
The CIRS spectra have also determined the stratosphericdistributions of photochemical compounds (Teanby et al., 2009b;
Coustenis et al., 2010). Latitudinal proﬁles in the stratosphere are
retrieved from nadir observations (e.g. Coustenis et al., 2007;
Teanby et al., 2006, 2009a), as well as vertical proﬁles at several
latitudes from limb observations (e.g. Teanby et al., 2007; Vinatieret al., 2007, 2010). These observations have determined the708

S. Lebonnois et al. / Icarus 218 (2012) 707–722

latitude–pressure distribution of many chemical compounds in Titan’s troposphere, conﬁrming the enrichment in the northern(winter) polar region compared to other latitudes, already observed at a different season (just after northern spring equinox)
with the IRIS instrument on board Voyager 1 (Coustenis andBézard, 1995). This characteristic feature is consistent with thesubsidence within the winter polar vortex, as discussed by Teanby
et al. (2008a). Comparing IRIS/Voyager 1 and CIRS/Cassini datasets
enables to investigate seasonal variations in the stratospheric composition (Teanby et al., 2008b; Coustenis et al., 2010).
On January 14th, 2005, the Huygens probe landed on Titan’ssurface. During its descent, the data acquired by the onboard
instruments determined vertical proﬁles of atmospheric temperature, pressure (Fulchignoni et al., 2005), wind speed (Bird et al.,2005; Folkner et al., 2006), haze (Tomasko et al., 2005) and composition (Niemann et al., 2005) at 10°S latitude. Properties of theboundary layer have also been inferred from Huygens data (Tokano
et al., 2006). Huygens results are reviewed by Lebreton et al.
(2009).
Stratospheric winds and composition have also been observed
from Earth using different techniques: with the Infrared Space
Observatory (Coustenis et al., 2003), with submillimeterradiotelescopes (Marten et al., 2002; Moreno et al., 2005), withhigh-resolution infrared heterodyne interferometry (Kostiuk
et al., 2001, 2005, 2006), with Doppler shift in reﬂected solar lines(Luz et al., 2005), and also through analysis of Titan stellar occultations (Hubbard et al., 1993; Bouchez, 2004; Sicardy et al., 2006).
These Earth-based observations, though providing less spatial
resolution, allow for a larger sampling of the seasonal variations.
Starting with Del Genio et al. (1993) and Hourdin et al. (1995),
several Titan atmosphere General Circulation Models have beendeveloped, based on different Earth GCMs: the Köln model (Tokano
et al., 1999), which has been used for many studies on Titan’s troposphere (for a review, see Tokano, 2009); a 2-dimensional model
mainly used to study the tropospheric methane cycle (Mitchell
et al., 2006); a non-hydrostatic GCM developed in Russia (Mingalev
et al., 2006); a Titan version of the PlanetWRF GCM (Richardson
et al., 2007); a Titan version of the CAM model (Friedson et al.,2009). Most of these models have difﬁculties reproducing the highzonal winds observed in Titan’s stratosphere. A recent version of
the TitanWRF model (Newman et al., 2011) obtains strong superrotation comparable to observations, and the authors demonstrate
how sensitive it is to assumed values of horizontal dissipation.
Development of the LMD Titan atmospheric model started inearly 1990s, based on the GCM dynamical core developed for the
Earth’s atmosphere. Hourdin et al. (1995) managed to spin up
the atmosphere of Titan from rest up to superrotation with zonal
winds comparable to observed values. The vertical range of thisGCM was from the surface up to 250 km. Radiative transfer wascomputed using a prescribed vertical proﬁle of haze and gases(inheritated from McKay et al. (1989)) horizontally uniform on
the planet. In order to include coupling with haze, Rannou et al.
(2002, 2004) developed a 2-dimensional axisymmetric (latitude–
altitude) version of the model that demonstrated strong interaction between haze distribution and dynamics. This model covered
altitudes up to 500 km and it was able to interpret many featuresof the haze layer, in particular the detached haze observed at
350 km altitude at the time of Voyager 1. This IPSL 2-dimensional
climate model (2D-CM) was further developed to include the photochemical model of Lebonnois et al. (2001), and the effects of thecoupling between chemical composition and dynamics were studied (Lebonnois et al., 2003; Hourdin et al., 2004). A recent comparison of the CIRS/Cassini results with the latest version of the IPSL
2D-CM was published by Crespin et al. (2008). The model reproduces the thermal structure observed in the troposphere andstratosphere, the abundance and vertical proﬁles of most chemicalcompounds in the stratosphere, and their enrichment in the winterpolar region. A reference simulation may be found on the web(Rannou et al., 2005), and a review of the interactions between
haze, dynamics and composition may be found in Lebonnois
et al. (2009).
This 2-dimensional climate model uses a speciﬁc parametrization to take into account the horizontal mixing by barotropic
waves, which were obtained in the 3-dimensional GCM of Hourdin
et al. (1995) but could not be present in the 2-dimensional reduction of the model. This parametrization was developed by Luz and
Hourdin (2003) and Luz et al. (2003). It uses diagnosed barotropic
instabilities to evaluate the horizontal mixing coefﬁcient associated with the corresponding waves. This mixing is part of theangular momentum budget and plays a key role in the superrotation, but no ad hoc tuning was done to get the successful simulation of observed features detailed previously.
However, the best way to conﬁrm the conclusions from the 2dimensional climate model is to take advantage of the increasein computing power available to develop a new version of the 3dimensional GCM. This model is now operational, simulations havebeen run and analysis is presented below. The model is detailed in
Section 2. In Section 3, the results for stratospheric circulation arecompared to observations and to the previous 2D-CM results. Some
aspects of the tropospheric circulation are described in Section 4.

2. Description of the GCM
2.1. Dynamical core
The LMDZ4 dynamical core, currently used for the Earth(Hourdin et al., 2006) and for Venus (Lebonnois et al., 2010) has
been adapted for Titan’s atmosphere. This dynamical core is basedon a ﬁnite-difference discretization scheme that conserves bothpotential enstrophy for barotropic nondivergent ﬂows, and total
angular momentum for axisymmetric ﬂows. A longitudinal ﬁlter
is applied in polar regions (poleward of 60° latitude) to limit theeffective resolution to that at 60°. The model is currently used witha resolution of 32 longitudes by 48 latitudes (11.25°  3.75°), but aresolution of 64  48 has also been tested to estimate the impact of
the longitudinal resolution. This increased resolution did not signiﬁcantly affect the results discussed in this paper. The dynamical
time-step is 3 min (reduced to 1.5 min for the higher longitudinal
resolution). The model uses a leapfrog time integration scheme,with a periodic predictor–corrector time-step. Horizontal
dissipation is done similarly to Earth and Venus, using an iterated
laplacian with a time constant of 2  105 s below roughly 400 hPa(25 km altitude), 1  105 s between 400 hPa and 4 hPa (25 and130 km), and 2  104 s above, with smoothed transitions over a
few layers.
The vertical grid of the model is based on hybrid coordinates,
though for the moment no topography is included, so model levels
basically coincide with constant-pressure levels, with a surfacepressure of 1.467 bar (Fulchignoni et al., 2005). We use 55 verticallevels, from the surface to roughly 550 km altitude. This grid isclose to the one used in the 2-dimensional studies, listed in Table
1 of Rannou et al. (2005). Note that the altitudes listed in that tablewere computed based on a constant gravity g. When taking into account the thickness of the atmosphere compared to Titan’s radius
and the corresponding variation of g with altitude, as it is the casein this GCM, the real altitudes of the uppermost levels are greaterthan those indicated in Rannou et al. (2005). Pressures and altitudes of layer mid-points are displayed in Table 1. In the top threelayers of the GCM, a sponge layer is included, with Rayleigh friction
damping horizontal winds to zero. Time constants are 2.5  104 s(0.3 Earth days) for the top layer, then 5  104 and 1  105 s forS. Lebonnois et al. / Icarus 218 (2012) 707–722

709

Table 1
Globally averaged pressure and altitude for the 55 vertical layers of the GCM.
Layer

Pressure(Pa)

Altitude
(km)

Layer

Pressure(Pa)

Altitude
(km)

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28

1.464  105
1.456  105
1.439  105
1.401  105
1.329  105
1.216  105
1.080  105
9.48  104
8.21  104
6.97  104
5.80  104
4.75  104
3.83  104
3.05  104
2.40  104
1.87  104
1.44  104
1.10  104
8.39  103
6.34  103
4.77  103
3.58  103
2.69  103
2.02  103
1.52  103
1.14  103
8.50  102
6.36  102

0.035
0.147
0.39
0.93
1.98
3.74
6.05
8.55
11.2
14.3
17.5
21.0
24.6
28.4
32.3
36.3
40.4
44.7
49.2
53.9
58.9
64.6
71.2
78.6
86.7
95.3
104
114

29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55

4.75  102
3.54  102
2.64  102
1.97  102
1.47  102
1.10  102
82
61
45
34
25
19
14
10
7.8
5.7
4.2
3.1
2.2
1.6
1.2
1.0
0.74
0.46
0.29
0.20
0.14

124
134
145
156
167
179
191
203
215
228
241
255
268
282
297
313
330
348
367
387
405
418
437
467
497
522
546

Fig. 1. Modeled latitude-pressure distribution of the zonally averaged temperature
after the Cassini arrival (Ls  305°). The altitude scale on the right axis isapproximate.

the two layers below. These values are identical to those used in
the IPSL 2D-CM.
Though gravitational tides have recently been added to ourmodel (Charnay and Lebonnois, 2011), the simulations presented
here do not include them. The impact on the results discussed hereis small.
2.2. Physical parametrizations
The physical parametrizations have been mostly taken from the
IPSL 2D-CM (see Rannou et al., 2004, 2005; Crespin et al., 2008).
We will discuss here the most important differences.
Though the radiative transfer is similar, the 3-dimensional
geometry allows inclusion of the diurnal cycle. The radiative transfer is computed frequently, 10 times per Titan day in these simulations. Though a higher frequency has been used for a study of thediurnal cycle near the surface (Charnay and Lebonnois, 2011), theinﬂuence of this parameter on the results discussed in this paper
is negligeable. Even without any diurnal cycle, results are close
to those presented here.
For the boundary layer scheme, we use a ‘‘Mellor and Yamada’’parameterization (Mellor and Yamada, 1982), taken from the Earthversion of the LMD GCM. This parameterization is fully describedin the Appendix B of Hourdin et al. (2002). The surface drag coefﬁcient is Cd = (0.4/ln(1 + z1/z0))2, where z1 is the altitude of the center of the ﬁrst layer (equal to 10 m), and z0 is the roughness
coefﬁcient, taken equal to 0.5 cm, as estimated from Huygens databy Tokano et al. (2006). The value of Cd is then equal to the value
deduced from Huygens data, 2.8  103 (Tokano et al., 2006).
The soil model used here for thermal conduction in the soil isthe same 11-layer model as for the other GCMs developed at
LMD (Hourdin et al., 1993). The model depends on the thermal
inertia of the surface, chosen here to be I = 2000 J m2 s0.5 K1. A
simulation was also performed with I = 340 J m2 s0.5 K1. The
inﬂuence on the results discussed in this paper is small. However,thermal inertia affects the temporal evolution of the surface tem-

Fig. 2. Modeled vertical proﬁles of temperature at latitudes 60S, equator, 60°N and
80°N (Ls  305°). The HASI/Huygens vertical proﬁle (10S) is plotted for comparison(bold line).

perature, as well as its latitudinal proﬁle at a given season. Other
surface parameters include the albedo (0.15) and the emissivity(0.95).
2.3. Haze microphysics and photochemistry
Due to computational cost, it is difﬁcult to compute photochemical and haze microphysical tendencies at every longitude.
However, the relative variation of the haze layer in the longitudinal
direction is extremely weak. No variation have been reported inobservations at the scale of the model grid (500 km). Therefore,
zonally averaging the microphysics is not anticipated to produce
signiﬁcant loss of information in the model, though this approximation precludes study of longitudinal variations of atmosphericcomposition related to the diurnal variability. We adopted the following technique: at each time-step, the 3-dimensional ﬁelds ofmicrophysical and photochemical tracers are zonally averaged, before computing the tendencies as diurnal averages in the same
microphysical and photochemical modules as in the IPSL 2D-CM(Rannou et al., 2004; Crespin et al., 2008). The tendencies are thenconverted back to 3-dimensional ﬁelds proportionally to the corre-

710

S. Lebonnois et al. / Icarus 218 (2012) 707–722

Fig. 3. Latitudinal proﬁles of modeled temperatures: (a) at 1, 2 and 5 hPa, Ls 305° to compare to CIRS/Cassini data at 1.8 hPa (crosses with error bars, Flasar et al., 2005) and (b)at 0.4, 1 and 2 hPa, Ls 9° to compare to IRIS/Voyager 1 data at 1 hPa (circles with error bars, Coustenis and Bézard, 1995).

Fig. 4. Modeled latitude-pressure distribution of the zonally averaged opacity ﬁeld at 700 nm, for (a) northern spring equinox and (b) northern winter solstice. The altitude
scale on the right axis is approximate.

No cloud microphysics was included in this GCM conﬁguration.
This work is underway for methane, ethane, and acetylene (basedon Rannou et al. (2006)) and will be the subject of a future publication. Inclusion of a cloud model would affect the results in the
troposphere, but essentially at second order. From preliminary results including the clouds (out of the scope of this paper), themethane cycle does not alter signiﬁcantly the results described
here.
As a starting point, we used the haze and composition distributions taken from Crespin et al. (2008), at northern spring equinox(Ls = 0°). However, haze opacity has been artiﬁcially reduced everywhere at the beginning of the simulation, to let the 3-dimensional
GCM reach its own equilibrium regime. Production has been
maintained at the value used in the 2D simulations, 1.0  1012
kg m2 s1, peaking at 1 Pa (roughly 400 km altitude).
Fig. 5. Annual variations of the modeled vertical proﬁle of the zonally averaged
equatorial opacity at 700 nm. The altitude scale on the right axis is approximate.

sponding tracer ﬁeld. The total mass is conserved in this scheme.
Only the microphysical and photochemical modules are affectedby this zonal average; haze and chemical compounds are advectedby the full 3-dimensional circulation. Photochemistry is computed
only once per Titan day (as in 2D), while the microphysics is computed 200 times per day.

2.4. Initial conditions and simulation parameters
The initial ﬁelds were all taken from the two-dimensional simulations described in Crespin et al. (2008) at northern spring equinox, including all components of the wind, and the temperature.
These 2-dimensional ﬁelds were used as uniform in longitude. Thismeans the model is not started from rest. The surface pressure wasadjusted to the Huygens value, 1.467 bar (Fulchignoni et al., 2005).
The simulations were run for 12 Titan years. During the ﬁrst
year, the total angular momentum decreased by 40%, due to the711

S. Lebonnois et al. / Icarus 218 (2012) 707–722

signiﬁcant adjustment of the wind proﬁle in the troposphere. After
that, its relative variations are limited to roughly 1.5% the second
year, 0.8% the third year, and less than 0.5% the following years.
The stability from 1 year to the other is good, though longer simulations may still yield some small variations. Test variations of input parameters have some impact on the total angular momentum(up to a few percent) depending on the inﬂuence of the modiﬁed
parameter on tropospheric wind speeds. Output is averaged over
1 Titan ‘‘week’’ of 11 Titan days, giving 61 weeks/Titan year (i.e.
673 Titan days, the latest 2 weeks having 12 days instead of 11).
A test simulation started from rest has been run with decoupledhaze, i.e. a ﬁxed vertical proﬁle of haze opacity, uniform in the horizontal. This haze distribution was taken from an averaged distribution obtained with the 2-dimensional model. It run for
15 Titan years, which may not be long enough (Newman et al.,2011), with the evolution of the total angular momentum around
0.7% per year after this time. The superrotation was much smaller
than in the simulation started from the 2-dimensional simulations,with maximum zonal wind speed limited to 40 m/s. Further analysis will be done in future work, in particular the inﬂuence of thealtitude distribution of the haze (strongly affecting the temperature structure), as well as the horizontal dissipation scheme, asNewman et al. (2011) showed the extremely large inﬂuence thisparameter may have on the spin-up. This is also a question understudy in Venus atmosphere modeling (Lebonnois et al., 2010,2011).
3. Results in the stratosphere
The modeled ﬁelds are compared here to available observations.
Emphasis is also given to differences with the ISPL 2D-CM results.
Discussions concerning this model refer to the work of Crespin
et al. (2008) (hereafter C08).
3.1. Temperature structure
The zonally averaged temperature structure in northern winter(Ls  300°) is shown in Figs. 1, 2 and 3a. As in the IPSL 2D-CM simulations, the temperature vertical proﬁles below 2 hPa (160 km)
compare quite well with observations (see Figs. 1–3, C08). However, above this level, modeled stratospheric temperatures are signiﬁcantly colder than observations. In addition, the absence of astratopause and mesosphere is still present, pointing to the same(a)radiative transfer problem as discussed in C08. The latitudinal features are similar to observations (Achterberg et al., 2008b) and to
the 2D simulations, but the polar contrasts and structures in thewinter (northern) stratosphere appear lowered in altitude in the
10–102 hPa region compared to observations. This difference
may be related to the haze distribution over the winter polar
region and to the misrepresentation of the atmospheric structure
above 0.3 hPa (250 km). As discussed in C08, one reason for thismisrepresentation is the modeled altitude of the detached hazelayer (as shown in the next section), which is much too low
compared to observations due to the model limitations in thealtitude domain.
The thermal structure after spring equinox, compared to
IRIS/Voyager 1 data, is not as accurate as in previous versions ofour model, as shown in Fig. 3b (to be compared to Fig. 14, Rannou
et al., 2004). This disagreement with observations must be related
to the haze distribution and associated circulation in this GCM
which are not fully satisfying, as discussed in more detail in the following sections. The temperature contrasts are smaller than in the2D simulations, especially around equinox, with consequences forwind speeds.
3.2. Haze distribution
The formation of the detached haze layer has been discussed inconnection with the meridional circulation and production region
by Rannou et al. (2002, 2004). The mechanism is the following.
In the production zone (around 1 Pa/400 km altitude), the particles
grow and are blown towards the winter pole by the mean meridional circulation. In the winter hemisphere, they settle down, leave
the production region and grow as fractal aggregates. Though they
accumulate dominantly over the winter polar region, they are
mixed in the stratosphere by the stratospheric dynamics. In thesummer hemisphere, ascending motions help maintain a population of particles in the detached haze layer, corresponding to the
top of the dominant meridional circulation, that includes both
small spheres from the formation region and fractal aggregatesfrom the main haze layer.
In the 3-dimensional simulations presented here, the resulting
detached layer is located around 6–3 Pa (300–350 km altitude)
most of the time, slightly lower than in the 2D-CM simulations(C08). This is consistent with smaller latitudinal temperature contrasts and weaker winds (see next section). The meridional shapeof the haze opacity at 700 nm is shown in Fig. 4 for northern spring(b)

Fig. 6. (a) Approximate altitude of the maximum opacity at 700 nm in the detached haze layer as a function of the solar longitude. (b) Vertical proﬁle of the zonally averaged
equatorial opacity, during northern winter and beginning of spring.

712

S. Lebonnois et al. / Icarus 218 (2012) 707–722

equinox and northern winter solstice. This shape is very close to
the results obtained with the 2D-CM. The position of the detachedhaze layer is closely related to the temperature increase above the
10 Pa level, which also affects the meridional circulation in a feedback loop. This increase of temperature right at the detached hazelayer position is also seen in the Huygens observations (Fulchignoni et al., 2005).
In Fig. 5, the evolution of the vertical proﬁle of the opacity isshown at the equator during 1 Titan year (year 12 of our simulation). The detached haze layer disappears around Ls  30° and210°, before building up again just after the solstices. As illustrated in Fig. 6a, it is mostly located in a constant altitude range(310–340 km in this simulation). The evolution of this haze layer
during northern winter and spring is plotted in Fig. 6b. It illustrates
the fact that quickly after equinox, the detached haze layer moves
downward before disappearing. Qualitatively, this might explain
the difference in altitude of the detached haze layer between the
Cassini mission (roughly 510 km, Ls  300° to 360°) and the Voyager 1 observation (around 350 km, Ls = 9°), despite a wrong altitude. The timing of this descent of the detached haze layer
corresponds to the position of the ascending motions around the
10 Pa level during this season. The ascending branch joining thetwo Hadley cells is located over the summer polar region, thenshifts equatorward to reach the equator at equinox. It gradually

Fig. 7. Mean zonal wind ﬁeld (colors) obtained with the simulations during year 6 at several seasons: (a) Ls = 9°, Voyager 1; (b) Ls = 124°, close to the period of the 1989 stellar
occultation; (c) Ls = 180°, northern autumn equinox. (d) Ls = 255°, close to the period of the 2001 stellar occultation; (e) Ls = 300°, beginning of Cassini observations and close
to the 2003 stellar occultation; (f) Ls = 330°, during Cassini observations. The zonally averaged meridional stream function is shown as black contours (solid line is clockwise
rotation, dashed line is anti-clockwise, unit is 109 kg/s). The altitude scale on the right axis is approximate.

S. Lebonnois et al. / Icarus 218 (2012) 707–722

shifts further toward higher latitudes, reaching the future summer
polar region after Ls = 30°. During this shift, the detached haze layersinks down into the main haze layer over both poles. At the equator, it also goes down, and even disappears everywhere well into
the equinox season, before building up again over the spring polar
regions when the solstice approaches. This behavior of thedetached haze layer is model dependent, of course, but it helps explain why this layer was observed lower after equinox (Voyager 1)than during the solstice season (Cassini). Since the Cassini mission
is monitoring the detached haze layer during equinox and the following season, it should provide extremely valuable informations
on the variations of the detached haze layer altitude during thiscrucial period.
The discrepancy between the modeled and observed positionsof the detached haze layer must certainly have an impact on theposition and amplitude of the zonal wind jet. We need to extend
the model vertically to improve this aspect of the simulations.
However, this is not straightforward to do, because of the atmospheric thickness (already far beyond the acceptable limit for theshallow water approximation (ZC  RT) where ZC is the verticalextension of the model and RT is the radius of Titan) and becauseof the radiative transfer (HCN needs to be taken into account in arealistic way, as demonstrated in Lebonnois et al. (2003)).
3.3. Zonal winds
The zonal wind ﬁeld is presented at several seasons in Fig. 7,
together with the meridional stream function. For Ls = 300–330°,
this circulation may be compared to the thermal winds retrieved
from CIRS/Cassini data (Achterberg et al., 2008b). The jet is both
too low and weaker compared to the retrieved thermal wind andto the 2D modeled zonal wind (C08), which is consistent withthe disagreement obtained in the thermal structure at the same
season (Fig. 1).
The seasonal variations of the superrotation have been observed
at least partially thanks to Earth-based observations, see Table 1 in
Kostiuk et al. (2010) and Fig. 13 in Sicardy et al. (2006). The modeled seasonal variations of the zonal winds at altitudes around
145–200 km (300–60 hPa) and 250–310 km (20–6 hPa) are shownin Fig. 8. Though the higher region should be considered for comparison to most observations, the lower position of the modeled
jet at the Cassini season indicates that the lower region of themodel may be more relevant for comparison.

713

Fig. 9. Vertical proﬁle of the zonal wind at 10°S, at a season close to the Huygens
probe landing. Huygens observed proﬁle is indicated with black dots.

The amplitude of the zonal winds are under-estimated in themodel, for all solar longitudes. During solstice seasons, modeledwinds are above 120 m/s around 200 km altitude (where the peakis located), while observations indicate a peak around 180–220 m/
s: (Hubbard et al., 1993), Ls = 128°; (Sicardy et al., 2006; Kostiuk
et al., 2005), Ls = 290°; (Achterberg et al., 2008b), Ls = 300–330°.
Right before northern summer solstice (Ls = 255°), Bouchez(2004) have retrieved a northern peak up to 220 m/s, but also indication of winds that are still around 160–200 m/s in the southernhemisphere. For this season, the model predicts the highest zonal
winds in the autumn hemisphere, but the winds in the springhemisphere have decreased around mid-season, before the observations reported in Bouchez (2004). This disagreement may be related to the seasonal variations in the position of the jet and to theamplitude question, as there are hints in Figs. 7d and 8a that strong
winds are still present in mid-southern latitudes at this season below 200 km altitude, though decreasing rapidly with time. Aroundequinox (Ls = 160–190°), (Kostiuk et al., 2001) reported high zonal
wind speeds (u = 210 ± 150 m/s). Again, modeled values are lower
than observed, but the winds obtained in the GCM at this season
are slightly higher than the peak of the solstice season, consistent
with the observed value compared to Cassini season retrievals.
As shown in Fig. 9, the wind proﬁle in the troposphere and lower stratosphere at the Huygens probe landing site position and sea-

Fig. 8. Seasonal variations of the latitudinal proﬁle of the mean zonal wind averaged in altitude regions (a) 145–200 km and (b) 250–300 km.

714

S. Lebonnois et al. / Icarus 218 (2012) 707–722

son is in reasonable agreement with the observed proﬁle. The
shape of the tropospheric zonal wind is signiﬁcantly improved
compared to the 2D simulations shown in C08. A minimum of zonal wind, as observed by Huygens around 20–40 hPa (60–75 km
altitude), is obtained in the model though at lower altitude (35–75 hPa; 50–65 km), as it was already the case in the 2D simulations(C08). This local minimum is not as deep as in the Huygens data.
For the radiative relaxation time, this region is a transition(Flasar et al., 1981). Above 10 hPa, it is less than half a Titan year,
so the seasonal variations of the solar forcing drives signiﬁcant seasonal variations in the temperature. Below 50 hPa, the radiative
relaxation time is more than 2 Titan years, so the temperature
structure is more sensitive to the annual mean of solar forcing.
This region corresponds also to the maximum of stability, withthe strongest temperature increase with altitude, just above thetropopause region. There, the meridional winds reverse compared
to dominant winds above, inducing a partial division in the Hadley-type cells (both for the pole-to-pole cell present most of theyear and for the two equator-to-pole cells present aroundequinoxes).
As shown in Hourdin et al. (1995), the temperature is very uniform in latitude in this region, both in radiative–convective simulations and in GCM simulations. In these ﬁrst GCM simulations of
the circulation in Titan’s atmosphere, a minimum of zonal wind
was already present in the 30–40 hPa region. In the simulations
presented here, the latitudinal contrasts of temperature are also

Fig. 11. Annual average of the latitudinal transport of angular momentum by mean
meridional circulation (MMC, dashed line), transient waves (dotted line) andstationary waves (dash-dotted line). Total is shown in solid line. Unit is 103 m3/s2,
positive values are northward.

very small, and they are even reversed at low-latitudes, with theequator temperature being colder than mid-latitude temperatures,at all seasons. This reversal DT (over a latitude region D/) is correlated with the decrease of zonal wind Du with altitude (Dz). They
are linked through the thermal wind equation:

Fig. 10. Annual average of (a) the solar heating rates, (b) the dynamical heating rates and (c) the infrared cooling rates in the lower stratosphere. Units are K/Titan day. The
blue contours (a and b) show the difference between the annual average of the temperature and its latitudinal mean (latitudinal contrast, in K). The black contours (b and c)
indicate the annual average of the zonal wind (in m/s). (For interpretation of the references to color in this ﬁgure legend, the reader is referred to the web version of thisarticle.)715

S. Lebonnois et al. / Icarus 218 (2012) 707–722

(a)

(b)

Fig. 12. Seasonal variations of the latitudinal transport of angular momentum by (a) mean meridional circulation and (b) transient waves. Unit is 103 m3/s2, positive values
are northward.

Fig. 13. Zonally averaged amplitude frequency spectrum for the zonal and meridional wind ﬁelds as a function of latitude, at pressure levels of 10 and 100 Pa. Chosen season
is Ls = 295°. The solid curve is the mean zonal rotation speed of the ﬂow (same unit as the frequency).

716

S. Lebonnois et al. / Icarus 218 (2012) 707–722
@
u2 tan /g @T
2Xu sin / þ
¼
;
@z
a
Ta @/

ð1Þ

where a = 2.575  106 m is the radius of Titan, X = 4.56  106 s1
is the rotation rate of Titan, / is the latitude, g = 1.35 m s2 is thegravity, z (m) is the altitude, T (K) is the temperature and u(m s1) is the zonal wind speed. This equation is presented here
in the thin atmosphere approximation. It is further approximatedto link DD/T to DDuz :Du
2u tan /g DT¼
2X sin / þ
:
a
Ta D/
Dz

ð2Þ

In the GCM simulation presented in this work, the following
values are estimated: u  16 m s1, T  90 K, /  p/12 (averaged
value for the region between the equator and 30°),
Du  4 m s1, Dz  1  104 m, DT  0.2 K and D/  p/6. The
ﬁrst term evaluates toDu
2u tan /¼ 2:3  109 s2 :2X sin / þaDz
The second term evaluates tog DT¼ 2:2  109 s2 :Ta D/
Given the approximations done, the agreement is quite good,
conﬁrming that the zonal wind and temperature ﬁelds are related
through the thermal wind equation. If we apply the same calcula-tion to the observed proﬁle, with u  20 m s1 (averaged value inthe decrease region), T  90 K, Du  40 m s1, Dz  1.5  104 m,
and /  sin /  tan /  D2/, DT may be evaluated as

DT ¼

Du 
u Ta
ðD/Þ2 ;
Xþ
a g
Dz

ð3Þ

i.e.

DT ¼ 5:65  ðD/Þ2 :If the meridional temperature inversion is located between
equator and mid-latitudes (D/ between 30 and 45°), this relation
should yield a latitudinal contrast in temperature around 1.5–
3.5 K. This may be possible to observe with the Cassini radio-occultation data (Schinder et al., 2011).
Fig. 10 shows the annual average of the solar heating rates,
heating due to dynamics, and cooling by infrared radiation (thesethree tendencies are balanced). Between the tropopause and
70 km, the solar heating at equator is compensated by the cooling
due to ascending motions, resulting in a slightly lower coolingneeded around equator than at mid-latitudes. Therefore, our model
suggests that this relative distribution of solar heating versus
dynamical tendencies forces a local minimum of temperature inthis region. Through the thermal wind equation, this local minimum induces a decrease in zonal wind with altitude, therefore
driving the local minimum in zonal winds.

Fig. 14. Zonally averaged amplitude frequency spectrum for the zonal and meridional wind ﬁelds as a function of pressure, at equator and 60°N. Chosen season is Ls = 295°.
The solid curve is the mean zonal rotation speed of the ﬂow (same unit as the frequency). The altitude scale on the right axis is approximate.

S. Lebonnois et al. / Icarus 218 (2012) 707–722

3.4. Transport of angular momentum and waves
Averaged over 1 Titan year, the transport of angular momentum
is quite similar to previous models (such as Hourdin et al., 1995).
Here, the angular momentum transport is computed for each latitude / as the integral over longitude k and pressure p of the quantity vua cos /dm, where dm is the element of mass (dm = a2cos /d/
dkdp/g), normalized by the mass:Rv uaR cos /dmdm3ðexpressed in m =s2 Þ:

As shown in Fig. 11, the overall poleward transport of angularmomentum by the mean meridional circulation (MMC) is partiallycompensated through transport by waves in the low- to mid-latitudes. Poleward of roughly 40°, angular momentum is transported
toward poles both by transient waves and by the mean meridional
circulation. The total poleward transport is balanced by exchangesof angular momentum between the surface and the atmosphere.
The zonally and annually averaged zonal wind at the surface isnegative below 30° latitude (transferring momentum from the surface into the atmosphere) and positive above (giving back momentum to the surface), as shown in Fig. 7. Fig. 12 shows the seasonalvariations of the transport by the MMC and waves. The evolution of
the MMC transport is consistent with the seasonal variations of the
Hadley cells. The transport by transient waves follows a more complex pattern: at mid- to high-latitudes, the transport is continuously poleward, while between 40°N and 40°S, it is directed717

towards the spring/summer pole, stronger in the autumn/winterhemisphere where it counteracts the transport by the MMC, andreversing a third of season after the equinox. Signiﬁcant variability
is apparent in this horizontal transport by transients. This is consistent with the results of Newman et al. (2011), who ﬁnd that stratospheric waves transporting angular momentum are particularly
active during some episodes in the fall and winter seasons.
To analyze the waves present in our simulation, high-frequency
output was recorded for 16 Titan days, 50 points per day, atLs  295° where strong wave activity is visible. The eddy components of u and v at a given point in the atmosphere (u0 (t) andv0 (t)) were computed as the difference between the ﬁelds u(t)and v(t) and their time average at this point. For each point in
the atmosphere, a fast fourier transform was then performed onu0 (t) and v0 (t) after multiplication by a triangle function to limit
the impact of the time window. The distribution of waves is illustrated in Fig. 13 at ﬁxed pressures (10 and 100 Pa), and in Fig. 14 at
ﬁxed latitudes (equator and 60°N).
The most visible groups of waves in the u0 ﬁelds are locatedaround frequencies of 6 cycles/Titan day (thereafter f1) and11 cycles/Titan day (thereafter f2), though above 100 Pa, thewaves are more widely distributed in frequency. For v0 , the f1 wave
is the most visible feature. The amplitudes of u0 and v0 around these
two frequencies are shown in the meridional plane in Fig. 15. Most
activity is located in the 103–10 Pa pressure range, where the zonal
wind jet is present, and at latitudes from 30°S to the northern(winter) pole, with u0 peaking near the equator while v0 is peaking
in mid-latitudes. The latitudinal gradient of the absolute vorticity

Fig. 15. Latitude–altitude distribution of the maximum amplitude of the transients in the frequency ranges 5–7 cycles/Titan day (f1) and 10–12 cycles/Titan day (f2) at
Ls = 295°. (a and c) Amplitude on zonal wind, the solid contours are the latitudinal gradient of the absolute vorticity (in 1011 m1s1). (b and d) Amplitude on meridional
wind, the solid contours are the mean zonal wind ﬁeld (in m/s).

718

S. Lebonnois et al. / Icarus 218 (2012) 707–722





@q
@
@u
¼
2X sin / @y a@/a@/ð4Þ

(where q is the zonal average of q, the absolute vorticity) is used as acriterion for barotropic instabilities (Rayleigh-Kuo criterion, Kuo,
1949; Vallis, 2006): it must change sign in the region of instability.
Fig. 15 conﬁrms that the region where the waves are the mostactive is also a region where @q
is low and often changes sign. There@yfore it is probable that these waves are barotropic in nature. Fig. 15
recalls work by Luz and Hourdin (2003) and Luz et al. (2003) (e.g.
Fig. 15 of Luz et al. (2003)) in which the latitudinal gradient of potential vorticity was used to parameterize the eddy momentum
transport in the IPSL 2D-CM. This conﬁrms that instabilities in the
2D-CM and 3D GCM have similar distributions (and similar consequences in terms of angular momentum transport). The distributionof the f1 wave shows a signiﬁcant activity on the poleward side of
the jet, compared to f2.
To better investigate the role of each wave in the transport ofangular momentum in this simulation, the transients have been
separated into three frequency regions, using ﬁlters applied on
the spectra: waves with frequencies lower than 4 cycles/Titan day(f0), waves with frequencies between 4 and 8 cycles/Titan day (i.e.
mainly f1), and waves with larger frequencies (i.e. mainly f2). Instantaneous horizontal ﬁeld snapshots of ﬁltered u0 are plotted inFig. 16. The chosen pressures correspond to regions of high amplitudes for each of these waves: around 200 Pa for f2, around 103 Pa
for f1, and very close to the surface for f0 (1.3  105 Pa). The stratospheric waves have wavenumber 1, while the f0 tropospheric wave(a)has higher wavenumbers. All these waves are eastward propagating. The dominant f0 wave seen at southern mid-latitudes inFig. 16a has a very low-frequency, around 0.2 cycles/Titan day.
The stratospheric waves present very similar characteristics tothe waves displayed in Newman et al. (2011), which had frequencies from around 5 cycles/Titan day during spin-up to around16 cycles/Titan day at the end of their simulation.
In the regions where the f0 wave is active, the zonal wind u issteadily increasing with altitude. Here, close to the surface, thedynamical regime is more baroclinic. A criterion similar to the Rayleigh-Kuo criterion for barotropic instabilities may be used to see
whether baroclinic instabilities can grow: the latitudinal gradientof the absolute vorticity @q
must change sign in the vertical direc@ytion (Holton, 1992; Vallis, 2006). This is visible in the structures
present in the troposphere in Fig. 15, close to the surface and close
to the pole. This reinforces the hypothesis that the tropospheric f0
waves are baroclinic in nature. Note that baroclinic instabilities
were seen at high latitudes in simulations of superrotating
atmospheres done by Mitchell and Vallis (2010), for intermediate
thermal Rossby numbers. Near-surface waves have also recently
been analyzed in a tropospheric Titan GCM simulation (Mitchell
et al., 2011), in which they are connected to the precipitation patterns observed near equinox (Turtle et al., 2011). Mitchell et al.
(2011) obtained mainly two wave modes, equatorial Kelvin waves
and mid-latitude waves propagating westward, conﬁned near thesurface and presenting similarities with the waves we obtain here
(though the f0 waves presented in Fig. 16a are slowly propagating
eastward, not westward).

(b)

(c)

Fig. 16. Snapshots of the horizontal ﬁeld of ﬁltered perturbations on the zonal wind, at Ls = 295°. (a) u0 (f0) in the lower troposphere, at 1.3  105 Pa; (b) u0 (f1) in thestratosphere, at 103 Pa; (c) u0 (f2) in the stratosphere, at 200 Pa.

719

S. Lebonnois et al. / Icarus 218 (2012) 707–722

The transport of angular momentum by each family of waves iscomputed using the ﬁltered u0 and v0 , and presented in Fig. 17,
where it is split between upper atmosphere, stratosphere and troposphere. This transport is composed of two main components.
The stratospheric waves (both f1 and f2 equally) transport a significant amount of angular momentum in the direction opposing themean meridional circulation, especially during the solstice season.
This transport also appears at low- to mid-latitudes in Figs. 11 and12b. However, another component of the transport appears in
Fig. 17: this is done by the low-frequency waves in the troposphere, at mid-to high-latitudes, despite their very low amplitude.
This transport is directed towards the pole and is present even inperiods of low barotropic wave activity. This component explains
most of the transport seen poleward of 40° latitude in Figs. 11
and 12b.
To verify whether the high-latitude tropospheric component
was related to the solar thermal tide, we have run a 1-year simulation without the diurnal cycle. During this year, the transportof angular momentum is not much affected; the same two components are found in the transport by waves, and the circulation isvery similar to the nominal simulation with a diurnal cycle. This
result is consistent with the small impact of the diurnal cycle obtained in Hourdin et al. (1995) and conﬁrms that the high-latitude
tropospheric transport in our simulations is not related to the solar
thermal tide.

4. Results in the troposphere
The annual cycles of temperature at the surface and at the tropopause are shown in Fig. 18. The surface temperature is roughly
1 K below the value observed by Huygens (93.65 K, Fulchignoni
et al., 2005), and the latitudinal contrast between winter polar regions and equator (around 2 K), as well as between summer polar
regions and equator (around 0.6 K) are both under-estimated compared to the values retrieved by Jennings et al. (2009) (respectively
3 K and 2 K). The contrasts seen after northern spring equinox(1.5–2 K, almost symmetric) are consistent with IRIS/Voyager 1
observations. The modeled temperature at tropopause is very close
to the Huygens observed value (70.4 K, Fulchignoni et al., 2005)
and does not vary appreciably with either latitude or season. Whenthe surface thermal inertia is reduced to I = 340 J m2 s0.5 K1, thetemporal variations of both surface and tropopause temperaturesare slightly enhanced, and the temporal evolution of the surface
maximum temperature is modiﬁed, with the maximum at
Ls  300° displaced towards southern mid-latitudes.
The vertical proﬁle of potential temperature at the Huygens
landing site (Tokano et al., 2006, Fig. 2) shows an increase ofroughly 5 K in the last 10 km, with a well-mixed boundary layer
conﬁned to 300 m above the surface. Shoulders in the proﬁle arealso seen around 800 m and 3 km altitude. A 2–3 km boundary
layer is also consistent with dune spacing, as shown in Lorenz

(a)

(b)

(c)

(d)

Fig. 17. Latitudinal transport of angular momentum in the simulation at Ls = 295°, for mean meridional circulation (solid lines), low-frequency transients (0–4 cycles/Titan
day, dotted lines), transients with frequencies between 4 and 8 cycle/Titan day (i.e. mainly f1, dashed lines) and high-frequency transients (above 8 cycle/Titan day, i.e. mainly
f2, dash-dotted lines). In panel (a), the relative angular momentum latitudinal transport is integrated for the entire column at each latitude, and is divided by the total mass ofthe column. In the other panels, it is integrated only for a fraction of the column, but it is still divided by the total mass of the entire column, so that the three plots add to giveplot (a). Unit is 103 m3/s2, positive values are northward.

720

S. Lebonnois et al. / Icarus 218 (2012) 707–722

Fig. 19. Potential temperature vertical proﬁle at 10°S (zonal average) and Ls = 300°,
corresponding to the Huygens descent. The HASI/Huygens potential temperature
proﬁle is also shown (bold line).

is slightly positive (eastward), though it is mostly negative in the
deepest 5 km, as in the observations. However, the amplitude ofthis westward wind is approximately half of the observed value.
When reduced thermal inertia is used, the amplitude of wind variations are enhanced, and the retrograde region (negative u) persists up to 7–8 km.
The comparison between model and observations is less satisfactory for the meridional wind. Southward wind (negative) is indeed obtained just above the surface, but again the amplitude is
only half the observed value. Between 1 and 5 km, both observed
and modeled winds are northward, though modeled values aremuch less than observations. Above 3 km, the modeled meridional
wind is almost zero up to the tropopause. The northward andsouthward winds present in the observations above 5 km are notobtained in the simulation. However, this observed proﬁle corresponds to an instantaneous proﬁle where waves and topographic
effects may have been present.
Speciﬁc studies dedicated to the troposphere are underway andwill be presented in future papers.
5. Conclusion
Fig. 18. Annual variations of the zonally averaged latitudinal temperature proﬁle at(a) surface and (b) tropopause (roughly 40 km altitude).

et al. (2010). In Fig. 19, the diurnally averaged potential temperature proﬁle in the GCM is shown for the deepest 10 km, at the Huygens landing site and season, together with the Huygens proﬁle.
The simulated boundary layer is located in the deepest 2 km, withan increase of 7 K between the surface and 10 km altitude. A detailed analysis of the diurnal and seasonal evolution of this nearsurface potential temperature proﬁle appears in a dedicated study(Charnay and Lebonnois, 2011). The modeled surface temperature
at this point (and also its latitudinal proﬁle) is sensitive to the
choices of thermal inertia and albedo, and is also sensitive to thevertical distribution of haze. Clouds (not included here) may also
affect the surface temperature.
The proﬁles of diurnally averaged zonal and meridional winds
in the troposphere at the Huygens landing site have been retrieved
from the probe’s trajectory (Tokano, 2009, Fig. 2). These are shown
together with the modeled wind proﬁles at the same location andseason in Fig. 20. As seen in Fig. 9, the vertical proﬁle of the zonal
wind in the troposphere is in very good agreement with the measurements of the Huygens probe. Just above the surface, the wind

Based on the current version of the LMDZ Earth General
Circulation Model dynamical core, we have developed a new model of Titan’s atmosphere, including physic parametrizations fromthe IPSL 2-dimensional Titan Climate Model (2D-CM). However,
due to computational limitations, microphysics of the haze andphotochemistry are still computed as zonal averages. The new
GCM simulation presented in this work has been initialized withthe 2D-CM reference state at northern spring equinox and run
for 12 Titan years.
Comparing this new simulation to previous 2-dimensional
simulations and to observations, the following conclusions may
be drawn:
 The modeled temperature structure is fairly similar between
both models, with the same problems in the upper stratosphere
where the temperatures are too cold, and above the 10 hPa level(approximately 300 km altitude) where no stratopause is
obtained and the temperature is not realistic compared toobservations. As in the 2D-CM, this is presumably related to
the limitations of the radiative transfer and to the upper boundary position, which limits the altitude of the haze production
region.

721

S. Lebonnois et al. / Icarus 218 (2012) 707–722

(a)

(b)

Fig. 20. Zonal (a) and meridional (b) wind vertical proﬁles at 10°S (zonal average) and Ls = 300°, corresponding to the Huygens descent. Values computed from the Huygens
descent trajectory dataset (Tokano, 2009) are also plotted (bold lines).

 The comparison to observed latitudinal proﬁles of temperaturein the lower stratosphere shows that the overall modeled temperature structure appears compressed in altitude. In spring,the comparison is not as good as for the 2D-CM, probably
because of the haze distribution.
 The detached haze layer is present most of the year, lower thanobserved as in the 2D-CM. It shrinks and disappears after equinox, due to the shift of the ascending region from one pole tothe other. The altitude difference between Cassini and Voyager
1 observations is a seasonal effect, now conﬁrmed by the Cassini monitoring of the detached haze layer over the equinox(West et al., 2011). Comparing simulations with observationsof this transition period should help to clarify the interaction
between the haze and the atmospheric circulation.
 The zonal wind peak modeled at Cassini season is lower inamplitude as well as in altitude in this new simulation compared to the 2D-CM. However, the temporal evolution of thezonal wind ﬁeld is qualitatively consistent with available
observations.
 Changing the planetary boundary layer scheme has greatly
improved agreement between the model and the tropospheric
proﬁle of the zonal wind from the Huygens descent.
 As in previous LMD models, a minimum in the zonal wind ﬁeld
is obtained close to the region where Huygens observed a stronger minimum in the equatorial lower stratosphere. In the simulation, this minimum is correlated with a local minimum in thetemperature ﬁeld, due to the relative distributions of solar forcing and dynamical transport of energy. According to the thermal
wind equation, the observed wind minimum could correspondto an equatorial temperature around 2 K lower than mid-latitudes at 20–40 hPa (60–75 km altitude). This region is also correlated with a reversal (change of sign) of the meridional wind.
 In this 3D simulation, the diurnal cycle has been taken into
account, but its role is minor. Barotropic instabilities, diagnosed
and parameterized in the 2D-CM, are now explicitely simulated
and many waves are visible in the stratosphere, mostly in theregion around the zonal wind maximum. The horizontal transport of angular momentum by the mean meridional circulation
is similar to the 2D simulation. However, the transport bytransient waves is more complicated than in previous work.
Barotropic waves transport momentum opposite to mean
meridional circulation in mid- to low-latitudes, but a second
component of the momentum transport appears at mid- tohigh-latitudes, due to low-frequency waves located close to
the surface. Instability criterion as well as the more baroclinic
characteristics of this region indicate that these waves are baroclinic waves. They transport momentum poleward at high latitudes all year long.
 In the lower troposphere, simulated potential temperature,zonal and meridional winds are compared to Huygens descentproﬁles. Though the amplitudes of wind variations are smaller
than observed, the vertical proﬁles are very similar, except forthe meridional wind above 3 km. Near-surface behavior of themodel is sensitive to uncertain input parameters such as albedo
and thermal inertia.
Due to the complex coupling between dynamics, haze particles
and radiative transfer, it is difﬁcult to adjust model parameters to
get a circulation that ﬁts observations everywhere. The modelneeds to be improved, especially above 300 km altitude where
the temperature structure is far from observations. This affects
the top of the haze layer as well as its time variability, and theresulting zonal wind structure and temperature distributions.
The upgrade of the radiative transfer module in the GCM and thevertical extension of the upper boundary are underway and should
help improve this situation.
Acknowledgements
The authors thank Frederic Hourdin for very useful discussions
on the model and on this paper. This work has been supported by
the project Exoclimats ﬁnanced by the Agence Nationale de laRecherche (ANR), and by the computation facilities of both theInstitut du Développement et des Ressources en Informatique Scientiﬁque (IDRIS) and the University Pierre and Marie Curie(UPMC).
References
Achterberg, R.K., Conrath, B.J., Gierasch, P.J., Flasar, F.M., Nixon, C.A., 2008a.
Observation of a tilt of Titan’s middle-atmospheric superrotation. Icarus 197,
549–555.
Achterberg, R.K., Conrath, B.J., Gierasch, P.J., Flasar, F.M., Nixon, C.A., 2008b. Titan’s
middle-atmospheric temperatures and dynamics observed by the Cassini
Composite Infrared Spectrometer. Icarus 194, 263–277.
Bird, M.K. et al., 2005. The vertical proﬁle of winds on Titan. Nature 438, 1–3.
Bouchez, A., 2004. Seasonal Trends in Titan’s Atmosphere: Haze, Winds, and Clouds.
Ph.D. thesis. California Institute of Technology.

722

S. Lebonnois et al. / Icarus 218 (2012) 707–722

Charnay, B., Lebonnois, S., 2011. Thermal structure and dynamics of Titan’s lower
troposphere. Nat. Geosci., in press.
Coustenis, A., Bézard, B., 1995. Titan’s atmosphere from Voyager infraredobservations. IV. Latitudinal variations of temperature and composition.
Icarus 115, 126–140.
Coustenis, A. et al., 2003. Titan’s atmosphere from ISO mid-infrared spectroscopy.
Icarus 161, 383–403.
Coustenis, A. et al., 2007. The composition of Titan’s stratosphere from Cassini/CIRS
mid-infrared spectra. Icarus 189, 35–62.
Coustenis, A., Jennings, D., Nixon, C., Achterberg, R.K., Lavvas, P., Vinatier, S., Teanby,
N.A., Bjoraker, G.L., Carlson, R.C., Piani, L., Bampasidis, G., Flasar, F.M., Romani,
P.N., 2010. Titan trace gaseous composition from CIRS at the end of the Cassini–
Huygens prime mission. Icarus 207, 461–476.
Crespin, A., Lebonnois, S., Vinatier, S., Bézard, B., Coustenis, A., Teanby, N.A.,
Achterberg, R.K., Rannou, P., Hourdin, F., 2008. Diagnostics of Titan’s
stratospheric dynamics using Cassini/CIRS data and the IPSL General
Circulation Model. Icarus 197, 556–571.
Del Genio, A.D., Zhou, W., Eichler, T.P., 1993. Equatorial superrotation in a slowly
rotating GCM: Implications for Titan and Venus. Icarus 101, 1–17.
Flasar, F.M. et al., 2005. Titan’s atmospheric temperatures, winds, and composition.
Science 308, 975–978.
Flasar, F.M., Achterberg, R.K., 2009. The structure and dynamics of Titan’s middle
atmosphere. Philos. Trans. Roy. Soc. A 367, 649–664.
Flasar, F.M., Samuelson, R.E., Conrath, B.J., 1981. Titan’s atmosphere: Temperatureand dynamics. Nature 292, 693–698.
Folkner, W.M. et al., 2006. Winds on Titan from ground-based tracking of the
Huygens probe. J. Geophys. Res. 111, E07S02.
Friedson, A.J., West, R.A., Wilson, E.H., Oyafuso, F., Orton, G.S., 2009. A global climate
model of Titan’s atmosphere and surface. Planet. Space Sci. 57, 1931–1949.
Fulchignoni, M. et al., 2005. In situ measurements of the physical characteristics ofTitan’s environment. Nature 438, 1–7.
Holton, J.R., 1992. An introduction to dynamic meteorology. In: International
Geophysics Series, third ed. Academic Press, San Diego, New York.
Hourdin, F. et al., 2006. The LMDZ4 general circulation model: Climate performanceand sensitivity to parameterized physics with emphasis on tropical convection.
Clim. Dyn. 27, 787–813.
Hourdin, F., Le Van, P., Forget, F., Talagrand, O., 1993. Meteorological variability andthe annual surface pressure cycle on Mars. J. Atmos. Sci. 50, 3625–3640.
Hourdin, F., Talagrand, O., Sadourny, R., Courtin, R., Gautier, D., McKay, C.P., 1995.
Numerical simulation of the general circulation of the atmosphere of Titan.
Icarus 117, 358–374.
Hourdin, F., Couvreux, F., Menut, L., 2002. Parameterization of the dry convective
boundary layer based on a mass ﬂux representation of thermals. J. Atmos. Sci.
59, 1105–1123.
Hourdin, F., Lebonnois, S., Luz, D., Rannou, P., 2004. Titan’s stratosphericcomposition driven by condensation and dynamics. J. Geophys. Res. 109,
E12005.
Hubbard, W.B. et al., 1993. The occultation of 28 Sgr by Titan. Astron. Astrophys.
269, 541–563.
Jennings, D.E. et al., 2009. Titan’s surface brightness temperatures. Astrophys. J. 691,
L103–L105.
Kostiuk, T. et al., 2001. Direct measurements of winds on Titan. Geophys. Res. Lett.
28, 2361–2364.
Kostiuk, T. et al., 2005. Titan’s stratospheric zonal wind, temperature, and ethane
abundance a year prior to Huygens insertion. Geophys. Res. Lett. 32, L22205.
Kostiuk, T. et al., 2006. Stratospheric global winds on Titan at the time of Huygens
descent. J. Geophys. Res. 111, E07S03.
Kostiuk, T. et al., 2010. High spectral resolution infrared studies of Titan: Winds,
temperature, and composition. Planet. Space Sci. 58, 1715–1723.
Kuo, H.I., 1949. Dynamic instability of 2-dimensional nondivergent ﬂow in abarotropic atmosphere. J. Meteor. 6, 105–122.
Lebonnois, S., Toublanc, D., Hourdin, F., Rannou, P., 2001. Seasonal variations inTitan’s atmospheric composition. Icarus 152, 384–406.
Lebonnois, S., Hourdin, F., Rannou, P., Luz, D., Toublanc, D., 2003. Impact of theseasonal variations of ethane and acetylene distributions on the temperatureﬁeld of Titan’s stratosphere. Icarus 163, 164–174.
Lebonnois, S., Hourdin, F., Rannou, P., 2009. The coupling of winds, aerosols andphotochemistry in Titan’s atmosphere. Philos. Trans. Roy. Soc. A 367, 665–682.
Lebonnois, S., Hourdin, F., Eymet, V., Crespin, A., Fournier, R., Forget, F., 2010.
Superrotation of Venus’ atmosphere analysed with a full General Circulation
Model. J. Geophys. Res. 115, E06006.
Lebonnois, S. et al., 2011. A comparative analysis of simpliﬁed general circulation
models of Venus atmosphere. In: Towards Understanding the Climate of Venus:
Application of Terrestrial Models to Our Sister Planet. Springer, Netherlands, inpress.
Lebreton, J.P., Coustenis, A., Lunine, J., Raulin, F., Owen, T., Strobel, D., 2009. Resultsfrom the Huygens probe on Titan. Astron. Astrophys. Rev. 17, 149–179.
Lorenz, R.D., Claudin, P., Andreotti, B., Radebaugh, J., Tokano, T., 2010. A 3 km
atmospheric boundary layer on Titan indicated by dune spacing and Huygens
data. Icarus 205, 719–721.
Luz, D., Hourdin, F., 2003. Latitudinal transport by barotropic waves in Titan’s
stratosphere. I. General properties from a horizontal shallow-water model.
Icarus 166, 328–342.
Luz, D., Hourdin, F., Rannou, P., Lebonnois, S., 2003. Latitudinal transport bybarotropic waves in Titan’s stratosphere. II. Results from a coupled dynamics–
microphysics–photochemistry GCM. Icarus 166, 343–358.

Luz, D., Civeit, T.R.C., Lebreton, J.P., Gautier, D., Rannou, P., Kaufer, A., Witasse, O.,
Lara, L., Ferri, F., 2005. Characterization of zonal winds in the stratosphere of
Titan with UVES. Icarus 179, 497–510.
Marten, A., Hidayat, T., Biraud, Y., Moreno, R., 2002. New millimeter heterodyne
observations of Titan: Vertical distributions of nitriles HCN, HC3N, CH3CN, andthe isotopic ratio 15N/14N in its atmosphere. Icarus 158, 532–544.
McKay, C.P., Pollack, J.B., Courtin, R., 1989. The thermal structure of Titan’s
atmosphere. Icarus 80, 23–53.
Mellor, G.L., Yamada, T., 1982. Development of a turbulent closure model for
geophysical ﬂuid problems. Rev. Geophys. Space Phys. 20, 851–875.
Mingalev, I.V. et al., 2006. First simulation results of Titan’s atmospheric dynamics
with a global 3-D non-hydrostatic circulation model. Ann. Geophys. 24, 1–15.
Mitchell, J.L., Vallis, G.K., 2010. The transition to superrotation in terrestrial
atmospheres. J. Geophys. Res. 115, E12008.
Mitchell, J.L., Pierrehumbert, R.T., Frierson, D., Caballero, R., 2006. The dynamics
behind Titan’s methane cloud. Proc. Nat. Acad. Sci. 103, 18421–18426.
Mitchell, J.L., Ádámkovics, M., Caballero, R., Turtle, E.P., 2011. Locally enhanced
precipitation organized by planetary-scale waves on Titan. Nat. Geosci. 4, 589–
592.
Moreno, R., Marten, A., Hidayat, T., 2005. Interferometric measurements of zonal
winds on Titan. Astron. Astrophys. 437, 319–328.
Newman, C.E., Lee, C., Lian, Y., Richardson, M.I., Toigo, A.D., 2011. Stratosphericsuperrotation in the TitanWRF model. Icarus 213, 636–654.
Niemann, H.B. et al., 2005. The abundances of constituents of Titan’s atmospherefrom the GCMS instrument on the Huygens probe. Nature 438, 1–6.
Rannou, P., Hourdin, F., McKay, C.P., 2002. A wind origin for Titan’s haze structure.
Nature 418, 853–856.
Rannou, P., Hourdin, F., McKay, C.P., Luz, D., 2004. A coupled dynamicsmicrophysics model of Titan’s atmosphere. Icarus 170, 443–462.
Rannou, P., Lebonnois, S., Hourdin, F., Luz, D., 2005. Titan atmosphere database. Adv.
Space Res. 36, 2194–2198.
Rannou, P., Montmessin, F., Hourdin, F., Lebonnois, S., 2006. The latitudinal
distribution of clouds on Titan. Science 311, 201–205.
Richardson, M.I., Toigo, A.D., Newman, C.E., 2007. PlanetWRF: A general purpose,
local to global numerical model for planetary atmospheric and climate
dynamics. J. Geophys. Res. 112, E09001.
Schinder, P.J., Flasar, F.M., Marouf, E.A., French, R.G., McGhee, C.A., Kliore, A.J.,
Rappaport, N.J., Barbinis, E., Fleischman, D., Anabtawi, A., 2011. The structure of
Titan’s atmosphere from Cassini radio occultations. Icarus 215, 460–474.
Sicardy, B. et al., 2006. The two Titan stellar occultations of 14 November 2003. J.
Geophys. Res. 111, E11S91.
Teanby, N.A. et al., 2008a. Titan’s winter polar vortex structure revealed by chemical
tracers. J. Geophys. Res. 113, E12003.
Teanby, N.A., Irwin, P.G.J., de Kok, R., Nixon, C.A., Coustenis, A., Bézard, B., Calcutt,
S.B., Bowles, N.E., Flasar, F.M., Fletcher, L., Howett, C., Taylor, F.W., 2006.
Latitudinal variations of HCN, HC3N and C2N2 in Titan’s stratosphere derived
from Cassini CIRS data. Icarus 181, 243–255.
Teanby, N.A., Irwin, P.G.J., de Kok, R., Vinatier, S., Bézard, B., Nixon, C.A., Flasar, F.M.,
Calcutt, S.B., Bowles, N.E., Fletcher, L., Howett, C., Taylor, F.W., 2007. Vertical
proﬁles of HCN, HC3N and C2H2 in Titan’s atmosphere derived from Cassini/CIRS
data. Icarus 186, 364–384.
Teanby, N.A., Irwin, P.G.J., de Kok, R., Nixon, C.A., Coustenis, A., Royer, E., Calcutt, S.B.,
Bowles, N.E., Fletcher, L., Howett, C., Taylor, F.W., 2008b. Global and temporalvariations in hydrocarbons and nitriles in Titan’s stratosphere for northern
winter observed by Cassini/CIRS. Icarus 193, 595–611.
Teanby, N.A., Irwin, P.G.J., de Kok, R., Jolly, A., Bézard, B., Nixon, C.A., Calcutt, S.B.,
2009a. Titan’s stratospheric C2N2, C3H4 and C4H2 abundances from Cassini/CIRS
far-infrared spectra. Icarus 202, 620–631.
Teanby, N.A., Irwin, P.G.J., de Kok, R., Nixon, C.A., 2009b. Dynamical implications of
seasonal and spatial variations in Titan’s stratospheric composition. Philos.
Trans. Roy. Soc. A 367, 697–711.
Tokano, T., 2009. The dynamics of Titan’s troposphere. Philos. Trans. Roy. Soc. A 367,
633–648.
Tokano, T., Neubauer, F.M., Laube, M., McKay, C.P., 1999. Seasonal variation of
Titan’s atmospheric structure simulated by a general circulation model. Planet.
Space Sci. 47, 493–520.
Tokano, T., Ferri, F., Colombatti, G., Mäkinen, T., Fulchignoni, M., 2006. Titan’splanetary boundary layer structure at the Huygens landing site. J. Geophys. Res.
111, E08007.
Tomasko, M.G. et al., 2005. Rain, winds and haze during the Huygens probe’s
descent to Titan’s surface. Nature 438, 1–14.
Turtle, E.P. et al., 2011. Seasonal changes in Titan’s meteorology. Geophys. Res. Lett.
38, L03203.
Vallis, G.K., 2006. Atmospheric and Oceanic Fluid Dynamics. Cambridge University
Press, Cambridge, UK.
Vinatier, S., Bézard, B., Fouchet, T., Teanby, N.A., de Kok, R., Irwin, P., Conrath, B.J.,
Nixon, C.A., Romani, P.N., Flasar, F.M., Coustenis, A., 2007. Vertical abundance
proﬁles of hydrocarbons in Titan’s atmosphere at 15°S and 80°N retrieved from
Cassini/CIRS spectra. Icarus 188, 120–138.
Vinatier, S., Bézard, B., Nixon, C.A., Mamoutkine, A., Carlson, R.C., Jennings, E.D.,
Guandique, E.A., Teanby, N.A., Bjoraker, G.L., Flasar, F.M., Kunde, V.G., 2010.
Analysis of Cassini/CIRS spectra of Titan acquired during the nominal mission. I.
Hydrocarbons, nitriles and CO2 vertical mixing ratio proﬁles. Icarus 205, 559–
570.
West, R.A. et al., 2011. The evolution of Titan’s detached haze layer near equinox in2009. Geophys. Res. Lett. 380, L06204.