Titan’s Atmospheric Albedo Asymmetry and Seasonal Variability Observed through the Cassini Imaging Science Subsystem

Using images from Cassini, we analyzed the north–south albedo asymmetry that has been observed in the atmosphere of Saturn’s moon, Titan. Suitable images from the Cassini Imaging Science Subsystem taken at 889 nm spanned from 2004 to 2017—around half of a Titan year—and revealed seasonal changes in the characteristics and orientation of the north–south asymmetry boundary. Such circumglobal features provide insight into the dynamics and circulation of the atmosphere more broadly. The albedo asymmetry has been observed to reverse for part of the Titan year, inverting the brighter and darker hemispheres; we also observed this inversion, along with the formation of additional banding briefly during the transition (around 2014–2016). A tilt in the rotation axis of Titan’s atmosphere with respect to the solid body rotation has previously been noted. Using robust edge-detection techniques, we likewise identified a tilt offset of a few degrees in the albedo transition boundaries. The azimuth of this tilt axis remained roughly fixed in inertial space, with some smaller possible seasonal fluctuations around the fixed direction noted.


Introduction
The atmosphere of Titan is a compelling area of research, not only for its uniqueness as the most substantial atmosphere of any moon in the solar system and as a rare example of atmospheric super-rotation, but also for the many similarities it shares with Earth's atmosphere.Titan is shrouded in a thick layer of haze made up of organic species produced through photochemical processes (Yung et al. 1984).There are a few notable features that have been observed in the atmosphere, including the detached haze layer, the north-south albedo asymmetry (NSA), polar hoods, and a number of more subtle bands at varying latitudes throughout the Titan year (Smith et al. 1981(Smith et al. , 1982;;Sromovsky et al. 1981;Tomasko & Smith 1982;Lorenz et al. 1997).The NSA, polar hoods, and other banding in particular have been observed to exhibit seasonal variations (Caldwell et al. 1992;Lorenz et al. 1997;Kutsop et al. 2022), and are crucial to understanding Titan's atmospheric dynamics and photochemistry.Here, we use the Cassini Imaging Science Subsystem (ISS) to characterize the NSA boundary and its variation throughout the Cassini mission.We also discuss its implications for Titan's atmospheric circulation more generally.
The atmosphere of Titan has a variation in albedo between the northern and southern hemispheres.The difference is most clearly visible in the Cassini ISS methane band filter at 889 nm (MT3), but it can also be seen in images using green, blue, violet, and other methane band filters.The feature was first observed in the 1980s in images captured by Voyager 1 and 2 (Sromovsky et al. 1981;Smith et al. 1982) and Pioneer 11 (Tomasko & Smith 1982).Through subsequent observations from the Hubble Space Telescope and Cassini, seasonal dependence has been observed, including in the north/south albedo ratio (Lorenz et al. 1997), the overall disk brightness (Lockwood & Thompson 2009), and the reversal of the bright and dark hemispheres for half of the orbital period (Caldwell et al. 1992).
Since the NSA and other band-like atmospheric features likely trace the atmospheric circulation, the characteristics of such features are a useful tool to inform our understanding of Titanʼs atmospheric dynamics more broadly.Previous studies have noted an offset or tilt in the rotation axes of several atmospheric features relative to the rotation axis of the solid body of Titan: the north polar zone (Sromovsky et al. 1981), the NSA boundary (Roman et al. 2009), middle atmosphere isotherms (Achterberg et al. 2008a(Achterberg et al. , 2011)), clouds (West et al. 2016), and gas abundances (Teanby et al. 2010;Sharkey et al. 2020).One recent study (Kutsop et al. 2022) also investigated seasonal trends in the orientation of polar and equatorial annuli using Cassini Visible and Infrared Mapping Spectrometer (VIMS) data.These previous studies have mostly indicated a tilt of around a few degrees for their respective features.In studies with sufficient temporal coverage (e.g., Achterberg et al. 2011 andKutsop et al. 2022), the azimuth of the tilt axis appears to be roughly fixed in inertial space, but may exhibit some seasonal fluctuations around some inertially fixed point offset from the north pole.General circulation models (GCMs) have shown similar results for seasonal fluctuations in tilt amplitude, but differ somewhat from existing observations of how the tilt azimuth changes throughout the Titan year (e.g., Tokano 2010).With relatively few studies of the atmospheric tilt over sufficiently long time spans, the seasonal dependence of the tilt orientation is not well constrained.Thus, further constraining the temporal trends in the tilt axis was a main objective in this study.
Overall, the tilt itself and the mechanism behind it are not yet well understood.Titanʼs atmosphere is unique with its superrotation and varying insolation throughout the year (see, e.g., Flasar & Achterberg 2009).Characterizing the tilt and seasonal dependence of Titanʼs atmosphere is an important step toward more accurate models and a more complete understanding of Titanʼs atmospheric dynamics.Thus, this work aims to better constrain the characteristics of the tilt and of the NSA and other atmospheric features visible at times throughout the Cassini mission using images that span roughly 14 yr, or nearly half of a Titan year.The expectation is that this careful analysis of the orientation of these features presumably advecting with the stratospheric winds can help constrain models of Titan's atmospheric dynamics.Additionally, this analysis expands on ideas and findings of previous research to provide a more comprehensive description of the seasonal variability and the semiannual hemispheric albedo reversal exhibited by Titanʼs atmosphere.This tracking of the temporal progression of these contrast features in Titan's atmosphere should also serve as a constraint on modeling (both theoretical and numerical) aimed at deciphering Titan's atmospheric circulation.Such work is especially imperative as major upcoming missions focus in on Titan, and it may have important implications for other superrotating planetary atmospheres as well (e.g., Venus).

