Constraining the Temporal Variability of Neutral Winds in Saturn's Low‐Latitude Ionosphere Using Magnetic Field Measurements

The Cassini spacecraft completed 22 orbits around Saturn known as the “Grand Finale” over a 5 months interval, during which time the spacecraft traversed the previously unexplored region between Saturn and its equatorial rings near periapsis. The magnetic field observations reveal the presence of temporally variable low‐latitude field‐aligned currents which are thought to be driven by velocity shears in the neutral zonal winds at magnetically conjugate thermospheric latitudes. We consider atmospheric waves as a plausible driver of temporal variability in the low‐latitude thermosphere, and empirically constrain the region in which they perturb the zonal flows to be between ±25° latitude. By investigating an extensive range of hypothetical wind profiles, we present and analyze a timeseries of the modeled velocity shears in thermospheric zonal flows, with direct comparisons to empirically inferred angular velocity shears from the Bϕ observations. We determine the maximum temporal variability in the peak neutral zonal winds over the Grand Finale interval to be ∼350 m/s assuming steady‐state ionospheric Pedersen conductances. We further show that the ionospheric currents measured must be in steady‐state on ∼10 min timescales, and axisymmetric over ∼2 h of local time in the near‐equatorial ionosphere. Our study illustrates the potential to use of magnetospheric datasets to constrain atmospheric variability in the thermosphere region.

