Longitudinal Filtering, Sponge Layers, and Equatorial Jet Formation in a General Circulation Model of Gaseous Exoplanets

General circulation models are a useful tool in understanding the three dimensional structure of hot Jupiter and sub-Neptune atmospheres; however, understanding the validity of the results from these simulations requires an understanding the artificial dissipation required for numerical stability. In this paper, we investigate the impact of the longitudinal filter and vertical ``sponge'' used in the Met Office's {\sc Unified Model} when simulating gaseous exoplanets. We demonstrate that excessive dissipation can result in counter-rotating jets and a catastrophic failure to conserve angular momentum. Once the dissipation is reduced to a level where a super-rotating jet forms, however, the jet and thermal structure are relatively insensitive to the dissipation, except in the nightside gyres where temperatures can vary by $\sim 100\,\mathrm{K}$. We do find, however, that flattening the latitudinal profile of the longitudinal filtering alters the results more than a reduction in the strength of the filtering itself. We also show that even in situations where the temperatures are relatively insensitive to the dissipation, the vertical velocities can still vary with the dissipation, potentially impacting physical processes that depend on the local vertical transport.


INTRODUCTION
General circulation models (GCMs) are an important tool in understanding the three-dimensional nature of hot Jupiter atmospheres given their potential for fast winds and large day-night temperature gradients.GCMs are faced, however, with the issue that simulations often require the inclusion of artificial dissipation to suppress gridscale fluctuations in order to maintain numerical stability, potentially altering the results, limiting our ability to reproduce observations, and taking the simulations in unphysical directions.Furthermore, there are physical dissipation mechanisms within hot Jupiter atmospheres that are ignored (e.g., gravity waves, shear instabilities; Watkins & Cho 2010;Fromang et al. 2016;Carone et al. 2020) or only approximated (e.g., magnetic drag; Perna et al. 2010;Rauscher & Menou 2013;Beltz et al. 2022) leading to a situation where the inclusion of a parametrised dissipation mechanism, representing the entirety of the unmodelled physics, may well be required to reproduce observed wind speeds even in high resolution simulations, as discussed in Heng et al. (2011) and Li & Goodman (2010).
Understanding the intrinsic dissipation in our models and the impact of dissipation included for numerical stability therefore becomes important as we seek to disentangle the dissipation required by unmodelled physics from the dissipation added due to numerical requirements.Intercomparisons of GCMs comparing idealised test ★ E-mail: christie@mpia.decases, such as Polichtchouk et al. (2014), can provide robust tests of individual components of a model while intercomparisons of specific planetary targets, such as THAI (Turbet et al. 2022;Sergeev et al. 2022;Fauchez et al. 2022), CAMEMBERT (Christie et al. 2022a), and MOCHA (Iro et al., in prep), can provide insights about the model variability across the range of GCMs in the presence of many coupled physical processes.These sorts of investigations, however, are often done for fixed, specific choices of the dissipation, as their intention is to look for differences between models, not within the models themselves.Some investigations of dissipation with these models do exist, however.The impact of dissipation order and effective resolution was investigated by Skinner & Cho (2021) within their pseudospectral code, concluding that an effective grid of ∼ 2000 × 1000 with a ∇ 16 diffusion operator is required to achieve numerical convergence within their model.Heng et al. (2011) also investigated the role of dissipation as a part of a larger investigation of multiple dynamical cores used within the dual-grey model Flexible Modelling System (fms) GCM.They found that variation between solution method and level of dissipation resulted in uncertainties in the wind speeds on the order of tens of percent, concluding that observational data are required to constrain the appropriate level of dissipation.While not varying the dissipation required for numerical stability, Koll & Komacek (2018) did investigate the role of explicit drag, approximating unmodelled physics, in the determination of the jet speed and examined the transition from Rayleigh drag dominating the dissipation of kinetic energy to the regime where numerical dissipation dominates.Hammond & Abbot (2022) examined the stabilising hyperdiffusion in the thor GCM and found that changes in diffusion order can result in zonal velocities changing by 50% with lower-order hyperdiffusion resulting in a slower equatorial jet.
In this paper, we examine the role of the longitudinal filter and the vertical sponge used within Met Office's Unified Model (UM) to stabilize simulations of hot Jupiters and sub-Neptunes, using a model of a 10× solar metallicity WASP-96b previously investigated by Zamyatina et al. (2024) as a test case.This paper is structured as follows: in Section 2, we outline the numerical setup used with the results presented in Section 3. A summary of the conclusions is presented in Section 4. In Appendix A, we present the eddy momentum flux gradients for each of the simulations presented in the paper, and in Appendix B, we present results from a short simulation with an increased horizontal resolution.

NUMERICAL SETUP
In this section we outline the details of the Met Office's unified model (UM) GCM as used in simulations of gaseous exoplanets, the details of the modeled planet WASP-96b, and the outline of the parameter studies that have been conducted.