Observations and Processing
The Cassini ISS collected thousands of images of Titan over the course of the mission.The ISS was made up of two cameras, a wide-angle camera and a narrow-angle camera, and used a variety of wavelength filters.For this work, images were chosen based on filter, spatial resolution, and viewing geometry.Only images from the 889 nm MT3 filter were used, as the NSA boundary was most prominent at this wavelength.We also required that the entire disk of Titan was visible in the image to allow for accurate navigation and calibration, and that the diameter spanned at least 100 pixels.Images were further constrained to those with phase angles below 70°to ensure the majority of the disk was visible.While images with phase angles greater than 30°and/or diameters less than 300 pixels had higher uncertainties in their measurements, we chose to include them to provide better temporal coverage, particularly later in the mission when more suitable images were not always available.A small number of additional images were removed from the data set that contained other bodies or anomalies that interfered with the analysis.A set of around 200 images remained that were suitable for our analysis, spanning 2004 to 2017 (see Appendix B for full listing and details of images used).For more details about the Cassini instrumentation, see Porco et al. (2004).
Images were identified and downloaded from the Planetary Data System (PDS) using the PDS Ring-Moon Systems Node's OPUS search service. 3The images had already been calibrated by the Cassini team using CISSCAL software, which removed many of the imaging artifacts and bad pixels (Knowles et al. 2020).Significant noise and both horizontal and vertical banding were still noted, however, so outliers were removed with median filters, and banding was reduced by subtracting out the vertical and horizontal means of the dark space around the planet disk in each image and median-filtering with a long, narrow boxcar.
USGS Integrated Software for Imagers and Spectrometers (ISIS; Sucharski et al. 2020) has convenient methods for navigating SPICE kernels and generating image files with the geometric information needed for image analysis.While the position and angle of the camera with respect to Titan was generally accurate enough for the analysis performed here, the center pointing was often visibly inaccurate, needing adjustments of a few to several pixels for most images.Similarly to Roman et al. (2009), the center-pointing adjustments were determined by fitting a circle to the disk formed by Titan's outermost haze layers.This was based on the findings from Achterberg et al. (2008b) that atmospheric rotation decays above 250 km and is very slow near altitudes of 500 km, resulting in the atmosphere at this height taking on a nearly spherical shape that is centered on the solid body.The lower haze layers that are subjected to faster, variable rotation and other dynamics result in a more complex, asymmetric oblate shape without an easily identified center point.However, while Roman et al. (2009) performed their fit to the detached haze layer near 500 km, which was clearly visible and relatively invariable over the time period their study covered, later in the Cassini mission the detached haze layer was found to both decrease in altitude and become more oblate before disappearing into the main haze layer entirely in the years surrounding the equinox-it eventually reappeared in 2016 (West et al. 2011(West et al. , 2018;;Seignovert et al. 2021).To work around this, we allowed the circle radius to vary and fit to the outermost visible haze layer in each image.While this led to some images in the latter half of the mission getting fit at lower altitudes where the atmosphere is less spherical, these were typically still well above 250 km, and close enough to spherical to achieve a good fit.All fits were also checked by eye to confirm reasonable results.Due to this variability and complex structure of Titan's outermost haze layers, even after these corrections the center pointing is still a source of uncertainty of up to a few pixels, which is taken into account in our error bars later on.
Assuming Titan's atmosphere at the altitudes of interest exhibits a slightly oblate spheroid shape, using the corrected pointing information, each pixel within the disk was assigned its proper coordinates and incidence and emission angles.Rather than mapping the image to a different projection, we chose to keep the data in image space to avoid the complications associated with stretching data to different resolutions.
The polar and equatorial radii used for mapping were given by the radii of the solid body of Titan plus the estimated altitude of the NSA feature.Due to spherical geometry and atmospheric opacity, data closer to the center of the disk probe deeper in the atmosphere than data near the limb.At 889 nm, measurements at the limb only reach to around 200 km altitude, while methane absorption in the lower altitudes of Titanʼs atmosphere means that any bright features near the center of the disk stem from haze above 60 km (Lorenz et al. 2001).The altitude of the NSA boundary itself is not well constrained and spans a range or ranges of altitudes; however, it was necessary to map the features to a discrete altitude for this analysis.Thus, similar to Roman et al. (2009), a Monte Carlo-like approach was used to generate results from around 50 to 200 km above the surface to cover this range of most likely altitudes.The variance in the results was later used to inform our error bars.
Images of Titan in the methane bands exhibit significant limb brightening, which was found to interfere with our analysis of the atmospheric tilt.Thus, it was necessary to flatten the brightness curve of the images as much as possible.While some attempts have been made to describe this brightening as a function of incidence and emission angles (e.g., Young et al. 2002), none of these yielded a flat enough image for our purposes, particularly very near the limb.
By binning the pixels from images in our data set according to their corresponding incidence and emission angles, it was possible to determine the average brightness of a pixel at a given incidence and emission angle.To isolate the effects of limb brightening/darkening from the albedo variations due to the NSA boundary and other atmospheric features, we selected images taken prior to the onset of the NSA reversal around 2014 and only included pixels north of the NSA boundary region (latitudes > −5°).The albedo is relatively flat north of the NSA boundary in this time period.The selected images spanned nearly a decade and a variety of lighting geometries, which further reduced correlations between any specific coordinates and incidence/emission angles (thus decoupling albedo and limb brightening) when the pixels were binned and averaged.We noted that phase angle did not seem to have a significant impact on limb brightening over the limited range of phase angles (0°-70°) included in our image set, so it was sufficient in this case to base the model on incidence and emission angle only.The resulting brightness values were tested on numerous images throughout the entire data set to ensure that this model consistently flattened limb effects at all phase angles included in our study, in both hemispheres, and in images beyond 2014 when the NSA reversed.Some smoothing and manual fine-tuning was performed as needed.This empirical determination of brightness as a function of emission and incidence angle enabled sufficient flattening of the images so that even subtle albedo variations near the limb could be identified.Higher phase angles and more of each image could be used here than in the previous similar study by Roman et al. (2009), resulting in more data points with less uncertainty.See Figure 1 depicting brightness based on incidence and emission angle, and Figure 2 for an example of a calibrated image.See Appendix A for a table of the resulting values.

