Temporal Evolution of Titan’s Stratospheric Temperatures and Trace Gases from a Two-dimensional Retrieval of Cassini Composite Infrared Spectrometer Data

We use a two-dimensional (2D) radiative transfer model of Titan, which allows the atmospheric structure to vary in both altitude and latitude, to retrieve the spatial distribution of temperature, haze extinction, and C2H2, C2H6, C3H8, CH3C2H, C4H2, and HCN gases, from Cassini Composite Infrared Spectrometer (CIRS) limb-mapping observations over the duration of the Cassini mission. We compare our results with previous analyses of CIRS limb observations using radiative models that only allow the atmosphere to vary in altitude. The temperature, haze, and gas composition retrieved with the 2D model mostly show the same broad spatial and temporal trends as previously published results from 1D models. However, there are some significant differences in the retrieved structure at the fall and winter poles poleward of 60°. Most noticeably, the HCN abundance in the depleted region near 65°N at 350 km in northern winter is stronger in the 2D retrievals than in previous 1D retrievals, and the 2D retrievals show very different structure from earlier 1D retrievals in the north polar C2H2 structure during early northern spring, with a strong depletion around 70°N at 0.02 mbar.


Introduction
The 13 yr long Cassini mission to the Saturn system allowed monitoring of Titan's atmosphere for nearly two full seasons, from northern midwinter (L S = 293°, where L S is the solar longitude relative to northern spring equinox) through the end of northern spring (L S = 93°). The Composite Infrared Spectrometer (CIRS) on board the Cassini orbiter performed regular thermal IR observations of Titan designed for monitoring of the spatial and seasonal variation of the temperature and composition of Titan's stratosphere and lower mesosphere from both nadir-sounding and limb-sounding data (Nixon et al. 2019). Early temperature determinations from CIRS data before the 2009 equinox showed that the highest temperatures were at the northern (winter) polar stratopause (Flasar et al. 2005;Teanby et al. 2007Teanby et al. , 2008aAchterberg et al. 2008;Vinatier et al. 2010a;Achterberg et al. 2011), attributed to adiabatic heating in the subsiding branch of the meridional circulation. Subsequent observations after equinox showed cooling of the northern polar stratopause and warming of the south polar stratopause (Teanby et al. 2012(Teanby et al. , 2017Vinatier et al. 2015Vinatier et al. , 2020Sharkey et al. 2021) consistent with the reversal of the meridional circulation. The post-equinox warming of the south polar stratopause was interrupted from 2012 to 2015 with a period of cooling, which has been explained by increased radiative cooling during midfall as the circulation transports radiatively active trace gases to lower altitudes (Teanby et al. 2017). These observations from CIRS are consistent with general circulation models (GCMs) of Titan, which predict that for most of the Titan year the mean meridional circulation in the stratosphere and mesosphere consists of a single pole-to-pole Hadley cell with rising motion at the summer pole and subsidence at the winter pole, and a short transition period around the equinoxes with two cells where the rising branch of the circulation moves through the low latitudes to the other pole Newman et al. 2011;Lebonnois et al. 2012;Lora et al. 2015;Lombardo & Lora 2023). A detailed description of Titan's stratospheric circulation can be found in the reviews by Flasar & Achterberg (2009) and Lebonnois et al. (2014).
Effects of the meridional circulation can also be seen in the distribution of trace hydrocarbons and nitriles whose distributions are monitorable by CIRS (C 2 H 2 , C 2 H 4 , C 2 H 6 , C 3 H 8 , CH 3 C 2 H, C 4 H 2 , C 6 H 6 , HCN, and HC 3 N). The trace gases are produced in the mesosphere up through the ionosphere through photochemistry and ion chemistry, and have sinks in the stratosphere through photochemistry or condensation, leading to vertical gradients in the volume mixing ratio that depend on the production and loss rates and the vertical mixing timescales (Wilson & Atreya 2004;Lavvas et al. 2008;Krasnopolsky 2009;Vuitton et al. 2019). Pre-equinox CIRS observations showed enhanced abundances of many trace gases at the north pole in both nadir (Coustenis et al. 2007(Coustenis et al. , 2010Teanby et al. 2008bTeanby et al. , 2009) and limb (Teanby et al. 2007;Vinatier et al. 2007;Teanby et al. 2008a;Vinatier et al. 2010a) observations. An unexpected feature of the limb observations was a local minimum in the vertical profiles of C 2 H 2 , CH 3 C 2 H, C 4 H 2 , HCN, and HC 3 N around 0.05 mbar near 60°N (Teanby et al. 2008a;Vinatier et al. 2010a), with sharp composition gradients on the poleward side, which Teanby et al. (2008a) dubbed the "mesospheric depleted zone," and attributed to inhibition of meridional mixing in the polar mesosphere by the winter polar vortex, analagous to terrestrial polar regions. Subsequent observations after equinox showed the north polar enhancements increasing through 2011 before decreasing, although the north polar enhancement persisted through late northern spring; enhancements at the south pole started after equinox at high altitudes (Bampasidis et al. 2012;Coustenis et al. 2013Coustenis et al. , 2018Coustenis et al. , 2020Vinatier et al. 2015Vinatier et al. , 2020Teanby et al. 2019;Mathé et al. 2020). As with the temperatures, the trace gas composition is consistent with a pole-to-pole meridional circulation reversing direction around equinox.
One aspect of Titan's atmosphere that is potentially important for modeling limb-viewing observations is related to its large vertical extent relative to Titan's solid body. As discussed by Achterberg et al. (2008), limb-viewing raypaths with tangent heights in the middle stratosphere can traverse a horizontal extent of up to 20°of great circle arc. While Titan's stratosphere shows little zonal variation, there are strong meridional gradients of temperature and trace gas composition at mid-to-high latitudes in the winter hemisphere, as discussed above. Limb observations of polar latitudes are of necessity taken with a near-equatorial sub-spacecraft latitude, and, as a consequence, the raypaths for limb-viewing observations at polar latitudes can traverse up to 15°-20°of latitude, as shown in Figure 1. As a result, a radiative transfer model that assumes that the atmosphere varies only in altitude may incorrectly calculate the contribution of emission from higher altitudes to observations in the lower atmosphere if there are significant latitudinal variations in temperature or opacity, resulting in errors in retrieved atmospheric structure.
In this paper, we use a two-dimensional (2D) radiative transfer model of Titan, which allows the atmosphere to vary in both altitude and latitude, to examine the effect of the latitudinal variations on previous published retrievals of temperature, haze, and gas abundance from CIRS limbmapping observations, which assumed a spherically symmetric atmosphere that varies only with altitude (a one-dimensional, or 1D, model). Section 2 describes the CIRS data to be used. The forward radiative transfer and retrieval models used are discussed in Section 3. Section 4 compares the results from a 2D and 1D retrieval from a limb observation during the T4 flyby. Section 5 presents the results of 2D retrievals from several pairs of observations with nearly complete latitude coverage, and discusses the differences from published retrievals that used a 1D model. The implications of the differences between the 2D and 1D retrievals are discussed in Section 6.

Observations
We use data fromCassini CIRS, a Fourier transform spectrometer covering a spectral range from approximately 10 to 1500 cm −1 , with a variable spectral resolution between 0.5 and 15.5 cm −1 . The spectral coverage was split into three ranges with separate focal planes. Focal plane 1 (FP1), covering the far-IR between roughly 10 and 600 cm −1 (16.7-1000 μm), was a single pixel with a circular field of view (FOV) of nominally 3.9 mrad radius and a Gaussian spatial response with a FWHM of 2.5 mrad. Focal planes 3 and 4 (FP3/FP4) cover the mid-IR between roughly 600 and 1100 cm −1 (9-16.6 μm) and 1150-1500 cm −1 (6.7-9 μm), respectively. FP3/FP4 were 1-by-10 arrays, with a pixel size of 0.273 mrad and approximately uniform spatial response over each pixel ), aligned parallel to each other and separated by 0.94 mrad. A detailed description of CIRS and its operation and calibration is given by Jennings et al. (2017).
The data we use is a set of limb-mapping observations, known as MIRLMBMAPs, generally taken between 5 and 9 hr of closest approach to Titan when the projected size of the FP3/4 FOVs on Titan's limb was approximately 40 km, comparable to the pressure scale height. In these observations, the FP3/4 arrays were placed perpendicular to Titan's limb at two altitude positions, covering between roughly 100 and 600 km altitude, at 5°latitude intervals over one quarter of the visible limb, integrating for about 4 minutes at each position. Full details of the observation design and a table of the observations are given in Nixon et al. (2019). The resulting latitude and time coverage of all MIRLMPMAP observations is shown in Figure 2. Due to viewing geometry constraints, limb observations of the polar latitudes were possible only when the sub-spacecraft latitude was near the equator, limiting polar limb observations to those times when Cassini was in a lowinclination orbit about Saturn. Consequently, there were no near-polar limb observations in late northern winter or northern midspring (southern midfall). Because these observations were intended primarily for temperature mapping, the data were taken at the lowest spectral resolution of 15.5 cm −1 to maximize the signal-to-noise ratio. Despite the low spectral resolution, these observations have also been used to retrieve trace gas abundances (Teanby et al. 2007(Teanby et al. , 2017Vinatier et al. 2015Vinatier et al. , 2020, although, as discussed later, the blending of emission bands may potentially cause errors in the retrieved abundances that are correlated between the overlapping gases. The data used were taken from the "local calibration" version (see Jennings et al. 2017 for details), which uses deepspace calibration spectra taken as close in time to the observations as possible. Spectra taken at spacecraft event times in which any of the FP4 spectra have a radiance at 1100 cm −1 less than -´----2. 10 W cm str cm 9 2 1 1 (approximately 4σ below 0) were discarded; this eliminates spectra whose calibration was affected by jitter in the scan mechanism. The jitter leads to misidentification of the zero path difference point in the interferogram, resulting in a nonphysical slope of the spectral continuum and negative radiances at the lowwavelength end of FP4. The remaining spectra were averaged for each detector pixel at each targeted observing position, giving typically 20-25 individual spectra in each average.

Retrieval Model
The 2D radiative transfer model used is based on the model from Achterberg et al. (2008Achterberg et al. ( , 2011 for retrieval of temperatures from CIRS limb data, updated to allow retrieval of gas abundances as well as temperature and haze extinction. The model allows the temperature, haze, and gas abundances of Titan's atmosphere to vary in both altitude and latitude, using 200 altitude layers spaced equally between the surface and 600 km, and 37 latitude points equally spaced between the poles, giving 5°latitude spacing, comparable to the latitude spacing of the observations. The atmosphere is assumed to be zonally uniform, which is consistent with CIRS nadir observations of the middle stratosphere (Sharkey et al. 2020).
The retrieval algorithm is an updated version of the constrained linear inversion algorithm used by Achterberg et al. (2008), generalized to allow for an arbitrary set of retrieved quantities selected from temperature, haze abundance, gas abundance, and vertical pointing offset of each observation pointing. The haze abundance is represented by the natural logarithm of the volume extinction coefficient at 100 cm −1 ; the trace gas abundances are represented by the natural logarithm of the volume mixing ratio. Assume we have N j observation positions, with radiance measurements I j (ν i ) at N i wavenumbers ν i for each observation, yielding an observation vector of size We wish to retrieve N k quantities, each represented by a state vector x (k) . The forward model is linearized around a set of reference profiles ( ) x k 0 : where I 0 is the modeled radiance for the reference profiles, the modeled radiance with respect to the model parameters: The retrieved model parameters can include any of the temperature, gas volume mixing ratio, haze volume extinction coefficient, or a vertical offset to the observation pointing.
We find a solution for Δx by minimizing the quadratic form: Here, E is the measurement covariance matrix, S (k) are constraint matrices that act to smooth the solutions by enforcing correlation between elements of Δx (k) , and α (k) determines the strength of the constraint. The left-hand term in Equation (3) is the usual function minimized in least-squares fitting. The additional constraint term is needed because the retrieval problem is mathematically ill-posed in that there are a wide range of solutions which fit the data to within the measurement uncertainty. The constraint term acts to filter out contributions to the solutions from short spatial scales, which are poorly constrained by the data. The constraint matrices S (k) assume a Gaussian distribution for the spatial correlation of the atmospheric parameters. We use an e-folding scale of 60 km (roughly 1.5 pressure scale heights in the stratosphere) in the vertical and 5°in latitude (corresponding to the observation spacing). Minimization of Equation (3) gives us the solution As the forward model is nonlinear, Equation (4) is applied iteratively, until the rms value of ΔI decreases by less than 1% per iteration. The sensitivity of the retrieval to instrument noise is given by Equation (6) is an estimate of the amount by which the retrieved solution would change if the measure radiance changed by 1σ, and as such is low both where the retrieval is well constrained by the data and where the retrieval is not sensitive to the data and reverts to the initial guess. The local maxima of ( ) s k 2 in altitude or latitude provides a rough estimate of the altitude and latitude range over which the retrievals are constrained by the data.
Our retrievals proceed in two parts. In the first part, temperature, haze abundance, and a vertical pointing correction for each integration position of the detector arrays are retrieved from FP4 limb spectra, using the spectral ranges from 1100 to 1140 cm −1 and 1210 to 1315 cm −1 , where the absorption is from the haze continuum and the ν 4 band of CH 4, respectively. The initial guess temperatures are taken from retrievals from nadir viewing mapping observations, averaged over five different time periods: flybys T0 to T35 (2004 January to 2007 August), T37 to T68 (2007November to 2010, T74 to T95 (2011 February to 2013 October), T97 to T114 (2014 January to 2015 November), and T117 to T126 (2016 February to 2017 April). These nadir retrievals were performed on individual spectra from MIDRTMAP observations (Nixon et al. 2019, their Section 7.1 and Table 9) using the retrieval method described in Achterberg et al. (2008). The initial guess haze extinction uses a vertical profile retrieved from CIRS FP1 data at 15°S by Anderson & Samuelson (2011), assumed uniform in latitude. The vertical profile of CH 4 is taken from Huygens GCMS data (Niemann et al. 2010), also assumed uniform in latitude. We then use FP3 limb spectra, from 600 to 850 cm −1 , to simultaneously retrieve the abundances of C 2 H 6 , C 2 H 2 , C 3 H 8 , CH 3 C 2 H, C 4 H 2 , C 6 H 6 , HCN, HC 3 N, and CO 2 , the haze extinction, and the vertical pointing correction of each observation location. The temperature, initial guess haze extinction, and initial guess pointing correction are taken from the FP4 retrievals. The initial guess trace gas abundances are taken from the CIRS nadir gas retrievals of Teanby et al. (2019), assumed uniform in altitude with saturation vapor pressure curves applied. It should be noted that, as discussed by Gille & House (1971), the effect of the pointing correction is to minimize the error in the pressure of the observation tangent points, not the absolute altitude. As such, the pointing correction adjusts not only for errors in the pointing information but also for errors in the calculation of the vertical pressure profiles, due to errors in the assumed surface pressure or in the temperature profile below the altitude for which the temperatures can be reliably retrieved. The retrieved pointing offsets thus cannot be interpreted strictly as corrections to the pointing information from the data set.
A potential caveat on these retrievals is that, except for the C 2 H 6 band at 822 cm −1 , the emission peaks for the gases are blended together at the 15.5 cm −1 resolution of the observations used. The primary emission features of HCN (713 cm −1 ) and C 3 H 8 (748 cm −1 ) are blended with the broad wing of the C 2 H 2 feature (729 cm −1 ); the emission lines of HC 3 N (663 cm −1 ), CO 2 (667 cm −1 ), and C 6 H 6 (673 cm −1 ) are blended into a single peak overlaying the far wing of the C 2 H 2 emission; and the C 4 H 2 (628 cm −1 ) and CH 3 C 2 H (633 cm −1 ) features are mixed. The retrievals show nonphysically large vertical abundance variations of CO 2 , which is expected to be nearly uniform due to its long photochemical timescale (Vinatier et al. 2015;Mathé et al. 2020); since the CO 2 emission feature is blended with those of C 6 H 6 and HC 3 N, we do not present results for those gases as they are unreliable.

Forward Model at 80°N
To examine the effect of including the meridional variation of temperature and composition upon the modeled limb radiances at polar latitudes, we calculated forward models using three assumptions: first, using the full 2D model with temperature and composition varying with latitude; second, a 1D model with the temperature and composition held uniform in latitude, with the values from the latitude of the tangent point of the modeled observations; and, finally, a hybrid model where the temperature and haze are allowed to vary in latitude, but the gas abundances are constant in latitude at the values for the tangent point of the observations. The model atmosphere uses temperature and haze extinctions retrieved from the MIRLMBMAP observation taken during the T4 flyby on 2005 April 1. The gas abundances use the latitude dependence from Coustenis et al. (2010), assumed uniform in altitude with mixing ratios constrained to remain below saturation.
The modeled spectra for a limb observation at 80°N at 150 km and 400 km altitude are shown in Figure 3 for FP3 and Figure 4 for FP4, at both 0.5 cm −1 and 15 cm −1 spectral resolutions, along with plots of the difference between the 1D and hybrid models and the full 2D model. For this model of a winter polar observation, comparison of the 1D and hybrid models at 150 km tangent height shows that ignoring the latitude variation of the temperature and haze extinction results in an overestimate of the level of the haze continuum by about 15%-20%, and an overestimate of the emission in the ν 4 methane band, caused by the 1D model assuming stratopause temperatures that are too warm. The difference between the 1D and 2D forward models becomes weaker with increasing altitude as the latitude range traversed by the raypaths decreases, becoming less than the noise level between 300 and 400 km. The effect of ignoring the latitude variation in abundance of the trace gases is unsurprisingly dependent on the strength of the latitude variation. For C 2 H 6 , C 2 H 2 , C 3 H 8 , and HCN, which have only a weak dependence on latitude, the errors in the 1D model are mainly caused by the latitude variation of temperature and haze, with the errors due to the latitude variation of gas abundance on the order of the noise level of the observations. However, for the gases with strong latitude variation, the errors from ignoring the latitude variation of the gas abundance are comparable to those from ignoring the latitude variation of temperature and haze. In particular, at the 15.5 cm −1 resolution of the MIRLMBMAPs we are using, ignoring the latitude dependence of the atmosphere produces errors in the spectral features due to C 4 H 2 , CH 3 C 2 H, HC 3 N, and C 6 H 6 on the order of 20% or larger and significantly greater than the noise level of the observations.

Retrievals from T4 Limb Map
To look at the effect of accounting for the latitudinal variation of the atmosphere on retrievals, we performed retrievals on the T4 MIRLMBMAP (2005 April 1) using both our 2D retrieval model and a 1D model. The 1D model used is similar to the 2D model with two differences: the atmosphere is assumed to be spherically symmetric and vary only in altitude, and log-pressure is used as the vertical coordinate. The retrieval algorithm is otherwise the same, and the same opacity sources are used. A subset of the fits to the 85°N spectra for both retrievals is shown in Figure 5. The 1D and 2D retrievals show generally similar fits to within the noise level, with the exception of the lowest tangent heights in the FP4 data where the 1D fit is slightly too warm. The resulting retrieved cross sections are shown in Figures 6 and 7. The plots are masked out at latitudes and pressures where the retrieval loses sensitivity to the data as estimated by Equation (6). In addition, the haze and gas retrievals are masked where the retrieved temperature is not sensitive to the data, and the gas retrievals are also masked where the retrieved haze abundance loses sensitivity to the data. . Top: modeled Titan spectra between 1100 and 1400 cm −1 for a limb-viewing geometry with a tangent point at 80°N latitude, 150 km altitude, at spectral resolutions of 0.5 cm −1 (left) and 15.5 cm −1 (right), using a full 2D radiative model with latitude variations of all atmospheric parameters (black lines), and a 1D radiative model with atmospheric parameters for 80°N (red lines). Second row: difference between the 2D forward model and the 1D model (red lines) at 0.5 cm −1 (left) and 15.5 cm −1 (right). The black dashed line is the estimated noise equivalent radiance for an average of 25 CIRS spectra. Third row: as top row, for 400 km altitude. Bottom: as second row, for 400 km altitude.
Overall, the results from the 2D and the 1D retrievals are qualitatively similar, but with some moderate differences in detail. Most importantly, the 1D retrievals are significantly colder in the mid-to-low stratosphere in the winter pole, with the retrieved temperatures up to 20 K colder, as a consequence of the 1D model overestimating the contribution to the thermal emission from the stratopause level due to not accounting for the latitude variations of temperature. The colder temperatures also cause the retrieved temperatures to lose sensitivity to the data at higher altitudes (lower pressures); the lower temperatures push the contribution functions to higher altitudes due to the extreme nonlinearity of the Planck function in the methane ν 4 band used for the temperature retrievals.
As a consequence of the lower retrieved temperatures in the winter polar mid-stratosphere, the 1D model also produces somewhat larger retrieved abundances at altitudes below the ∼0.5 mbar pressure level poleward of ∼60°N in most gases. This is most noticeable in the C 2 H 6 abundance, which is roughly uniform in both latitude and altitude in the 2D retrieval apart from a weak decrease at pressures greater than ∼2 mbar, while the 1D retrievals show an increase in abundance at pressures greater than ∼0.5 mbar with a maximum centered near 1 mbar north of 50°N.
There are other observable difference between the 1D and 2D retrievals near the polar stratopause. The most prominent difference is in the HCN abundance. The 2D retrieval shows a strong local minimum in the HCN mixing ratio near 60°N and 0.01 mbar and extends more weakly southward, which is not seen in the 1D retrievals. A similar but weaker local minimum is seen in the 2D retrievals for C 4 H 2 . Although not seen in our 1D retrievals, a local minimum in HCN and C 4 H 2 at 60°N and 0.01 mbar, but weaker than in our 2D retrievals, was seen in 1D retrievals from 0.5 cm −1 resolution CIRS limb data by Teanby et al. (2008a, their Figure 5). We note, however, that Vinatier et al. (2010a) also found similar local minima in C 2 H 6 and C 2 H 2 that are not seen in our 2D retrievals.
We note a couple of caveats in the interpretation of our gas retrievals. First, the retrieved C 2 H 2 abundances show local maxima in pressure near 0.02 mbar and 0.5 mbar that are anticorrelated with structure in the retrieved HCN and C 3 H 8 abundances. Since the HCN and C 3 H 8 emission features are blended with the C 2 H 2 emission at the low spectral resolution of the data used, it is possible that some of the observed structure is not real but is caused by the data lacking the necessary information to cleanly distinguish between variations in C 2 H 2 , HCN, and C 3 H 8 . Second, the HCN retrievals show a local maximum near 1 mbar at the north pole, just above the altitude at which condensation occurs. Since our model enforces that abundances remain below the saturation vapor pressure, errors in the temperature profile or the saturation vapor pressure curve may cause the abundance below the condensation level to be too low, which is compensated for by increasing the abundance just above the condensation level to fit the observations. Both of these effects are sometimes seen in the retrievals from synthetic data used to tune the retrieval parameters. Additionally, it should be kept in mind when comparing our results to previous 1D retrievals that many of the earlier retrievals used earlier versions of the CIRS data set, and thus that some of the differences in the results may be caused by changes in the calibration. In particular, the haze and temperature retrievals are sensitive to the absolute calibration Figure 6. Retrievals of temperature (top row), haze volume extinction at 100 cm −1 (second row), C 2 H 6 volume mixing ratio (third row), and C 2 H 2 volume mixing ratio (bottom row) from the T4 MIRLMBMAP observation on 2005 April 1. The left column is retrievals using a 2D retrieval model, the right column is from a 1D retrieval model. of the continuum level, and differences in the temperature retrieval propagate into the gas abundance retrievals.

Results
Meridional cross sections of our retrieved temperature, haze, and gas abundances for C 2 H 6 , C 2 H 2 , HCN, C 4 H 2 , CH 3 C 2 H, and C 3 H 8 are shown in Figures 8 through 12 for eight pairs of observations, each of which provides nearly complete latitude coverage. The observations used to construct these cross sections are listed in Table 1. There is one cross section for mid-northern winter (L S = 310°), one for late northern winter (L S ≈ 340°), three for early southern fall (L S ≈ 10°-30°), and three for late southern fall (L S ≈ 65°-75°).  HCN (top row), C 4 H 2 volume mixing ratio (second row), CH 3 C 2 H volume mixing ratio (third row), and C 3 H 8 volume mixing ratio (bottom row) from the T4 MIRLMBMAP observation on 2005 April 1. The left column is retrievals using a 2D retrieval model, the right column is from a 1D retrieval model.

Temperature
Figure 8 shows our retrieved cross sections of temperature profiles. The temperatures in mid-to-late northern winter are consistent with the retrievals from Achterberg et al. (2008Achterberg et al. ( , 2011 using both limb and nadir data. The warmest temperatures, just below 210 K, are seen at the winter polar stratopause near 0.01 mbar, nearly a decade in pressure lower and ∼20 K warmer than the stratopause at lower latitudes. Postequinox, our retrievals show the same general temperature structure as the retrievals from the same data by Vinatier et al. (2015Vinatier et al. ( , 2020 and Teanby et al. (2017). In the period just after equinox, the north polar stratopause cools by around 20 K and becomes lower in altitude, while the low-latitude stratopause becomes slightly colder and drops in altitude; by 2012 there is a weak double stratopause. The south polar stratopause becomes elevated with only a decrease in temperature between 2010 and 2012, and south polar temperatures in the lower stratosphere become slightly colder. In the later parts of northern spring, the north polar stratopause lowers to 0.1 mbar, the same altitude as the stratopause at lower latitudes, while remaining about 10 K warmer than the low-latitude stratopause. At south polar latitudes, the stratopause becomes warmer than earlier, while temperatures in the middle and lower stratosphere become up to 40 K colder than 3 yr earlier, consistent with nadir observations (Teanby et al. 2017).
The observed temperatures at the polar stratopause are broadly consistent with the variations in the polar stratopause being caused by adiabatic heating (or cooling) from the meridional circulation. GCM models of Titan predict a global meridional circulation in the stratosphere and mesosphere during much of the Titan year, with rising motion at the summer pole and subsidence at the winter pole producing heating of the winter stratopause Newman et al. 2011;Lebonnois et al. 2012;Lora et al. 2015), although the models do not produce the elevated stratopause. The circulation then reverses near the equinox, with a period during which there are two cells with rising at low latitudes and subsidence at both poles. However, the behavior at the south polar stratopause after equinox cannot be explained solely by the modeled circulation, which does not explain the cooling seen between 2011 and 2012. Teanby et al. (2017) has suggested the observed cooling of the south polar mesosphere is caused by enhanced radiative cooling at the pole due to the observed increase in radiatively active trace gases. Figure 9 shows the retrieved abundances of C 2 H 6 . From 2006 to 2012, C 2 H 6 shows only a weak latitudinal variation, with roughly a factor of 2 higher abundance at north polar latitudes than southern, at abundances consistent with CIRS nadir retrievals (Coustenis et al. 2010;Teanby et al. 2019). The abundance is roughly constant with altitude at pressures between 0.01 and 1 mbar, except at high southern latitudes in southern fall, at abundances within a factor of 2 (and usually closer than that) of the nadir retrievals by Teanby et al. (2019). During the later part of southern summer, the meridional gradient reverses and there is a maximum in abundance at the south pole. A stronger vertical gradient also develops in the southern hemisphere in late fall. Our retrievals show the same general pattern as the limb retrievals by Vinatier et al. (2010aVinatier et al. ( , 2015Vinatier et al. ( , 2020, although our retrieved abundances show less vertical variation. In particular, we do not see the local minimum in C 2 H 6 near 70°N and 0.02 mbar noted by Vinatier et al. (2010a), and we find weaker vertical gradients in the period just after equinox.

C 2 H 2 , HCN, and C 3 H 8
Figures 10-12 show retrieved cross sections of C 2 H 2 , HCN, and C 3 H 8 , respectively. As discussed in Section 4.2, there are noticeable correlations between the results of these molecules, which may indicate that the low spectral resolution of the data may not allow completely independent retrieval of all three gases such that true abundance variations may lead to spurious apparent variations in the others.
C 2 H 2 shows a generally weak meridional variation in the middle stratosphere, except for a roughly factor of 2 enhancement northward of about 60°N that persists through 2015, and a similar enhancement in the south in 2015-2016. Outside the polar regions, the seasonal variation is minimal. The vertical structure of C 2 H 2 at most latitudes shows a broad, weak maximum between about 0.1 and 1 mbar, a minimum around 0.5 to 0.1 mbar, and a higher maximum near 0.02 mbar. In 2015 and 2016 there is also a rapid drop above the ∼0.01 mbar level except within about 30°of the south pole, also seen in the same data by Vinatier et al. (2020), which does not show clearly in Figure 10 as it occurs above the altitude at which the temperature retrievals become unreliable. The broad, weak maximum at 0.1 to 1 mbar is also seen in the retrievals of low-spectral-resolution limb data by Vinatier et al. (2015Vinatier et al. ( , 2020, as well as the high-resolution data (Vinatier et al. 2010a;Mathé et al. 2020). We also see a strong localized minimum around 0.05 mbar and 60°N-70°N in the data from 2010 to 2012; the high-altitude north polar C 2 H 2 profiles in these early northern spring data are very different from the retrievals by Vinatier et al. (2015). However, at pressures less than about 0.2 mbar, the variations in the C 2 H 2 become roughly anticorrelated with the variations in HCN and C 3 H 8 , in contrast to limb retrievals at higher spectral resolution where the vertical variations in C 2 H 2 and HCN show similar structure (Vinatier et al. 2015;Mathé et al. 2020), suggesting that our retrievals may not be cleanly separating the effects of C 2 H 2 and HCN on the spectra at higher altitudes. HCN shows a larger meridional gradient with abundances near the north pole roughly an order of magnitude larger than the south pole through 2012. The north polar enhancement increases through 2011, and disappears at pressures less than about 0.05 mbar by 2015 while weakening at higher pressures, consistent with retrievals from nadir observations (Bampasidis et al. 2012;Coustenis et al. 2020). HCN shows a strong enhancement at high altitudes at southern polar latitudes in 2015-2016. At pressures greater than ∼0.1 mbar, the south polar temperatures at this time are low enough that HCN condenses out. Our retrieved structure at high altitudes near the winter pole is very different from previous retrievals from 1D models (Teanby et al. 2007;Vinatier et al. 2010aVinatier et al. , 2015Vinatier et al. , 2020, similar to the differences in our 2D and 1D retrievals of T4 data in Section 4.2-we see a strong depletion of HCN near 0.01 mbar and 60°latitude that is considerably stronger in the pre-equinox than seen by Teanby et al. (2008a) and Vinatier et al. (2010a), and persists through 2016 despite not being seen in the retrievals from the same 2010-2012 maps by Vinatier et al. (2015), nor in the near-equinox high northern latitude retrievals from high-resolution data by Mathé et al. (2020). A prominent feature in the HCN cross sections is a maximum at around 0.1 mbar at nearly all latitudes in all maps, with a sharp drop in the HCN abundance at lower pressures, except at the winter pole. A similar maximum in the HCN vertical profiles is seen in retrievals of the T4 and T6 MIRLMBMAPs by Teanby et al. (2007), in MIRLMBMAP retrievals by Vinatier et al. (2015Vinatier et al. ( , 2020, and in retrievals from 0.5 cm −1 limb observations pre-equinox by Vinatier et al. (2010a) and throughout the Cassini mission by Mathé et al. (2020), albeit typically about one scale height lower in the atmosphere.
C 3 H 8 shows generally similar trends to C 2 H 2 , with a roughly factor of 2 enhancement within 30°of the winter pole and little seasonal variation outside the poles, but with a slightly larger vertical gradient with roughly factor of 2 increase in altitude over the middle and upper stratosphere in the southern hemisphere. Between ±60°our results show the same trends seen in the 1D limb retrievals from 0.5 cm −1 CIRS data by Mathé et al. (2020), although they did not see the winter polar enhancements. On the other hand, the south polar enhancement in 2015-2016 is also seen in the 1D retrievals from the same MIRLMBMAPs by Vinatier et al. (2020); however, their C 3 H 8 retrievals also show a localized enhancement northward of ∼60°N at ∼0.5 mbar that is weaker to nonexistent in our retrievals. Figures 13 and 14 show retrieved cross sections of C 4 H 2 and CH 3 C 2 H, respectively. Both the meridional and vertical gradients are stronger than for the gases with longer chemical lifetimes (C 2 H 2 , C 2 H 6 , C 3 H 8 , and HCN) with C 4 H 2 showing a roughly order-of-magnitude increase in abundance at the north pole through 2012 and at the south pole in 2015-2016. The polar enhancement in CH 3 C 2 H is somewhat weaker, roughly a factor of 3. Both gases show an order-of-magnitude increase in abundance at high altitudes (p < ∼0.2 mbar) poleward of 60°N in the period after equinox, and the development of a south polar enhancement in 2015, consistent with earlier analysis of CIRS nadir (Bampasidis et al. 2012;Coustenis et al. 2013Coustenis et al. , 2020Teanby et al. 2019) and limb (Vinatier et al. 2010a(Vinatier et al. , 2015Mathé et al. 2020;Vinatier et al. 2020) data. Figure 15 shows the retrieved cross sections of haze extinction from FP4 spectra. Outside the polar regions, the haze extinction is roughly constant in latitude within the uncertainties, with a slightly lower extinction in the southern hemisphere prior to equinox. There is a roughly factor of 2 enhancement poleward of 60°N after equinox at pressures less than 0.5 mbar; a similar enhancement at the south pole is visible in the late 2015 map. Haze retrievals from FP3 show the same general trends as the FP4 haze retrievals, but are considerably noisier. Our haze retrievals are generally consistent with previous 1D retrievals of haze from 0.5 cm −1 data (Vinatier et al. 2010b) and of retrievals from the 15.5 cm −1 limb maps (Vinatier et al. 2015, although the north polar to the rotation axis (Flasar et al. 2005). To estimate the zonal winds, the temperatures were first smoothed meridionally using smoothing B-splines. Equation (7) was then integrated along cylinders parallel to Titan's rotation axis assuming a boundary condition of winds in solid body rotation at 100 km at 4 times Titan's rotation rate, consistent with the wind measurements from Huygens Doppler Wind Experiment (Folkner et al. 2006).

Haze
The resulting winds are shown in Figure 16. Due to the necessary smoothing of the temperatures and the poorly constrained boundary condition, these winds should be taken only as a rough indication of the actual wind fields. Nonetheless, the results are broadly consistent with the winds estimates of Vinatier et al. (2020) from CIRS limb data and Sharkey et al. (2021) from CIRS nadir data, showing the north polar vortex, with wind speeds of ∼200 m s −1 persisting until late 2010, and a south polar vortex with similar wind speeds forming sometime between 2012 and 2015; both the limb and nadir temperature data indicate that the south polar vortex forms at high latitudes in 2013 and moves to midlatitudes by 2016.

Conclusions
Overall, the composition of Titan's atmosphere retrieved from CIRS limb data when using a 2D radiative model, accounting for the latitudinal variation of the atmosphere, follows the same broad spatial and temporal trends as previous 1D retrievals, albeit with some nontrivial differences at the fall and winter poles where there are significant meridional gradients in temperature and composition poleward of about 60°. The primary differences between our 2D and previous 1D retrievals are most noticeable in the HCN retrievals at higher Figure 16. Cross sections of zonal mean winds, calculated from the temperature data of Figure 8 using the gradient wind approximation, assuming a lower boundary condition assuming winds at 100 km in solid body rotation at 4 times Titan's rotation rate. altitudes near the north pole, and in the north polar C 2 H 2 retrievals in early northern fall. The "mesospheric depleted zone" near 65°N at 350 km noted by Teanby et al. (2008a) in C 2 H 2 , HCN, C 4 H 2 , and HC 3 N and seen by Vinatier et al. (2010a) in C 2 H 6 as well, is much stronger in HCN in our 2D retrievals than in previous 1D retrievals, and persists through the latest MIRLMBMAP in early 2016. We see a similar strong depletion of HCN at the south pole in 2016, which is not seen in the retrievals by Vinatier et al. (2020) or Mathé et al. (2020). Our 2D retrievals of C 2 H 2 at the north pole in early northern spring also show very different structure from Vinatier et al. (2015) and Mathé et al. (2020), with a strong localized depletion centered around 70°N near 0.02 mbar.
Despite the differences in detail, the broad trends in our retrievals are similar to previous results from 1D models, and our results are consistent with the generally accepted picture of the stratospheric and mesospheric circulation. Our temperature retrievals show the same general trends as earlier limb retrievals (Teanby et al. 2008a;Achterberg et al. 2008Achterberg et al. , 2011Vinatier et al. 2010aVinatier et al. , 2015Vinatier et al. , 2020Teanby et al. 2017). In midto-late northern winter, there is an elevated, hot north polar stratopause that cools by around 20 K by early northern spring, and becomes less elevated by later northern spring. The temperature behavior at the south pole is more complex, with the south polar stratopause warming slightly and becoming elevated soon after equinox, then cooling sometime between 2012 and 2015 while the lower stratosphere cools significantly (∼50 K at the south pole); the south polar stratopause then warms again in early 2016. The majority of the trace gases, except C 2 H 6 , are enhanced at the northern pole from winter through early northern spring, with the north polar abundances peaking in early northern spring, and the enhancement generally being stronger for those gases with shorter chemical lifetimes, with the notable exception of HCN, which is strongly enhanced while having the longest chemical lifetime of the observed gases. All of the observed gases become enhanced at high altitudes (pressures less than ∼0.5 mbar) at the south pole sometime in southern midfall between 2012 and 2015. The observed changes in temperature and trace gas abundance are broadly consistent with the meridional circulation predicted by GCMs Newman et al. 2011;Lebonnois et al. 2012;Lora et al. 2015), with a single meridional circulation cell rising at the summer pole and subsiding at the winter pole, and reversing direction around equinox. As pointed out by Teanby et al. (2017), the temporary cooling of the south polar stratopause can be explained by an increase in radiative cooling as the pole becomes enriched in radiatively active gases. The winter polar enhancements are generally confined to poleward of 60°, particularly for the strongly enhanced species, with sharp meridional gradients in C 4 H 2 , CH 3 C 2 H, and HCN, consistent with inhibition of meridional mixing at the edge of the polar vortex as discussed by Teanby et al. (2008aTeanby et al. ( , 2017 and Sharkey et al. (2021). Figure 6. Retrievals of temperature (top row), haze volume extinction at 100 cm −1 (second row), C 2 H 6 volume mixing ratio (third row), and C 2 H 2 volume mixing ratio (bottom row) from the T4 MIRLMBMAP observation on 2005 April 1. The left column shows retrievals using a 2D retrieval model; the right column shows retrievals using a 1D retrieval model.