The Unified Model
The UM, through its Endgame dynamical core, solves the full, deepatmosphere, non-hydrostatic Navier Stokes equations (Mayne et al. 2014b;Wood et al. 2014) and has been adapted to model hot Jupiters (e.g., Mayne et al. 2014a;Amundsen et al. 2016) and sub-Neptunes (e.g., Drummond et al. 2018;Mayne et al. 2019).The radiative transfer is handled by the Socrates radiative transfer code based on Edwards & Slingo (1996) which was adapted to the H 2 and He dominated atmosphere appropriate to hot Jupiters and sub-Neptunes in Amundsen et al. (2014Amundsen et al. ( , 2016Amundsen et al. ( , 2017)).Within the configuration adopted here, previously used in Christie et al. (2021), the included opacity sources are H 2 O, CO, CH 4 , NH 3 , Li, Na, K, Rb, Cs and H 2 -H 2 and H 2 -He collision-induced absorption (CIA).Opacities are computed using the correlated-k method and ExoMol line lists (Tennyson & Yurchenko 2012;Tennyson et al. 2016), and the chemical abundances are computed using the analytic model of Burrows & Sharp (1999) which was extended by Amundsen et al. (2016) to include alkali species.While the UM is capable of more sophisticated chemical modelling (e.g., Drummond et al. 2018;Zamyatina et al. 2023Zamyatina et al. , 2024)), as the purpose here is to investigate the impact of dissipation choices, the analytic model is sufficient.
To achieve stability, the UM implements two methods of dissipation and diffusion in the hot Jupiter/sub-Neptune setup.To limit reflections at the upper computational boundary, a vertical "sponge" is employed which damps vertical velocity  near the boundary, The height-dependent dissipation coefficient  sp is given by where  is the vertical coordinate with  = 0 at the base of the computational domain,  = / top ,   is the dimensionless height at which the sponge layer begins relative to the top of the top of the computational domain  top , and  is a coefficient regulating the strength of the vertical sponge 1 .As the height at which the sponge begins is measured relative to the location of the upper boundary, the effect of changing   can be accomplished by raising or lowering the upper boundary, subject to issues of numerical stability.The use of a sponge to limit reflections at the upper boundary is common in studies of both planetary and exoplanetary atmospheres, although the details vary between specific implementations (see e.g., Carone et al. 2020;Deitrick et al. 2020).
The second form of dissipation is a longitudinal filtering of the horizontal wind u h = (, ), where , , and  index the longitudinal, latitudinal, and vertical zones, respectively, and the filtering constant  is given by and is parametrised by   which is the number of applications of the operator needed for a grid-scale perturbation to decay by a factor of .The first order Shapiro filter corresponds to the limit of   → 0, or, equivalently,  = 1/4.For many hot Jupiter simulations performed using the UM (e.g., Mayne et al. 2017;Lines et al. 2018;Christie et al. 2021), the value of   is taken to be 1; however, it has been found that this level of dissipation can, in certain cases greatly influence the dynamics (e.g., Christie et al. 2022b;Zamyatina et al. 2024).This choice of filtering is equivalent to a longitudinal diffusion operator applied to the horizontal wind, where  is the longitude.The effective longitudinal diffusion coefficient  eff is given by where  is the latitude with  = 0°being the equator, Δ is the longitudinal grid spacing, and Δ dyn is the dynamical timestep.We note that as this filter operates exclusively on the horizontal velocity without a corresponding transport of gas density, neither the linear nor the angular momentum are conserved throughout this process, except in the exceptional case where there are not longitudinal gradients in gas density.While applying the filter instead to the horizontal momentum u ℎ would improve conservation of axial angular momentum, it would not necessarily achieve the desired result of suppressing grid-scale perturbations in horizontal velocities as only perturbations in momentum would be filtered.Higher order 1 Drummond et al. (2018) implemented a polar profile to the sponge wherein the vertical coordinate  used in determining  sp has a latitudinal dependence, where  is the latitude with the equator located at  = 0.This change has the effect of causing the sponge layer to descend to lower altitudes near the poles, with   only representing the height at which the sponge layer begins at the equator.Both the horizontally uniform and the polar sponge profile have been used in simulations of hot Jupiters and mini-Neptunes.
approaches can avoid this issue (e.g., Deitrick et al. 2020); however, this is beyond the scope of this paper.This choice of stabilizing dissipation differs from other GCMs in that it is a lower-order operator and applied only to the horizontal winds, and only in the longitudinal direction.While the focus of the operator is to suppress grid-scale perturbations and noise, it also has the effect of introducing a latitudinally-dependent diffusion acting on the horizontal winds due to the filter acting on a latitude-longitude grid.
Unlike other studies (e.g., Liu & Showman 2013;Carone et al. 2020;Schneider et al. 2022), studies using the UM have not included a "Rayleigh" drag at the inner boundary.As with the other types of dissipation, this paramtrisation can be included purely for reasons of numerical stability (Carone et al. 2020;Schneider et al. 2022), or it could represent a real physical mechanism such as magnetic drag (e.g" Liu & Showman 2013), although these parametrisations remain poorly constrained.Studies employing this interior drag could remove spurious accelerations in the deep atmosphere, but could also artificially impact the angular momentum exchange between the deep and upper atmosphere (e.g., Liu & Showman 2013;Carone et al. 2020, also see Schneider et al. 2022, Appendix A).

WASP-96b
For the purposes of investigating the impact of the longitudinal filtering and damping on the dynamics, we opt to model the atmosphere of WASP-96b.While the polar filtering can be seen to impact models of other targets resulting in counter-rotating jets, such as with the GJ 1214b 100× solar metallicity,  sed = 0.1 case in Christie et al. (2022b), WASP-96b was seen to show this behaviour (see Zamyatina et al. 2024) that can be easily reproduced without requiring additional computationally-intensive physics (clouds, non-equilibrium chemistry, etc.).
WASP-96b's metallicity is taken to be 10× solar metallicity, with the abundances of all species other than H 2 and He having their abundances relative to H increased by a factor of 10, with the elemental abundances taken from Lodders et al. (2009).The analytic chemistry formulation of Burrows & Sharp (1999) is found to be sufficient to reproduce the issue of counter-rotating jets and a breakdown in the conservation of angular momentum, so while Zamyatina et al. (2024) employs equilibrium and non-equilibrium chemistry schemes, we deem those unnecessary for this study and opt for computational simplicity.
The simulation parameters are presented in Table 1.The simulations are initialized in a wind-free hydrostatic state using a 1D profile generated by the radiation-convection code atmo (Tremblin et al. 2015;Drummond et al. 2016).As the flows found in GCMs can heat the deep atmosphere through transport processes unmodelled in 1D (Tremblin et al. 2017), and as it is prohibitive to run these models until convergence in the deep atmosphere is achieved, we note that there will be a temperature inversion in the final state.This can be avoided by using a "hot start" (i.e., using a profile with the tempreature profile artifically shifted to hotter temperatures; Amundsen et al. 2016;Tremblin et al. 2017); however, as our interest is in investigating behaviour seen in Zamyatina et al. (2024), we opt for the initialization used there.

The Parameter Study
To investigate the impact of the various dissipation mechanisms, we perform simulations of WASP-96b with values from   = 1 to   = 12 (equivalent to  = 0.158 to  = 0.019).This range encompasses the default value used in previous studies using the UM (  = 1; see, e.g., Mayne et al. 2019;Christie et al. 2021Christie et al. , 2022b;;Zamyatina et al. 2023), the value used in Zamyatina et al. (2024) (  = 6), and the lowest dissipation rate for which the model ran stably for the desired run time of 1000 Earth days (  = 12; for brevity going forward days will refer to Earth days).While we focus on a sponge depth of   = 0.75, as has been used in the previously mentioned papers, we include cases with a raised sponge where   = 0.9, with   = 1, 3, and 6.For the case of   = 12 and   = 0.9 the simulation did not run stably.The impact of the sponge could also be investigated through varying the sponge constant .Qualitatively, we would expect similar results, as both parameter studies serve to reduce the vertical damping in parts of the atmosphere; however, reducing  explicitly reduces the vertical damping at the outer boundary, unlike altering   , which potentially impacts the effectiveness of the sponge in preventing reflections at the upper boundary.Because of this, we opt to focus on an investigation varying   .We also note that resolution is expected to influence the dynamics of the jet, both in terms of the implicit dependence within the filtering operator and the numerical limitations of the dynamical core, but also due to the dynamics of smaller scales being resolved with increased resolution impacting the overall results.While this provides another avenue of investigation and has been explored by Skinner & Cho (2021) and Sainsbury-Martinez et al. ( 2019), we opt to investigate the impact of other parameters instead.

RESULTS
We present here the results from the suite of simulations investigating the impact of different aspects of the dissipation within the UM hot Jupiter model.Unless otherwise specified, the analysis is done after the pressure, temperature, and velocities have been averaged over the final 100 days of the simulation.In cases where the vertical coordinate is mapped from a height-based grid to a pressure-based grid, the temporal averaging occurs before the mapping.

The Longitudinal Filter
The motivation for this investigation is a failure in the model to conserve axial angular momentum, resulting in the formation of counter-rotating jets, as seen in Christie et al. (2022b) andZamyatina et al. (2024).As the primary mechanism for suppressing grid-scale perturbations within hot Jupiter simulations done with the UM is the longitudinal filtering of the horizontal velocities and given that that filtering does not necessarily conserve angular momentum, we first examine modifications to the application of the filter and their impact on the jet.

Reducing the Filter Strength
The most direct change that can be made to the filtering is through the e-folding parameter   .For the default value   = 1 as well as   = 2, a super-rotating jet forms in the first 100 days; however, this jet quickly erodes, resulting in the majority of the atmosphere counter-rotating, with an equatorial counter-rotating jet (see Figures 1 and 2).The flow in these cases forms a pair of dayside vortices on either side of the counter-rotating jet (see Figure 2, two leftmost panels), very different from the expected formation of a super-rotating equatorial jet expected from theoretical analysis (Showman & Polvani 2011;Tsai et al. 2014) and seen in GCM simulations of close-in gas giants (e.g., Showman et al. 2008;Menou 2012;Mayne et al. 2017).The dayside flow morphology is reminiscent of the linear flow in Showman & Polvani (2011); however, unlike in that analysis, the nightside is not symmetric.These two counter-rotating cases show significantly worse conservation of axial angular momentum relative to the remaining cases (see Figures 11 and 12).For the   = 1 case specifically, the simulation loses approximately 3% of its initial axial angular momentum during the 1000 day run.A similar breakdown in angular momentum conservation was seen in the 100× solar metallicity,  sed = 0.1 model of GJ 1214b in Christie et al. (2022b).The rate at which angular momentum is lost during the simulation is not linear in   (or ), as seen in Figure 12, hinting that the loss process might be self reinforcing.These unphysical cases also show discrepant behaviour in their kinetic energy evolution (Figure 13) with the kinetic energy increasing much faster than the cases with super-rotating jets, resulting in the kinetic energy being a factor of ∼ 2 larger for these cases compared to the super-rotating cases.Unsurprisingly, these cases show very different thermal structures (see Figures 8 and 9) with less heating deep in the atmosphere and, due to the reversal in the jet, warmer morning terminators and cooler evening terminators.As these cases exhibit obvious numerical breakdowns and unphysical results, we feel they are easily dismissed and do not dwell on them further. 2 2 The   = 1 and 2 cases do exhibit a pair of dayside vortices reminiscent of the "modons" found in the simulations of Skinner et al. (2023) (cf., their Figure 6c; see also ? and Skinner & Cho 2021); however, as the formation of the vortices in their simulations is robust to changes in dissipation, and given In contrast, for the more weakly damped cases (  = 3, 4, 6, and 12), a zonal jet is seen to form (Figure 1, rightmost three panels) with the expected flow morphology (Figure 2, rightmost three panels).These cases also show improved conservation of the axial angular momentum (Figures 11 and 12).For comparison, the   = 12 case loses 0.03% of its initial angular momentum over the 1000 day run, a factor of 100 improvement over the   = 1 case.The peak zonal velocities approach 6 km s −1 , with the   = 3 and 4 cases having the slowest zonal wind at 5.8 km s −1 .The counter-rotating regions at low pressures near the poles seen in Figure 1 are not counter-rotating jets, but are instead due to the retrograde velocities in the nightside gyres being larger than the prograde velocities on the dayside, as can be seen in Figure 2. Similarly, in the   = 6 and   = 12 cases where the super-rotating region extends to the poles at  ∼ 0.1 bar, this is due to the slower velocities in the gyres and not a jet that extends around the planet.The centers of the gyres also move towards the poles as suppression of longitudinal variations is reduced.
At the equator, the vertical profiles of the local zonal wind (Figure 3) do not show significant variation with   with the peak of the profiles roughly correlated with the bottom of the sponge.At a latitude of 45°(Figure 4), however, there is a larger variability in the local zonal wind with   , typically ∼ 0.5 km s −1 , with the largest dependency on   seen in the anti-stellar and morning profiles.At the mid-latitudes, unlike at the equator, the peak of the profiles do not correspond to the bottom of the sponge and show a more complex, longitudinally-dependent vertical behaviour 3 .The meridional velocities (Figure 5) show some variability at a latitude of 45°, especially in the anti-stellar profile where, for example, at a pressure of 1 mbar the   = 3 and   = 12 velocities differ by 0.5 km s −1 , in part due to changes in the direction -northward to southward -occurring at different pressures.The local dependence of the meridional velocity on   , however, varies greatly with the choice of latitude.A similar examination of the meridional velocities as in Figure 5, except at a latitude of 60°(not shown), finds improved agreement in the antistellar profile but larger variability with   in the profile along the morning terminator.In general, greater variability with   is seen in the horizontal winds outside of the equatorial jet compared to variability of the zonal velocities within the jet.
The equatorial vertical velocities for these cases are shown in Figure 6.The velocities within the equatorial jet do not show indication of agreement between 0.1 and 1 bar with the vertical velocities in the   = 12 case being up to ∼ 2× larger than the vertical velocities in the   = 6 case.This region is just above the temperature peak seen in Figure 8 and shows the impact of horizontal filtering on the vertical transport of potential temperature (Tremblin et al. 2017;Sainsbury-Martinez et al. 2019).Were the initial pressure-temperature profile used within the simulations to be closer to the final, converged state, or if a "hot start" were used (Amundsen et al. 2016;Tremblin et al. 2017), it is unclear as to whether these differences in this region would have been seen.The vertical velocities at the equator also show indications of not being converged in the morning profiles between 10 −2 and 10 −3 bar with differences in this region between the   = 6 and   = 12 cases being up to ∼ 20%.
The pressure-temperature profiles at the equator and at a latitude of 45°are shown in Figures 8 and 9, respectively.For the cases where   = 3, 4, 6, and 12, the temperature profiles show good agreement  For the more heavily damped cases (  = 1 and 2), a super-rotating jet does not develop.
Horizontal wind speed at 1 mbar averaged between days 900 and 1000 for the   = 0.75 cases.The substellar point is located at a longitude of 180°, the center of each panel.The grey streamlines indicate the direction of the flow.For the more heavily damped cases (  = 1 and 2), a super-rotating jet does not develop, instead forming a fast counter-rotating jet.
except at pressures ≲ 10 −2 bar on the nightside and morning terminator.This region corresponds to the nightside gyres which, due to the lack of external radiative forcing and the longer dynamical timescales, are impacted more by the small velocity perturbation allowed by the weaker filtering.As discussed above, the temperatures at ∼ 1 bar also increase with reduced filtering as the transport of potential temperature into the deeper atmosphere is enhanced.While we focus on   as the natural parametrisation for the dissipation from a numerical point of view, this approach is resolution dependent (see Eqn. 7) with  eff being the dynamically relevant quantity.In Appendix B we test that a sufficient increase the horizontal resolution, corresponding to a reduction of  eff , does result in the formation of a super-rotating jet for case of   = 1.

Flattening the diffusion profile
As, by design, the longitudinal filtering acts on the grid without consideration of the metric, the effective diffusion coefficient varies with latitude (see Equation 7) which has the potential for impacting the equatorial jet.To investigate this possibility, we run a simulation with   = 0.75 and change the filter coefficient K to partially account for the implicit latitudinal gradient: where  0 = 57.92°.This choice corresponds to   = 3 in the polar caps || >  0 while holding the effective diffusion coefficient ( eff ) constant in the equatorial region − 0 <  <  0 , resulting in   = 12 at  = 0°.While it may be possible to run the simulation with a larger region of constant  eff , this choice restricts the filtering to values that have been investigated and shown to form a super-rotating jet, facilitating the interpretation of results as it is less likely that any differences will be due to extreme choices in the local filtering.
The zonal velocities are shown in Figure 14 (rightmost panel).The peak zonal velocity is 6.79 km s −1 , faster than for the   = 3 (5.87km s −1 ) and   = 12 (6.09km s −1 ) filtering profiles which bracket the flattened  eff profile; moreover, the vertical profiles of the zonal velocity (Figure 3) show the zonal velocities to be consistently larger than in any of the other   values examined with   = 0.75.This indicates that while there is an increase in zonal velocity with the reduced filter strength, the profile itself impacts the velocity to a Substellar Anti-stellar greater degree than the strength, at least within the range of values required for numerical stability.
The temperature structure is impacted in the same regions as seen in the   parameter study investigated in the previous section: the equatorial temperatures are relatively insensitive to the filtering, except at  ∼ 1 bar where the downward transport of entropy results in deep atmosphere heating and the previously discussed temperature inversion.At the mid-latitudes (Figure 9), the temperatures are hotter at  ≲ 10 −2 bar around the nightside and morning terminator -the region of the nightside gyres -compared to the   = 3 case, where it was seen in the previous section that weaker filtering resulted in cooler mid-latitude temperatures in the same region.This reinforces the idea that relative to the jet, the temperatures in the gyres are more sensitive to the local, smaller-scale transport processes.

The Sponge
The fiducial suite of simulations were performed with the sponge parameter of   = 0.75, i.e., the sponge covers the upper quarter of the computational domain.As the sponge strength is parametrised based on height (see Equation 2), the sponge strength is not constant along isobaric surfaces, and for some pressures the isobaric surface only partially intersects the sponge.This can be seen in Figure 1 where the lower dotted lines indicate where the surface begins to intersect the sponge and the upper dotted line indicates the lowest point where the isobaric surface is entirely within the sponge.Due to the sponge extending into the jet, reaching, at least partially, pressures of 10 −2 bar, we opt to investigate the impact of moving the sponge upward to   = 0.9. 4he zonal velocities for the simulations with   = 0.9 are shown in Figure 14 and the horizontal velocities at 1 mbar are shown in Figure 15.The case of   = 1 still forms a counter-rotating jet with a breakdown in the conservation of axial angular momentum (Figure 11) and an excess of kinetic energy relative to the less dissipative cases (Figure 13).The rate of axial angular momentum loss and the growth rate of the kinetic energy show larger differences between the two sponge cases for   = 1 compared to the super-rotating cases where we see a greater degree of agreement between the sponge cases; however, the resulting flow morphology and temperature structure are largely the same for both   = 1 cases.The cases which form a super-rotating jet (  = 3 and   = 6), a similar level of dissipation (Figure 12) is seen compared to the deep sponge (  = 0.75) cases, indicating that the sponge is not indirectly influencing the rate of axial angular momentum loss to a significant degree.The   = 3 case shows a faster peak zonal velocity than the   = 6 case (6.18 km s −1 compared to 5.92 km s −1 ; see Figure 10), while in the deep sponge case the   = 3 and   = 6 cases have essentially the same peak zonal windspeed (5.87 km s −1 and 5.85 km s −1 , respectively).The zonal velocities also begin to decrease with pressure higher in the atmosphere compared to the deep sponge case, indicating that the sponge serves to suppress the jet.Although not investigated here, Substellar the zonal wind speeds are likely also impacted by the upper boundary as the the maximum equatorial pressures at the upper boundary are between 0.3 mbar and 0.45 mbar, depending on the simulation.Equatorial zonal wind speeds at pressures lower than this are certainly altered by the upper boundary; however, the exact impact at larger pressures is unclear.While varying the location of the upper boundary is beyond the scope of this parameter study investigated here, we remain cognizant of its potential impact on the results.The temperature differences between the two sponge cases are small, with the largest differences occurring around the anti-stellar point at the equator where the temperatures differ by ≲ 50 K for pressures ≲ 8 × 10 −2 bar.This is expected since as the nightside is cooler and thus has smaller scale heights, the sponge does not reach the high pressures it does on the dayside due to the sponge layer being parametrised based on height, reducing its impact on the flow.The sponge can be seen to have a larger influence, perhaps unsurprisingly, on the vertical velocity .At the terminators along the equator (Figure 6) , the velocities at pressures ≲ 10 −2 bar can differ by up to 7 m s −1 , comparable to the peak vertical velocities seen in those regions.Moreover, sharp transitions to the sponge and the damping of the vertical velocities to zero can be seen at low pressures.This is especially apparent at mid-latitudes along the evening terminator (see Figure 7, lower right panel).While the temperatures are not significantly altered by the sponge, physical processes not considered here which depend on the local transport properties (e.g., disequilibrium chemistry, aerosols) may have their results altered by its presence.Similarly, post-processed models that make use of the local vertical velocities extracted from GCMs (e.g., Tsai et al. 2023) may also include impacts from the nature of the sponge applied in the original simulations.

Other Impacts of Damping
While the cases of   = 1 and   = 2 can safely be discounted as unphysical, the differing levels of damping within the subset of simulations that do form a super-rotating jet allows for the impact of the damping to be better understood.In this subsection we look at a few other aspects of the simulations affected by dissipation.

Spin Up
Since the choice is made to run the simulation for a fixed time period as running to convergence might not be feasible, the spin-up period cannot be discounted in determining the state at the end of the simulation.Figure 16 shows the evolution of the instantaneous peak zonal wind, sampled every 10 days, over the length of each simulation.Within the first 200 days, the zonal winds oscillate in peak amplitude, with the less damped simulations having larger oscillations.Within this early phase, up until day 400, the deep sponge simulations (  = 0.75) spin up faster; however, at day 1000 the shallow sponge simulations (  = 0.9) have faster zonal winds than those simulations with the same   .There is not, unfortunately, an obvious trend with   with the   = 12 simulation having the fastest zonal wind for the deep sponge cases and   = 3 having the fastest zonal wind for the shallow sponge cases.The simulation with a flattened  eff stands out from the rest of the simulations as spinning up faster than any model with a fixed   and having the fastest jet after 1000 days.In all cases, the jets continue to accelerate, with rates on the order of a few ×10 −4 km s −1 day −1 , based on an estimate of the rate of change over the final 100 days of each simulation.As with the jet speeds, there is not an obvious trend in the jet acceleration rate with Substellar   ; however, the simulations with a shallow sponge do have higher zonal wind accelerations in the last 100 days than those with a deep sponge.For the shallow sponge, the   = 3 and   = 6 cases are accelerating at 5.7 × 10 −4 km s −1 day −1 and 4.6 × 10 −4 km s −1 day −1 while for the deep sponge the   = 3 case accelerates the fastest at the end at 2.9×10 −4 km s −1 day −1 .For comparison, the simulation with the flattened  eff had an acceleration of 4.8 × 10 −4 km s −1 day −1 .This indicates that while the shallow sponge may result in a faster continued acceleration of the jet, the integration times required to see a significant increase in the jet speed based on a linear extrapolation are prohibitive.
We additionally look at the evolution of the axial angular momentum at fixed pressure intervals over the duration of the simulation (see Figure 17).The lowest pressure interval is 10 −2 bar <  < 10 −4 bar.At even lower pressures, isobaric surfaces in some of the simulations intersect the upper boundary, complicating the interpretation of the results.The   = 1 and   = 2 cases are also excluded as loss of angular momentum makes comparison difficult5 .As is shown in the Figure 17 (left panel), the deep ( > 1 bar) atmosphere loses axial angular momentum as the jet above forms, resulting in the deep atmosphere to slowly counter-rotate (see Figures 1 and 14 for the zonal winds).While some of this can be attributed to the global loss of angular momentum, the trend is seen in all the simulations, including the   = 6 and   = 12 cases which have improved conservation properties, indicating that this effect is predominantly due to the transport of angular momentum to lower pressures.These latter cases exhibit a leveling off of the axial angular momentum by day 600 both in the deep atmosphere and in the intermediate atmosphere (1 bar <  < 10 −2 bar).While this could be interpreted as the upper atmosphere approaching a dynamically steady state, the atmospheric kinetic energy (Figure 18), although showing some slowing, continues to evolve.We thus conclude that even with the shorter evolutionary timescales of the upper atmosphere relative to the deep atmosphere, the upper atmosphere in the simulation has not fully converged.

Vertical mixing
Extracting an effective vertical mixing parameter   from GCM simulations has become an important input to one-and twodimensional models of higher microphysical complexity (e.g., Powell et al. 2019;Tsai et al. 2024).While tracer-based methods (e.g., Parmentier et al. 2013) may provide more precise estimates of the global mixing, the analysis of Komacek et al. (2019) provides a straightforward estimate of   , where  and  are the root-mean-squared zonal and vertical velocities, computed on isobaric surfaces, and the advection timescale is estimated as  adv =  p /.The radius  p is taken to be the inner boundary radius (see Table 1).For the purposes of this estimate, we ignore the chemical timescale  chem and only consider extremely long-lived tracers in the atmosphere.We also note that this reduces Substellar to the traditional estimate of vertical mixing,   ∼ , where  is the scale height, when / p ∼ /.This has been shown to overestimate the global mixing (Parmentier et al. 2013); however it has been used, with an overall scaling applied, when an estimate of the local mixing rate is needed, as opposed to a global value (Tsai et al. 2023).
Figure 19 shows the estimates for   based on Equation 9. We find that reduced vertical damping by the sponge leads to a 60% increase in   above 10 −2 bar, which, given that the parameter space of mixing rates can span many orders of magnitude, is a relatively modest increase.At pressures greater than 10 bar, larger differences are present, with the simulations with increased damping showing order of magnitude decreases in vertical mixing.This region, however, is not converged in the simulation, as shown in Figure 8 where the region below 10 bar has yet to be heated by the downward advection of entropy (Tremblin et al. 2017).We further note that the   profiles computed from instantaneous outputs as opposed to a time-average show even larger variations in the mixing at pressures larger than 10 bar, while at pressures less than 10 bar the   estimates show little variation relative to their time-averaged values.

CONCLUSIONS
In this paper, we demonstrate that sufficiently strong longitudinal filtering of the horizontal velocity can lead to a suppression of a zonal equatorial jet and a breakdown in the conservation of axial angular momentum; however, for the case of WASP-96b investigated here, we find that once the filter is weak enough for a super-rotating jet to form, the speed of the jet and the thermal structure are relatively insensitive to the filter coefficient .An exception to this is within the lower pressure (≲ 10 −2 bar) nightside gyres where the temperatures can vary with  up to ∼ 100 K.A larger influence on the jet is seen from the latitudinal dependence of the filtering.As a latitudelongitude grid is used within the model, and since the filtering acts on the grid directly without involving the metric, filtering with a constant value of  results in stronger diffusion at latitudes near the equator.By altering  to partially account for this, resulting in a constant effective diffusion parameter  eff near the equator, we show that the latitudinal profile of  can impact the zonal jet to a greater degree than the choice of   , at least within the range of values that maintain numerical stability.We also find that the choice of dissipation can impact the vertical velocities, even in regions where only minor temperatures differences are seen.This may affect chemical transport processes, although that was not investigated in this paper.Given the impact of such dissipation on the simulations, observations are required to constrain the magnitude of the combined numerical and physical dissipative mechanisms.Additionally, as there are several physical processes likely to impart dissipation which are not explicitly included in GCM simulations (e.g., magnetic drag, gravity waves, shocks), although instructive, higher spatial and temporal resolution simulations do not provide an ultimate 'ground truth' nor remove the requirement for observational constraint to determine the dynamical structures of gas giant atmospheres.Substellar The simulations in this work were performed using the DiRAC Data Intensive service at Leicester, operated by the University of Leicester IT Services, which forms part of the STFC DiRAC HPC Facility (www.dirac.ac.uk).The equipment was funded by BEIS capital funding via STFC capital grants ST/K000373/1 and ST/R002363/1 and STFC DiRAC Operations grant ST/R001014/1.DiRAC is part of the National e-Infrastructure.As supercomputing is an energyintensive endeavor, we note that the production runs used in this paper required 241496 CPU-hours, resulting in 500 kg of CO 2 being emitted, based on estimates for UK electricity generation6 .This excludes emissions due to testing and analysis and is certainly an underestimation of the carbon impact of this research.
The analysis of the simulation data and the production of plots for this paper made use of the following Python packages: aeolus (Sergeev et al. 2022), iris (Hattersley et al. 2023), matplotlib (Hunter 2007), and numpy (Harris et al. 2020).
The authors would also like to thank the anonymous referee for their comments which greatly improved the quality of the paper.Substellar Axial Angular Momentum the eddy momentum flux gradients for interested readers.Generally, with the exception of the unphysical cases where a counter-rotating jet forms, the flux gradients show the same qualitative morphology for all cases, which is also similar to results found for the radiative cases in the hot Jupiter simulations performed in Mayne et al. (2017) (cf., their Figure 14, right column).The decreases in dissipation and damping do result in changes in the small-scale structure, but we do not investigate the issue any further.
To compute the relevant flux gradients, we follow Hardiman et al. (2010) and Mayne et al. (2017) in decomposing the flow into the sum of a zonally-averaged component and a perturbation from that mean (e.g., the zonal velocity  becomes  =  +  ′ where  and  ′ are the Kinetic Energy [×10 27 J] zonal mean and perturbations, respectively).The time evolution of the zonally-averaged momentum   can then be written as

APPENDIX B: RESOLUTION
As the effective dissipation coefficient  eff scales as  eff ∝ (Δ) 2 , a doubling of the longitudinal resolution results in a quartering of  eff .At the fiducial resolution used in this study, a change from   = 1 to   = 3 -transitioning from the counter-rotating to the super-rotating cases -reduces  eff by a factor of 2.23.Thus, one may expect that a   = 1 simulation with the longitudinal resolution doubled would result in a super-rotating jet.As discussed in Heng et al. (2011), it is not always possible when investigating the behaviour of a dissipation operator to decouple the effect of the operator from the intrinsic effects of changing the grid.We nonetheless opt to verify explicitly that the increase in resolution results in the formation of a superrotating jet.
To this end, we perform a test simulation with   = 1 and the horizontal grid resolution doubled in both the latitude and longitude, resulting in a grid of 288 longitude points and 180 latitude points.The simulation is able to run with the same dynamical timestep as the fiducial simulations, Δ dyn = 30 s, avoiding complications in interpreting the results due to the timestep changing as well (see Equation 7).Further refinements in the grid would, unfortunately, require a reduction in Δ dyn , increasing  eff . 7As the intention is to only demonstrate that a jet forms, we run the simulation for 200 days, sufficient for the model to spin-up.
The results from the simulation are shown in Figure B1.As a quartering of  eff relative to the   = 1 case is equivalent to increasing   to ∼ 5.81, we compare the results of the higher resolution simulation to both the   = 1 case at the fiducial resolution as well as   = 6 at the fiducial resolution.As is shown in Figure B1, the increased resolution simulation with   = 1 forms a super-rotating jet qualitatively similar to the   = 6 result at the fiducial resolution, with peak mean zonal wind speeds of 5.04 km s −1 and 5.01 km s −1 , respectively, further supporting the idea that the physical breakdown is due to the  eff .We see little effect from the resolution in the mean zonal wind or in other quantities (not plotted here) beyond the impact on  eff .Resolution dependent effects may become apparent with a reduction in   in the increased resolution simulation; however, that is beyond the scope of this paper.This paper has been typeset from a T E X/L A T E X file prepared by the author.
7 While a resolution study would be valuable for the understanding of the behaviour of the um generally, especially given the resolution dependence found in the hot Jupiter simulations of Skinner & Cho (2021) and the relatively high resolution requirements needed for convergence within their model, such a study should be designed based on the dissipation requirements to run at the highest resolutions, not as an off-shoot of a study that looks at low-resolution behaviour.t K = 1, s = 0.9 60°30°0°+30°+60°t K = 3, s = 0.9 60°30°0°+30°+60°t K = 6, s = 0.9 60°30°0°+30°+60°F lattened K eff .   estimated for each simulation using Eqn.9, computed using the average velocities and pressures over the final 100 days of the simulations.The simulations which resulted in counter-rotating jets have been omitted as they represent unphysical results.Vertical Eddy Momentum Flux Gradient (×10 6 kg m 2 s 2 ) Figure A1.The vertical eddy momentum flux gradient for the average atmospheric state between days 900 and 1000 for each value of   in the   = 0.75 parameter study (see Equation A1).For the two over-damped cases,   = 1 and 2, the vertical transport of eddy momentum can be seen to be vastly different from the cases in which a super-rotating equatorial jet develops,   = 3, 6, and 12.Note that the vertical coordinate in these plots is the height , as opposed to the pressure as used in other plots and that the plots show only the central region of the computational domain.A1) for the average atmospheric state between days 900 and 1000 for the   = 0.9 cases (three leftmost panels) and the case with a flattened  eff and   = 0.75 (rightmost panel).Note that the vertical coordinate in these plots is the height , as opposed to the pressure as used in other plots and that the plots show only the central region of the computational domain.Meridional Eddy Momentum Flux Gradient (×10 6 kg m 2 s 2 ) Figure A3.The meridional eddy momentum flux gradient for the average atmospheric state between days 900 and 1000 for each value of   in the   = 0.75 parameter study (see Equation A1).For the two over-damped cases,   = 1 and 2, the vertical transport of eddy momentum can be seen to be vastly different from the cases in which a super-rotating equatorial jet develops,   = 3, 6, and 12.Note that the vertical coordinate in these plots is the height , as opposed to the pressure as used in other plots and that the plots show only the central region of the computational domain.A1) for the average atmospheric state between days 900 and 1000 for the   = 0.9 cases (three leftmost panels) and the case with a flattened  eff and   = 0.75 (rightmost panel).Note that the vertical coordinate in these plots is the height , as opposed to the pressure as used in other plots and that the plots show only the central region of the computational domain.Mean Meridional Momentum Flux Gradient (×10 6 kg m 2 s 2 ) Figure A7.The mean meridional momentum flux gradient for the average atmospheric state between days 900 and 1000 for each value of   in the   = 0.75 parameter study (see Equation A1).For the two over-damped cases,   = 1 and 2, the vertical transport of eddy momentum can be seen to be vastly different from the cases in which a super-rotating equatorial jet develops,   = 3, 6, and 12.Note that the vertical coordinate in these plots is the height , as opposed to the pressure as used in other plots and that the plots show only the central region of the computational domain.A1) for the average atmospheric state between days 900 and 1000 for the   = 0.9 cases (three leftmost panels) and the case with a flattened  eff and   = 0.75 (rightmost panel).Note that the vertical coordinate in these plots is the height , as opposed to the pressure as used in other plots and that the plots show only the central region of the computational domain.

