Simulating Jupiter’s weather layer. Part II: Passive ammonia and water cycles

We examine the ammonia and water cycles in Jupiter’s upper troposphere and lower stratosphere during spin-up of a multiple zonal jet global circulation using Oxford’s Jupiter General Circulation Model. Jupiter’s atmosphere is simulated at 512×256 horizontal resolution with 33 vertical levels between 0.01 and 18bar, putting the lowest level well below the expected water cloud base. Simulations with and without a 5.7 Wm 2 interior heat source were run for 130000–150000 d to allow the deep atmosphere to come into radiative-convective-dyna- mical equilibrium, with variants on the interior heating case including varying the initial tracer distribution, particle condensate diameter, and cloud process timescales. The cloud scheme includes simple representations of the ammonia and water cycles. Ammonia vapour changes phase to ice, and reacts with hydrogen sulphide to produce ammonium hydrosulphide. Water changes phases between vapour, liquid, and ice depending on local environmental conditions, and all condensates sediment at their respective Stokes velocities. With interior heating, clouds of ammonia ice, ammonium hydrosulphide ice, and water ice form with cloud bases around 0.4bar, 1.5bar, and 3bar, respectively. Without interior heating the ammonia cloud base forms in the same way, but the ammonium hydrosulphide and water clouds sediment to the bottom of the domain. The liquid water cloud is either absent or extremely sparse. Zonal structures form that correlate regions of strong latitudinal shear with regions of constant condensate concentration, implying that jets act as barriers to the mixing. Regions with locally high and low cloud concentrations also correlated with regions of upwelling and downwelling, respectively. Shortly after initialisation, the ammonia vapour distribution up to the cloud base resembles the enhanced concentration seen in Juno observations, due to strong meridional mean circulation at the equator. The re-semblance decays rapidly over time, but suggests that at least some of the relevant physics is captured by the model. The comparison should improve with additional microphysics and better representation of the deep ammonia reservoir. called Part I hereafter), we examined the spin-up of a global banded jet structure with super-rotating equatorial flow in simulations of Jupiter’s weather layer using a new generation of Oxford’s Jupiter General Circulation Model (GCM) based on the


Introduction
Jupiter's vertical cloud structure has been deduced primarily from 1D equilibrium cloud condensation models (Lewis, 1969;Weidenschilling and Lewis, 1973;Sánchez-Lavega et al., 2004;Wong et al., 2015) combined with ground-based and spacecraft observations (Banfield et al., 1998;Gierasch et al., 2000;Matcheva et al., 2005;Reuter et al., 2007;de Pater et al., 2016;Li et al., 2017), and the Galileo entry probe (Ragent et al., 1998). While CH 4 is the most abundant molecule in Jupiter's atmosphere after H 2 and He, it does not condense at Jupiter's temperatures and pressures. Instead, condensates of ammonia (NH 3 , as solid and aqueous solution), water (H 2 O, as liquid and solid) are expected, along with solid ammonium hydrosulphide (NH 4 SH) formed by the reaction between gaseous NH 3 and H 2 S.
Pre-Juno, the uppermost observed clouds were solid NH 3 around 0.75-0.25 bar (Banfield et al., 1998), while the presence of H 2 O clouds at depth ( > 3 bar) was inferred by observations of highly convective storms , and lightning on Jupiter's night side (Little et al., 1999), although direct evidence for both H 2 O and NH 4 SH has been sparse (Carlson et al., 1996;Simon-Miller et al., 2000;Matcheva et al., 2005). The Galileo probe measured a much lower H 2 O concentration than expected, although the probe entered a region of the atmosphere that was anomalously bright around 5 μm, indicative of a region devoid of clouds in the upper troposphere.
In Young et al. (2018b, called Part I hereafter), we examined the spin-up of a global banded jet structure with super-rotating equatorial flow in simulations of Jupiter's weather layer using a new generation of Oxford's Jupiter General Circulation Model (GCM) based on the MITgcm. This model of Jupiter's upper troposphere and lower stratosphere has been developed over the last 15-20 years, originally based on the dynamical core of the UK Met Office External Unified Model (ExtUM) (Yamazaki et al., 2004;Yamazaki and Read, 2006). Zuchowski et al. (2009b,c) used it to develop a simple cloud parameterisation over a limited area of the southern hemisphere, and a single-column moist convection parameterisation (Zuchowski et al., 2009a).
The aim of this paper is to study the distribution of condensate clouds under various conditions using the new generation of the Jupiter model. This paper focuses on passive tracers, which are advected by the flow and may be transformed into other tracers, but which do not feed back onto the flow's evolution via radiative effects, latent heating, or moist convective processes. In a forthcoming paper we will include an active water cycle with feedbacks from latent heat release and moist convection; the developments and experiments presented in this paper and in Part I will put those results into context.
The cloud scheme in the new generation of our Jupiter model is described in Section 2, followed by a description of the simulations in Section 3. We investigate the cloud dynamics in 4, and conclude in Section 5. Appendix A contains a complete description of our Jupiter cloud scheme, and Appendix B defines the buoyancy frequency accounting for molar mass gradients.

