Herschel map of Saturn's stratospheric water, delivered by the plumes of Enceladus

Context. The origin of water in the stratospheres of Giant Planets has been an outstanding question ever since its first detection by ISO some 20 years ago. Water can originate from interplanetary dust particles, icy rings and satellites and large comet impacts. Analysis of Herschel Space Observatory observations have proven that the bulk of Jupiter's stratospheric water was delivered by the Shoemaker-Levy 9 impacts in 1994. In 2006, the Cassini mission detected water plumes at the South Pole of Enceladus, placing the moon as a serious candidate for Saturn's stratospheric water. Further evidence was found in 2011, when Herschel demonstrated the presence of a water torus at the orbital distance of Enceladus, fed by the moon's plumes. Finally, water falling from the rings onto Saturn's uppermost atmospheric layers at low latitudes was detected during the final orbits of Cassini's end-of-mission plunge into the atmosphere. Aims. In this paper, we use Herschel mapping observations of water in Saturn's stratosphere to identify its source. Methods. Several empirical models are tested against the Herschel-HIFI and -PACS observations, which were collected on December 30, 2010, and January 2nd, 2011 (respectively). Results. We demonstrate that Saturn's stratospheric water is not uniformly mixed as a function of latitude, but peaking at the equator and decreasing poleward with a Gaussian distribution. We obtain our best fit with an equatorial mole fraction 1.1 ppb and a half-width at half-maximum of 25{\deg}, when accounting for a temperature increase in the two warm stratospheric vortices produced by Saturn's Great Storm of 2010-2011. Conclusions. This work demonstrates that Enceladus is the main source of Saturn's stratospheric water.


Introduction
The interiors of Giant Planets are supposedly rich in oxygen (Owen & Encrenaz 2003;Gautier et al. 2001;Hersant et al. 2004). In the deep hot layers, where thermochemical equilibrium prevails, H 2 O is the most abundant oxygen species. Oxygen species are transported from deep down towards higher lev-‹ Herschel is an ESA space observatory with science instruments provided by European-led Principal Investigator consortia and with important participation from NASA. els, but only CO (and CO 2 in Jupiter and Saturn) can reach the stratosphere because of the tropopause cold trap (Lodders & Fegley 1994;Wang et al. 2015;Cavalié et al. 2017). H 2 O was thus only observed in the tropospheric saturation layers of Jupiter and Saturn (Larson et al. 1975;de Graauw et al. 1997). The detection of oxygen species (H 2 O, CO and CO 2 ) in the stratospheres of the Giant Planets and Titan Bézard et al. 2002;Lellouch et al. 2002;Coustenis et al. 1998;Samuelson et al. 1983;Noll et al. 1986;Marten et al. 1993;Burgdorf et al. 2006) thus demonstrated that the upper atmospheres of the Article number, page 1 of 16 arXiv:1908.07399v1 [astro-ph.EP] 20 Aug 2019 A&A proofs: manuscript no. Saturn_H2O_map_observations giant planets are contaminated by external sources from their close and/or more distant environments.
While oxygen-rich interplanetary dust particles (IDP) produced from asteroid collisions and from comet activity are a ubiquitous source (Landgraf et al. 2002;Poppe 2016), it seems to be the main H 2 O source only at Uranus and Neptune (Moses & Poppe 2017) and the emerging overall picture however looks more complex. Other sources can actually be at work, like local sources from planetary icy environments (rings, satellites) (Strobel & Yung 1979;Connerney 1986;Prangé et al. 2006;Waite et al. 2018;Mitchell et al. 2018;Hsu et al. 2018), or cometary "Shoemaker-Levy 9 (SL9) type" impacts (Lellouch et al. 1995). In Jupiter's stratosphere, H 2 O, CO, and CO 2 come from the SL9 comet fragments Lellouch et al. 2002Lellouch et al. , 2006Cavalié et al. 2008aCavalié et al. , 2012Cavalié et al. , 2013). An older comet-impact component was even proposed to explain its the CO in the lower stratosphere , even if IDP are another option (Moses & Poppe 2017). Comets are also the probable source of CO beyond Jupiter (see review in Mandt et al. 2015 andPoppe 2017), as seen in Saturn (Cavalié et al. 2009(Cavalié et al. , 2010, Uranus , and Neptune (Lellouch et al. 2005(Lellouch et al. , 2010Luszcz-Cook & de Pater 2013;Moreno et al. 2017).
In Saturn, neither ISO nor SWAS disk-averaged observations Bergin et al. 2000) had sufficient spectral resolution and/or high enough signal-to-noise ratio to identify the source of external H 2 O. Now, a comet impact is probably not the cause of the observed stratospheric H 2 O (Moses & Poppe 2017) because of the contradiction between (i) the relatively ancient impacts required to fit the CO observations of Cavalié et al. (2010) and (ii) the relatively short diffusion timescale from the deposition level of cometary material in such impacts (Lellouch et al. 1995;Moreno et al. 2003) down to the H 2 O condensation level (Ollivier et al. 2000;Moses et al. 2000Moses et al. , 2005, which is located between the 1 and a few mbar level, depending on the latitude. In addition, the H 2 O/CO ratio in Saturn's stratosphere ("0.15 -Moses & Poppe 2017) is too large to be characteristic of a cometary impact, that delivers mostly CO (H 2 O/COă0.01 in Neptune). This tends to indicate a source that is steadier than a discrete comet impact. Previous observations of Saturn with Herschel-HIFI led to the detection of an H 2 O torus at the orbit of Enceladus (Hartogh et al. 2011), fed by the plumes of this moon (Hansen et al. 2006;Porco et al. 2006;Waite et al. 2006). The fate of H 2 O from this torus is eventually to spread in Saturn's system (Cassidy & Johnson 2010) and a fraction is predicted to fall into Saturn's stratosphere with a distribution centered around the equator. Based on these models, and comparing with the Herschel observations, Hartogh et al. (2011) tentatively concluded that Enceladus was the source of Saturn's stratospheric water. It should be noted that a fraction is also expected to feed Titan's atmosphere as well, and models show that the flux expected at Titan and originating from the Enceladus water torus could explain the observations (Moreno et al. 2012;Lara et al. 2014;Dobrijevic et al. 2014;Rengel et al. 2014;Hickson et al. 2014).
Recent in situ measurements of Cassini during the proximal orbits of its end-of-mission demonstrated the existence of a flux of material from the rings to Saturn's atmosphere. The Magnetospheric Imaging Instrument (MIMI), the Ion and Neutral Mass Spectrometer (INMS), and the Cosmic Dust Analyzer (CDA) instruments have measured an infall of material originating from the D-ring: (i) neutral icy grains within 2˝around the equator, (ii) various gases including H 2 O (concentrated at the 24% level) in a latitudinal band of 8˝centered on the equator ), (iii) and charged grains transported to higher latitudes along the magnetic field lines Hsu et al. 2018). This latter phenomenon is the long anticipated ring rain phenomenon (Connerney & Waite 1984;Connerney 1986;Prangé et al. 2006;Moore et al. 2010Moore et al. , 2015O'Donoghue et al. 2017). The mass flux of infalling dust and gas are estimated to be "5 kg/s  and "10 4 kg/s Perry et al. 2018), respectively. The sources that can explain the presence of H 2 O in Saturn's stratosphere are thus in principle IDP, Enceladus plumes and Saturn's rings.
Disentangling the various sources of externally supplied water in Giant Planet stratospheres, and thus in Saturn, was a key objective of the Herschel key program HssO (Herschel Solar System Observations) (Hartogh et al. 2009). In section 2 of this paper, we present the first disk-resolved mapping observations of H 2 O in Saturn's stratosphere and a disk-averaged observation, obtained with the Photodetector Array Camera and Spectrometer (PACS) (Poglitsch et al. 2010) and the Heterodyne Instrument for the Far Infrared (HIFI) (de Graauw et al. 2010) (respectively), two instruments onboard the ESA Herschel Space Observatory (Pilbratt et al. 2010). Using a combination of models presented in section 3, we derive the H 2 O meridional distribution in Saturn's stratosphere in section 4. We discuss its origin according to our results in section 5 and give our conclusion in section 6.

