# Tracking Short-term Variations in the Haze Distribution of Titan's Atmosphere with SINFONI VLT

Fiona Nichols-Fleming, Paul Corlies, Alexander G. Hayes, Máté Ádámkovics, Patricio Rojo, Sebastien Rodriguez, Elizabeth P. Turtle, Juan M. Lora, Jason M. Soderblom

► To cite this version:

Fiona Nichols-Fleming, Paul Corlies, Alexander G. Hayes, Máté Ádámkovics, Patricio Rojo, et al. Tracking Short-term Variations in the Haze Distribution of Titan's Atmosphere with SINFONI VLT. The Planetary Science Journal, 2021, 2, 12 pp. 10.3847/PSJ/abffd7. insu-03589887

HAL Id: insu-03589887

https://insu.hal.science/insu-03589887

Submitted on 26 Feb 2022

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

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

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

Distributed under a Creative Commons Attribution 4.0 International License
<figure><img src="image_2.png" /></figure>

# Tracking Short-term Variations in the Haze Distribution of Titan's Atmosphere with SINFONI VLT

Fiona Nichols-Fleming<fnref n="1" />, Paul Corlies<fnref n="2" />,<fnref n="3" />, Alexander G. Hayes<fnref n="2" />, Máté Ádámkovics<fnref n="4" />, Patricio Rojo<fnref n="5" />, Sebastien Rodriguez<fnref n="6" />, Elizabeth P. Turtle<fnref n="7" />, Juan M. Lora<fnref n="8" />, and Jason M. Soderblom<fnref n="3" />

Received 2020 August 3; revised 2021 May 3; accepted 2021 May 6; published 2021 September 6

## Abstract

While it has long been known that Titan's haze and atmosphere are dynamic on seasonal timescales, recent results have revealed that they also exhibit significant subseasonal variations. Here, we report on observations of Titan acquired over an eight-month period between 2014 April and 2015 March with the Spectrograph for Integral Field Observations in the Near Infrared instrument on the Very Large Telescope using adaptive optics. These observations have an average five-day cadence, permitting interrogation of the short-period variability of Titan's atmosphere. Disk-resolved spectra in the *H* and *K* bands (1.4–2.4 µm) were analyzed with the PyDISORT radiative transfer model to determine the spatial distribution and variation of stratospheric haze opacity over subseasonal timescales. We observed a uniform decrease in haze opacity at 20°N and an increase in haze opacity at 250–300°E and ~40°N over the span of our observations. Globally, we found variations on the order of 5%–10% on timescales of weeks, as well as a steady, global increase in the amount of haze over timescales of months. The observed variations in haze opacity over the short timescales of our observations were of similar magnitude to long-period variations attributed to seasonal variation, suggesting rapid dynamical processes that may take part in the distribution of hazes in Titan's atmosphere.

*Unified Astronomy Thesaurus concepts: Ground-based astronomy (686); Atmospheric variability (2119);*
*Saturnian satellites (1427)*

*Supporting material: data behind figures*

## 1. Introduction

Titan's dense atmosphere and photochemically produced hydrocarbon haze have been prominent features of study, as the first resolved images of the moon were acquired by the Voyager missions four decades ago (Hanel et al. 1981; Smith et al. 1981, 1982). The haze distribution has been used to trace the global circulation of the atmosphere (Smith et al. 1981; Lorenz et al. 2001; Rannou et al. 2004; Larson et al. 2015; West et al. 2018), while the properties of sunlight scattered by the haze particles have been used to constrain their shape, size, and abundance distributions (Rages et al. 1983; McKay et al. 1989; Tomasko et al. 2008; Mishchenko et al. 2016), which vary spatially and with altitude in Titan's atmosphere (Rages & Pollack 1983; Lorenz et al. 2001; Tomasko et al. 2005, 2008; Anderson et al. 2008; Doose et al. 2016). Thus, the study of the physicochemical properties of Titan's haze and its evolution, along with trace gaseous species, are fundamental for better comprehension of the atmosphere's photochemistry, structure, and dynamics (Strobel 1974; Yung et al. 1984; McKay et al. 1989, 2001; Rannou et al. 2002; Wilson & Atreya 2004; Lavvas et al. 2008a; Tomasko et al. 2008; Krasnopolsky 2009).

Voyager observations showed a hemispheric asymmetry in the haze opacity, which was later linked to vertical and meridional motions above 300 km (Toon et al. 1992; Rannou et al. 2002; West et al. 2018). This circulation, as well as varying rates of haze production, were also invoked to reconcile the discrepancy between predicted and observed geometric albedo and brightness variations (Hutzell et al. 1993, 1996). After the descent of the Huygens probe, a coupled photochemical–microphysical model of haze, incorporating various chemical pathways of haze formation (Lavvas et al. 2008b), self-consistently reproduced the local distribution implied by measurements from the Descent Imager/Spectral Radiometer (DISR) instrument on board the probe (Tomasko et al. 2005, 2008; Lavvas et al. 2008b, 2011). Over the years, observations of atmospheric haze have identified significant structure and variability (Teanby et al. 2009; Penteado et al. 2010; Vinatier et al. 2015, 2020; Karkoschka 2016; Ádámkovics & de Pater 2017; Seignovert et al. 2017; West et al. 2018) and a wide variety of atmospheric models have been developed to investigate its formation, evolution, and interaction with the atmosphere's circulation (Toon et al. 1992; Rannou et al. 2004; Lavvas et al. 2008a; Lebonnois et al. 2012; Larson et al. 2015; Lora et al. 2015). However, the impacts of seasonal changes and the large-scale dynamics on the global distribution and seasonal variations of haze remain relatively unconstrained, so information about the temporal evolution of the haze on a variety of

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