• Magnetic field data reveal lowlatitude field-aligned currents at Saturn and are used to infer variability of neutral flows in ionosphere • We test many modeled wind profiles and find the peak neutral flows are likely variable by ∼350 m/s near the equator over a 5 months interval • We present empirically inferred timeseries of inter-hemispheric angular velocity shears in neutral flows for the low-latitude ionosphere axisymmetric internally generated planetary field (Burton et al., 2010;Cao et al., 2011Cao et al., , 2019, and no significant source of magnetospheric plasma in the inner magnetosphere (within the orbit of Enceladus, < 4R S ).
A plausible mechanism for generating these field-aligned currents was first suggested by Khurana et al. (2018), who proposed that if neutral zonal winds were present in Saturn's thermosphere and asymmetric about the magnetic equator, then a velocity shear in the neutral zonal flows between the two magnetic hemispheres would drive ionospheric currents which close via the magnetosphere as field-aligned currents. By projecting an averaged 1-bar tropospheric zonal wind profile inferred from cloud-tracking measurements at Saturn (Read et al., 2009) into the thermosphere, Khurana et al. (2018) estimated a ∼20 nT B ϕ signature on the magnetic field line connected to the inner-edge of the D-ring, which is consistent with the typical peak in B ϕ observations within this region.
Field-aligned currents can also be generated by dispersive Alfvén waves provided that the wavelength of the perturbation is comparable to the local kinetic length scales in the magnetosphere. The ion densities in the inner magnetosphere were observed to be of the order of 10 2 -10 3 cm −3 during the Cassini Grand Finale orbits (Morooka et al., 2009), which correspond to ion inertial lengths of 0.38R S and 0.03R S respectively. The central B ϕ feature is observed over a length scale of ∼0.6R S which is much larger than the upper limit of the ion inertial length estimated in the inner magnetosphere, therefore it is unlikely that the large scale B ϕ signatures observed on Revs 271-292 are generated by dispersive Alfvén waves. However, it is plausible that the small scale fluctuations in the measurements could be generated by this mechanism Southwood et al. (2020). Provan et al. (2019b) find that that the orbit to orbit variability in the large scale B ϕ structure is overall uncorrelated with the residual poloidal magnetic field components, and Hunt et al. (2019) show that the B ϕ observations are generally symmetric across the magnetic equator on timescales of <30 min. These results are overall inconsistent with what one might expect if the observed B ϕ feature was generated by the spacecraft encountering standing Alfvén waves in the inner magnetosphere, where the Alfvén travel time is on the order of a few seconds. Therefore, the inter-hemispheric ionospheric current system proposed by Khurana et al. (2018) is the most probable source of the large scale B ϕ feature observed on the Grand Finale orbits.
A more rigorous consideration of the ionospheric currents was later presented by Vriesema et al. (2019), who used a fully resolved ionosphere and five different modeled wind profiles to evaluate a temporally static, but spatially variable B ϕ profile. One of their five wind profiles, loosely based on the averaged zonal flow profile evaluated by Read et al. (2009), showed good overall agreement with the typical B ϕ profile observed during the Grand Finale orbits. If inter-hemispheric atmospheric coupling is indeed present at Saturn, then the B ϕ observations from the Grand Finale orbits provide a unique opportunity to study the temporal evolution of netural zonal flows in Saturn's thermosphere on week-long timescales over a ∼5 months interval.
The discussions by Khurana et al. (2018) and Vriesema et al. (2019) conclude that there must be neutral zonal flows present at thermospheric altitudes to explain the typical B ϕ signature observed during the Grand Finale orbits. This discussion is supported by Provan et al. (2019b) who show that the sense of measured B ϕ in the magnetosphere is directly proportional to the magnitude and sense of the inter-hemispheric angular velocity neutral wind shears. The typical observation from Revs 271-292 shows a positive B ϕ signature in the intra D-ring region, however, a reversal in the sense of B ϕ is also observed on three Revs, which would require a significant change in the thermospheric neutral winds according to the derivations from Provan et al. (2019b), or significant perturbations from unknown magnetospheric sources.
While the near-equatorial 1-bar winds at Saturn have been observed to be more or less in steady state between the Voyager era and the Cassini Mission (García-Melendo et al., 2011;Read et al., 2009), Müller-Wodarg et al. (2019) and Brown et al. (2020) have recently presented in-situ observations of atmospheric waves in Saturn's thermosphere from measurements taken by the Ion Neutral Mass Spectrometer (INMS) and the Ultraviolet Imaging Spectrograph (UVIS) instruments on the Cassini spacecraft during the Grand Finale. Atmospheric waves are a well-known driver of variability in the upper atmospheres of Earth and Mars, and Müller-Wodarg et al. (2019) suggest that atmospheric waves could assist in the vertical transport of the 1-bar wind structure to higher altitudes, where they could also introduce temporal variability in the zonal flows.
In this study, we create an extensive selection of hypothetical zonal wind profiles and use them to model B ϕ along the spacecraft trajectory of the Grand Finale orbits using the theoretical approach from Provan AGIWAL ET AL.  (2019b). Through comparison of the modeled B ϕ profiles with the B ϕ observations, we infer constraints on the temporal variability of neutral zonal winds in the low-latitude ionosphere at Saturn, assuming atmospheric waves to be the most likely perturbation source.

The Data
The B ϕ observations from Revs 271-292 at 1 s resolution are presented in Figure 1a, except for Rev 277 which is omitted due to a large data gap at periapsis. A running median of the observations is shown by a solid black line, which represents the "typical" B ϕ profile modeled by Khurana et al. (2018) and Vriesema et al. (2019). The data are shown as a function of the magnetically mapped ionospheric latitude θ i which indicates the ionospheric latitude (in the spacecraft's local hemisphere) that is magnetically connected to the position of the spacecraft in the magnetosphere, and should not be confused with the latitudinal position AGIWAL ET AL.
10.1029/2020JE006578 3 of 18 Measurements taken in the northern kronographic hemisphere (z > 0 R S ) are shown relative to their northern footpoint θ i > 0° and measurements from the southern hemisphere (z < 0 R S ) are shown relative to their southern footpoints θ i < 0°. The conjugate hemisphere footpoints θ c are shown along the x-axis. The innermost field line traversed by the spacecraft is indicated by dashed gray lines at (14°, −4°), thus the gap between the data shows magnetic field lines that the spacecraft did not traverse. The median B ϕ observation is shown by a solid black line. The detrended B ϕ observations are then shown in panel (b) where a fourth order polynomial is subtracted from each Rev, to isolate the intra D-ring signature Provan et al., 2019b). Panel (c) shows an averaged angular velocity Ω profile of the 1-bar zonal winds at Saturn relative to an assumed planetary rotation rate of ∼10 h 34 min (Read et al., 2009) as a function of latitude. The ring-plane mapping ρ eq of θ i is shown on all panels by shaded gray regions for the A-D rings.
of the spacecraft itself. The magnetic mappings are determined using a degree-3 internal field model for Saturn , which is sufficient in this region.
The spacecraft trajectory was inbound from the northern hemisphere, and outbound in the southern hemisphere. For ease of discussion in future sections, we shall refer to measurements from the northern kronographic hemisphere (shown with respect to the northern ionospheric latitude θ i > 0°) as the "inbound" part of the orbit, and the measurements from southern kronographic hemisphere (θ i < 0°) as the "outbound" part of the orbit. However, this is not strictly true as orbit periapsis occurred at approximately −6° spacecraft latitude, which is in the southern kronographic hemisphere and magnetically connected to θ i between −9° and −15° for Revs 271-292 (the spread in θ i values is due to slight orbit by orbit variability in the spacecraft trajectory).
Shaded gray panels in Figure 1 indicate the latitudes which are magnetically connected to the main rings in the equatorial plane and can be used as a visual guide to locate magnetically conjugate latitudes. The magnetically conjugate latitude in the opposite hemisphere to the spacecraft θ c , and the perpendicular distance between Saturn's spin axis and the position of a magnetic field line in the equatorial plane ρ eq are also given along the x-axis. The innermost magnetic field line traversed by the spacecraft maps to (14°, −4°) and is indicated on Figure 1 by vertical dashed gray lines, therefore the "gap" in the data between these lines represents field lines that were not traversed by the spacecraft. It should be noted that due to a slight northward offset of Saturn's magnetic axis by ∼0.05R S from the rotational axis (Burton et al., 2010;Dougherty et al., 2018), magnetically conjugate latitudes are asymmetric about the rotation equator.
Strong spatial gradients in B ϕ indicate the flow of field-aligned currents in the magnetosphere, which constrain the observed field-aligned current system in Figure 1a to lie inside of 1.24R S in the equatorial plane (D ring outer-edge), which corresponds to latitudes between + 30° and − 25° in Saturn's ionosphere. While significant gradients in B ϕ are not measured outside this region (with the exception of the magnetic field lines which map to the B-ring during the outbound part of each Rev) there is a ∼ 20 nT spread in the observed B ϕ values in both hemispheres associated with the A-C rings. Provan et al. (2019b) show that this spread is in part due to the planetary period oscillations (PPO), and thus subtract a baseline azimuthal field from each Rev by fitting a fourth order polynomial to the B ϕ observations outside the D ring outer-boundary in order to focus their analysis on the intra D-ring signature. We employ their baseline subtraction method on the observations presented in Figure 1a, the results of which are presented in Figure 1b. It can be seen that the subtraction acts to reduce B ϕ to almost zero over the A-C rings, and causes minor adjustments to the intra D-ring signature on the order of ∼5 nT which is consistent with the PPO amplitudes of ∼5−8 nT estimated in this region by Provan et al. (2019a). Figure 1c shows the averaged profile of the angular velocity of the tropospheric zonal winds Ω relative to planetary rotation evaluated by Read et al. (2009) at the 1-bar pressure level from cloud-tracking measurements. The assumed planetary rotation period for the profile shown is 10 h 34 min, which is in close agreement with the most recent rotation period value of 10 h 33 min and 38 s estimated by Mankovich et al. (2019) using Saturn ring seismology. The latitudinal profile of the azimuthal velocity of the tropospheric zonal winds U are presented in Supplementary Information for comparison.
The magnitude and sense of the inter-hemispheric angular velocity shear on magnetically conjugate footpoints presented in Figure 1c is near-zero on magnetic field lines in the region between the A-C rings (with the exception of the magnetic field lines mapping to the middle of the B-ring, where a nonzero angular velocity shear is present, and will be discussed more later), and positive on magnetic field lines in the intra D-ring region. These trends are overall consistent with the near-zero B ϕ observations on A-C ring magnetic field lines, and positive median B ϕ observations on the intra D-ring region magnetic field lines. Vriesema et al. (2019) show that the typical B ϕ observations are best modeled by a zonal wind profile based on the Read et al. (2009) profile, but with the ionospheric values scaled to half those in the equatorial troposphere. Their study employs modeled Pedersen conductivities from the "Saturn Thermosphere Ionosphere Model" (STIM; Müller-Wodarg et al., 2006, 2012Moore et al., 2004Moore et al., , 2010, which has previously shown good agreement with observed Pedersen conductivities in Saturn's ionosphere, and will be discussed a little more in Section 2.1 as it is also used in our modeling. In the absence of any direct observations of thermospheric neutral winds at Saturn, we use the Ω profile presented in Figure 1c scaled by half as the baseline hypothetical wind profile in our study. The angular velocity shears from the Read et al. (2009) profile show an overall positive correlation with the median B ϕ observations, and a half-scaled version of the Read et al. (2009) profile was shown to be a good fit for certain Grand Finale orbits by Vriesema et al. (2019). This baseline wind profile gives us a starting point from which we then create the aforementioned range of hypothetical zonal wind profiles that are used to investigate the variability in the B ϕ observations, in an attempt to constrain the true underlying thermospheric neutral wind profiles. Provan et al. (2019b) show that in the presence of neutral atmospheric flow shear between the two ends of a given field line, the horizontal ionospheric meridional current per radian of azimuth I m , positive northward and equal at the two ends of the field line is given by

Modeling the Temporal Variability in the Azimuthal Magnetic Field
where, Ω N,S are the neutral flow angular velocities at the two ends of the field line, Σ PN,S are the Pedersen conductances (or the height-integrated Pedersen conductivities) at the two ends of the field line, ρ iN,S are the perpendicular distances of the footpoints of the field lines in the ionosphere from the magnetic/spin axes, B iN,S are the ionospheric field strengths, and |b iN,S | are the magnitudes of the ionospheric field strength normal to the current layer. The geometric factors B iN,S /|b iN,S | account for the tilt of the magnetic field lines relative to the normal to the layer, and arise from the requirement that the current flows entirely horizontally in the ionospheric layer (Cowley & Bunce, 2003). Provan et al. (2019b) assume that the parameters defined in Equation 1 are in steady state and azimuthally symmetric in the ionosphere, but differ between the northern and southern hemispheres. Their model is equivalent to the fully resolved ionosphere model presented by Vriesema et al. (2019) in a thin sheet approximation of the ionosphere (see Appendix A). The field-aligned current flowing above the ionospheric layer is then given by the variation of I m with latitude. From the integral form of Ampère's law, the azimuthal magnetic field generated by the current system above the ionosphere on each field is given by where ρ is the perpendicular distance of the point of measurement from the magnetic/spin axis, and μ 0 is the permeability of free space.

Steady-State, Axisymmetric, Modeling Parameters
The magnetic field and geometric parameters in Equation 1: B iN,S , b iN,S , and ρ iN,S are evaluated using a degree-3 internal field model  at ∼1,000 km above the surface of Saturn. They are assumed to be in steady state for the duration of the Grand Finale orbits since any variability in B iN,S , which is on the order of ∼ 20,000 nT, due to the variability in the B ϕ component, which is on the order of ∼20 nT, is negligible.
In the absence of charged dust, Pedersen conductivities scale linearly with the ion densities, which are equal to the electron densities in a quasi-neutral ionosphere. Radio occultations of the mid-and low-latitude ionosphere at Saturn have shown temporally variable electron densities on week-long timescales (Hadid et al., 2018b; Kliore et al., 2009;Nagy et al., 2006), and recent in-situ measurements from the Cassini Grand Finale in the near-equatorial upper atmosphere (between − 15° and 5°) have shown the Pedersen conductivities at ∼1,000 km above the peak conducting layer of the ionosphere to be variable by up to an order of magnitude between the ∼6.5 days periapsis passes . However, there are no observations of the magnitudes or extent of variability in Σ P in the peak conducting layer of the low-latitude ionosphere during Revs 271-292.
Equations 1 and 2 show that the magnitude of I m can be scaled by varying either the Σ P parameter, or the angular velocity shear ΔΩ = Ω S −Ω N . However, the sense of I m , and therefore B ϕ (which is both positive and negative in observations), is dependent only on the sense of ΔΩ. Since both parameters are able to scale the magnitude of I m , we assume Σ P to be variable in latitude but otherwise in steady state for this study. This assumption allows us to investigate an upper limit on the extent of variability in the angular velocity of the zonal flows, however we acknowledge that in reality, the ionospheric currents will also be sensitive to the variability in the Pedersen conductances. Hunt et al. (2019) calculate the intra D-ring field-aligned currents from the B ϕ measurements during the Grand Finale orbits, and show that the large scale structure of the currents appear to be in steady state (on timescales < 30 min) between the inbound and outbound legs of each Rev, but are variable from orbit to orbit. They also find that the small scale structure in the currents does not show the same level of conjugacy between the inbound and outbound legs of each Rev. Their findings overall support the presence of a large scale current system which is mostly in steady state, and thus we assume that both the Pedersen conductances and the neutral flow angular velocities are also in steady state during any given Cassini pass, but that the latter can be variable from orbit to orbit. Due to our assumptions, we are not expecting to model the small-scale structure observed on the B ϕ profiles, as understanding those fluctuations and their source is beyond the scope of this study.

Ionospheric Pedersen Conductances From STIM
The STIM , 2012Moore et al., 2004Moore et al., , 2010) is a global time-dependent general circulation model of Saturn's coupled thermosphere and ionosphere. The ionospheric plasma densities and conductivities are calculated by considering solar and particle ionization sources and ionospheric chemical reactions, and have shown good agreement with observations from southern hemisphere summer at Saturn in the low-and mid-latitude regions (Moore et al., 2006(Moore et al., , 2010. As there are no observations of Σ P in the region of interest during Revs 271-292, we use modeled Pedersen conductances from STIM that are evaluated under northern summer conditions at Saturn. The trajectory of the spacecraft is variable by ∼ 5 h in local time (LT) between inbound (pre-noon sector) and outbound (post-noon sector), with periapsis near local noon, and a minor drift of ∼ 0.1 h in LT between each Rev. However, the derivation of Equations 1 and 2 from Provan et al. (2019b) require Σ P to be azimuthally symmetric. Therefore, we choose the conductivities evaluated at local noon in STIM to use when evaluating our B ϕ models, since orbit periapsis occurs between 13 h (Rev 271) and 11 h (Rev 292) in LT for the Grand Finale orbits.
We calculate the Pedersen conductance, Σ P by integrating the Pedersen conductivities over an ionospheric layer of thickness ∼350 km centered around the peak conducting layer (∼1,000 km above the 1-bar pressure level) at 1° latitude resolution. The modeled Σ P values are presented in Figure 2a, where circle markers show Pedersen conductances evaluated at local noon, and vertical bars show the range of Σ P values evaluated at the LT position of the spacecraft between Revs 271-292 as a function of θ i . The blue and orange colors indicate the Σ P values evaluated at the northern and southern ends of each field line, respectively. The dip in the southern ionosphere Σ P values between −15 and −35° in latitude is due to a shadow cast by the rings on the ionosphere (Hadid et al., 2018a;Moore et al., 2004).
The local noon Σ P values are largely consistent with the Σ P modeled in the intra D-ring region at the LT position of the spacecraft, with some minor differences in the southern hemisphere Σ P for the inbound spacecraft trajectories which were in the pre-noon sector. The conductances are variable by ∼1.5S between Revs 271 and 292, which on the D-ring field lines corresponds to a maximum uncertainty in the modeled values of ∼4 nT which is at the limit of the B ϕ resolution during the periapsis passes of the Grand Finale orbits, and thus negligible. Shebanits et al. (2020) determined the Pedersen conductivity close to the peak conducting layer between − 15° and 10° in latitude to be 10 −4 − 10 −5 Sm −1 from in situ measurements taken during Revs 288-293 (with Rev 293 being Cassini's final orbit), which integrated over an ionospheric layer of ∼350 km corresponds to a Pedersen conductance of 3.5 − 35S. These conductances are significantly larger than previous observations have suggested, and this is primarily due to the presence of charged dust in the near-equatorial (±5°) ionosphere, with some dust still being present at slightly higher (between 10° and −15°) latitudes. The modeled Σ P values presented in Figure 2a are ∼5S in the region between −5° and −15° in latitude, which is consistent with the lower limit of the observed range but could be variable by up to an order of magnitude on certain orbits if the in-falling dust from the rings is able to survive all the way down to the peak conducting region.
Since the modeled Σ P profiles from STIM fall within the observed range of Σ P values by Shebanits et al. (2020), we find it suitable to use the modeled Σ P profile in this work given our assumption of steadystate Pedersen conductances. However, we keep in mind that any ΔΩ values constrained from our study at the near-equatorial latitudes in the southern ionosphere are subject to being scaled by up to a factor of 10 due to ionospheric variability in the Pedersen conductances. Although, we might expect the extent of AGIWAL ET AL.
10.1029/2020JE006578 7 of 18 variability to be overall smaller as the presence of charged dust is primarily constrained to ±5° in latitude, which corresponds to magnetic field lines significantly below the spacecraft periapsis altitude (as shown in Figure 1). It is more difficult to constrain the extent of the variability in Σ P for the higher latitudes in our study, as we do not have observations of ionospheric densities or conductivities in this region at local noon. However, Moore et al. (2010) show Σ P to vary between 1S and 10S in the low-to-mid latitude ionosphere from dawn and dusk radio occultations of Saturn's ionosphere (Kliore et al., 2009), thus we might expect the local noon ionospheric conductivities in the low-mid-latitude regions to be temporally variable as well.

Creating Hypothetical Neutral Angular Velocity Wind Profiles
The baseline angular velocity Ω profile mentioned in Section 1.1 is presented in Figure 2b as a solid black line. The 14° and − 4° latitudes are treated as boundaries at which we scale the baseline Ω profile, and we do not attempt to model in detail of how the zonal flows may behave within these boundaries since the dynamics of this region cannot be inferred from the B ϕ observations. We explore 92 hypothetical angular velocity profiles in this study, which are comprised of profiles that are variable from the baseline Ω profile in the southern hemisphere within the yellow envelope shown in Figure 2b, and are then combined with one of two configurations in the northern hemisphere, henceforth N1 or N2, shown in Figure 2 by a solid black line and dashed orange line respectively. The N1 profile is simply the baseline profile, and the N2 profile is mirrored from N1 such that the zonal flows are westward (Ω < 0 rad s −1 ) instead of eastward (Ω > 0 rad s −1 ).
While atmospheric waves have been identified observationally in Saturn's thermosphere (Brown et al., 2020;Müller-Wodarg et al., 2019), their latitudinal distribution is presently unknown. It can be seen from Figure 1b that the temporal variability in B ϕ observations is largely confined between magnetically conjugate latitudes of 33° and − 25° which are magnetically connected to the outer-edge of the D-ring in the equatorial plane. The variability being confined to low latitudes could suggest the presence of an equatorially confined wave in Saturn's thermosphere, such as a planetary Rossby wave, which arises due to the rotation of a planet. Planetary waves can be symmetric or antisymmetric about the rotational equator, however the region in which they are present is expected to be confined to rotationally conjugate latitudes; therefore, we constrain the region of variability in our study to be between ± 25° (indicated by dot-dashed blue lines in Figure 2b).
Because of the offset of the magnetic equator and the assumption of rotationally conjugate boundaries for wave activity, any putative region of wave activity will have more of its area occupied by the southern footpoints of magnetic flux tubes. Thus, we vary the baseline Ω profile more significantly in the southern hemisphere than the northern hemisphere to create our Ω m profiles, to ensure that the different Ω m profiles do not simply scale the magnitude of I m , but also allow for variable gradients in I m evaluated along the spacecraft trajectory.
We introduce variability in the southern hemisphere by scaling the baseline Ω at the − 4° boundary by a scale factor s, between s = −1 and s = 2 at 0.2 intervals, while creating a fixed point at − 25°, − 20°, − 15°, and − 10° respectively. We then use Gaussians with peaks at the scaled values and tails at the fixed points to replace the baseline profile in the nominal region, and thus create our hypothetical profiles. For the fixed points between − 10° and − 20°, the s > 1 scale factors would not be significantly different to the profiles evaluated with the fixed point at − 25°, therefore those profiles are omitted. The extent of the scaling is shown by the yellow envelope in Figure 2b, and the Ω m corresponding to s = −1 and s = 2 for the fixed point at − 25° are shown by dashed black lines. The Ω m corresponding to s = −1 for the remaining fixed points are shown by dashed purple lines, which introduce turning points in the wind profile at these fixed points, allowing for a more variable selection of hypothetical wind profiles.
Using Equations 1 and 2, a modeled B ϕ field is evaluated for each of the hypothetical Ω m profiles. To demonstrate that the Ω m profiles sufficiently capture the extent of variability observed in B ϕ during the Grand Finale orbits, in Figure 2c, we present the modeled B ϕ output for the s = −1 and s = 2 hypothetical Ω m profiles with fixed points at − 25° combined with the N1 (blue dashed lines) and N2 (orange dashes lines) profiles respectively. The darker shaded region enclosed between N1(s = 2) and N2(s = −1) lines shows that the more symmetric models, or models with a very slight asymmetry across the equatorial boundary, show better overall agreement with the data. The highly asymmetric profiles such as N1(s = −1) and N2(s = 2) appear to largely over-estimate the zonal flows near the inner equatorial boundary. The near-zero B ϕ output associated with the A-and C-rings, where a nonzero velocity component exists in Ω but the evaluated shears are near-zero, highlights the key dependence of the model on angular velocity shears.

Example Cases: Revs 271, 278, 289, and 292
The full selection of modeled B ϕ profiles were compared to the observed B ϕ for each Rev by a simple rootmean-square-error (RMSE) analysis in the intra D-ring region. For all 21 Grand Finale orbits considered, the best 10% of the models show good agreement with the observed B ϕ , provided that their RMSE was less than ∼ 8 nT, which is on the order of deviations from mean-field behavior seen in the measurements. The best-fit model outputs and their associated Ω m profiles are presented in Figure 3 for four representative cases with differing B ϕ profiles: Revs 271, 278, 289 and 292, where the lowest RMSE best-fit models * B  , and their associated * Ω m profiles are shown using red lines. The model output for the s = 1 case, 1 s B  corresponding to the baseline Ω profile are also presented for comparison.
Rev 271 represents the typical B ϕ observed on 15 of the Grand Finale orbits, where the inbound and outbound observations are largely symmetric on magnetically conjugate latitudes, and a strong positive gradient in B ϕ is seen on the inbound leg of the intra D-ring region. This behavior is well modeled by the 1 s B  model, with some disparity in the magnitude of B ϕ for |θ i | < 20°. The best-fit model * B  , shows much better AGIWAL ET AL.
10.1029/2020JE006578 9 of 18 agreement with the gradients and magnitude of the observed B ϕ profile due to a small increase in the angular velocity of the southern hemisphere winds (by s = 1.4) from the baseline Ω profile.
It is, however, impossible for us to separate the effects of variable angular velocity shears from variable Pedersen conductances for cases such as Rev 271 using our simple approach. We know from Equation 1 that a simple global scaling of the Σ P profile would also be able to improve agreement between 1 s B  and the observed B ϕ . As both sources of variability are plausible, the * Ω m values shown in Figure 3 are subject to being scaled relative to variability in the Σ P values (up to a factor of 10 for |θ i | < 15°), with no overall difference in the modeled B ϕ . For orbits such as Rev 292, however, where the B ϕ observations are negative on field lines mapping to the near-equatorial latitudes, we can be certain that the observed variability is at least in part dependent on a reversal in sense of the angular velocity shears, and then the magnitude of the negative angular velocity shear would scale relative to the Σ P values.
The Rev 271 observations are also well described by Ω m profiles with an N2 profile in the northern hemisphere and a smaller s factor profile with a fixed point at − 15°. We are unable to comment on the likelihood of equatorially asymmetric zonal flows in Saturn's thermosphere, as it would require a significant deviation from the zonal flows observed at the 1-bar pressure level presented in Figure 1c. However, this result does indicate a need for caution when inferring thermospheric behavior from magnetic field measurements in the context of absolute zonal flow values, as the magnetic field observations are only directly indicative of the angular velocity shear.
The other three Revs in Figure 3 all show similar results, where the best-fit models correspond to Ω m profiles with both N1 and N2 profiles in the northern hemisphere. The small scale structure of the observed B ϕ profiles is largely captured in the red envelope which indicates the region in which the models are able to acceptably emulate the observed behavior, however it should be noted that the small scale structure in B ϕ could be unrelated to the inter-hemispheric current system modeled in this study. While we do not think that the central B ϕ feature observed on Revs 271-292 is due to dispersive Alfvén waves, it is entirely plausible that the small scale structure observed during Rev could be generated by this mechanism. For Revs 271, 289 and 292, the best-fit * B  shows excellent agreement with the large scale structure of the observed B ϕ . For Rev 278, however, the best-fit model appears to only capture the mean-field behavior of the observations. The deviation from the mean field behavior for Rev 278 is observed to be on the order of ∼ 5 − 10 nT, which is on the order of the small-scale fluctuations observed in B ϕ during the outbound leg of Rev 271 and inbound leg of Rev 292. Fluctuations of smaller magnitude are also seen on Rev 289, 292, and most of the other Grand Finale orbits (shown in Figure 1). The multiple reversals of B ϕ gradients observed during Rev 278 could be an indication of highly complex zonal flow structures that are not captured by our hypothetical Ω m profiles. Alternatively, they could be related to other sources of spatial or temporal variability, such as LT variability in spacecraft position; in-falling charged ring material from the planetary rings in the near-equatorial latitudes (Hsu et al., 2018;Shebanits et al., 2020); or dispersive Alfvén waves in the magnetosphere (Southwood et al., 2020). The solar wind may also have some influence over the low-latitude ionospheric currents if equatorward transport of heat from the auroral latitudes can be assisted by mechanisms such as gravity waves breaking in the thermosphere (modeled by Müller-Wodarg et al. (2019) at Saturn), or the time-dependent response of the giant planet atmosphere to solar-wind shocks (modeled by Yates et al. (2014) at Jupiter).
As briefly discussed in Section 1.1, we also observe nonzero gradients in B ϕ on magnetic field lines that are connected to the B-ring during the outbound leg of most of the Grand Finale orbits. A nonzero B ϕ signature is also seen at the latitudes corresponding to these magnetic field lines (47°, −41°) in the modeled B ϕ profiles presented in Figure 3 due to the negative angular velocity shear in the baseline Ω profile employed. However, the modeled signature is observed for both the inbound and outbound legs of the orbit due to the axisymmetric assumptions made in our model, and the gradients in the modeled B ϕ for θ i < 0° are in the opposite sense to the gradients observed in the data.
The asymmetry between the inbound and outbound leg of the observation could perhaps be understood by asymmetries in the Σ P values modeled at the LT position of the spacecraft, shown in Figure 2a. For the inbound leg, the Pedersen conductances used in our model are larger than the modeled Σ P at the LT position of the spacecraft for the A-and B-ring latitudes, however for the outbound leg, the Σ P used in our model lie in the middle of the range of Σ P values evaluated at the spacecraft LT. If we assume a nonzero steady-state angular velocity shear for these latitudes, then the modeled B ϕ is likely over-estimated for the inbound part of the orbit due to overestimating the Σ P in the southern ionosphere.
The difference in the sense of the observed and modeled gradients cannot be understood from our analysis, however, since a neutral thermospheric zonal wind profile presented by Brown et al. (2020) for the mid-to high-latitudes inferred from stellar occultations show the zonal flows to be ∼−200 ms −1 and −300 ms −1 relative to planetary rotation, and at the B-ring latitudes, 47° and −41° respectively. These zonal flows correspond to an angular velocity shear of ∼ 1.9 rad s −1 , which is equal in magnitude and sense to the ΔΩ value evaluated at these latitudes from our baseline Ω profile.
The gradients in the observed B ϕ on the B-ring field lines are largely repeatable from orbit to orbit, therefore we think it is unlikely that they are driven by inter-hemispheric shears in angular velocities of neutral zonal flows. Instead, the currents associated with these gradients would likely have to close in the magnetosphere nearer the equatorial plane, presumably in the ring material. Sulaiman et al. (2019) present evidence of whistler mode waves being launched in the magnetosphere due to electrons traveling away from the ring plane in both hemispheres along magnetic field lines mapping to the B-ring, where the opposite flow of currents in both hemispheres may be indicative of ring-planet coupling in the Saturnian magnetosphere. Their reported direction of electron flow is also consistent with positive deflection in the B ϕ observations in the southern hemisphere, which corresponds to the flow of field-aligned currents into the ring plane. Further work is needed to determine to precise mechanism associated with the B-ring signature, however that is beyond the scope of this study.
The Ω m profiles have been successful in describing the mean field behavior for all 21 Grand Finale orbits in the intra D-ring region using profiles scaled between s = −1 and s = 1.6 for the models which satisfy the broader criteria for a good fit, and between s = −0.2 and s = 1.4 if only considering the lowest RMSE bestfit models from each Rev. This constrains the temporal variability in ΔΩ to ∼ 5.5 × 10 −6 rad s −1 at the − 4° boundary, which is where the neutral flows generally peak in our sampled interval. The temporal variability in angular velocity flows is equivalent to Δu ϕ ∼350 ms −1 for the 5 months interval of the Grand Finale orbits, where u ϕ = Ωρ and ρ is taken to be ∼1R S .

Timeseries of Angular Velocity Shears in Neutral Zonal Winds
We re-arrange Equations 1 and 2 in order to empirically determine the angular velocity shears ΔΩ(S − N) along the spacecraft trajectory for each Rev. We use the detrended B ϕ presented in Figure 1b at 1 min resolution and assume the same steady-state, axisymmetric Pedersen conductivities used in our modeling to solve Equations 1 and 2 for ΔΩ(S − N).
In Figure 4, we present a timeseries of these empirically constrained angular velocity shears, henceforth ΔΩ′, as a function of Cassini Rev number using light green and dark green lines for the inbound and outbound parts of each Rev. The ΔΩ′ values are presented for five magnetically conjugate latitudes in our inferred region of wave activity, with the southern latitudes at − 23°, − 17°, − 13°, − 10°, and − 5° in Figures 4a−4e, respectively. We also show the angular velocity shears from to the lowest RMSE * Ω m profiles for each Rev as black circle markers connected by a dashed black line for comparison. The corresponding local * Ω m values at the northern and southern end of the field line are shown using blue and orange triangle markers, with corresponding vertical bars showing the range of Ω m values from the wider range of models that also fit the observed B ϕ well. We again note that the range of velocities suggested are subject to being scaled relative to Σ P .
Additionally, we present new modeling results from STIM, where Müller-Wodarg et al. (2019) have forced the simulation at its lower boundary with an atmospheric Kelvin wave inferred from the model of Medvedev et al. (2013) with a period of six Saturn rotations. The inclusion of this wave leads to temporal variability in the modeled zonal flows, which is consistent with our hypothesis. We calculate the range between the minimum and maximum angular velocity shear values ΔΩ(Δt) observed in an interval of 18 Saturn days at 1 day resolution for each of the magnetically conjugate latitude pairs presented in Figures 4a−4e, and present this range as a shaded purple region in each panel. The corresponding Ω N,S ranges at the local northern and southern latitudes are shown in blue and orange respectively, outside each panel.
The ΔΩ′ profiles presented in Figure 4 show excellent conjugacy between the inbound and outbound profiles for all Revs in Figures 4d−4e, where the spacecraft traverses the same magnetic flux tube ∼8-10 min apart. This indicates that the ionospheric currents must be in steady-state on timescales of ∼10 min and azimuthally symmetric between 11 and 13 h in LT between −5° and −10° in latitude in the southern hemisphere, and between 15° and 19° in latitude in the northern hemisphere. However, the inbound and outbound ΔΩ′ profiles presented in Figures 4a−4c are largely nonconjugate, which is mainly due to the polynomial subtraction process outlined in Section 1.1, where the detrended B ϕ profiles presented in Figure 2c show some asymmetries between the inbound and outbound passes for each Rev. If such asymmetries are physical, the nonconjugacy could indicate that the Pedersen conductances or angular velocity shears must either be spatially variable across ∼2 h of LT at the latitudes shown, or temporally variable on timescales longer than 10 min, which is consistent with results of Hunt et al. (2019)    We investigated the possibility of this nonconjugacy being a spatial effect due to the LT asymmetries in the modeled Σ PS values between the pre-noon (inbound) and post-noon (outbound) sectors shown in Figure 2a. The ΔΩ′ profiles evaluated using the azimuthally variable Pedersen conductances along the spacecraft trajectory showed little difference from the ΔΩ′ profiles presented in Figure 4 for observations from the intra D-ring region. While this does implies that the nonconjugacy is more likely due to temporal variability in I m , it does not necessarily rule out spatial variability as a source for the nonconjugacy, since the Σ P values used in our investigation are modeled values, and we assume azimuthally symmetric zonal flows, which may be an oversimplification.
The modeled * ΔΩ m profiles show very good agreement with the ΔΩ′ profiles presented in Figure 4 corresponding to the −5°, −10° latitude panels where the ΔΩ′ profiles show excellent conjugacy, while in the remaining three panels, the * ΔΩ m profiles seem to be consistent with the average behavior of the inbound and outbound ΔΩ′ profiles. Given that our modeled B ϕ profiles, which are perfectly conjugate between the two hemispheres due to the steady-state and axisymmetric assumptions made when modeling I m , are optimized using RMSE, the models which return the lowest RMSE are generally consistent with the statistical mean of the behavior between the two hemispheres. If we had fit to the inbound and outbound data separately, we may have been able to improve the agreement between the ΔΩ′ and * ΔΩ m profiles presented for the higher latitudes. However, in the absence of observational constraints on Σ P we find it suitable to limit our investigation to orbit by orbit variability in angular velocity shears.
The overall agreement between our modeled results and the observations indicates that the range of * Ω m models presented in Figure 4 in both the northern and southern hemisphere could potentially be used to impose constraints on the magnitude and sense of neutral zonal flows in future ionospheric modeling. The Ω N,S ranges from STIM presented outside each panel in Figure 4 additionally support some of the assumptions made when modeling the Ω m profiles in this study. The range of northern hemisphere Ω values are generally much smaller than the range of southern hemisphere values, which is consistent with our treatment of varying the baseline Ω profile more significantly in the southern hemisphere than the northern for a given latitude. The STIM ranges also show that it is possible to have Ω profiles that are asymmetric about the magnetic equator since the range of northern Ω values indicates mostly westward flows for all panels in Figure 4, but the range of southern Ω values corresponds to both eastward and westward flows.
The range of ΔΩ(Δt) values observed in the STIM output show that the simulation is able to evaluate both positive and negative angular velocity shears during an interval of 18 Saturn rotations (timescales of ∼1 Cassini Rev) for the latitude pair corresponding to − 5° in the southern hemisphere, however for all other latitude pairs, the evaluated shears are either entirely positive (Figures 4b−4d), or entirely negative (Figure 4a). The entirely positive/negative angular velocity shears are inconsistent with ΔΩ′ profiles presented in Figure 4, and the extent of variability in the ΔΩ(Δt) values is much smaller than the extent of variability seen in the ΔΩ′ and * ΔΩ m profiles. However, the simulations presented here illustrate the principle of atmospheric waves causing ionospheric perturbations which can be observed in the magnetosphere, which are overall consistent with the observations from the Cassini Grand Finale.

Conclusions
We investigate variable ionospheric zonal neutral winds as a source of temporally variable low latitude field-aligned currents observed during the Cassini Grand Finale orbits. Through consideration of magnetic field geometry at Saturn, we constrain the region of atmospheric variability to be between ±25° and create an extensive selection of hypothetical angular velocity Ω m profiles that are used to model B ϕ in a thin-sheet approximation of Saturn's ionosphere as derived by Provan et al. (2019b). The modeled B ϕ are compared to 21 Revs from the Cassini Grand Finale orbits, and a diverse range of hypothetical Ω m profiles are found to show comparably good agreement with the observed B ϕ for each Rev.
We infer the angular velocity shears along the spacecraft trajectory directly from the Grand Finale B ϕ observations using the Provan et al. (2019b) approach, and use them to present the first empirically constrained timeseries of angular velocity shears in neutral ionospheric winds at Saturn. These results are well matched by the best-fit modeled angular velocity shears in our study, particularly at near-equatorial latitudes. The overall agreement between our modeled results and the observations indicates that the range of angular velocities corresponding to the best-fit models could be used to impose constraints on the expected angular velocity of neutral zonal flows in future ionospheric modeling works.
Our analysis shows that the ionospheric currents must be in steady-state on timescales of ∼ 10 min and azimuthally symmetric between 11 and 13 h in LT between −5° and −10° in latitude in the southern hemisphere, and between 15° and 19° in latitude in the northern hemisphere. We constrain the temporal variability of the peak zonal flows in the southern hemisphere to ∼350 ms −1 over an interval of ∼5 months, which is subject to being scaled by variability in Pedersen conductances, which could be variable by up to a factor of ∼10 at near equatorial latitudes .
Our study has illustrated the use of magnetospheric datasets to potentially constrain atmospheric variability in the ionospheric layer. However, our approach is kinematic in the sense that detailed physical processes responsible for the wind variability were not explicitly modeled. Using our results to infer constraints on global Giant Planet thermospheric models such as STIM would be the natural next step toward a dynamical understanding of Saturn's low latitude thermosphere and ionosphere.

Appendix A: Equivalence of the Provan et al. and Vriesema et al. Formulations in the Thin Ionospheric Current Layer Approximation
In the following Appendix, we will refer to all equations from Provan et al. (2019b) with the prefix P, and all equations from Vriesema et al. (2019) with the prefix V.
The formulas describing the electric current system and magnetic perturbations generated by Saturn's near equatorial ionospheric wind dynamo employed in the present study (Equations 1 and 2) are those derived by Provan et al. (2019b), based on the initial simplifying assumptions of steady state conditions, azimuthal symmetry, and equipotential field lines. The first two assumptions require the electric field to be poloidal, with zero azimuthal component, such that the ionospheric Pedersen current is also poloidal, while the Hall current is wholly azimuthal and axisymmetric. The third assumption, equipotential field lines, further requires the angular velocity of plasma, ω, associated with the E × B drift, to be constant along each field line. Provan et al. (2019b) also approximated the ionosphere as a thin conducting layer, in which a field-aligned current flows together with the Pedersen current perpendicular to the inclined magnetic field lines, such that the total ionospheric current in the meridian is directed horizontally as required. Current continuity then requires the angular velocity of the plasma, in terms of the angular velocities of the neutral atmospheric flow in the current layer at the two ends of the field lines (N) and south (S), Ω nN,S , to be given by The corresponding horizontal meridional ionospheric current per radian of azimuth, positive northwards and equal at the two ends of each field line, is given by Here, Σ PN,S are the height-integrated Pedersen conductivities at the two ends of the field line, ρ iN,S are the perpendicular distances from the magnetic/spin axes, B iN,S are the ionospheric field strengths, and |b iN,S are the magnitudes of the ionospheric field strength normal to the current layer. Equations A1 and A2 are equations PA7 and PA8. The field-aligned current flowing above the ionospheric layer is then given by the variation of I m with latitude. From the integral form of Ampère's law, the azimuthal magnetic field generated by the current system above the ionosphere on each field is given by where ρ is again the perpendicular distance from the magnetic/spin axis, and μ 0 is the permeability of free space.
The subsequent modeling paper by Vriesema et al. (2019) is based on the same set of initial assumptions as indicated above, that is, steady state, axisymmetry, and equipotential field lines, but now includes ionospheric and thermospheric models that are resolved in height (pressure level) as well as in latitude. They also specialize to the case of an axisymmetric dipole planetary magnetic field, not necessarily implied in the above analysis, which in general is offset along the spin axis from the planet's center. The analysis then proceeds using an orthogonal curvilinear dipole coordinate system (α, β, ϕ), in which coordinate α is constant on a dipole flux shell, β defines position along a field line varying from −∞ at the southern hemisphere dipole pole to + ∞ at the northern hemisphere dipole pole, and ϕ is the usual azimuthal angle. In terms of usual spherical polar coordinates (r, θ, ϕ) defined relative to the dipole center, α and β are given by where R S is Saturn's radius and L the McIlwain parameter (the largest distance of a field line from the magnetic/rotation axis in R S , and 2 cos . r    (A5) Similar to Provan et al. (2019b), the current continuity condition is then applied to the meridional Pedersen and field-aligned currents, the Hall current again being azimuthal and axisymmetric, and thus not contributing to continuity. This condition then yields the following solution for the field-transverse meridional electric field E α associated with the plasma flow in and above the ionosphere (equation V23) (  , σ P is the ionospheric Pedersen conductivity, u ϕ the thermospheric azimuthal wind speed, and the integrals are taken along individual field lines from just under the ionosphere (where σ P drops to zero) in the southern hemisphere, point β 1 (α), to just under the ionosphere in the northern hemisphere, point β 2 (α). We note from the definition of the coordinates in Equations A4 and A5 that the planetary magnetic field B has only a β component, given by where B is the planetary field strength given by 0 3 2 1/ 2 1 3 (1 3cos ) . in the (shifted) dipole coefficient for Saturn, and the negative sign in Equation A7 is required by the definition of β, which increases along field lines from south to north, opposite to the direction of Saturn's planetary field. The h factors in Equation A6 are scale factors that relate increments in the coordinates to spatial displacements (i.e., the Lamé coefficients for dipole coordinates), which can be written as In particular, the relationship between an increment dβ and the distance ds along a field line, taken positive in the direction of the field, is, using Equation A10 The constancy of h α E α along field lines is therefore equivalent to the constancy of the plasma angular velocity along field lines, as may be expected. Writing the thermospheric azimuthal velocity as u ϕ = ρΩ, where Ω is the neutral atmospheric angular velocity, and substituting Equations A7, A9, A10, A11, and A14 into Equation A6 then yields the following expression from the plasma angular velocity where the integrals are now taken along a field line from just under the ionosphere in the north (position s N where σ P goes to zero) to just under the ionosphere in the south (position s S where σ P similarly goes to zero). The plasma angular velocity on a field line is thus given by an integral of the thermospheric angular velocity on the field line weighted principally by the Pedersen conductivity σ P .
Integration of Ampère's law gives the azimuthal magnetic field as (equation V28) where the Pedersen current density j α is given by To compare with the results of the analysis of Provan et al. (2019b) in Equations A1 and A3 above, we must consider Equations A15 and A6 in the limit that the ionosphere is taken to be a thin conducting layer. We consider field lines such as those traversed during the Cassini proximal periapsis passes in which the integrals along field lines in Equation A15 can be taken to consist of two contributions, one from the northern ionosphere, the other from the southern, separated by a field segment above the ionosphere where σ P ∼ 0, which thus makes no contribution to the integral. We further assume that on the two field segments that pass through the ionosphere, ρ i , Ω n , and B i in Equations A1 and A2 may be taken to be approximate constants, though not in general the same in the northern and southern hemispheres. We also write