Analysis
We needed a method to identify the location of the NSA boundary that was robust enough to handle the variety of image geometries, variations in the boundaryʼs shape within and between images, and other changes and features that developed over time.When methods previously used by Roman et al. (2009) yielded unsatisfactory results, particularly for images later in the mission, edge-detection and computer-vision techniques were used to develop a more sophisticated algorithm.The method we settled on resembles a specialized Canny edge detector, following a similar multistep process of smoothing, finding the intensity gradient of the image, thresholding/identifying possible edges, and edge tracking.
In addition to the noise removal described in the calibration steps, a low-pass blurring filter (which calculated mean pixel values with a 5 × 5 to 11 × 11 kernel, depending on image resolution) was applied to reduce any remaining noise, artifacts, and smaller atmospheric features that were not of interest.After smoothing, the value of the intensity gradient was determined at each usable pixel in the visible disk, where larger gradients indicate stronger albedo transitions.Since the boundaries we were interested in lie roughly along latitude lines, a directional derivative of intensity was calculated for every pixel in the local northward direction using a 3 × 3 kernel, which weighted the surrounding pixels according to the direction.This maximized the signal from the latitudinal albedo variations, while any lingering brightness variations due to limb effects were suppressed since they ran roughly perpendicular to the direction vectors for most pixels in most image geometries.
The array of directional derivatives was limited to pixels between −40°and 40°latitude, and divided longitudinally into 100 slices of equal pixel area.The values of the derivatives in each slice were then fit to a degree-10 polynomial as a function of latitude, where local maxima and minima were considered possible boundaries.The high degree of polynomial allowed the several most significant albedo transitions in each slice to be identified, and it provided a robust determination of the approximate midpoint of each transition despite some variation in the shape of the curves across the image.An example of this polynomial fitting is shown at the bottom of Figure 3.
With the latitudes of all possible major edges identified for each slice, the final step was to track these edge candidates across the image.Each maxima and minima was sorted into a group with other maxima and minima, respectively, that were at similar latitudes in nearby slices.Any group with more than 90 related edge candidates (out of the 100 slices) was considered a strong edge and used in the subsequent analysis.Groups with less than 90 were thrown out as false edges or weak edges with too few points for accurate analysis.See Figure 3 for an example of strong versus weak edge points in an image.
Finally, the location and orientation of the boundaries were determined by fitting the points in each group to a plane.Leastsquares fitting was performed in two steps: the normal vector of the plane was determined first, followed by the mean latitude.This reduced the correlation between the tilt of the boundary and the mean latitude, leading to more consistent results from the fitting functions.