Original content from this work may be used under the terms of the Creative Commons Attribution 4.0 licence. Any further distribution of this work must maintain attribution to the author(s) and the title of the work, journal citation and DOI.

<footnote n="1">¹Department of Earth, Environmental and Planetary Sciences, Brown University, Providence, RI 02912, USA; pcorlies@mit.edu</footnote>

<footnote n="2">²Department of Astronomy, Cornell University, Ithaca, NY 14853, USA</footnote>

<footnote n="3">³Department of Earth, Atmospheric, and Planetary Science, Massachusetts Institute of Technology, Cambridge, MA 02139, USA</footnote>

<footnote n="4">⁴Advanced Technology Center, Lockheed Martin, Palo Alto, CA 94304, USA</footnote>

<footnote n="5">⁵Departamento de Astronomia, Universidad de Chile, Santiago, Chile</footnote>

<footnote n="6">⁶Université de Paris, Institut de Physique du Globe de Paris, CNRS, Paris, France</footnote>

<footnote n="7">⁷Johns Hopkins Applied Physics Lab, Laurel, MD 20723, USA</footnote>

<footnote n="8">⁸Department of Earth and Planetary Sciences, Yale University, New Haven, CT 06520, USA</footnote>
timescales is helpful for better understanding the dominant processes in Titan's atmosphere.

Although the in situ measurements from the Huygens probe are the most detailed available, the haze is neither uniform spatially nor static temporally, as documented by a variety of observations. For example, Hubble Space Telescope (HST) data show that aerosol size varies with altitude (Karkoschka & Lorenz 1997). Subsequent HST observations, acquired two years later on the opposite hemisphere, found that these measurements differed significantly (Lorenz et al. 1999). Additional variations in haze brightness were observed on seasonal timescales as well as with latitude (especially near the pole) (Lorenz et al. 2001, 2004). Young et al. (2002) compiled a three-dimensional map of Titan's haze variations in altitude using six HST filters and showed that the atmospheric aerosols above 16 km had a peak opacity near the equator and the atmospheric aerosols below 16 km had peak optical depths in both north and south mid-latitudes. Meier et al. (2000) tentatively identified banded haze structures near the surface with HST observations from 1997 and 1998, and Karkoschka (2016) analyzed the 1997 through 2004 HST observations to show that there are two separately varying haze opacity components found above and below 100 km.

Similarly, analysis of Cassini Visual and Infrared Mapping Spectrometer (VIMS; Brown et al. 2004) acquired between 2004 and 2008 revealed both gradual gradients and discontinuities in haze distributions (Penteado et al. 2010; Rannou et al. 2010). VIMS also observed features such as the north polar hood (Rannou et al. 2012; Hirtzig et al. 2013), and a combination of VIMS, Cassini Imaging Science Subsystem (ISS), and Cassini Composite Infrared Spectrometer (CIRS) data were used to identify and study a tropical haze band (de Kok et al. 2010).

Haze asymmetries and evolution can also be observed from ground-based telescopes in the near-infrared (Gibbard et al. 1999; Roe et al. 2002; Ádámkovics et al. 2004, 2006). Gibbard et al. (2004) measured the 2.0 µm haze opacity with speckle and adaptive optics imaging at the W. M. Keck Observatory and found that the haze opacity in the southern hemisphere decreased by a factor of two between 1996 and 2004, while the equatorial haze opacity stayed approximately constant. Ádámkovics & de Pater (2017) used observations from the OH Suppressing Infrared Imaging Spectrograph (OSIRIS) instrument at the Keck observatory acquired between 2006 and 2015 to investigate seasonal and meridional variations in haze opacity, and found that the haze above 20 km showed significant nonlinear meridional variations, while the haze below 20 km varied nonmonotonically over seasonal timescales.

While it has been demonstrated that seasonally driven dynamics can redistribute hazes on global and seasonal timescales, little is known about the distribution of hazes at smaller spatial scales and short temporal intervals. Recent work with VIMS and ISS data has shown variations in the haze on timescales of hours to days (Carrasco et al. 2018; Rodriguez et al. 2018; West et al. 2018; Seignovert et al. 2021) during Cassini flybys of Titan, and ground-based telescopes can complement these observations by offering similar observation cadence over months-long temporal baselines.

Here, we present ground-based observations of haze opacity variations in the near-infrared that are complementary to the recent works of Ádámkovics & de Pater (2017), West et al. (2018), Carrasco et al. (2018) and Rodriguez et al. (2018). Our

44 observations, acquired from 2014 April through 2015 March, have higher temporal, spectral, and spatial resolution than previous ground-based near-infrared work, and cover a different season. We interpret these new haze evolution observations with a radiative transfer model, as presented in Section 3. We present the results of fitting the haze variations in Section 4, and discuss the results and conclude in Section 5.

## 2. Observations

Titan was observed from 2014 April to 2015 March using the Spectrograph for INtegral Field Observations in the Near Infrared (SINFONI) on the Very Large Telescope (VLT) at the Paranal Observatory. SINFONI is an integral field spectrograph that uses adaptive optics to allow for diffraction-limited spatial resolution of Titan's disk. As an integral field spectrometer, each pixel is fiber fed to a spectrograph, which has a spectral resolution of $\Delta\lambda \approx 1$ nm and a corresponding resolving power of $R = \lambda/\Delta\lambda \approx 1500$ (Eisenhauer et al. 2003; Bonnet et al. 2004).