Model description
The model is based on the atmospheric mode of the MITgcm (Marshall et al., 1997a;1997b;Adcroft et al., 2004;, with several Jupiter-specific parameterisations added. It is a new generation of the model initially developed by Yamazaki et al. (2004Yamazaki et al. ( , 2005; Yamazaki and Read (2006); Zuchowski et al. (2009aZuchowski et al. ( , 2009bZuchowski et al. ( , 2009c, which used the UK Met Office Unified Model, and is fully described (apart from the cloud scheme, presented below) in Part I.
The MITgcm solves the primitive equations for an incompressible fluid on a rotating sphere. In the Jupiter model we use the MITgcm's hydrostatic atmospheric mode with an ideal gas equation of state, longitude-latitude Arakawa-C grid, vector-invariant form of the momentum equation, centred second order advection-diffusion for potential temperature, with a 4th-order Shapiro filter, zonal filter poleward of 45°latitude, omitting the vertical Coriolis terms, linear free surface, and with no harmonic or bi-harmonic viscosity.
The model's horizontal domain covers the whole globe with resolution 512 × 256 (called resolution M hereafter) with a = t 300 s time step. There are 33 vertical levels, approximately equally spaced in log pressure such that the lower face of the lowest level is at there. The lower boundary is free-slip. The Jupiter parameterisations are fully described in Part I. In summary, we include a two-stream radiation scheme that simulates incoming shortwave solar heating (as a function of latitude based on the annual mean solar flux), emission and absorption of longwave heating from the atmosphere, and an optional interior heat source. A vertical diffusion term parameterises small scale turbulent mixing such as Kelvin-Helmholtz and symmetric instabilities. Columns unstable to dry convection are adjusted towards a neutrally stable profile over a prescribed timescale. We include linear drag at the bottom model level to ensure numerical stability, and a sponge layer acting on velocity and temperature eddies in the top three levels, to damp wave reflection from the top of the domain.
The cloud scheme that is the focus of this paper adds a 4D tracer advection equation is the total derivative, τ is a tracer mass mixing ratio, t is time, G tracer is a tracer tendency from sources and sinks, = u v v ( , , 0) h is the horizontal velocity on pressure surfaces, where u and v are the zonal and meridional velocities respectively, ∇ p is the gradient operator in pressure coordinates, and = Dp Dt / is the vertical velocity in pressure coordinates (positive downwards). Tracer advection is performed using a second order flux limiter scheme.
The model carries seven substances as passive tracers transported by model winds: solid and gaseous (vapour) ammonia (NH 3 ), gaseous hydrogen sulphide (H 2 S), solid ammonium hydrosulphide (NH 4 SH), and gaseous, liquid, and solid water (H 2 O). These represent the two main condensible cycles in Jupiter's atmosphere: the ammonia cycle and the water cycle. The tracers are passive in the sense that the mass of the condensate does not contribute to the equation of state, and the various transformation processes do not have any thermodynamic feedback onto the flow. An active cloud scheme that parameterises latent heat release and moist convective processes will be the focus of an upcoming paper.
The cloud scheme parameterises several processes, which are resolved by the model as tendencies on the right hand side of the tracer advection equation Eq. (1). Each component of the cloud scheme is described in full in Appendix A. NH 3 can change phase between vapour and solid (ice) depending on the local saturation mass mixing ratio (SMMR) (Sánchez-Lavega et al., 2004). Its gaseous form can also react with H 2 S to form solid NH 4 SH ice (Lewis, 1969;Weidenschilling and Lewis, 1973). H 2 O undergoes phase changes between its vapour, liquid, and solid forms, depending on the local saturation mass mixing ratio (to determine whether any condensate is formed) and the local melting point (to determine whether liquid or solid condensate form). All four condensates also fall through the atmosphere following Stokes' law, representing condensate sedimentation (Böttger, 2003), but there is currently no representation of precipitation. This scheme updates the cloud parameterisation in Zuchowski et al. (2009c) for the new generation of Oxford's Jupiter GCM, reconfiguring it for the global context, recasting the phase changes as finite-time tendencies rather than instantaneous adjustments, and adding the reaction to form NH 4 SH.

Simulations run
We ran several simulations with different configurations of interior heating and cloud parameters. Table 1 lists the parameters for each run.
The two base runs A and B are the same runs as described in Part I. They were initialised from rest, with no clouds, and with a temperature field generated by a radiative-convective model (RCM) containing only the radiation and dry convection part of the code. The only difference was that Run B had a 5.7 W m 2 interior heat flux and Run A did not. The flow was initially spun up at half the horizontal resolution in each direction (called resolution L hereafter) with a 600 s time step, until equilibrated (when the total kinetic energy and top-of-atmosphere net radiative flux had stabilised). Then the resolution was doubled to resolution M and the time step halved to 300 s, and it was run to equilibrium again, which took until 147500 d and 128700 d for Runs A and B, respectively.
Only at this point was the cloud scheme turned on. It was not turned on from the start in order to minimise the computational time required for the flow to reach equilibrium, because advecting seven tracers reduces the model speed by about 70% while not affecting the flow as the tracers are passive. We ran the model for an additional several thousand days. We found during testing that only several thousand days were required for a tracer distribution to reach equilibrium; they were considered equilibrated once the average rate of change of total tracer mass for each tracer approached zero. Several versions of the more realistic interior heated case B were run, all started at the same point in Run B, with different configurations of the cloud scheme. These variants are listed in Table 1; they altered the initial tracer profile (Runs B2-B5), the condensate particle diameter (Run B6), and the ratio between the NH 3 phase change and the NH 4 SH reaction timescales (Runs B7 and B8).
The tracers were initially set either to a fixed multiple of the solar abundance (Runs B, B3-B8), or to abundances measured by various spacecraft and ground-based observations (Run B2). Fig. 1 shows these initial profiles. Only the gaseous tracers were initialised; all the condensate fields began empty. In all cases the horizontal symmetry of the initial vapour concentrations was broken by applying uniformly distributed multiplicative noise with amplitude 0.1%. If a multiple of the solar composition was used, the mass mixing ratio was calculated from the solar abundance, and was the same mass mixing ratio at all levels unless a cut-off pressure was specified. Most of the observational abundances we used are from Irwin (2009, Table 4.6), except H 2 S above 0.1 bar (Owen et al., 1980), and H 2 O and H 2 S below 3.6 bar (Niemann et al., 1998). The observed H 2 O profile we used is the depleted profile measured by the Galileo probe. At the time the simulations were run the new observations of ammonia abundances observed by Juno were not yet available .
We used 64 and 128 cores respectively for the resolution L and M parts of the runs. Running on 128 cores, where about 50% of linear speedup was possible, the full model with tracers runs about 0.24 simulated days per core-hour, or 42000 core-hours per 10000 d.

Results
Part I described the flow in both runs in some detail. The distribution of condensate clouds primarily depends on the vertical temperature structure, which differs considerably between Runs A and B, but also to some extent on the choices made in the cloud scheme. In Fig. 2 we reproduce some of Part I Fig. 4, which shows the properties of the flow most important for our discussion of the clouds. In both cases global bands of zonal jets form. Run A has a wide sub-rotating equatorial jet, and Run B a thinner super-rotating equatorial jet with regular off-equatorial jets that migrate towards the equator. In both cases the tropopause is around 0.1 bar. In Run A the atmosphere is strongly convectively stable at the base of the model, while in Run B the atmosphere below about 0.3 bar is close to a state of neutral stability with respect to the dry adiabat.
For continuity with Part I, we analyse the clouds using the same times as the analysis in that paper, i.e., at the end of the simulation for Run A, and at 133000 d for Run B and its variants, unless indicated otherwise. Fig. 3 shows saturation vapour pressure curves for ammonia and water phase changes, and the equilibrium reaction constant for ammonium hydrosulphide formation, using the global mean temperature at the point during Runs A and B when clouds were added, assuming solar mixing ratios for all trace species. This predicts that NH 3 (s) will form around 500 mbar when there is interior heating and slightly lower without interior heating, NH 4 SH(s) should form around 1.5 bar and H 2 O(s) around 3-4 bar with interior heating, but that the water condensate is likely to form too high for liquid water. It also predicts that without interior heating both NH 4 SH(s) and H 2 O(s) should form throughout the model.  The total tracer mass tendencies are shown in Fig. 7. Both ammonia and water cycles equilibrate within several thousand days, with the ammonia cycle equilibrating faster, after about 3000 d compared with 6000 d for water. In the ammonia cycle the tendency for NH 3 (s) is much smoother than for the other three species. This is because in Run B the NH 4 SH(s) reaction formation timescale (2400 s) is shorter than the ammonia phase change timescale (10000 s).  . Saturation vapour pressure curves for ammonia and water (left) and equilibrium constant for the NH 4 SH formation reaction (right). Grey and black lines show global mean temperature profiles (left) and equilibrium constants for the formation of NH 4 SH at those temperatures (right, ln (K reac ), Eq. (A.7)) for Runs A and B, respectively. Where these are to the left of the coloured lines, condensate can form. In the left panel, the red and blue lines show the saturation vapour pressure curves ϵp*/q (Eq. (A.4) divided by the mole fraction) for well-mixed solar abundances of nitrogen and oxygen, at the start of Runs A and B. The line changes from blue to cyan where the temperature favours a liquid rather than a solid water condensate. In the right panel, the green line shows the sum of log partial pressures of NH 3 (g) and H 2 S(g) at the beginning of the runs (Eq. (A.7)). The temperature profiles used in these calculations are from the start of the period with tracers. (For interpretation of the references to colour in this figure legend, the reader is referred to the web version of this article.) NH 3 forms a cloud base around 500 mbar, NH 4 SH cloud forms around 1.5-2 bar, and a H 2 O cloud forms around 4 bar (Figs. 4 and 6a). These are at approximately the altitudes predicted by the saturation vapour pressure curves (Fig. 3). The condensate is largely concentrated at the equator and near the poles, there is a latitudinal variation related to the zonal jets, and the cloud deck at the equator is quite deep for both NH 4 SH and H 2 O, extending almost a decade in pressure.
Away from the poles, the cloud concentrations are generally highest where the temperature is highest (Fig. 2b), with a maximum concentration at the equator. The cloud base is higher at warmer latitudes (i.e., nearer the equator), which is the expected behaviour based on the saturation vapour pressure curves (consider moving the B lines in Fig. 3 to the left as temperature decreases-the cloud base falls). Near the equator both NH 4 SH(s) and H 2 O(s) clouds are quite deep, extending up to around 300 mbar at the equator itself; all three cloud decks have their tops at the equator at around the same pressure. This is likely due to the meridional circulation at and near the equator, which is about two orders of magnitude stronger at the equator in the troposphere than at ± 10°latitude, and extends up to around 300 mbar, so the flow will lift much of the condensate up to this altitude.
The zonal mean condensate structure is modulated in latitude by the zonal jets. Fig. 5 shows horizontal sections through the three cloud decks at the pressures where their concentrations are highest. This shows a latitudinal modulation of the cloud concentrations according to the position of the zonal jets. There is a correlation between the latitudes where there is large zonal wind shear, i.e., where |∂⟨u⟩/∂y| is large, and latitudes where the cloud concentration is approximately constant. Latitudes with local maxima in concentration ( Fig. 4) are typically in regions of anticyclonic shear flow (∂⟨u⟩/∂y > 0 in the northern hemisphere), and latitudes with local minima in concentration are typically in regions of cyclonic shear. Conversely, at the peaks of the zonal jets the latitudinal gradients in cloud concentration are highest.
This correlation is strong in all three condensates, but breaks down near the poles and in the equatorial jet. This pattern is consistent with a region of strong mixing where the zonal shear is largest, homogenising the cloud condensate concentration, while preventing mixing across the latitudes where the zonal shear is small, i.e., at the jet peaks, which act as a barrier to mixing.
In Fig. 8 we zoom into a small region of the atmosphere (a 15°by 15°horizontal square). On a more local scale than the zonal jets (i.e., on the scale of individual eddies), there are correlations between condensate concentrations and the prognostic variables. First, as at the jet scale, at a local scale there is also a correlation between regions of strong vorticity (either negative or positive) and regions where the condensate concentration is approximately constant, due to turbulent mixing. Conversely, where vorticity is close to zero there are gradients in the condensate concentration.
Second, there is also a correlation between high condensate concentrations and regions where the vertical velocity is upwards, and vice versa for regions of low condensate concentration. In equilibrium, air above the condensate clouds will be strongly depleted in both vapour and condensate tracer, because any condensate will sediment downwards into the cloud layer, and the amount of vapour that remains decreases exponentially fast with falling pressure. Bulk upwards motion will bring vapour up from below, which subsequently condenses, while regions where the bulk motion is downwards bring air depleted in condensate down from above.
Finally, there is a positive correlation between the eddy temperature and the condensate concentration. This is clearer for NH 3 (s) and NH 4 SH (s) than for H 2 O(s). This might appear paradoxical given that higher temperatures decrease condensate production by increasing the saturation mass mixing ratio and the NH 4 SH equilibrium constant. However, note that this correlation is in the eddy fields, where the variability is very weak ( ∼ 0.1 K) compared with the latitudinal gradient in (full) temperature across the panel, which is some ten times larger. Note also that there is a correlation between the vertical velocity and the eddy temperatures. Therefore the correlation between eddy temperature and concentration is because both are correlated with vertical motions. Warm, condensate-rich air is brought up from below, and cool, condensate-depleted air is brought down from above.
There is a peak in tracer concentration near the poles. This is particularly clear in the NH 3 (s) case (Fig. 4b), where the maximum in the tracer concentration is around ± 85°, and a high concentration of NH 3 (s) extends right to the top of the domain in the southern hemisphere. The atmosphere is coldest in this region (Fig. 2b), so we might expect some condensate in the polar stratosphere, but probably not as much as is observed. The polar stratosphere has a much weaker meridional circulation than the tropical troposphere, and the weak temperature inversion in the polar stratosphere will suppress vertical motion. Therefore any condensate (and vapour) entering this region may not be able to escape.
However, we expect this result is partly due to numerics, although the exact reason is uncertain. There is a strong temperature inversion at the model's uppermost level (Fig. 3), which was described in Part I. This causes this region to become a tracer sink, as temperatures there are typically above the solidus and hence no condensate will form. As tracer accumulates there it will not be advected out (as the region is strongly convectively stable) nor will it sediment out (as the temperature is too high for condensate to form). However, this does not explain why this occurs near the poles only. We were unable to test whether the strong temperature inversion is responsible as we could not remove the strong temperature inversion at the top of the model in the MITgcm's current configuration.
It is also unclear why the same does not occur for NH 4 SH(s) and H 2 O(s), although note both do have slightly deeper clouds right at the pole ( Fig. 4c and d). It is possible that because the cloud base is lower for these two species, within the convectively neutral part of the atmosphere, that it is not prevented from circulating by a strongly convectively stable temperature profile, as NH 3 (s) is in the upper atmosphere.

Distribution of condensates with no interior heating
In the case with no interior heating, Run A, the temperature inversion in the lowest parts of the domain (Fig. 2b) means that H 2 O(s) clouds condense throughout the domain, and at the lower temperature than Run B the equilibrium between NH 4 SH(s) and NH 3 (g)/H 2 S(g) is strongly in favour of the condensate (Fig. 3).
Both these condensates subsequently sediment out over the several thousand days of the simulation towards the bottom of the model, and do not re-sublime as they fall, unlike in Run B. As a result they both accumulate below 10 bar (we do not include this cross-section in Fig. 4). Fig. 9 shows the Stokes velocity (Eq. (A.12)) and the estimated time for each condensate to fall to the bottom of the model domain, for both Runs A and B. The Stokes velocity is dominated by the Knudsen number term above about 0.1 bar, and by Stokes drag below. The lower temperatures in Run A produce faster Stokes velocities in the lowest part of the domain. As a result the time for condensate to fall from the top of the model domain to the bottom is about 5000 d for NH 4 SH(s) and 7000 d for H 2 O(s), and there is enough time for both to sediment out completely. The condensates in Run B, however, re-sublime before they reach the bottom of the model.
The saturation vapour pressure curve for NH 3 (s) means it condenses in the upper troposphere (Fig. 3), as in Run B. As predicted by the saturation vapour pressure curves, the NH 3 (s) clouds in Run A are at a lower level than Run B, due to the lower temperatures. Fig. 6a shows the global mean condensate profiles in both Runs A and B, where these differences can be seen clearly.

Varying the initial tracer distribution
Runs with different variations on the default cloud scheme are presented in Figs. 6, 10, and 11. These variations generally have a second-order effect on the distribution of condensates. Fig. 10 shows cases varying the initial conditions (Runs B2-B5). This has a noticeable effect on the equilibrium state (Figs. 6b and 10), although advection of the tracers by the flow does remove a lot of these initial differences. Much of the latitudinal banding remains.
When the simulation starts with the observed profiles (Fig. 10a-c), the NH 3 (s) cloud deck is some 10 times denser than the control case Run B, at all latitudes (although it does not extend to the top of the model as in Run B). The NH 4 SH(s) cloud is similar to Run B. This is due to the increased initial NH 3 (g):H 2 S(g) ratio in the observed case over the 1 × solar case at the pressures where NH 3 (s) is expected to There was no liquid water condensate. Supplementary Data S1 shows an animation of the water ice column density (i.e., integrated over the depth of the model) over 200 d during Run B. condense (Fig. 1). In the 1 × solar case the initial ratio of mass mixing ratios is 2.2:1, but in the observed case in the upper troposphere the ratio is almost 400:1. Hence the trade-off between the NH 3 (s) phase change and the NH 4 SH(s) formation reaction favours NH 3 (s) formation. Note also that NH 3 (s) and NH 4 SH(s) are both depleted in the stratosphere relative to Run B; this is due to the much reduced initial concentrations of the gaseous species compared with the 1 × solar case, so there is much less vapour to turn into condensate.
The water condensate profile is similar for all initial conditions except for 2 × solar. Most of the condensate is near the equator, with some at higher latitudes. As the amount of water decreases the total amount of condensate decreases proportionally (see dot-dashed blue line in Fig. 6b relative to the solid line), but the fraction at high latitudes also decreases. When 2 × solar water abundance is used there is a sink of water condensates (both solid and liquid in this case) at the water cloud base around 50°S. We think this is probably a numerical artefact, as it is also the only case in which there is substantial liquid water condensation. In all our test runs of this scheme we found that liquid water very seldom condensed (very occasionally at the model level below the solid water cloud). However, it was very intermittent in space and occasionally caused the water ice cloud to become unstable, with a large amount of water going into one of the condensed phases and not being able to re-evaporate. As the liquid water cloud only occurred very rarely this was not a major problem, but in this particular case it seems to have collected all the water condensate in one place.
For the 0.5 and 1 × solar (control) cases, the H 2 O(s) cloud is thicker than the NH 3 (s) and NH 4 SH(s) clouds, even at the peak of the NH 3 (s) and NH 4 SH(s) condensate distributions (Fig. 6a). Given that Jupiter's atmosphere at the observed cloud level is dry, this implies that our thick water cloud is over-estimating the amount of water in the upper troposphere. In fact only the observed profile (based on the Galileo entry probe and generally considered to underestimate the global water abundance in Jupiter's atmosphere) yields a NH 3 (s) cloud thicker than the H 2 O(s) cloud at the altitude where ammonia condenses. Even the case where we restricted the initial water abundance to be solar below 4 bar and zero elsewhere, the model advects enough water to higher altitudes to make the water cloud thicker than the ammonia cloud at the ammonia cloud base. This suggests that the model overestimates the meridional circulation near the equator, bringing more water to higher altitudes than occurs on the real planet.
Restricting the initial condition to tracer below a specific pressure level (Fig. 10f-h) does not change the final distributions very much, as most of the initial tracer is in the deep atmosphere anyway, and the tracers are advected throughout the atmosphere well enough for the lack of initial tracer above the predicted cloud base to be overcome. The overall concentrations are slightly lower (Fig. 6b, dashed line) as there is just not as much tracer in the atmosphere. Note, however, that the strong peak in NH 3 (s) near the poles is not present in this case (Fig. 10f), suggesting further that this has something to do with the temperature inversion in the uppermost model layer.

Varying the condensate particle diameter
Figs. 6c and 11 show the case with 10 μm particles (Run B6). Varying the condensate particle diameter only affects the sedimentation process.
We experimented with particle sizes to see how large a particle was required for condensate to sediment out of the atmosphere faster than it was replenished, depleting the condensate cloud decks, and a radius of 10 μm was large enough. Fig. 9 shows the Stokes velocity and time to fall to 18 bar for this run. The Stokes velocity is some 100 times larger than the other runs, and hence the time to sediment out is only around 100 d. This is reflected in the zonal cross-sections (Fig. 11), where the cloud decks are depleted everywhere, and in the global mean (Fig. 6c), which is depleted by 1-2 orders of magnitude.
Near the poles there is still some condensate, albeit only 10% of the concentration in Run B. The cloud base is slightly lower than in the nominal case Run B, as the Stokes velocity varies as d p 2 in the lower atmosphere, so particles will fall 100 times faster and hence will fall farther before they are re-sublimed.
Note that there is a sparse but non-zero H 2 O(l) cloud formed in Run B6 (Fig. 6c). This may occur if the H 2 O(s) condensate falls far enough that the temperature goes above the melting point before re-subliming. As the condensates are able to fall farther when the particle diameter is larger, this is able to occur in this case.

Varying the ammonia process timescales
Varying the ratio between the NH 4 SH reaction timescale τ reac and the NH 3 (s) phase change timescale τ pchg (Runs B7-B8) had little visible effect on the condensate fields in zonal cross-section, and so these are not shown.
There are trends in the global mean profile (Fig. 6d) at altitudes where the condensate concentrations are small, however. As τ reac increases, there is a small increase in NH 3 (s) in the stratosphere at the expense of NH 4 SH(s), as the NH 3 (s) formation process occurs faster and so uses up excess vapour first. NH 4 SH(s) also increases below 2 bar, because as τ reac increases, sedimenting NH 4 SH(s) is able to fall farther before it decomposes into NH 3 (g) and H 2 S(g).
We noted earlier that in Run B the NH 3 (s) tracer mass tendency was less noisy than the other species in the ammonia cycle (Fig. 7), because τ pchg < τ reac . As τ reac increases in Runs B7 and B8 (figures not shown), the time series for the other species involved in the NH 4 SH(s) formation reaction become less noisy, although parity between the two processes is only reached in Run B8.

Comparison with Juno vapour mixing ratios
Recent microwave radiometer observations of Jupiter's weather layer by Juno  have shown a heterogeneous distribution of ammonia vapour in both latitude and altitude, and while reasons have been suggested , a definitive explanation will require further study. The equatorial jet between 0-5°N is strongly ammonia-rich relative to the surrounding latitudes, and ammonia-poor between 2-10 bar. We do not expect our model to match these observations precisely, but can make some general comparisons with the new data.  Fig. 4), shows concentrations below the ammonia condensation level for the three vapour species in our model. It shows Run B2, where simulations were initialised with the (pre-Juno) observed profiles. For all three vapour species the equatorial regions have enhanced mixing ratios relative to other latitudes, and relative to the initial concentrations in Fig. 1. In this region the vapour is well mixed (i.e., the mixing ratio is approximately the same at all altitudes up to the cloud base), unlike in the initial profiles. Poleward of this region there is a latitude band where the vapour is depleted by over 50%, and the vapour is not well Fig. 7. Total tracer mass tendencies (summed over the whole atmosphere) for each species during Run B. There is one point every (Earth) day, and the grey line shows the 100 d running mean. H 2 O(l) is omitted as it did not form at any point during this run. Note that all the condensate concentrations (lower row) start at zero but their tendencies are generally negative; this is because during the first few timesteps (before the first point plotted in these figures) the tendency is positive and 3-4 orders of magnitude larger for a short period (and correspondingly negative for the gaseous tracers), as all the vapour in the upper atmosphere rapidly converts to condensate.
mixed, more like the initial state. The region of enhancement at the equator is wider than in the Juno observations, extending latitudinally over more than just the equatorial zonal jet.
Examining the vapour evolution with time reveals what is happening. The meridional circulation at the equator is strong (Part I, Fig. 4) compared with other latitudes. During the first few hundred days after the clouds are initialised, shown in Fig. 13, a well-mixed region of enhanced mixing ratio develops in a thin band at the equator and, at this point, appears to be similar to the Juno observations in some respects.
Over time, more vapour is lifted at the equator, the vapour already aloft spreads latitudinally by turbulent mixing, and the weaker midlatitude meridional circulation cells (Part I, Fig. 12) slowly lift vapour from the deep atmosphere at other latitudes. So what we see is probably a transient towards a well-mixed vapour distribution. Some of the enhanced mixing ratios at higher latitudes are also likely to be due to resublimation of condensate. Note, for example, how the maximum in H 2 S at 60°N/S (Fig. 12) matches up with the edge of the high latitude NH 4 SH cloud in Fig. 10b.
Over much longer periods, a permanent enhanced vapour mixing ratio at lower latitudes and depleted mixing ratio at higher latitudes might establish itself. Fig. 4 in Part I showed that, while the meridional mass circulation is primarily on the scale of individual jets, there is also a weaker, hemisphere-scale latitudinal gradient in the meridional circulation, which may generate a latitudinal contrast over time, but not over the small length scales seen by Juno.
Nevertheless, this shows that the model's strong equatorial  Fig. 2. The colour scales in the three cases using observations are stretched relative to the Run B plots in Fig. 4.   Fig. 11. Zonal-time mean cloud concentrations for 10 μm particle size (Run B6). Concentrations are averaged over the same time period as for Run B in Fig. 2. The colour scales for NH 4 SH and H 2 O have been scaled by a factor 0.1 relative to the Run B plots in Fig. 4. meridional circulation can generate, at least temporarily, enhanced mixing ratios at the equator, leading to latitudinal contrasts in vapour mixing ratios. Over short time periods this can generate a remarkably Juno-like distribution, although we do not claim that this is the only mechanism at work in those observations. The observed depletion at higher latitudes is not something that the model is expected to sustain.
The new observations show that there are some (as yet unknown) relevant physical or chemical processes missing from the model that create this heterogeneity in the vapour mixing ratios. This is not surprising, as we do not claim to include all of Jupiter's relevant physical processes in this model. Yet it is encouraging that the model can spontaneously produce some aspects of the observed tracer distribution.

Conclusions
In this paper we have shown how the water and ammonia cycles respond to the spin-up of a global jet structure using the new generation of Oxford's Jupiter GCM described in Young et al. (2018b, Part I), a new generation of the model using the MITgcm. We used it to investigate the distribution of cloud condensate with and without an interior heat flux, and with a range of cloud scheme parameters. The cloud scheme was updated to work with the global model, and the reaction between NH 3 (g) and H 2 S(g) to produce NH 4 SH(s) was added. This version does not include the moist convective parameterisation developed by Zuchowski et al. (2009a), whose effects on the global dynamics will be presented in a forthcoming paper.
In the interior heating case, condensate clouds of NH 3 (s), NH 4 SH(s), and H 2 O(s) were found near altitudes predicted by the saturation vapour pressure curves. H 2 O(l) was very sparse, only appearing when the initial water abundance was high or the particle diameter was large. The vertical extent of the H 2 O(s) cloud in particular was large, to such an extent that would be clearly detectable were it to occur on the real planet. In the case without heating from below there is a temperature inversion near the bottom of the model, and so NH 4 SH(s) and H 2 O(s) condense throughout the domain and fall to the bottom of the model. Zonal bands and more local regions of horizontal shear were correlated with regions of approximately constant cloud concentration, implying that jets (where the zonal shear is weakest) are barriers to mixing. In general, regions of high and low cloud concentrations also correlated with regions of upwelling and downwelling, respectively.
Certain aspects of the cloud scheme remain to be improved in future work. The model's vertical resolution under-resolves the cloud scale height by about a factor 1.5 (based on data from Sánchez-Lavega et al., 2004, Table III). However, we also exclude precipitation, so have deeper cloud decks than in reality, and so our modelled clouds are better resolved than this. There is a trade-off between affordability and physical reality given the simplicity of the scheme. To resolve, say, three levels per cloud scale height would require a 4-5 increase in runtime, or less for a targeted increase in vertical resolution at predicted cloud levels.
There is currently no contribution from the tracer species to the equation of state via the air density. The effect of this on static stability, which is close to zero in the deep atmosphere [see Part I Fig. 6], can be estimated using the expression in Appendix B. Fig. 14 shows the predicted contribution of molar weight gradients to the buoyancy frequency and static stability. The static stability increases by 0.01-0.02 K km 1 between 5 and 0.2 bar. Most of the difference is due to water, even the peak at 400 mbar near the ammonia cloud level, as this is due to the drop off in water ice at that pressure. This increase is enough to stabilise the deep atmosphere above the water cloud base, but it is not as stable as the Galileo probe measurements (Magalhães et al., 2002). We note that other authors (e.g., Lian and Showman, 2010) have transported vapour species only, which may be more stable, but introduces other concerns such as replenishment of the water reservoir via a deep source.
Condensate can become trapped in certain regions of the atmosphere in both the NH 3 (s) and H 2 O(l) tracers. Some 1D experiments we conducted showed that if the model top is not at the vacuum then overheating of the uppermost layer may be avoided (Mendonça et al., 2015, also note this). The effects of cloud opacity on solar and infrared absorption are not currently accounted for, and are a necessary advance for more realistically simulating radiative heating in the cloud layer. They would be critical if, for example, the model were extended to simulate Saturn's stratosphere, with its dense haze.
Since these simulations were run, the Juno spacecraft  has challenged our assumptions about the ammonia distribution Fig. 12. Vapour volume mixing ratios (q/ϵ) during Run B2 for (from top) ammonia, hydrogen sulphide, and water. Run B2 was initialised from (pre-Juno) observations (solid lines in Fig. 1). Plots show zonal means averaged over the same period as for Run B in Fig. 2. The figures have been drawn to be comparable with Li et al. (2017, Fig . 4), but note deep ammonia concentrations are slightly higher here due to the higher initial deep concentration in the pre-Juno data (600 ppm below 8.6 bar). We also show a wider latitude and smaller pressure range. in Jupiter's weather layer Ingersoll et al., 2017). Evidently, some relevant physical processes are currently missing from the model. These observations, along with future measurements of the deep water abundance, will allow us to better constrain tracer abundances below the visible clouds. The relevant dynamics could then be better examined by a focused study where the deep volatile abundance is used to initialise the model, for example, with the observed zonal jet structure imposed. This would also allow a deeper domain to be used closer to the ∼ 100 bar reached by the Juno microwave radiometer. The time to reach radiative equilibrium for the deeper domain would be a problem to overcome, but by running with only a handful of zonal grid points to radiative equilibrium, or by increasing both the solar and internal heating by a factor 10, it might be possible to reduce it.

Accompanying data
Data from 154760-154860 d in Run A and 132900-133000 d in Run B can be obtained from the Oxford University Research Archive-Data (https://ora.ox.ac.uk) (Young et al., 2018a).
where the three terms on the right hand side refer to phase changes (A.1), the reaction to form NH 4 SH(s) (A.2), and sedimentation (A.3), respectively.

A1. Phase changes
In each grid cell the saturation mass mixing ratio (SMMR) is calculated for NH 3 (gas-solid) and H 2 O (gas-solid and gas-liquid) transitions  Bulk composition of Jupiter's atmosphere as well as the physical properties of NH 4 SH (not used in the bulk atmosphere calculation). Some properties of the tracer species used in the cloud scheme are also listed. The bulk molar mass M b =2.33 g mol 1 is a weighted average (by bulk mole fraction) of the six most abundant bulk atmosphere species (excluding NH 4 SH), and so should be within 0.01% of the true value. The bulk molecular radius = r 119 bulk × 10 12 m is a weighted average of the H 2 and He Van der Waals radii. All quantities are from Lide (1995), except X and f (Irwin, 2009, Table 2.2a and b), ρ p and C p for NH 3 (s) (Blum, 1975;Overstreet and Giauque, 1937), and other heat capacities (Lide, 2004 Young, et al. Icarus 326 (2019) 253-268 (Sánchez-Lavega et al., 2004).
where p* is the saturation vapour pressure (SVP) of the gas, a function of latent heat (and hence temperature) derived from the Clausius-Clapeyron equation where R v is the specific gas constant for the tracer, and L is the latent heat of the phase change, which is constant or linear in T. In general, p* is given in Pa by using the parameters in Table A.2 (this general form is given here due to the various definitions of the SVP used in the source material).
The total mass mixing ratio for NH 3 or H 2 O (vapour plus any condensates) is calculated: If this exceeds the relevant SMMR then the model forces the vapour concentration towards the SMMR and the condensate concentration towards the excess. If the total mass mixing ratio is less than the SMMR then the vapour is forced towards the total and the condensate towards zero. The timescale τ pchg for these processes is a tuneable parameter, set to 10000 s in all the runs (following Lian and Showman, 2010). For water phase changes, the relevant SMMR curve depends on whether the temperature of the cell is above or below the melting point of water. This is calculated as a linear function of log p between -9°C at 1000 bar and 0°C at 1 bar (Lide, 1995, p.6-53), and between there and the triple point of water at 611.73 Pa (Lide, 1995, p.6-11). Below the triple point only water ice can form. In cells where the temperature is above the melting point of water, water ice is always forced towards zero regardless of the SMMR, and vice versa for liquid in cells below the melting point.
In summary, the tracer tendency is where q eqm is given in Table A.3 for each process.

A2. NH 4 SH formation reaction
In Jupiter's atmosphere ammonia and hydrogen sulphide gases may react to form ammonium hydrosulphide ice: This process is also parameterised by the model. This process and the NH 3 phase change are in competition for gaseous NH 3 . Tendencies from both these processes are calculated once per time step, but the timescales are set separately. We assume that the reaction timescale τ reac is shorter than the phase change timescale, so set the timescale for the NH 4 SH reaction as short as possible. However, we must make sure that the ammonia gas concentration does not go negative when both processes are extracting gas from the vapour reservoir, and after testing 2400 s turned out to be the shortest practical timescale. Runs B7 and B8 change the relative timescales for these two processes. The reaction scheme calculates the equilibrium between the three species, then uses the difference between the current composition and the equilibrium to determine a forcing tendency for the tracer advection equation.   Young, et al. Icarus 326 (2019) 253-268 In each grid cell, first we calculate the potential partial pressures of NH 3 and H 2 S that would occur if all NH 4 SH currently in the cell were converted back into vapour. NH 4 SH will form when the combined partial pressures of NH 3 and H 2 S gases exceed the equilibrium constant (Lewis, 1969, Eq. ( 25) (Lewis, 1969, to 5 s.f., and converted to correct units). 1 In terms of the mass mixing ratio, the partial pressure = pp qp/ , where ϵ is the ratio between the tracer and bulk molar mass ( In grid cells where the partial pressures of the gases are not high enough to form NH 4 SH, the equilibrium partial pressure of NH 4 SH is set to zero and the equilibrium partial pressures of the gases are set to the potential partial pressures calculated at the start of the process. The equilibrium partial pressures are converted back to equilibrium mass mixing ratios, and a tendency term to the tracer advection is then calculated by forcing the current concentrations towards the equilibrium, given by

A3. Sedimentation
This parameterisation calculates sedimentation tendencies for the liquid and solid tracers using Stokes' law and a Van Leer method for sedimentation. The scheme was first used in our planetary modelling by Böttger (2003) and Lee (2006). Stokes law for drag on a sphere is used to calculate the sedimentation rate of the liquid and solid species, using the Cunningham slip factor to take account of kinetic effects at low pressure. Sedimentation is performed using the Van Leer I method described by Hourdin and Armengaud (1999).
First we calculate the vertical sedimentation speed of the tracer particles using Stokes' law for resistance to a moving sphere to calculate the terminal sedimentation velocity. The vertical velocity due to Stokes' law is (Lee et al., 2010, Eq. (2)) = V d g C µ 18 p p n stokes 2 (A.12) where d p is the diameter of the condensate particle, ρ p is its density, μ is the dynamic viscosity of the bulk atmosphere, g is the gravitational acceleration, and C n is the Cunningham slip factor (Cunningham, 1910;Lee et al., 2010, Eq. (4)) = + C K 1 4 3 n n (A.13) K n is the Knudsen number, the ratio of the mean free path in the bulk atmosphere λ p (Chapman and Cowling, 1970, Eq. (4.21.3)) to the condensate particle radius: where r bulk is the radius of molecules in the bulk atmosphere and k B is the Boltzmann constant. Many of these quantities used the data listed in Table A.1. The dynamic viscosity of the bulk atmosphere was calculated using Brokaw (1968, Eq. (10)) for a mixture of 86.2% H 2 , 13.5% He, and 0.2% CH 4 by mole fraction, which accounts for 99.9% of Jupiter's bulk atmospheric composition. Using gas viscosities for these species from Lide (1995, p.6-243) we get These viscosities were fit to a 3rd order polynomial, assuming no variation with pressure. The equivalent column mass of bulk atmosphere transferred downwards by sedimentation during one time step is then where δt is the model time step and ρ is the atmospheric density. The Stokes velocity and mass transfer are calculated on the bottom and top of each cell face. Then, using a (mass conserving) Van Leer I scheme (Hourdin and Armengaud, 1999) with slopes limitation (maximum slope = 2), we calculated how tracers are redistributed among model levels by sedimentation. There is no sedimentation through the upper and lower boundaries of the model domain. Finally, as the model requires a tendency rather than an adjustment, this is finally converted to a tracer mass mixing ratio tendency for each condensate in terms of the tracer in each cell before and after sedimentation:

A4. Adjustments
Three additional steps are done as adjustments at the end of each time step, rather than as forcing of the tracer advection equation. First, each of the tracer fields has the MITgcm's zonal filter applied to it. Second, each field is searched for negative values, which appear after the zonal filter at points where there is no tracer, because of the Gibbs phenomenon. Where negative values are found, cells with positive tracer mass are scaled so they are reduced by exactly the right amount to fill in the global negative tracer mass: where M and + M are the global negative and positive tracer masses respectively. Finally, to conserve tracer mass over long simulations, every 24 h the global mass of each tracer is rescaled to its initial global mass calculated at the start of the run. Note, however, that for the ammonia cycle we want to conserve both the total amount of NH 3 and the total amount of H 2 S separately. But there are three rescaling factors (for NH 3 , H 2 S, and NH 4 SH separately), but only two constraints (the masses of NH 3 and H 2 S at the beginning of the run). This is therefore an under-constrained problem. Instead, we conserve the combined total mass of the three species as a whole, but this means the NH 3 :H 2 S ratio could drift during a run.
The model was programmed to stop should this rescaling be larger than 1 part in 100 over a single 24 h period. In practice, this was never triggered, and the true correction was typically 1 part in 10 10 , 7 8 showing that tracer mass is conserved well by the MITgcm. By the end of each run the NH 3 :H 2 S ratio had typically drifted by 1 part in 10 5 , and at most by 1 part in 1000.

Appendix B. Contribution to buoyancy frequency from molar mass gradients
Following Durran and Klemp (1982), the buoyancy frequency taking into account tracer mass, but excluding latent heat contributions from saturated air, is and other quantities are as defined in Part I

Supplementary material
Supplementary Data S1 shows an animation of the water ice column density in Run B in equilibrium over 200 d. Supplementary material associated with this article can be found, in the online version, at doi:10.1016/j.icarus.2018.12.002.