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

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.


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. 1981Smith et al. , 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 McKay et al. 1989;Tomasko et al. 2008;Mishchenko et al. 2016), which vary spatially and with altitude in Titan's atmosphere Lorenz et al. 2001;Tomasko et al. 2005Tomasko et al. , 2008Anderson 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. 1989McKay et al. , 2001Rannou 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(Hutzell et al. , 1996. After the descent of the Huygens probe, a coupled photochemicalmicrophysical model of haze, incorporating various chemical pathways of haze formation (Lavvas et al. 2008b), selfconsistently reproduced the local distribution implied by measurements from the Descent Imager/Spectral Radiometer (DISR) instrument on board the probe (Tomasko et al. 2005(Tomasko et al. , 2008Lavvas et al. 2008bLavvas et al. , 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. 2015Vinatier et al. , 2020Karkoschka 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 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(Lorenz et al. , 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 1997and 1998 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 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. 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) 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.

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 Δλ ≈ 1 nm and a corresponding resolving power of R = λ/ Δλ ≈ 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.

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 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. 2017Teanby et al. , 2019Sylvestre 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).

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 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.

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.

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  , 5, 10, 15, 20, 25, 30, 35, 40, 45, 50, 55, 60, 65, 70, 75, 80 parameters. We then run a Levenberg-Marquardt (LM) minimization 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 variability, taken to be over the 2.20-2.35 μm region. 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 H with three sigma error bars for the same selection of observations as displayed in Figure 1. Early in the  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. 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.

Results
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 variations 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 distributions of haze scaling factor for each observation are provided in as the data behind Figure 4.

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. , 2004Lebonnois 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. , 2004Larson 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 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 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. 2017Teanby et al. , 2019Sylvestre 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 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.  Á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 (δτ 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 the to 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.   . 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. 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.