Herschel^††thanks: Herschel is an ESA space observatory with science instruments provided by European-led Principal Investigator consortia and with important participation from NASA. map of Saturn’s stratospheric water, delivered by the plumes of Enceladus

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₂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$\times$5 spatial pixels ($\sim$50″$\times$50″ on sky), each covering a short instantaneous wavelength range sampled by 16 spectral pixels. Our observation of Saturn consists of a 3$\times$3 raster map with 3″ separation at 66.44 $\mu$m. We thus pointed the 25-pixel array over 9 different positions to record an oversampled 225-point map of Saturn at 66.44 $\mu$m. For the observation of such an intense far-infrared 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″ 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 $\sim$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″ in RA and 0.3″ 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 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₂O line is detected in all positions within the disk and on positions within $\sim$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-$\sigma$ rms of 0.0023 (in units of microns $\times$ % of the continuum) when adding quadratically the baseline removal uncertainty and the spectral noise. We use this value in our $\chi^{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 $\sim$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.

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₂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″ (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$\cdot$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 $\sim$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″ 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_{USB}(H+V)=0.469\pm 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 $\chi^{2}/N$ calculations are based on this rms value, which does not include absolute calibration uncertainties that are estimated to be 5%.

The pressure-latitude thermal field we use in our radiative-transfer 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 $\sim$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. (2009, 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}$.

In this paper, we test two types of H₂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₂O. The mole fraction of H₂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₂O, centered on the equator, with a peak H₂O mole fraction at the equator $y_{\mathrm{eq}}$, and a half-width at half-maximum (HWHM) $\sigma$, such that:

where $y_{\mathrm{eq}}$ is the H₂O mole fraction above the condensation level at the equator, and $\phi$ is the planetocentric latitude. Here also, latitudinally-dependent condensation is computed from our thermal field. The recently observed ring source (Waite et al., 2018; 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₂, not H₂O), another potential source for Saturn’s stratospheric H₂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₂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.

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₂-H₂, H₂-He, and H₂-CH₄ pairs are included following Borysow et al. (1985, 1988) and Borysow & Frommhold (1986). More recent updates to the H₂-H₂ 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₄ respectively, according to Conrath & Gautier (2000) and Fletcher et al. (2009b). We account for the absorption caused by the H₂O line, as well as surrounding NH₃ and PH₃ broad lines. We take the NH₃ and PH₃ distributions from Davis et al. (1996) and Fletcher et al. (2009a), respectively. We compute the relevant H₂/He broadening parameters from Levy et al. (1993, 1994) for PH₃, from Fletcher et al. (2007) for NH₃, and from Brown & Plymate (1996) and Dick et al. (2009) for H₂O, for the H₂/He mixture described above and the appropriate temperature range. The 1097 GHz line has an opacity of $\sim$0.3 over the HIFI beam in disk-centered conditions, and the 66.44 $\mu$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₂O profiles given in Fig. 15.

We note that, contrary to the 66.4 $\mu$m line seen by PACS, the 1097 GHz line seen by HIFI lies on the far wings of a PH₃ line. A change in the PH₃ abundance therefore influences the line-to-continuum ratio of the H₂O line and subsequently on the derived H₂O abundance. Increasing the PH₃ abundance by a factor of 2 decreases the H₂O abundance required to fit the line by $\sim$15%.

This model significantly improves on the model already described in Cavalié et al. (2008b, 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 line-of-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 $\mu$m by Spilker et al. (2003, 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 lines-of-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 disk-averaged 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¹¹1In the case of the HIFI spectrum, which was not absolutely calibrated, we have rescaled the spectrum to match our model continuum. The analysis is then done in terms of line-to-continuum ratio..

An initial approach was to try to fit the PACS map with a uniform distribution of H₂O which would be representative of an IDP source. As shown in Fig. 5, the H₂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 $\chi^{2}/N$ is obtained for an H₂O mole fraction of 6$\times 10^{-10}$ above the condensation level (roughly a few mbar, varying with latitude), but the $\chi^{2}/N$ is as high as 18.9 (i.e., more than 4-$\sigma$ 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₂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 $\sim$3.5$\times 10^{-9}$, which is consistent with Fletcher et al. (2012) but incompatible with the PACS data. When we simulate the H₂O line at 556.936 GHz with this distribution, we obtain a line-to-continuum ratio that is consistent with the previous observation of Bergin et al. (2000). This shows that the disk-averaged H₂O influx responsible for our observations of $\sim$6$\times 10^{5}$ cm${}^{-2}\cdot$s^-1 (Moses & Poppe, 2017), which translates into a disk-integrated mass flux of $\sim$8 kg/s, is surpassed by orders of magnitude by the extraordinarily high flux recently measured by INMS of $\sim$10⁴ 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₂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., $\sim$ 10-15 yr) to diffuse down to observable levels ($\sim$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.

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₂O in Saturn’s stratosphere.

In what follows, we test several cases of H₂O meridional distributions as parametrized in eq. 1. We first explore the $y_{\mathrm{eq}}$–$\sigma$ parameter space. Fig. 6 shows the $\chi^{2}/N$ obtained as a function of these two parameters. There is a series of marginally acceptable solutions ($\chi^{2}/N$$<$9). The best fit is obtained for $\sigma$$=$25^∘ and $y_{\mathrm{eq}}$$=$1.4 ppb, with $\chi^{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.

As in the uniform distribution case, the HIFI data requires a $\sim$5 times higher $y_{\mathrm{eq}}$ than the PACS data for $\sigma$$=$25^∘ to minimize $\chi^{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₂O vertical profile above the condensation layer should thus (i) reduce the width of the HIFI line by removing some H₂O in the levels just above the condensation layer, (ii) increase the HIFI line amplitude by adding some H₂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₂O vertical distribution above its condensation level (Section 4.4).

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 $\sim$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 $\sim$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₂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₂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₂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₂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 ($\sigma$,$y_{\mathrm{eq}}$) combinations, with a best fit ($\chi^{2}/N$$=$1.1) for $\sigma$$=$25^∘ and $y_{\mathrm{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₂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₂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₂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_{\mathrm{min}}$ that can be added to equation 1. We find that it cannot exceed 0.06 ppb, at the 2-$\sigma$ 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₂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., $\sim$5 times more H₂O than in the PACS best-fit model is needed to fit the HIFI line. It remains unchanged in terms of $y_{\mathrm{eq}}$ for a given $\sigma$, compared to Section 4.2, 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₂O mole fraction of 4$\times 10^{-10}$ above the condensation level. With a $\chi^{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₂O distribution and confirmed the invalidity of the meridionally uniform distribution.

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₂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₂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₂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₂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_{\mathrm{gradient}}$, above which the H₂O abundance is held constant, and $n$ which is the $\log(y)/\log(p)$ 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_{\mathrm{gradient}}$ and $n$, with values of $p_{\mathrm{gradient}}$ ranging from 0.01 mbar to 1 mbar ($p_{\mathrm{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₂O vertical profile above the condensation layer reduces the incompatibility between the HIFI and PACS $y_{\mathrm{eq}}$ values from a factor of $\sim$5 down to a factor of $\sim$2.4 for $n$$=$2–3 and $p_{\mathrm{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_{\mathrm{gradient}}$$=$0.1 mbar, $y_{\mathrm{eq}}$$=$9$\times 10^{-8}$ and $\sigma$$=$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_{\mathrm{eq}}$$=$2.2$\times 10^{-7}$ (see Fig. 17).

As it stands, we find no H₂O distribution that enables us to fully reconcile the HIFI and PACS data in terms of H₂O abundance (through the value of $y_{\mathrm{eq}}$), even if introducing a positive gradient 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₃ disk-averaged abundances derived from observations (e.g., Weisstein & Serabyn 1996 vs. 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₃ distribution influences the continuum of the 1097 GHz line (but not of the 66.4 $\mu$m one), and thus the derivation of the H₂O abundance from this line, quite significantly. The thermal field also bears some uncertainties ($\sim$3 K at the tropopause and $\sim$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₂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.