Figure 1 .
Figure 1.The mean zonal wind in the   = 0.75 models, averaged over the last 100 days.The dashed line indicates the boundary between super-rotating and counter-rotating flows.The dotted lines indicate the pressure levels at which the pressures partially (lower line) and fully (upper line) intersect the sponge layer.For the more heavily damped cases (  = 1 and 2), a super-rotating jet does not develop.

Figure 3 .
Figure 3.The zonal velocities  at the equator averaged over the final 100 days.The circles indicate where the profiles begin to intersect the sponge.As the counter-rotating simulations are unphysical, they have been omitted for clarity.

Figure 4 .
Figure 4.The zonal velocities  at a latitude of 45°averaged over the final 100 days.The circles indicate where the profiles begin to intersect the sponge.As the counter-rotating simulations are unphysical, they have been omitted for clarity.

Figure 5 .
Figure 5.The meridional velocities  at a latitude of 45°averaged over the final 100 days.The circles indicate where the profiles begin to intersect the sponge.As the counter-rotating simulations are unphysical, they have been omitted for clarity.

Figure 6 .
Figure 6.The vertical velocities  at the equator averaged over the final 100 days.The circles indicate where the profiles begin to intersect the sponge.As the counter-rotating simulations are unphysical, they have been omitted for clarity.