Throughout the campaign, the *H* and *K* bands which correspond to the wavelength range 1.45–1.7 µm and 1.9–2.25 µm, respectively, were used simultaneously for each observation. These observations used a 0.′′8 × 0.′′8 field of view with a pixel scale of 25 mas per pixel, corresponding to 180 km at Titan's equator. The integration time for the standard star was four seconds in all of the observations and for Titan was fifteen seconds initially, increasing to a few minutes for later observations. Four exposures were coadded to cover the entire disk and increase signal to noise; observations were then reduced with the SINFONI data reduction pipeline.

The average time between observations was five days. Frequently, target of opportunity observations allowed for multiple observations of Titan on a given night. The longest gap in the observations occurred between 2014 October and 2015 February, when Saturn was in conjunction and therefore not observable from Earth.

## 3. Methods

### 3.1. Radiative Transfer Model

To analyze the reflected sunlight through Titan's atmosphere, we use the PyDISORT radiative transfer model, a Python implementation of the widely used DISORT radiative transfer code (Stamnes et al. 1988), which has been developed specifically for Titan (Ádámkovics et al. 2016; Ádámkovics & de Pater 2017). The Huygens Gas Chromatograph Mass Spectrometer (GCMS) and Atmospheric Structure Instrument (HASI) data provided the methane mole fraction, altitude, pressure, and temperature from around 1400 km to the surface (Fulchignoni et al. 2005; Niemann et al. 2010) and the DISR data provide the aerosol opacity structure, single-scattering albedos, and phase functions (Tomasko et al. 2008).

PyDISORT utilizes a plane parallel discrete ordinate routine to model the radiative transfer of Titan's atmosphere. It includes contributions from gaseous absorption, collision induced absorptions, aerosol (multiple) scattering, and reflections from Titan's surface. Gaseous absorption is modeled using the correlated *k*-coefficients that are interpolated on a temperature and pressure grid to relevant Titan-like conditions as a function of altitude. Likewise, collision induced absorptions are also modeled for Titan-like temperature and pressures (McKellar 1989; Lafferty et al. 1996; Hartmann et al. 2017).
To account for scattering, we use the haze scattering phase function as derived in Tomasko et al. (2008), extrapolating the phase function for altitudes >80 km to the surface, as prescribed in Campargue et al. (2012). Though DISR provided altitude-dependent values for the single-scattering albedo, these only apply to wavelengths <1.6 µm, and so we utilize a constant single-scattering albedo in both altitude and spectral region for each of our two regions of interest. These values are 0.85 for 1.65–1.72 µm and 0.75 for 2.15–2.22 µm, as derived in Hirtzig et al. (2013). Finally, we use a constant surface albedo of 0.10, although the surface albedo does not affect our results, given as we are not fitting surface-sensitive channels.

Over the course of the Cassini mission, variations in temperature of a few degrees on the surface (Jennings et al. 2019) and 20°–30° in the stratosphere (Schinder et al. 2020) have been observed. CIRS was able to measure vertical temperature profiles of Titan's atmosphere over the course of the Cassini mission (e.g., Vinatier et al. 2015; Teanby et al. 2017, 2019; Sylvestre et al. 2018; Coustenis et al. 2020; Mathé et al. 2020). Here we use seasonally minimum/maximum profiles as measured by CIRS to measure the influences of temperature in the model. To determine if these variation contribute significantly to our analysis, we model reasonable ranges of the observed temperature profile as measured with CIRS over the latitudes/season of our observations. While temperature variations in the stratosphere may vary on the scale of ≈20 K, we do not see variations in the modeled spectra above the noise of the observations and conclude that the single Huygens temperature profile is sufficient for this analysis (see Appendix A).

## 3.2. Altitude Sensitivity Tests

Ádámkovics & de Pater (2017) used two scaling factors to constrain the temporal variations in the haze opacity, one for altitudes <20 km and one for altitudes >20 km. The altitude of 20 km was chosen as an ad hoc critical altitude at which to split the model atmosphere. They found relatively stable retrievals for altitudes >20 km, but generally poorer fits (i.e., large spatial or rapid temporal variations) for altitudes <20 km, impacting the interpretation of these data.

To simplify this analysis, and to more cleanly define the critical altitudes at the tropopause, we adopt regions of interest for this work that are only sensitive to vertical opacity of Titan's hazes for altitudes >40 km. We limit ourselves to these altitudes to eliminate sources of error and degeneracies related to surface variations and any variability in tropospheric methane abundance (Lora & Ádámkovics 2017). To determine which wavelengths would be sensitive to changes in the haze opacity at altitudes above 40 km, a sensitivity test similar to the one conducted in Ádámkovics & de Pater (2017) was performed. This new test accounts for potential differences in sensitivity resulting from the differing spectral resolutions between OSIRIS and SINFONI as well as determines the wavelengths of interest for the 1.6 µm window. Figures 8 and 9 in Appendix B plot the results of this analysis for the 1.6 and 2.0 µm windows, respectively. As expected, we find that wavelength regions with higher *I*/*F* are sensitive to lower altitudes and vice versa, and that the regions of the spectra sensitive to the stratosphere are the far redward wings of the two methane windows which correspond to 1.65–1.8 and 2.15–2.30 µm. To prevent fitting to regions in which strong

<table><caption>**Table 1**
Summary of the Nodes Used to Generate a Look-up Table for Fitting Observations to the Radiative Transfer Model</caption><thead><tr><td>Parameter</td><td>Nodes</td></tr></thead><tbody><tr><td>Haze scaling factor</td><td>0.3, 0.5, 0.7, 0.9, 1.1, 1.3, 1.5, 1.7</td></tr><tr><td>Incident angle</td><td>0, 10, 20, 30, 40, 50, 60</td></tr><tr><td>Emission angle</td><td>0, 10, 20, 30, 40, 50, 60</td></tr><tr><td>Azimuth angle</td><td>0, 5, 10, 15, 20, 25, 30, 35, 40, 45, 50, 55, 60, 65, 70, 75, 80</td></tr></tbody></table>