Error Analysis
Throughout the analysis, possible sources of error were identified, mitigated, and accounted for in the reporting of the results.Uncertainties in center pointing and altitude of features, image artifacts and noise, physical variations, and fitting errors affected the accuracy of the calculations of atmospheric tilt.To address the many, sometimes complex sources of error, a Monte Carlo-like method was applied.The analysis was performed repeatedly on each image while adjusting each of the image parameters over their ranges of expected uncertainty, for a total of 180 runs per image.The resulting distribution of results for each image allowed us to determine the most likely values for the atmospheric tilt and estimate the possible error.
Both the center-pointing uncertainties and the unconstrained altitude of the NSA boundary feature contributed to possible mapping errors.While center-pointing inaccuracies were reduced by fitting to the spherical upper haze layer (as described in Section 2.1 and in Roman et al. 2009), uncertainty of up to a few pixels remained in some images, particularly those at higher phase angles where part of the disk was not illuminated.The uncertain altitude of the features added further uncertainty to the mapping process.Both error sources become more pronounced near the planet limb where the curvature of the planet exacerbates any mapping errors, yet the points near the limb are extremely valuable for constraining the NSA boundary orientation.
The images also exhibited significant noise and artifacts, specifically in the form of both vertical and horizontal banding, which were strong enough to interfere with the NSA boundary identification process.As discussed above in Section 2.1, these artifacts persisted after the initial CISSCAL calibration by the Cassini team (Knowles et al. 2020), so median filters and boxcar means were used to reduce them further.We aimed to suppress the noise and artifacts as much as possible while minimizing the loss of real planetary features, but this is a lingering source of minor errors that we have accounted for.
We also noticed that throughout the mission, subtle "bumps" in the NSA boundary could be observed, warping parts of the boundary outside of the assumed flat ring/plane.In addition, the width and gradient of the boundary varied spatially and temporally, making it more complicated to consistently determine the center point of the boundary.The analysis methods described in Section 2.2 were developed largely through trial and error to ensure a consistent and robust pipeline that can handle such variations.However, some uncertainty is still assumed, and is accounted for by injecting artificial Gaussian noise into the boundary points in the various Monte Carlo runs.
Lastly, there is some uncertainty in the fitting functions used throughout the pipeline.However, these errors are more easily calculated, and thus were propagated through the analysis and incorporated into the final uncertainty calculations.