Figure 7 .
Figure 7.The vertical velocities  at a latitude of 45°averaged over the final 100 days.The circles indicate where the profiles begin to intersect the sponge.As the counter-rotating simulations are unphysical, they have been omitted for clarity.

Figure 8 .
Figure 8. Temperature profiles at the equator, averaged over the last 100 days.Colors indicate the value of   (see legend in the upper left panel) while linestyles indicate the sponge parameter (solid lines correspond to   = 0.75 while dashed lines correspond to   = 0.9).

Figure 9 .Figure 10 .
Figure 9. Temperature profiles at a latitude of 45°, averaged over the last 100 days.As with Figure 8, the colors indicate the value of   (see legend in the upper left panel) while linestyles indicate the sponge parameter (solid lines correspond to   = 0.75 while dashed lines correspond to   = 0.9.

Figure 11 .
Figure 11.The axial momentum evolution for each case, normalized to the initial axial angular momentum.

Figure 12 .
Figure 12.Rate of dissipation of the normalized axial angular momentum based on a linear fit to the axial angular momentum evolution in Figure 11.

Figure 13 .
Figure13.The evolution of the kinetic energy for each simulation.The simulations which resulted in counter-rotating jets show a significantly faster growth in their total kinetic energy.
various transport terms are labelled therein.The meridional and vertical eddy momentum flux gradients as well the mean meridional and vertical momentum flux gradients for all the simulations investigated in this paper are shown in Figures A1 through A8, inclusive.

Figure 14 .Figure 15 .
Figure14.The mean zonal wind in the   = 0.9 simulations (left three panels), and the simulation with a flattened  eff and   = 0.75 (rightmost panel) , averaged over the last 100 days.The dashed line indicates the boundary between super-rotating and counter-rotating flows.The dotted lines indicate the pressure levels at which the pressures partially (lower line) and fully (upper line) intersect the sponge layer.For the more heavily damped cases (  = 1 and 2), a super-rotating jet does not develop.

Figure 16 .Figure 17 .
Figure16.The evolution of the peak zonal wind speed for each case that resulted in a super-rotating jet.

Figure 18 .
Figure18.The evolution of kinetic energy within three regions of the atmosphere.Unlike the axial angular momentum in Figure17, no normalisation is applied.

Figure A2 .
Figure A2.The vertical eddy momentum flux gradient (see EquationA1) for the average atmospheric state between days 900 and 1000 for the   = 0.9 cases (three leftmost panels) and the case with a flattened  eff and   = 0.75 (rightmost panel).Note that the vertical coordinate in these plots is the height , as opposed to the pressure as used in other plots and that the plots show only the central region of the computational domain.

Figure A4 .
Figure A4.The meridional eddy momentum flux gradient (see EquationA1) for the average atmospheric state between days 900 and 1000 for the   = 0.9 cases (three leftmost panels) and the case with a flattened  eff and   = 0.75 (rightmost panel).Note that the vertical coordinate in these plots is the height , as opposed to the pressure as used in other plots and that the plots show only the central region of the computational domain.

Figure A8 .
Figure A8.The mean meridional momentum flux gradient (see EquationA1) for the average atmospheric state between days 900 and 1000 for the   = 0.9 cases (three leftmost panels) and the case with a flattened  eff and   = 0.75 (rightmost panel).Note that the vertical coordinate in these plots is the height , as opposed to the pressure as used in other plots and that the plots show only the central region of the computational domain.