methane opacity limits altitude sensitivity, we do not fit the full redward wings but rather fit 1.65–1.72 and 2.15–2.22 µm.

## 3.3. Data Reduction and Calibration

To remove telluric features from our data, we utilize the standard processing technique of observing a standard star to spectrally calibrate the observed flux from Titan. First, a sky template was created using the ESO SkyCalc tool for the time and pointing of the observation. Telluric features were then removed from the observed data by dividing the data by the sky template. Finally a photometric calibration was performed by measuring the observed flux of the standard star, following the process described in Ádámkovics et al. (2016), applying this correction independently to the *H* and *K* bands.

A small wavelength shift (1–2 spectral channels) was also found to result from the VLT calibration pipeline. To account for this shift we fit our observations (to get an approximate model spectrum), measured the offset between our observation and our model fit, shifted the observation the correct number of wavelength bins, and then finally refit the corrected observation to derive the best-fit haze scaling factor.

## 3.4. Spectrum Fitting

The process of fitting spectra with PyDISORT is computationally expensive. Each generation of a model spectrum takes about two minutes on a standard workstation and, as a result, the full minimization process requires ~20 minutes to find a best-fit spectrum, assuming only 10 spectral generations are required. To reduce the computation time of our fitting routine, we create a five dimensional look-up table over our spectral regions (1D), viewing geometry (3D), and haze scaling factor (1D); which we then use for the rapid inversion of the best-fit haze scaling parameter. Fitting our data in this way takes only seconds, dramatically improving computation time. For this analysis, the haze scaling factor is applied as a constant scaling in haze opacity relative to the DISR profile above 20 km and the wavelengths fit constrain this scaling further to the stratosphere as described in Section 3.2. The nodes of our look-up table are listed in Table 1. To validate this procedure, several fits were also performed using the full minimization process. The results for these test cases were identical to within error (1σ).

As mentioned previously, we only focus on the stratospheric haze layers, and therefore we fit the regions of the spectra sensitive only to the stratosphere. These regions are 1.65–1.72 and 2.15–2.22 µm and consist of 282 spectral channels.

To find best-fit parameters for our observations, we defined a function that will return a spectrum from our linearly interpolated look-up table for a given set of observational
<figure><img src="image_4.png" /><figcaption>**Figure 1.** Observed and best-fit spectra for one latitudinal band of the observations selected for Figure 5. The pixels of interest are shown in yellow on the inset image of the disk averaged between 1.56 and 1.61 µm, from which averaged spectra are shown in black. The best-fit model spectrum, determined from our minimization routine, is shown in green. The optimal haze scaling factor and latitude of each fit are listed on the associated panels along with the subobserver latitude and longitude, and date of the observation.</figcaption></figure>

parameters. We then run a Levenberg-Marquardt (LM) mini-
mization on this function with a constant error in the spectrum of
0.002 and an initial value for the haze scaling factor of 1.0. The
error is based on the spread of a histogram of the I/F values of
all the spectra used in this work in the wavelength range from
2.30 to 2.45 µm, following Ádámkovics et al. (2010). This
wavelength region is commonly considered a region of stable
signal, so it can be used to model the residual error in the
observations. We fit the histogram with a Gaussian, and the
square root of the variance of the best-fit histogram was used as
our spectral error. Our fitting routine returned errors on the
derived parameters based on the gradient of the steepest descent
in the minimization routine. For a more detailed explanation of
our consideration of errors see Appendix C.

To investigate the variations in Titan’s haze, we look at
averaged meridional variations across the center of Titan as
well as co-varying spatial variations across the entire disk. To
increase the signal to noise in our analysis, averages of several
pixels were used for these analyses. In the case of latitudinal
variations, we average 10 pixels centered on the horizontal
center of the disk for that observation in each vertical row.
When looking for two-dimensional spatial distributions, we
averaged pixels in 3 × 3 squares. In all cases, we require that

each pixel is inside of the disk and that the cosine of the
average incidence angle and the cosine of the average emission
angle were both greater than 0.80. This requirement removed
any nonphysical variations caused by the degradation of our
model at the extreme viewing geometries toward the edge of
the disk. Accounting for these restrictions we find a spatial
coverage of ~70° in latitude (from ~10°S to ~60°N) and
~100° in longitude. To normalize for inter-epoch variations,
each observation is scaled to a common value in a region of
high methane opacity, and therefore minimal spectral varia-
bility, taken to be over the 2.20–2.35 µm region.

4. Results

Figure 1 shows the observed and best-fit spectra for a selection of the 44 observations covering a variety of subobserver longitudes and temporal separations between observations. Overall, the model agrees well with the observations, with the exception of a slight overprediction of absorption in the ~1.7 µm region, suggesting an overestimation of the methane opacity in this spectral region.

Figure 2 shows the latitudinal distribution of the derived haze scaling factor *f*<sub>H</sub> with three sigma error bars for the same selection of observations as displayed in Figure 1. Early in the
<figure><img src="image_5.png" /><figcaption>**Figure 2.** Observed and best-fit spectra for the ~44°N latitude bin on 2015 March 8. The pixels of interest are shown in yellow on the inset image of the disk averaged between 1.56 and 1.61 µm, from which averaged spectra are shown in black. The best-fit model spectrum, determined from our minimization routine, is shown in green. The optimal haze scaling factor and latitude of the fit are listed on the figure along with the subobserver latitude and longitude, and date of the observation.</figcaption></figure>