Results
We analyzed a few hundred Cassini ISS images distributed from 2004 to 2017 to determine the orientation of the NSA boundary and thus the tilt axis of the atmospheric circulation more generally.At the 889 nm wavelength used here, the northern hemisphere appeared brighter at the start of the data set and reversed so that the southern hemisphere was brighter at the end.We found that the NSA boundary was tilted with respect to the solid body rotation axis by roughly 4°, fluctuating by a couple degrees apparently with the seasons.The azimuth angle of the tilt axis was roughly fixed in inertial space, its offset from the subsolar longitude moving westward throughout the year.When plotted in an inertial reference frame, there is some evidence of oscillation on a period much shorter than the Titan year.The NSA boundary had a mean latitude around 9°S (when the northern hemisphere was brighter) before migrating southward and eventually being replaced by a new inverted boundary (at the time of this writing the southern hemisphere is brighter) in the northern hemisphere.
The tilt amplitude of the NSA boundary with respect to the solid body north pole versus time is shown in the upper panel of Figure 4.The tilt amplitude started off around 5°and decreased to around 3°before the NSA reversal around 2014.The amplitude then seemed to increase before the original boundary faded away.The new boundary in the northern hemisphere had approximately a 3°amplitude.Uncertainty in the measurements was generally 1°-2°for the first half of the mission, but increased around the boundary transition period (2012 and later) to 2°-3°for most measurements.
The tilt amplitude started with around a 70°offset from the subsolar longitude in 2004.Throughout the course of the Cassini mission, the azimuth offset moved westward at nearly the same rate as the progression of the solar longitude (L s ), meaning the azimuth remains approximately fixed in inertial space (see the bottom panel in Figure 4).Though generally following a linear trend with L s , the results hint at some shorterterm fluctuations, suggesting the atmospheric tilt axis may oscillate around a fixed point on timescales of less than a Titan year.Data at the end of the mission for the new southward brightness boundary had high uncertainties for the azimuth calculations due to smaller tilt amplitudes (in some cases approaching 180°), and it is unclear if the trend continues or if there is actually a 180°change.Similar ambiguity in azimuth measurements from this time period appears in Kutsop et al. (2022).However, the small tilt amplitude here means that only a slight perturbation of the tilt axis can drastically change the azimuth, so it is unlikely that the large variation in azimuth calculations indicates any real drastic phenomenon is occurring.
The mean latitude of the NSA boundary is shown in the middle panel of Figure 4. Another figure, Figure 5, displays the latitudinal brightness over time to show more comprehensively the changes in brightness across the whole visible disk.We found that the mean latitude remained fairly steady at 9°S until the boundary began to migrate south starting around 2011 as the northern hemisphere began to darken.This original boundary faded away in 2016 at around 20°S when the southern hemisphere took over as the brighter hemisphere.The mean latitude measurements are generally expected to be accurate to within a few degrees.Meanwhile, a new inverted boundary arose in the northern hemisphere as a dark cap extended southward and eventually encompassed most of the northern hemisphere, completing the NSA reversal.Figure 6 shows example images from before, during, and after the NSA reversal.Numerous bands and brightness transitions can be seen in the middle image before settling back into the north-south asymmetry pattern.