Observations
In the following sections, we present a summary of the performed observations (Table 1) and the relevant geometry of Saturn (as obtained from the JPL/Horizons database)

PACS observations
We observed Saturn with Herschel/PACS in the framework of the HssO Key Program on January 2, 2011, with the aim of mapping the distribution of H 2 O in the stratosphere of the planet. We used PACS in its line spectrometry mode (Poglitsch et al. 2010). The integral field spectrometer consists of 5ˆ5 spatial pixels ("50 2ˆ5 0 2 on sky), each covering a short instantaneous wavelength range sampled by 16 spectral pixels. Our observation of Saturn consists of a 3ˆ3 raster map with 3 2 separation at 66.44 µm. We thus pointed the 25-pixel array over 9 different positions to record an oversampled 225-point map of Saturn at 66.44 µm. For the observation of such an intense farinfrared continuum source, we had to use a non-standard mode, in which the spectrometer readout electronics were configured to the shortest possible reset intervals of 1/32 s, to avoid detector saturation. The half-power beam width (HPBW) of Herschel is 9.4 2 for the PACS observations, which results in covering latitudes from´23˝to 45˝for a beam centered on the planet. The spectral resolving power is "2500-3000.
The basic data reduction was run with HIPE 8.0 (Herschel Interactive Processing Environment, Ott 2010) and additional processing (flat-fielding, outlier removal and rebinning) was performed with standard IDL tools, providing us with the 225 spectra presented in Fig. 1 for all raster positions of the map. After subtracting the coordinates of Saturn given by the JPL/Horizon ephemeris, we find from the continuum distribution that the map has a residual pointing offset of 2.5 2 in RA and 0.3 2 in dec. The position of the line peak presents some scatter, which results from a combination of two effects: (i) a Doppler shift caused by the rapid rotation of the planet, and (ii) wavelength shifts induced by the non-uniform illumination of the instrument slit on some of the raster positions (Poglitsch et al. 2010). A baseline, caused by standing waves excited by the strong continuum emission of the ote: OD means Herschel operational day, ∆t is the total integration time. Ang. Diam. is Saturn's equatorial angular diameter, λ obs and φ obs are the longitude and latitude (respectively) of the sub-observer point, λ @ and φ @ are the longitude and latitude (respectively) of the sub-solar point, a NP is the North Polar angle. The solar longitude of Saturn (L S ) is 17˝. All latitudes in this table are planetographic, and all longitudes are System III West longitudes. The physical parameters of Saturn (right column) have been obtained from the JPL/Horizons database. planet in the instrument, was then removed from the raw spectra by fitting a polynomial. This step introduces on average a 10-15% uncertainty on the line amplitude. The H 2 O line is detected in all positions within the disk and on positions within "one beam from the planetary limb. For the reasons detailed in Cavalié et al. (2013), reasonable uncertainties cannot be achieved on an absolute flux calibration and we have to express the spectra in terms of line-to-continuum ratio (l{c) to cancel out the absolute response. In addition, the line width is purely instrumental (and Gaussian) and can vary from one position to another because of varying beam filling factors on the detector array. Therefore, the map must be analyzed in terms of line area. We fit the lines with a Gaussian and compute their area from the Gaussian fit parameters. According to the noise on the line peak intensity and line width, we estimate that the line-line-area maparea map can be safely analyzed for positions within the planetary disk or in the vicinity of the limb (see Fig. 2). We find a 1-σ rms of 0.0023 (in units of micronsˆ% of the continuum) when adding quadratically the baseline removal uncertainty and the spectral noise. We use this value in our χ 2 {N calculations (N is the number of pixels in the map). The average S/N (signal-to-noise ratio) around the limb, i.e., where the lines have the biggest contrast, is "15-45 on the line area, with maximum value at the eastern and western limbs and a minimum at the northern limb. At a given latitude, the increased emission at the limb (compared to the disk center) is a geometrical effect, caused by limb-brightening in the line and limb-darkening in the continuum. There is enhanced emission at low latitudes with respect to mid-to-high latitudes, and a local minimum around the North Pole. The resulting continuum and line-area maps are presented in Fig. 2   The line area in micronsˆ% of the continuum is displayed. The S/N in this map varies from "15 to "45 (the 1-σ rms is 0.0023 micronsˆ% of the continuum), and enables us to demonstrate the difference of water emission between the poles and equatorial region despite the limited spatial resolution. The map shows a lack of emission at high northern latitudes, and to a lower extent in the high southern latitudes as well.

HIFI observations
In addition, and in an attempt to check the PACS results at the planetary scale, we use the disk-averaged observation of Saturn conducted with the HIFI instrument (de Graauw et al. 2010) at 1097 GHz on December 31, 2010, with a dual beam-switch mode (Roelfsema et al. 2012), which enables pointing to the deep sky for calibration by moving the secondary mirror instead of the primary. We have chosen this strong H 2 O line because it is not affected by the Enceladus torus absorption which is seen only in lines like the 557 GHz and 1113 GHz lines (Hartogh et al. 2011) that have a low lower state energy relative to the ground state (23.7944 and 0.0 cm´1, respectively); the lower state energy of the 1097 GHz line is 136.7617 cm´1. Given the Herschel HPBW is 19.3 2 (i.e., bigger than Saturn) and the spectral resolution of 1.1 MHz of the Wide Band Spectrometer (WBS), the line appears smeared because of the fast rotation of Saturn (9.9 km¨s´1 at the equator). We processed the data with the standard HIPE 8.2.0 pipeline (Ott 2010) up to level 2 for the H and V polarizations. We then remove standing waves by using a Lomb (1976) algorithm before averaging the two polarizations and smoothed the spectral resolution to 10 MHz. This last operation barely changes the line shape as the width of the line is already "35 MHz because of the aforementioned rotational smearing. The resulting line is displayed in Fig. 3, and we note a small asymmetry and the fact that the line-center frequency is not aligned with the line rest frequency. This asymmetry is caused by a pointing offset of 1.5 2 in the planetary western limb direction that we account for in what follows. We performed no absolute calibration and thus cannot constrain the North-South pointing offset, which must be of the order of a few arc seconds according to the observatory pointing uncertainty at this frequency (Roelfsema et al. 2012). We account for the double sideband (DSB) response of the instrument by assuming a normalized sideband ratio G US B pH`Vq " 0.469˘0.016 (Higgins et al. 2014) and a ratio between the continuum in the upper and lower sidebands of 1.037 (according to our continuum model), and produce a single sideband (SSB) spectrum that we use for our analyses in terms of line-to-continuum ratio (l{c) in this paper. The total uncertainty on the line-to-continuum ratio, when accounting for the uncertainties caused by the normalized sideband ratio, the baseline removal process, and the spectral noise, is 3%, i.e., 0.08 K. All our χ 2 {N calculations are based on this rms value, which does not include absolute calibration uncertainties that are estimated to be 5%.  Fig. 3. Water line at 1097.365 GHz observed by the HIFI spectrometer on December 31, 2010. The spectrum is expressed in terms of kelvins on the DSB scale. The asymmetry seen in the line shape is caused by a pointing offset of 1.5 2 in the planet western limb direction. The North-South pointing offset is of the order of a few arcsec, but cannot be constrained further because of calibration uncertainties. The total uncertainty on the line-to-continuum ratio is 3%.

Modeling
In this section, we present the thermal and radiative-transfer models that have been applied in our data analysis.

Thermal model
The pressure-latitude thermal field we use in our radiativetransfer modeling is extracted from the spline interpolation made over time to fill in missing parts of Cassini's Composite Infrared Spectrometer (CIRS) time series spanning the full mission (Fletcher et al. 2017). We can thus extract the zonally-averaged stratospheric temperatures pertinent to the dates of the Herschel observations. At the time of our observations, the CIRS dataset features significant temperature anomalies associated with the planetary storm that erupted at the end of 2010 (Sánchez-Lavega et al. 2011;Fischer et al. 2011;Fletcher et al. 2011). This storm significantly changed the northern stratospheric temperatures in a 20˝latitude band centered on 40˝N planetocentric. Two hot stratospheric vortices formed well above the tropospheric storm and were nicknamed "beacons" because of their bright appearance in the IR and the rapid rotation of Saturn. They eventually merged in April-May 2011 to form a huge vortex in which a temperature increase of up to 80 K was measured at the 2-mbar level. But, because our observations occur quite early in the time evolution of the storm and beacons, we first choose to produce a thermal field representative of storm-free conditions by removing any coverage between 20N-60N and 2011-2012 in the spline smoothing, and will include the beacons later in the paper (in Section 4.3). We take the thermal field produced for the date closest to our observations, i.e., December 28, 2010. The stratospheric temperatures in the 0.1-10 mbar range, i.e., roughly the range in which the HIFI and PACS lines are sensitive (see Fig. 8 right), are quite symmetric in latitude with respect to the equator up to "40˝, with a gradient from 40S to the South Pole ă5 K and a´10 K gradient from 40N to the North Pole, as already hinted for observations at similar L S by Sinclair et al. (2013). The CIRS observations used 600-1400 cm´1 in nadir geometry, therefore probing only the 0.2-5 mbar and 70-250 mbar ranges. Outside of these ranges, the profile goes back to an a priori profile built by averaging Guerlet et al. (2009Guerlet et al. ( , 2010 limb observations, which probe higher stratospheric altitudes than the study of Fletcher et al. (2017), between latitudes of 45N and 45S. However, whether this a priori is valid for all latitudes is questionable. Therefore, we choose to extrapolate isothermally the temperature profiles for pressures lower than 0.2 mbar. Between 20S and 20N, we use the thermal profiles of February 2010 derived by Guerlet et al. (2011) from Cassini/CIRS limb viewing geometry observations in order to benefit from the higher vertical resolution of such observations. With such data, the equatorial quasi-periodic oscillation (Fouchet et al. 2008;Orton et al. 2008;Guerlet et al. 2011) is better resolved than in the nadir data, and the vertical sensitivity in the temperature retrieval is extended to 10´2 mbar. Fig. 4 displays the thermal field. We note that the temperatures in the submillibar range are in agreement with temperatures derived from Voyager 2 occultation experiments at several latitudes by Lindal et al. (1985) for a similar L S .

H 2 O distribution empirical models
In this paper, we test two types of H 2 O spatial distributions that are representative of two different sources: IDPs and the Enceladus plumes.  The IDP source is modeled with a meridionally uniform distribution of H 2 O. The mole fraction of H 2 O is uniform above the local condensation level, which is computed following Fray & Schmitt (2009) at each latitude.
The Enceladus source is modeled with a 2-parameter Gaussian-shaped latitudinal distribution of H 2 O, centered on the equator, with a peak H 2 O mole fraction at the equator y eq , and a half-width at half-maximum (HWHM) σ, such that: where y eq is the H 2 O mole fraction above the condensation level at the equator, and φ is the planetocentric latitude. Here also, latitudinally-dependent condensation is computed from our thermal field. The recently observed ring source Hsu et al. 2018;Mitchell et al. 2018) shares common properties with the Enceladus source: it is centered on the equator and has, to the first order, a Gaussian shape.
The ring atmosphere (mostly O 2 , not H 2 O), another potential source for Saturn's stratospheric H 2 O, has lower densities than the Enceladus neutral torus (Johnson et al. 2006;Cassidy & Johnson 2010). Neutral-neutral collisions are therefore inefficient, and even ion-neutral charge exchange collisions are slow. Unlike elsewhere in the magnetosphere, the co-rotating plasma flow speed is similar to the orbital speed. The ring atmosphere is therefore lost at higher latitudes like ring rain: the atmosphere is lost via photoionization followed by precipitation along field lines (Moore et al. 2015). We note that the neutral component of their H 2 O distribution is actually tied to Enceladus. This is why we do not consider the ring atmosphere as a significant enough source to explain our observations, and will not include it in our models.

Radiative-transfer model
We model the observations with a line-by-line non-scattering model that solves the radiative-transfer equation in 1D on several thousand lines-of-sight, thus highly oversampling our observational spatial resolution. We use the spectral line parameters from the JPL Molecular Spectroscopy Database (Pickett et al. 1998). The spectra of collision-induced absorption caused by H 2 -H 2 , H 2 -He, and H 2 -CH 4 pairs are included following Borysow et al. (1985Borysow et al. ( , 1988 and Borysow & Frommhold (1986). More recent updates to the H 2 -H 2 absorption from ab initio quantum theory exist (Orton et al. 2007;Fletcher et al. 2018a), but they do not make any significant difference in this spectral region. We adopt mole fractions of 11.8% and 0.47% for He and CH 4 respectively, according to Conrath & Gautier (2000) and Fletcher et al. (2009b). We account for the absorption caused by the H 2 O line, as well as surrounding NH 3 and PH 3 broad lines. We take the NH 3 and PH 3 distributions from Davis et al. (1996) and Fletcher et al. (2009a), respectively. We compute the relevant H 2 /He broadening parameters from Levy et al. (1993Levy et al. ( , 1994 for PH 3 , from Fletcher et al. (2007) for NH 3 , and from Brown & Plymate (1996) and Dick et al. (2009) for H 2 O, for the H 2 /He mixture described above and the appropriate temperature range. The 1097 GHz line has an opacity of "0.3 over the HIFI beam in disk-centered conditions, and the 66.44 µm line has an opacity that ranges from 0.1 to 0.4 over the PACS beam for disk-centered and limb conditions, for the best fit model H 2 O profiles given in Fig. 15.
We note that, contrary to the 66.4 µm line seen by PACS, the 1097 GHz line seen by HIFI lies on the far wings of a PH 3 line. A change in the PH 3 abundance therefore influences the line-tocontinuum ratio of the H 2 O line and subsequently on the derived H 2 O abundance. Increasing the PH 3 abundance by a factor of 2 decreases the H 2 O abundance required to fit the line by "15%.
This model significantly improves on the model already described in Cavalié et al. (2008bCavalié et al. ( , 2013 regarding the handling of the geometry. The observation plane is sampled with an irregular grid so that the limb and rings, i.e., the regions with the highest variations in terms of continuum and line emission, are sufficiently sampled without hampering the computational efficiency. We account for the 3D ellipsoidal geometry of Saturn and rings at the time of the observations, using the parameters listed in Table 1. This means we compute the local latitudes, longitudes, and altitudes, at each sampled location on each lineof-sight, and are thus able to use thermal and abundance fields that are fully 3D, before applying the beam convolution. In January 2011, the rings had an inclination of 12.4˝, as seen from Herschel. We thus account for the contribution of the rings to the continuum emission. The A-B-C ring system brightness temperatures are computed in the wavelength range of our observations by interpolating the observational results obtained in the infrared with Cassini/CIRS (Flandes et al. 2010) and at 2 cm with the Very Large Array (van der Tak et al. 1999), accounting for the brightness temperature roll-off seen around 200 µm by Spilker et al. (2003Spilker et al. ( , 2005. We use the solar elevation as seen from the rings given in Table 1, and the optical depths of the different rings from Altobelli et al. (2014) and Guerlet et al. (2014). Results are in general agreement with Dunn et al. (2005), when considering averages for the whole ring system. In practice, we identify geometrically the lines-of-sight that cross the rings and we account for the ring absorption and thermal emission, once all the lines-of-sight on the planet and limb are computed. Beyond the planetary disk, we find the positions of the rings and add their thermal emission. We then create a regular grid of linesof-sight from which we compute the beam-averaged spectra, in such a way that the beam is highly oversampled by the grid. We produce maps of line integrated emission (for PACS) and diskaveraged spectral line (for HIFI) that can directly be compared with the observation. We present the simulations as well as the residuals resulting from the subtraction of the model from the observations 1 .

Meridionally uniform distribution
An initial approach was to try to fit the PACS map with a uniform distribution of H 2 O which would be representative of an IDP source. As shown in Fig. 5, the H 2 O line-area map is a direct translation of the thermal field and we find a strong gradient between the North and South poles. The minimum χ 2 {N is obtained for an H 2 O mole fraction of 6ˆ10´1 0 above the condensation level (roughly a few mbar, varying with latitude), but the χ 2 {N is as high as 18.9 (i.e., more than 4-σ away from the data). The observed decrease of the line area around the southern limb is not reproduced (see Fig. 5), as a results of the warmer stratospheric temperatures in the South compared to the North (see Fig. 4). On the contrary, the strongest emission comes from the south when we use a meridionally uniform distribution of H 2 O. In addition, there is not enough emission at low latitudes. Such a solution is thus not satisfactory and can be discarded at this stage.
We note that an acceptable fit to the disk-averaged HIFI data can be obtained for a uniform mole fraction of "3.5ˆ10´9, which is consistent with Fletcher et al. (2012) but incompatible with the PACS data. When we simulate the H 2 O line at 556.936 GHz with this distribution, we obtain a line-tocontinuum ratio that is consistent with the previous observation of Bergin et al. (2000). This shows that the diskaveraged H 2 O influx responsible for our observations of "6ˆ10 5 cm´2¨s´1 (Moses & Poppe 2017), which translates into a disk-integrated mass flux of "8 kg/s, is surpassed by orders of magnitude by the extraordinarily high flux recently measured by INMS of "10 4 kg/s coming from the inner ring system. This proves that the ring source observed by Cassini in 2017 cannot be the cause of the H 2 O observed by Herschel in 2010-2011. The ring source must therefore be more recent than our observations, or at least this source had not been active for long enough (i.e., " 10-15 yr) to diffuse down to observable levels ("0.1-1 mbar). According to Waite et al. (2018), such a large flux is probably linked to the appearance of clumps of material in the D68 ringlet in 2015 (Hedman & Showalter 2016), and is unsustainable for more than 1 Myr without depleting the rings. In conclusion, we will thus not consider the ring source in our subsequent modeling of the Herschel observations.

Meridionally variable distributions
Despite the low spatial resolution of the PACS observations, the contrast observed between the eastern/western limbs and the northern/southern ones (Fig. 2 bottom), combined with the knowledge of the altitude/latitude thermal field (Fig. 4), provides us with a good tool for deriving the meridional distribution of H 2 O in Saturn's stratosphere.
In what follows, we test several cases of H 2 O meridional distributions as parametrized in eq. 1. We first explore the y eq -σ parameter space. Fig. 6 shows the χ 2 {N obtained as a function of these two parameters. There is a series of marginally acceptable solutions (χ 2 {Nă9). The best fit is obtained for σ"25˝and y eq "1.4 ppb, with χ 2 {N"6.3, and the corresponding line area and residual maps are displayed in Fig. 7. We cannot obtain better overall fits with this model, because of a remaining excess of emission in the northwestern limb (and to a lower extent in the northeastern limb). Minimizing the residuals implies compensating partially for this excess, which results in too much emission in the southern hemisphere. Thus, we first try to improve the fit by reassessing our temperature field in the next section, since the temperature field at the time of the observations was influenced by the onset of Saturn's Great Storm of 2010-2011. Fig. 6. χ 2 {N as a function of the H 2 O meridional gaussian distribution parameters y eq and σ (gaussian HWHM in degrees). Acceptable solutions are found for σ ranging from "20˝to "40˝, with a minimum for σ"25˝and y eq "1.4 ppb. Corresponding line area and residual maps are shown in Fig. 7. As in the uniform distribution case, the HIFI data requires a "5 times higher y eq than the PACS data for σ"25˝to minimize χ 2 {N. The fit, shown in Fig. 8 left, is also only marginally acceptable. The line is too broad and not strong enough. This is an indication that models as simple as constant profiles above the condensation level are not optimal. The HIFI line-center contribution function peaks indeed at slightly lower pressures than the PACS line contribution function (see Fig. 8). Introducing a positive gradient in the H 2 O vertical profile above the condensation layer should thus (i) reduce the width of the HIFI line by removing some H 2 O in the levels just above the condensation layer, (ii) increase the HIFI line amplitude by adding some H 2 O at higher levels to which the HIFI line center is still sensitive, and (iii) help to reconcile the HIFI and PACS data as they probe different altitudes.
In an attempt to improve the fit to the PACS data and to improve the compatibility with the HIFI data, we apply a two-step approach, in which we first account for the temperature changes caused by Saturn's 2010-2011 Great Storm (Section 4.3), and then apply a parametrized gradient in the H 2 O vertical distribution above its condensation level (Section 4.4).   Fig. 5 for a gaussian distribution of H 2 O around the equator with y eq "1.4 ppb and σ"25˝. The overall fit is acceptable, even if the fit is not good at the northwest limb. A better fit in this region can be obtained by increasing y eq , but the model would then depart even more from the data in the southern hemisphere. It is already the case in this model, which is just the best balance between the northwest and southern emission.

Accounting for Saturn's Great Storm of 2010-2011
The Herschel observations were conducted a few weeks after the formation of a huge storm system around 40˝in the northern hemisphere of the planet (Sánchez-Lavega et al. 2011;Fischer et al. 2011;Fletcher et al. 2011). Between December 2010 and April 2011, it formed a cool region surrounded by two confined warm regions in Saturn's stratosphere, referred to as "beacons" (B1 and B2) given their appearance in the infrared combined with Saturn's fast rotation.
According to Cassini/CIRS observations made on the same day as the PACS data were recorded (and two days after the HIFI observations), the B1 and B2 beacons were located in the following longitude ranges: 300W-340W for B1, and 220W-250W for B2. They were thus both in the PACS field-of-view. Each one was located "halfway between the central meridian (289˝) and the planetary limb (west for B1 and east for B2). Only B2 was in the HIFI field-of-view, located around the central meridian (227˝). At that time, CIRS did only observe their northern edges and an increase of "5-10K had been measured there at the mbar level (Fletcher et al. 2012). Given the limited spatial coverage of the Cassini data (and the absence of ground-based thermal data) at that date, we have assumed a latitudinal extent of B1 and B2 from 30˝N to 50˝N (planetocentric), and tried several temperature increase combinations within B1 and B2 compared to adjacent longitudes. As the emission excess seen in the data is more obvious in the northwest than in the northeast, we assume a stronger temperature increase in B1 compared to B2. We have tested the following combinations:`6 K/`3 K, 10 K/`5 K,`15 K/`10 K, and`20 K/`15 K for B1/B2 (respectively), and applied these temperature increases uniformly above the 10-mbar level at all latitudes/longitudes in B1/B2, following the vertical trends in temperature found in these early development stages of B1 and B2 (see Fig. 8a in Fletcher et al. 2012). The choice of 10 mbar as a cut-off pressure for the temperature increase was guided by the following argument: this level is below the H 2 O condensation level, so that we do not probe this level, while it is still well above the level where the continuum is generated. It is therefore a pragmatic way to ensure that the temperature changes in the beacon will impact the H 2 O lines. For B1 and B2, we also account in the modeling of the thermal field for the longitudinal smearing of 15˝caused by the PACS integration time. Our thermal field thus becomes a 3D field in our subsequent modeling (see Fig. 9), and we also recompute the condensation level of H 2 O in the beacons according to this new thermal field (see Fig. 10).
We obtain the best-fit results for the Gaussian meridional distribution of H 2 O with uniform temperature increases in B1 and B2 of 10 K and 5 K, respectively. The resulting 3D temperature field is shown in several latitudinal cross-sections in Fig. 9. With this field, we significantly improve the fit to the PACS data and find a range of very good-fit (σ,y eq ) combinations, with a best fit (χ 2 {N"1.1) for σ"25˝and y eq "1.1 ppb. The line area and residual maps are displayed in Fig. 11 and the column density as a function of latitude is shown in Fig. 12. We stress that we did not change the H 2 O abundance within the beacons in this work (only the pressure of the condensation level). The results seem to be indicating some enhancement in H 2 O abundance from the excess of emission in the observations compared to our simulation. In addition, the temperature profile within the beacons was more complex than our idealized model, with a peak at 0.5 mbar and a drop at lower pressures (Fletcher et al. 2012). This should result in the need for more H 2 O in the mbar region to compensate for the fainter emission above the temperature peak. We leave this for future analysis of the beacon emissions.
With this model, we can even constrain a background and meridionally uniform flux (e.g., IDP) represented by a variable y min that can be added to equation 1. We find that it cannot exceed 0.06 ppb, at the 2-σ limit (Fig. 13), which is an order of magnitude lower than the disk-averaged contribution of the gaussian distribution. This confirms the small contribution from IDP, as predicted by Moses & Poppe (2017). New photochemical simulations are now required to determine the corresponding IDP H 2 O flux upper limit.
Accounting for the temperature increase in B2, which is in the HIFI field-of-view does not improve the situation regarding the PACS/HIFI incompatibility, i.e., "5 times more H 2 O than in the PACS best-fit model is needed to fit the HIFI line. It remains unchanged in terms of y eq for a given σ, compared to Section 4.2,

Contribution Function
HIFI best-fit model at 1097GHz -nadir PACS best-fit model at 66µm -nadir PACS best-fit model at 66µm -limb . Both functions are computed at their respective observed spectral resolution, and the contribution of the continuum emission has been subtracted. While the PACS line essentially probes near the condensation level, the HIFI line probes slightly lower pressures. We note that we have rescaled the HIFI contribution function by multiplying it by a factor of 1.45 for an easier comparison with the PACS one. and the fit to the HIFI data keeps the same flaws (wings too broad and line center too weak). Now that we have a good fit to the PACS data, we have to briefly turn back to the meridionally uniform distribution model to check whether the improvement in the temperature field induces any improvement in the fit to the data for such a thermal model. Fig. 14 shows the line area and residual map obtained for an H 2 O mole fraction of 4ˆ10´1 0 above the condensation level. With a χ 2 {N"16.2, such a model remains invalid. After leveling out the thermal effects from the line emission map, the uniform distribution still results in high latitudes being too bright and low latitudes being too faint compared to the observations. Accounting for the effect on the thermal field of Saturn's Great Storm of 2010-2011 has enabled us to significantly improve the quality of the fit to the PACS data with the Gaussian meridional H 2 O distribution and confirmed the invalidity of the meridionally uniform distribution.

Introducing a gradient in the H 2 O vertical profile above its condensation layer
In this section, we seek to reconcile the PACS and HIFI observations while building on the improvements in the PACS interpretation that resulted from including Saturn's Great Storm effect in our modeling, by introducing a gradient in the H 2 O vertical distribution above the condensation level.
As the HIFI line probes slightly higher than the PACS line, and as the HIFI line is too broad in the wings and too faint at line center, we introduce a parametrized vertical gradient in the H 2 O profile at all latitudes to make it less abundant above the condensation level and more abundant at higher altitude. This way, the empirical H 2 O vertical profiles will be closer to physical profiles computed with photochemical models (e.g., Ollivier et al. 2000;Moses et al. 2000). As the PACS data require less H 2 O than the HIFI data, we see this gradient as a means to reconcile the two datasets. To achieve this, we introduce two additional parameters: a cut-off pressure p gradient , above which the H 2 O abundance is held constant, and n which is the logpyq{ logppq slope of the profile between the cut-off pressure and the condensation level, as already used in previous works (Marten et al. 2005;Rezac et al. 2014). We test several combinations of p gradient and n, with values of p gradient ranging from 0.01 mbar to 1 mbar (p gradient is always smaller than the pressure of the condensation level, whatever the latitude) and values of n ranging from 0.5 to 3 (an example is shown in Fig. 15).
We find that adding a positive gradient in the H 2 O vertical profile above the condensation layer reduces the incompatibility between the HIFI and PACS y eq values from a factor of "5 down to a factor of "2.4 for n"2-3 and p gradient ă0.3 mbar in the range of values we have tested. An example for PACS is shown in Fig. 16 for n"2, p gradient "0.1 mbar, y eq "9ˆ10´8 and σ"25˝. The corresponding vertical profiles as a function of latitude are displayed in Fig. 15. The corresponding best fit of the HIFI line is obtained for the same set of parameters, except y eq "2.2ˆ10´7 (see Fig. 17).

Final remarks on the HIFI/PACS discrepancy
As it stands, we find no H 2 O distribution that enables us to fully reconcile the HIFI and PACS data in terms of H 2 O abundance (through the value of y eq ), even if introducing a positive gradient Article number, page 10 of 16 T. Cavalié et al.: Herschel map of Saturn's stratospheric water, delivered by the plumes of Enceladus Fig. 9. 3D thermal field used in simulations accounting for 10 K temperature increases in B1 and 5 K for B2. B1 is located between 300W-355W, and B2 between 220W-265W (when accounting for longitudinal smearing during the PACS integration), and both between 30 N and 50 N (Fletcher et al. 2012). The top panels are meridional cross-sections isolating B1 (left) and B2 (right). The bottom left panel is a zonal cross-section at 40 N, and the bottom right is a pressure cross-section at 2 mbar, in which both beacons can be identified. above the condensation level significantly improves the situation.
There are some systematic and random errors we have not considered that may explain at least in part the remaining discrepancy between the HIFI and PACS derived abundances. For instance, there is some scatter in the PH 3 disk-averaged abundances derived from observations (e.g., Weisstein & Serabyn 1996vs. Fletcher et al. 2012, its meridional distribution was only measured once by Fletcher et al. (2009a) and any temporal evolution is therefore unconstrained. As noted in Section 3.3, the PH 3 distribution influences the continuum of the 1097 GHz line (but not of the 66.4 µm one), and thus the derivation of the H 2 O abundance from this line, quite significantly. The thermal field also bears some uncertainties ("3 K at the tropopause and "4 K at 1-5 mbar, as detailed in Fletcher et al. (2018b)). Given that we implement a full 3D thermal field and that the opacities of the observed H 2 O lines are of the same order but not strictly identical, the uncertainties on the thermal field may be an additional cause of discrepancy between the HIFI and PACS lines.
Finally, we stress that the meridional distribution may be more complicated than the one given by Eq. 1, and the vertical gradient may also vary with latitude. This study can thus be seen as a stepping stone for future computations of physical profiles with more complete altitude-latitude photochemical models.

Discussion
In situ measurements carried out by several Cassini instruments during the end-of-mission orbits and final plunge into the atmosphere (all in 2017) have provided strong evidence for material infalling from the rings onto Saturn's atmosphere, either in gaseous or solid form, either neutral or ionized (Hsu et al. 2018;Mitchell et al. 2018;Waite et al. 2018). While small infalling icy grains, constrained to˘2˝around the equatorial plane, rain onto Saturn with a modest mass flux of 5 kg/s , H 2 O gas seems to be feeding Saturn's upper atmosphere at a much higher rate. From Cassini's Ion and Neutral Mass Spectrometer (INMS), Waite et al. (2018) found that a global integrated mass flux of 4800-45000 kg/s of gaseous material was raining onto the planet's atmosphere at the time of the measurements, with 24˘5% H 2 O, within a latitudinal band of 8˝centered on the equator. Using data from the same instrument, Perry et al. (2018) find an even larger influx of 2-20ˆ10 4 kg/s. This is 2 to 4 orders of magnitude higher than the influx of H 2 O of Article number, page 11 of 16 A&A proofs: manuscript no. Saturn_H2O_map_observations  Fig. 10. H 2 O vertical profiles of the PACS best fit model (see Fig. 11), in which the latitudinal distribution is a Gaussian centered around the equator with σ"25˝, and y eq "1.1 ppb. The effect of the temperature increases in the B1 and B2 beacons is accounted for when computing the condensation level, as shown by the yellow and cyan lines (corresponding to B1 and B2, respectively). The local equatorial minimum at 1 mbar results from condensation at a temperature minimum caused by the Saturn quasi-periodic oscillation (see Fig. 4). The corresponding line area and residual maps are shown in Fig. 11. 7-25 kg/s derived from previous disk-averaged observations and models Bergin et al. 2000;Moses et al. 2005;Hartogh et al. 2011). Therefore, this material infalling from the inner ring system, seen in 2017 by Cassini, cannot be the source of the H 2 O observed by Herschel in December 2010 and January 2011. An earlier Cassini finding provided another credible candidate: Enceladus and its H 2 O plumes (Hansen et al. 2006;Porco et al. 2006;Waite et al. 2006). Cassini's first observations of the plumes at the south pole of the moon venting out significant amounts of H 2 O triggered a modeling effort to evaluate the fate of the outgassed molecules (Jurac & Richardson 2007;Cassidy & Johnson 2010;Moore et al. 2010). For instance, Cassidy & Johnson (2010) found that neutral/neutral scattering was the main cause of spreading of the molecules in Saturn's system and predicted that Enceladus plumes formed an H 2 O neutral torus at the orbital distance of the satellite. This torus was confirmed by the observations of Hartogh et al. (2011). According to Cassidy & Johnson (2010), this torus is likely to be a significant source of H 2 O for Saturn's rings, the planet itself and even Titan. The latter was confirmed by Moreno et al. (2012). At Saturn, the meridional distribution of infalling neutral H 2 O originating from Enceladus plumes is predicted by Cassidy & Johnson (2010) to be Gaussian and centered on the planet equator with a HWHM σ"15˝. The disk-averaged mass flux was estimated to 6ˆ10 5 cm´2¨s´1, i.e., "8 kg.s´1, by Hartogh et al. (2011) by fitting the Enceladus source to the absorption caused by the Enceladus torus on the Saturn H 2 O line at 557 GHz. Comparing this 6ˆ10 5 cm´2¨s´1 flux to the required external flux invoked to explain Saturn's water (1˘0.5)ˆ10 6 (Moses et al. 2000), Hartogh et al. (2011) concluded that the Enceladus plumes were the likely source of Saturn's stratospheric H 2 O.
In this paper, we have analyzed a Herschel-PACS map of H 2 O emission coming from Saturn's stratosphere. We have presented clear evidence that H 2 O in Saturn's stratosphere is not meridionally uniform, and that a distribution in which its abundance peaks at the equator and decreases exponentially towards higher latitudes fits the data. This result enables us to identify Fig. 11. Same as Fig. 5 for a Gaussian distribution of H 2 O around the equator with y eq "1.1 ppb and σ"25˝, after accounting for the temperature increases in the B1 and B2 stratospheric beacons caused by Saturn's Great Storm of 2010-2011 at the time of the PACS observations. The overall fit is very good, even if some emission excess remains in the northwest limb. Increasing the temperature further in B2 would improve the fit in this region, but degrade the fit at the North Pole because of the large beam of PACS.
Enceladus as the main source of Saturn's stratospheric H 2 O, and to reject the IDP as a primary source. However, a faint meridionally uniform background flux cannot be ruled out by our observations. Our data indicate that this background is an order of magnitude lower than the disk-averaged contribution of the Enceladus source, in agreement with the prediction from Moses & Poppe (2017).
The difference we find between the predicted width of the influx meridional distribution with the observed one allow us to compute a rough estimate of meridional eddy mixing at lowto-mid latitudes. If we assume that K yy " L 2 {t, we then take L"10500 km as the distance between 15˝latitude and 25˝latitude. These latitudes correspond to the HWHMs of the input flux from Cassidy & Johnson (2010) and our best model. For the dif-Article number, page 12 of 16 T. Cavalié et al.: Herschel map of Saturn's stratospheric water, delivered by the plumes of Enceladus Fig. 12. Saturn's stratospheric H 2 O meridional distribution as derived from Herschel-PACS mapping observations on January 2, 2011: (left) mole fraction above the local condensation level as a function of latitude, and (right) corresponding column density as a function of latitude. The distribution is a gaussian centered on the equator with y eq "1.1 ppb and σ"25˝, corresponding to the PACS map best-fit model (see Fig. 11). Background image credits: NASA/JPL-Caltech/Space Science Institute. fusion timescale between those latitudes, we take t equal to the downward diffusion timescale from the top of the atmosphere, where the material is delivered, to the condensation level (150 y according to Moses & Poppe 2017), where H 2 O was mapped by Herschel. We find that K yy "2ˆ10 8 cm 2¨s´1 , and this value is nominally lower by a factor of 10 than early modeling results from Friedson & Moses (2011). A more accurate derivation of the input fluxes as a function of latitude at Saturn that are responsible for the observations will require additional simulation work with 2D photochemical models (Hue et al. 2015(Hue et al. , 2016(Hue et al. , 2018.

Conclusion
The findings of our paper can be summarized as follow:  (2017). -We can improve the compatibility between the results obtained from the HIFI and PACS observations by adding a positive gradient above the condensation level in the H 2 O vertical profiles, in an attempt to bring our profiles closer to physical profiles. However, we do not manage to fully reconcile the disk-resolved PACS data and the disk-averaged HIFI data with our empirical models, probably because of remaining systematic and random errors and because of the simplicity of our meridional distribution model. Full 2D photochemical modeling is now necessary to move to the next step.
It remains to be seen what meridional distribution of input fluxes are required to reproduce the Herschel observations, and 2D photochemical modeling is required at this stage, starting for example with those of Hue et al. (2015Hue et al. ( , 2016Hue et al. ( , 2018. Given their extraordinary sensitivity, bandwidth, and spectral and spatial resolutions, future complementary ALMA and JWST observations will certainly help to shed light on the recent extraordinary influx of exogenic material from the inner ring system to Saturn's stratosphere, as captured by Cassini in its final orbits Perry et al. 2018;Mitchell et al. 2018;Hsu et al. 2018). These future observations (Lellouch 2008;Norwood et al. 2016),   Fig. 16, in which the latitudinal distribution is a Gaussian centered around the equator with σ"25˝, and y eq "9ˆ10´8. The effect of the temperature increases in the B1 and B2 beacons is accounted for when computing the condensation level, and a positive gradient is added for pressures higher than p gradient "0.1 mbar and down to the condensation level, with a slope n"logpyq{ logppq"2. The value of y eq applies to păp gradient . Condensation between 0.1 and 1 mbar occurs in the equatorial region as a result of a local temperature minimum caused by the quasi-periodic oscillation (see Fig. 4). The corresponding line area and residual maps are shown in Fig. 16.
Chalmers University of Technology -MC2, RSS & GARD; Onsala Space Observatory; Swedish National Space Board, Stockholm University -Stockholm Observatory; Switzerland: ETH Zurich, FHNW; USA: Caltech, JPL, NHSC. The Herschel spacecraft was designed, built, tested, and launched under a contract to ESA managed by the Herschel/Planck Project team by an industrial consortium under the overall responsibility of the prime contractor Thales Alenia Space (Cannes), and including Astrium (Friedrichshafen) responsible for the payload module and for system testing at spacecraft level, Thales Alenia Space (Turin) responsible for the service module, and Astrium (Toulouse) responsible for the telescope, with in excess of a hundred subcontractors.
HCSS / HSpot / HIPE is a joint development (are joint developments) by the Herschel Science Ground Segment Consortium, consisting of ESA, the NASA Herschel Science Center, and the HIFI, PACS and SPIRE consortia.   . 17. Best fit to the HIFI data for a Gaussian distribution of H 2 O around the equator with y eq "2.2ˆ10´7 and σ"25˝, after accounting for the effect of the temperature increases in B1 and B2 and a positive gradient for pressures higher than p gradient "0.1 mbar and down to the condensation level with a slope n"logpyq{ logppq"2. The best-fit model to the PACS data from Fig. 16 is shown for comparison. Same layout as in Fig. 8 left.