campaign, a slight latitudinal gradient was observed from 0 to 20°N, with a remaining relatively flat scaling factor poleward of 20°N. In late 2014 May, we observed a sudden increase in the haze opacity at all latitudes, and the emergence of the 0–20°N enhancement. The following two months after this global enhancement show relatively stable haze distributions.

We find approximately flat distributions of the haze scaling factor in latitude for each observation, with the exception of a few interesting features; however, we found that the overall haze scaling factor could vary by up to 10% on timescales of weeks. One feature in particular that we saw in all of the observations, spanning a range of 12 months, was a dip in haze scaling factor at around 20°N latitude. Later observations showed an enhancement between this dip and the equator as well as an occasional swelling between 20 and 60°N latitude that flattens over time (see Figure 3).

The global variations from epoch-to-epoch suggest varia-
tions in the overall haze opacity on short timescales are similar
in magnitude to the variations in haze opacities observed on
seasonal timescales (Karkoschka 2016; Ádámkovics & de
Pater 2017). This suggests rapid dynamical variations in Titan's
stratosphere; however, exactly where in the atmospheric
column these variations are dominant is not constrained.

In addition to fitting constant latitude bins, we compiled a
two-dimensional map of haze scaling factor in latitude and
longitude using all of our observations, shown in Figure 4. This
plot shows the median haze scaling factor for each latitude and
longitude over the course of our observations, with a resolution
of one degree per pixel. The square areas that correspond to the
regions fit by our minimization routine were projected onto
latitude and longitude maps providing haze scaling factors for
those regions. There appears to be large-scale structure to the
spatial variation of the haze, including a local maximum in the
haze scaling factor at ~250–300°E and ~30–45°N, which
is present in all observations of the campaign, as well as a

minimum at ~20°N for all longitudes, the evolution of which is
plotted in Figure 3. The individual two-dimensional distribu-
tions of haze scaling factor for each observation are provided in
as the data behind Figure 4.

**5. Discussion and Conclusions**

Based on these 44 observations of Titan we observe global
changes in the haze on timescales of months that are similar in
magnitude to the previously observed seasonal changes
(Karkoschka 2016; Ádámkovics & de Pater 2017), as well as
smaller-scale variations in the spatial distribution of the haze on
weekly timescales.

Most notable in the observations is a 10% increase in the
haze scaling factor over the first few months of observations
(see Figure 5). Three-dimensional models that attempt to
simulate the distribution of Titan's haze suggest that an
increase in haze in the atmosphere above 200 km is expected
around the time of these observations, resulting from the
overturning circulation after northern vernal equinox (Rannou
et al. 2002, 2004; Lebonnois et al. 2012; Larson et al. 2015).
However, the speed of the enhancement is particularly stark;
while most models predict an increase in the haze during this
time, variations driven by circulation typically happen more
slowly, on the scale of two to three Earth years based on
changes to the altitude of the detached haze (see Figure 10 in
Larson et al. 2015). Thus, the observed rapid evolution of the
stratospheric haze over a period of a few months suggests that
rapid perturbations can be as important as seasonal variation in
determining the opacity of Titan's hazes (Rannou et al.
2002, 2004; Larson et al. 2015).

At latitudes greater than 50°N, we find a general decrease in
the haze opacity for all of our observations. While the limb
viewing geometry could produce similar effects, the asymmetry
observed between the equatorial and high northern latitudes in
<figure><img src="image_6.png" /><figcaption>**Figure 3.** Latitudinal distribution of haze scaling factor with three sigma error bars for a selection of observations spanning a variety of subobserver longitudes and temporal separations. The date and subobserver longitude of each observation is indicated in the top right of each panel.</figcaption></figure>

all our observations indicates that this decrease is real (see Figure 3) and is consistent with a Hadley cell circulation decreasing the abundance of hazes over the summer pole (Rannou et al. 2004; Lebonnois et al. 2012; Larson et al. 2015). It is believed that as Titan changes seasons, the mean meridional circulation clears the pole that is moving into summer of the hazes, driving them to the pole that is moving into winter (West et al. 2018). While we cannot observe southern latitudes, the profiles in the northern hemisphere are consistent with this interpretation. Furthermore, these observations provide a constraint on the timescale of the seasonal shifts of the overturning circulation of Titan's atmosphere, suggesting the meridional motions of hazes is well underway by several

years after northern vernal equinox (2009 August), but before northern summer solstice (2017 May). Finer constraints on the time, duration, and structure of the transition between seasons will require similarly cadenced observations as those presented here, but over a longer temporal baseline.

Another interesting feature is the evolution of hazes at equatorial to mid-latitudes. The campaign begins showing a general decrease in haze opacities at these latitudes, but the observations subsequently evolve to form a local enhancement in the hazes from 0 to 20°N for much of the campaign (see Figure 3). Comparison with observations from Ádámkovics & de Pater (2017) in 2015 July, which show no such enhancement, suggests that this feature may have dissipated on the timescale of
<figure><img src="image_7.png" /></figure>

**Figure 4.** Spatial distribution of the median haze scaling factor for all of our observations. Fits are projected onto a latitude and longitude map with a resolution of one degree per pixel. The medians of each pixel on the maps are compiled to create an overall spatial distribution and then averaged in either dimension to produce equatorial and meridional variations over the full length of our observations. The error bars on the equatorial and meridional variations represent the standard deviation of the values in the associated latitudinal or longitudinal strip. The data behind this figure is available in a numpy array in the .tar.gz package.
(The data used to create this figure are available.)

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

**Figure 5.** Temporal variability of Titan's haze scaling factor from our observations. Disk-averaged values of haze scaling factor are shown with $1\sigma$ error bars for each observation representing the propagation of errors for each spectra fit across the disk. Our observations show a monotonic global increase for the first few months of the campaign along with some short-term variability on the timescale of weeks. The data behind this figure is available in a numpy array and machine readable format.
(The data used to create this figure are available.)