Discussion
The purpose of this study was to analyze the tilt of the rotation axis of Titan's atmosphere, using the highest-  resolution imagery available (i.e., Cassini ISS), to gain further insight into the atmospheric dynamics of Titan's super-rotating atmosphere.Additionally, by characterizing the location and orientation of major banding features in Titan's atmosphere throughout the Cassini mission, we can determine the seasonal variations of the features themselves and the atmosphere as a whole.This enables a deeper understanding of the underlying dynamics and provides further constraints for Titan atmospheric models (both theoretical and numerical).
A seasonal dependence for the albedo asymmetry was already known, as past observations revealed the reversal in hemispheric brightness (e.g., Caldwell et al. 1992).However, the seasonal dependence of the atmospheric tilt the asymmetry reveals was less well constrained until a recent paper by Kutsop et al. (2022).Our results are complementary to those of Kutsop et al. (2022) in that we use an entirely different data set (Cassini ISS versus VIMS) and independent methods.Our data set included additional temporal coverage of the NSA, particularly beyond 2012 covering the NSA reversal (measurements from this time period in Kutsop et al. 2022 are primarily from the north polar annulus rather than the equatorial annulus/NSA boundary), and the spatial resolution of ISS images generally exceeds that available in VIMS.Along with the additional steps taken here during image processing and analysis to robustly and precisely navigate features and their tilts, our analysis may prove to be the most precise assessment of the atmospheric tilt and its seasonal behavior.Overall, our results are generally in accord with those of Kutsop et al. (2022), and show that the azimuthal angle of the atmospheric rotation axis is roughly fixed in inertial space (perhaps with some slight fluctuations around this angle on timescales shorter than seasonal), and there is a seasonal fluctuation in the amplitude of the tilt angle.
In the top panel of Figure 4, we see the amplitude decrease from around 5°to around 3°-4°near the time the NSA reverses.This is in agreement with the results from the Tokano (2010) GCM, which showed a similar decrease in amplitude before a more rapid increase back to the previous amplitude when the asymmetry reversed.This increase seems to be exhibited by the original northward brightness gradient boundary, but not by the new boundary that arises.However, the measurements of the new boundary are quite noisy and uncertain due to the more difficult viewing geometries and somewhat fuzzier boundary toward the end of the mission, so such a trend would not necessarily fall outside our error bars.It is also possible that the more drastic circulation changes occurring around the transition period perturbed the boundary to some extent.
The bottom panel of Figure 4 shows the azimuth angle offset from the subsolar longitude.Throughout the mission, we see this angle increase roughly along with the solar longitude (L s ), indicating the atmosphere was tilted at nearly the same angle in inertial space throughout the mission.This result is also generally in agreement with previous results such as Achterberg et al. (2011) and Kutsop et al. (2022), which similarly found that the azimuth angle is roughly fixed in inertial space, and Sharkey et al. (2020), which found that assuming the tilt was fixed in inertial space resulted in less zonal asymmetry in thermal emissions.Kutsop et al. (2022) also noted slight oscillations of the tilt axis and offered additional analysis and interpretations of some fine-grained trends.We did not deem such a fine estimation of oscillations feasible with the noise and uncertainty in our data, so will only compare the coarsest trends.However, while we cannot comment specifically on any finer trends, we do see what appear to be slight, perhaps sinusoidal fluctuations in the azimuth offset; note the data points in the bottom panel seemingly alternating between falling above and below a fixed azimuth offset, represented by the dashed line.Lastly, while Roman et al. (2009) and Achterberg et al. (2008a) did not have enough temporal coverage to determine any seasonal dependence, the average values they calculated for amplitude and azimuth angle offset support these results, as well.
Our observation of the NSA reversal is also in agreement with past observations (e.g., from Voyager and Hubble; Caldwell et al. 1992), where the same seasonal dependence was noted.New to this study is the more comprehensive analysis of what happens in between, at least from the point of view of the 889 nm methane band.Here, we saw how the original NSA boundary seemed to migrate south before fading away, while the new boundary gradually strengthened and moved equatorward from the north.Two major boundaries coexisted briefly, along with a number of other more subtle bands.
At times throughout the mission, the NSA boundary exhibited features that may be attributable to some weather phenomena affecting the region.Such features were more clearly evident when viewing the gradient of the image, and consisted of a partially split NSA boundary (i.e., two discrete latitudinal brightness changes) and bumps or waves where the boundary deviated significantly from a flat ring.While such features complicated our primary analysis somewhat as they deviated from our expected model of the boundary, they do suggest more interesting and complicated dynamics are occurring at smaller scales and may offer an opportunity for future deeper analysis.The phenomena occurred both during and outside of the NSA reversal transition period, but during the nontransition times such phenomena were much more isolated and distinguishable from their more predictable surroundings.
Cassini provided data for around half of a Saturn year.While this gave us some of the best temporal coverage applied to this problem so far and enabled some analysis of seasonal trends, it would be useful to extend the coverage further to cover the entire Titan year and beyond.This could potentially be done using existing data from Voyager or Hubble, or perhaps with future observations.The lack of usable images around 2014 was also particularly unfortunate since this was a key period during the NSA reversal.A set of data with more complete coverage of the transition period could provide a more complete look at the changes the atmosphere experiences.
Additionally, only images from the 889 nm filter were used, and the NSA boundary was assumed to have arisen in a certain region of the atmosphere.The NSA is present in a range of wavelengths from IR to UV.We know different wavelengths reveal different altitudes in Titan's atmosphere, and the NSA itself is wavelength dependent, with the bright and dark hemispheres inverting at wavelengths below 440 nm.Analyzing images taken through different filters would provide additional data points and probe different regions/features of the atmosphere.Finally, there are additional atmospheric features (i.e., the polar annuli discussed in Kutsop et al. 2022) that were not included in this work but are visible in many Cassini ISS images that could also be analyzed to offer more data points and compare atmospheric behavior at differing latitudes.