a few months, capping this local enhancement to approximately 12 months in duration. We cannot determine whether the variations in hazes in our observations are caused by a source or sink of the stratospheric hazes and more 3D modeling of the photochemical haze production and dynamics are required.

One possibility could be the result of rising motion in the stratosphere resulting from the local insolation maximum at these latitudes for this season (Lora et al. 2015). Indeed, upwelling in the equatorial region during this time period is predicted to explain the independent measurements of Bézard et al. (2018) and Vinatier et al. (2020) and 3D GCMs predict an enhancement in hazes above 200 km at equatorial latitudes (see Figure 5 of Lebonnois et al. 2012, Figure 10 of Larson et al. 2015). The results presented here also favor the prediction of Larson et al. (2015) of the reappearance of the detached haze layer at

equatorial latitudes between mid 2014 to early 2015. Although the reappearance of the detached haze layer was not observed until late 2015 or early 2016 (West et al. 2018; Seignovert et al. 2021), the haze extinction in the atmosphere above 350 km varied significantly at the time of our observations (see Figure 2 of West et al. 2018; Figure 15 of Seignovert et al. 2021). This upper atmosphere variability is consistent with our results, suggesting the high variability observed in the upper atmosphere by ISS is also observed in the stratosphere. Variations in haze opacity of this magnitude were also observed in the atmosphere above 70 km by VIMS between 2009 May and 2010 July (see Supplementary Information of Rodriguez et al. 2018).

Further, the local enhancement results in the apparent appearance of a dark band at ~20°N (see Figure 4), which could also be a similar feature to the dark bands observed by
the Voyager 1 mission at northern mid-latitudes (Smith et al. 1981). Our observations are approximately half-way between Titan's spring equinox and summer solstice, and the Voyager 1 encounter occurred around Titan's spring equinox, suggesting the possibility of the seasonal recurrence of this feature. Combined with the local enhancement at equatorial latitudes described previously, this would suggest a dynamic origin for this feature.

If diurnal variations in the haze exist (previously seen in Coustenis et al. 2001; Hirtzig et al. 2006), we do not see strong variations between 9 am and 3 pm local time in our observations. This lack of visible variation, therefore, places a limit on the times and magnitudes over which the diurnal differences in Titan's haze may exist.

Finally, we compare our observations to previous work, which has shown a long-term variability in Titan's hazes (Karkoschka 2016; Ádámkovics & de Pater 2017). Our work complements these longer baseline studies by observing at a higher temporal cadence. From this comparison we find short timescale variations in the haze of the same magnitude in addition to a consistent increase in global haze over the first few months of the observing campaign (see Figure 5). The data sets from this work as well as Karkoschka (2016) and Ádámkovics & de Pater (2017) were each constructed using a unique model for observations from a different instrument (SINFONI, HST, and OSIRIS, respectively), therefore we do not mean to compare the absolute agreement of the three models, but rather compare the results to show context for the magnitude of the variations seen in this work. Combined, these three sets of observations suggest variations in the hazes at multiple timescales, thus requiring high temporal cadence to fully understand the dynamical properties of Titan's hazes, and suggests that long-term trends must also be considered in the context of short-term variability to understand and interpret the complex dynamics of Titan's atmosphere from haze circulation. A more thorough analysis using multiple telescopes could better inform the relative contributions of short-term and seasonal variability of Titan's hazes.

These observed changes in haze opacity, both globally and on smaller scales, indicate that there is either more variability in the production of stratospheric haze, or more variability in the redistribution of hazes into the stratosphere, than previously appreciated. While we have proposed some explanations for observed features, more detailed work on GCMs is necessary beyond the scope of this work. These short timescale variations have been seen previously in the detached haze and atmosphere from Cassini data (Rodriguez et al. 2018; West et al. 2018;

Seignovert et al. 2021), but this work shows that they can be observed from the ground as well and that further ground-based monitoring campaigns can provide valuable insight into the types of variations and features described in this work.

This work was based on observations collected at the European Southern Observatory under ESO programs 093.C-0557 and 094.C-0422. Fiona Nichols-Fleming gratefully acknowledges the summer 2017 REU program in the Cornell Center for Astrophysics and Planetary Science at Cornell University and the NSF for financial support under NSF award AST-1659264. Paul Corlies gratefully acknowledges NASA for financial support under NASA Earth and Space Science Fellowship NNX14AO31H S03. We thank Peter Gierasch and Maryame El Moutamid for providing us with CIRS temperature data to constrain the model sensitivity to variations in temperature. Finally, we thank two anonymous reviewers for the careful and considerate reviews, which have helped to significantly improve the final manuscript.

# Appendix A
## The Effect of Temperature Variations on Spectra

CIRS was able to measure vertical temperature profiles of Titan's atmosphere over the course of the Cassini mission (e.g., Vinatier et al. 2015; Teanby et al. 2017, 2019; Sylvestre et al. 2018; Coustenis et al. 2020; Mathé et al. 2020). To quantify the importance of temperature variations on the derived haze scaling factor, we used three temperature profiles measured by CIRS for northern, southern, and equatorial latitudes as well as the temperature profile measured by Huygens (see Figure 6) to produce four different modeled spectra for the same viewing geometry. There are slight variations between these spectra, as shown in Figure 7, but these variations are well within our error of 0.002 indicating that there is not significant error introduced by using the Huygens temperature profile at all latitudes. The atmospheric temperatures can also vary throughout Titan's seasons, which may introduce other errors as our observations are from a different season than when the Huygens temperature profile was measured, but these variations should be on the same order as those in Figure 6, and therefore it should be the case that these effects are well within our error of 0.002 as well.

# Appendix B
## Altitude Sensitivity Analysis

To determine which SINFONI wavelengths are sensitive to which altitudes, we use a methodology similar to that of
<figure><img src="image_9.png" /><figcaption>**Figure 6.** Measured temperature profiles from Huygens and three additional latitudinal regions from CIRS. There are minimal temperature variations observed for the surface with larger variations up to about 40 K in the upper part of the stratosphere.</figcaption></figure>

<figure><img src="image_10.png" /><figcaption>**Figure 7.** Modeled spectra for the same viewing geometry using the four temperature profiles shown in Figure 6. Uniform error bars of 0.002 applied to the modeled spectrum using the Huygens profile are shown by the gray envelope. Modeled spectra using the temperature profiles from the northern, southern, and equatorial latitudes are shown in red, blue, and green, respectively. Differences between the modeled spectra become more pronounced for the longer wavelengths.</figcaption></figure>
Ádámkovics & de Pater (2017), in which we construct a high-resolution model that sampled Titan's atmosphere every ~2.5 km over the entire spectral range of our observations. For each altitude, we insert a small, homogeneous addition of haze opacity ($\delta\tau_H = 0.1$) for that layer and generate a new model for comparison to the un-altered reference. As the hazes are efficient scatterers, increasing their opacity acts to increase Titan's reflectivity, with the strongest variations corresponding to the altitude at which the wavelength is most sensitive. In this way, we empirically determine the altitude

each SINFONI channel is most sensitive to, similar to a contribution function. Figures 8 and 9 display the result of this analysis. In all, we find wavelengths in the methane windows, and therefore weaker methane absorption, to be sensitive all the way to Titan's surface. Conversely, wavelengths in the methane band, and therefore stronger methane absorption, are only sensitive to the highest altitudes. For our final selection, we chose a combination of wavelengths most sensitive to Titan's stratosphere from approximately 40 to 300 km.

<figure><img src="image_11.png" /><figcaption>**Figure 8.** Figure showing the sensitivity of wavelength channels in the 1.6 µm window to changes in haze opacity at particular altitudes. The channels between 1.65 and 1.8 µm are sensitive to changes in the stratosphere (above 40 km). The black line in the upper panel is the original *I/F* without the addition of any haze variation.</figcaption></figure>

<figure><img src="image_12.png" /><figcaption>**Figure 9.** Figure showing the sensitivity of wavelength channels in the 2 µm window to changes in haze opacity at particular altitudes. The channels between 2.15 and 2.3 µm are sensitive to changes in the stratosphere (above 40 km). The black line in the upper panel is the original *I/F* without the addition of any haze variation.</figcaption></figure>
<figure><img src="image_13.png" /><figcaption>**Figure 10.** Histograms of derived haze scaling factors for three MCMCs with the solid black line indicating the best-fit scaling factor. Each spectrum was offset uniformly by a random value sampled from a Gaussian distribution centered on zero with a standard deviation of one, two, or three times 0.002, the error of our observations, before being fit by our minimization routine.</figcaption></figure>

<figure><img src="image_14.png" /><figcaption>**Figure 11.** Histograms of derived haze scaling factors for an MCMC with the solid black line indicating the best-fit scaling factor. Each wavelength of the spectrum was offset by a random value sampled from a Gaussian distribution centered on zero with a standard deviation of 0.002, the error of our observations, before being fit by our minimization routine. A best-fit Gaussian distribution with a mean of 1.043 and a standard deviation of 0.0127 is shown in blue as well.</figcaption></figure>
## Appendix C
### Justification of Errors

To get a sense for an upper bound on the errors associated with our derived haze scaling factors, we ran three MCMCs in which a chosen spectrum is uniformly offset by a value randomly sampled from a normal distribution centered at zero with a standard deviation of one, two, or three times the error in the observations of 0.002 as used in the main text. Figure 10 contains the histograms of the results from these MCMCs where the solid black line in each is the best-fit scaling factor for the original, unchanged spectrum. As expected, with larger offsets allowed, the distribution of scaling factors becomes wider. A fit of the histogram in Figure 10(a) provides an estimate of the error for a systematic offset in the data of approximately 20%. This is likely an overestimate of the true errors as it would represent a systematic error source.

To get a more accurate representation of errors, we ran another MCMC where each channel of our selected spectrum was given a unique offset randomly sampled from the Gaussian distribution centered on zero with a standard deviation of 0.002. This produced the histogram shown in Figure 11 where the solid black line represents the best fit for the data without an offset. When this distribution is fit with a Gaussian (shown in blue on Figure 11), the estimate of the errors is approximately 1% which is on the same order as the errors reported by our minimization routine justifying our use of the minimization derived errors as our reported errors.

## ORCID iDs

Fiona Nichols-Fleming https://orcid.org/0000-0002-7700-5139

Paul Corlies https://orcid.org/0000-0002-6417-9316

Máté Ádámkovics https://orcid.org/0000-0003-1869-0938

Sebastien Rodriguez https://orcid.org/0000-0003-1219-0641

Juan M. Lora https://orcid.org/0000-0001-9925-1050

Jason M. Soderblom https://orcid.org/0000-0003-3715-6407

## References

* Ádámkovics, M., Barnes, J. W., Hartung, M., & de Pater, I. 2010, *Icar*, 208, 868

* Ádámkovics, M., & de Pater, I. 2017, *Icar*, 290, 134

* Ádámkovics, M., de Pater, I., Hartung, M., et al. 2006, *JGRE*, 111, E07S06