Conclusions
To accurately describe and model the dynamics and circulation of Titan's atmosphere, observational data collected over seasonal (or longer) timescales is an essential source of insight and model constraints.Here, we analyzed Cassini ISS data from approximately half of a Titan year to reveal both qualitative and quantitative seasonal trends in the appearance and orientation of Titan's atmosphere.We extended and improved upon previous work-first through the use of the now-complete Cassini data set, and additionally through more robust image-processing and feature-detection techniques.We developed an empirical photometric model to describe the unusual limb brightening exhibited by Titan in the 889 nm filter, which increased the amount of usable data near the limb.The NSA boundary was identified in a more robust way using image smoothing, directional gradients, and connecting edge candidates to capture the overall NSA boundary regardless of minor variations in its characteristics and without false detections of smaller features.
We found that the vector normal to the NSA boundary, and thus the rotation axis of Titan's atmospheric mean zonal flow, is offset from the poles of the solid body by 4°on average.The azimuth of the atmosphere's axis is roughly fixed in inertial space, its offset from the subsolar longitude increasing at about the same rate as L s , or alongside the progression of Titan's orbit and seasons.

Figure 1 .
Figure 1.The empirically determined photometric model for Titan in the 889 nm methane band filter, with I/F shown as a function of incidence and emission angles.

Figure 2 .
Figure 2. Top: image W1536388236_1 after reduction of noise and banding artifacts in units of I/F.Only the pixels considered useful to the analysis are shown, i.e., parts of the disk with emission angle <85°and incidence angle <87°.Middle: the surface brightness model (in I/F) produced when the empirical photometric model is interpolated for the incidence and emission angles of each pixel in the image.Bottom: the resulting flattened image after the calibrated image is divided by the brightness model.

Figure 3 .
Figure 3. Top: image N1496574260_1 after limb-brightening reduction.Middle: the brightness gradient in the northward direction (north is to the lower right in this image).The locations of gradient extrema, as determined from polynomial fits, are indicated on the images, where black points indicate extrema determined to belong to a strong edge (used for subsequent analysis) and white points indicate extrema belonging to weaker edges (not used).Bottom: a few slices from the above image gradient, demonstrating the degree-10 polynomial fitting.

Figure 5 .
Figure 5.The relative brightness by latitude after limb-brightening reduction is shown throughout the Cassini mission.Brightness data were interpolated through the time gaps between images.This figure shows the NSA boundary move south and eventually fade, while a new reversed boundary arises slightly to the north.Due to the image geometries used in this study, the data near the poles were often low quality or unavailable throughout the mission.

Figure 4 .
Figure 4.The location and orientation of the NSA boundary changes over time and shows a seasonal dependence.Results are shown for both the initial major NSA boundary where brightness increases in the northward direction (blue) and for the new reversed boundary that arises near the end of the Cassini mission (red).The top panel shows the amplitude of the tilt with respect to the solid body north pole, the middle panel shows the mean latitude of the NSA boundary, and the bottom panel shows the offset of the azimuth from the subsolar longitude.The dashed line in the bottom panel represents a best fit to the azimuth offset if it was fixed in inertial space.Deviations from this line in the data may indicate that seasonal or other atmospheric phenomena affect the tilt azimuth.

Figure 6 .
Figure 6.Sample images from before (N1567440863_1), during (N1786775340_1), and after (N1852202278_1) the NSA reversal show how the NSA evolves over time.The north angle for each image is oriented toward the top of the figure (though the actual north pole may point somewhat into/out of the page), and the tilt in the NSA boundary is clearly visible in the top panel.
(This table is available in its entirety in machine-readable form.)