* Ádámkovics, M., de Pater, I., Roe, H. G., Gibbard, S. G., & Griffith, C. A. 2004, *GeoRL*, 31, L17S05

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

* Anderson, C. M., Young, E. F., Chanover, N. J., & McKay, C. P. 2008, *Icar*, 194, 721

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

* Bonnet, H., Conzelmann, R., Delabre, B., et al. 2004, *Proc. SPIE*, 5490, 130

* Brown, R. H., Baines, K. H., Bellucci, G., et al. 2004, *SSRv*, 115, 111

* Campargue, A., Leshchishina, O., Wang, L., et al. 2012, *JQSRT*, 113, 1855

* Carrasco, N., Tigrine, S., Gavilan, L., Nahon, L., & Gudipati, M. S. 2018, *NatAs*, 2, 489

* Coustenis, A., Gendron, E., Lai, O., et al. 2001, *Icar*, 154, 501

* Coustenis, A., Jennings, D. E., Achterberg, R. K., et al. 2020, *Icar*, 344, 113413
de Kok, R., Irwin, P. G. J., Teanby, N. A., et al. 2010, *Icar*, 207, 485

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

* Eisenhauer, F., Abuter, R., Bickert, K., et al. 2003, *Proc. SPIE*, 4841, 1548

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

* Gibbard, S. G., de Pater, I., Macintosh, B. A., et al. 2004, *GeoRL*, 31, L17S02

* Gibbard, S. G., Macintosh, B., Gavel, D., et al. 1999, *Icar*, 139, 189

* Hanel, R., Conrath, B., Flasar, F. M., et al. 1981, *Sci*, 212, 192

* Hartmann, J.-M., Boulet, C., & Toon, G. C. 2017, *JGRD*, 122, 2419

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

* Hirtzig, M., Coustenis, A., Gendron, E., et al. 2006, *A&A*, 456, 761

* Hutzell, W. T., McKay, C. P., & Toon, O. B. 1993, *Icar*, 105, 162

* Hutzell, W. T., McKay, C. P., Toon, O. B., & Hourdin, F. 1996, *Icar*, 119, 112

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

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

* Karkoschka, E., & Lorenz, R. D. 1997, *Icar*, 125, 369

* Krasnopolsky, V. A. 2009, *Icar*, 201, 226
* Lafferty, W. J., Solodov, A. M., Weber, A., Olson, W. B., & Hartmann, J.-M. 1996, *ApOpt*, 35, 5911

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

* Lavvas, P., Griffith, C. A., & Yelle, R. V. 2011, *Icar*, 215, 732

* Lavvas, P. P., Coustenis, A., & Vardavas, I. M. 2008a, *P&SS*, 56, 67

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

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

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

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

* Lorenz, R. D., Lemmon, M. T., Smith, P. H., & Lockwood, G. W. 1999, *Icar*, 142, 391

* Lorenz, R. D., Smith, P. H., & Lemmon, M. T. 2004, *GeoRL*, 31, L10702

* Lorenz, R. D., Young, E. F., & Lemmon, M. T. 2001, *GeoRL*, 28, 4453

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

* McKay, C. P., Coustenis, A., Samuelson, R. E., et al. 2001, *P&SS*, 49, 79

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

* McKellar, A. R. W. 1989, *Icar*, 80, 361

* Meier, R., Smith, B. A., Owen, T. C., & Terrile, R. J. 2000, *Icar*, 145, 462

* Mishchenko, M. I., Dlugach, J. M., Yurkin, M. A., et al. 2016, *PhR*, 632, 1

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

* Penteado, P. F., Griffith, C. A., Tomasko, M. G., et al. 2010, *Icar*, 206, 352

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

* Rages, K., Pollack, J. B., & Smith, P. H. 1983, *JGR*, 88, 8721

* Rannou, P., Cours, T., Le Mouélic, S., et al. 2010, *Icar*, 208, 850

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

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

* Rannou, P., Le Mouélic, S., Sotin, C., & Brown, R. H. 2012, *ApJ*, 748, 4

* Rodriguez, S., Le Mouélic, S., Barnes, J. W., et al. 2018, *NatGe*, 11, 727

* Roe, H. G., de Pater, I., Macintosh, B. A., et al. 2002, *Icar*, 157, 254

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

* Seignovert, B., Rannou, P., Lavvas, P., Cours, T., & West, R. A. 2017, *Icar*, 292, 13

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

* Smith, B. A., Soderblom, L., Batson, R. M., et al. 1982, *Sci*, 215, 504
* Smith, B. A., Soderblom, L., Beebe, R. F., et al. 1981, *Sci*, 212, 163

* Stamnes, K., Tsay, S. C., Jayaweera, K., & Wiscombe, W. 1988, *ApOpt*, 27, 2502

* Strobel, D. F. 1974, *Icar*, 21, 466
* Sylvestre, M., Teanby, N. A., Vinatier, S., Lebonnois, S., & Irwin, P. G. J. 2018, *A&A*, 609, A64

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

* Teanby, N. A., de Kok, R., & Irwin, P. G. J. 2009, *Icar*, 204, 645

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

* Tomasko, M. G., Archinal, B., Becker, T., et al. 2005, *Natur*, 438, 765

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

* Toon, O. B., McKay, C. P., Griffith, C. A., & Turco, R. P. 1992, *Icar*, 95, 24

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

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

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

* Wilson, E. H., & Atreya, S. K. 2004, *JGRE*, 109, E06002

* Young, E. F., Rannou, P., McKay, C. P., Griffith, C. A., & Noll, K. 2002, *AJ*, 123, 3473

* Yung, Y. L., Allen, M., & Pinto, J. P. 1984, *ApJS*, 55, 465