Towards long-term simulations of planetary-scale vortices and storms on Jupiter and Saturn

Long-term simulations of planetary vortices and storms are essential for improving our understanding of the atmospheric dynamics on gas giants such as Jupiter and Saturn. These simulations allow investigating the physical mechanisms that govern the lifespan and stability of features like Jupiter’s Great Red Spot (GRS) [1] and Saturn’s giant storms [2]. Furthermore, long-term integration is necessary to reproduce and assess observed morphological changes in key atmospheric features—such as the gradual shrinking or 90-days longitudinal oscillations of the GRS [3, 4]—, recorded over decades by spacecraft missions as well as ground- and space-based telescope observations. By capturing these evolving processes, simulations can test the adequacy of proposed forcing and dissipation mechanisms, assess the role of convective injection [5], and inform planetary circulation models [6], among others. The objective of this work is to improve current planetary atmospheric models to reproduce relevant scales of motion on Jupiter and Saturn with very small dissipation over time scales of years. In order to understand the atmospheric dynamics of gas giants, we require models capable of reproducing and simulating the processes involved. Specifically, we employ the Shallow Water (SW) equations to reproduce fluid dynamics, focusing on a thin atmospheric layer influenced by rotational effects and zonal winds. Despite their relative simplicity, the SW equations capture essential dynamics that govern atmospheric phenomena like planetary-scale vortices and storms, e.g., Saturn’s Great White Spot [2]. To simulate these phenomena within the model, we introduce a localized Gaussian-shaped pulse perturbation representing a convective source. Model parameters are tuned to optimize agreement between simulated features and observational data from spacecraft instrumentation and telescope observations. SW equations are typically solved using a semi-discrete approach, where spatial discretization is applied first, followed by time integration. However, both conventional spatial and temporal schemes present limitations that hinder long-term simulations of planetary atmospheres. Large-scale perturbations introduced in the model trigger a dynamical response that involves substantial instabilities, such as gravity waves, which represent a challenge for the numerical methods. To address this, many SW models employ advanced spatial discretization techniques designed to suppress unphysical oscillations. A commonly used approach is the Total Variation Diminishing (TVD) finite volume scheme, which incorporates gradient-based limiters that dynamically switch between lower- and higher-order formulations. While these limiters are effective at managing steep gradients and instabilities, they can also introduce excessive numerical dissipation, which compromises the conservation of the flow properties. To address these limitations, the objective of this work is twofold: 1) Implementing entropy viscosity algorithms to mitigate numerical dissipation [7]; 2) Adopting higher-order strong stability preserving Runge-Kutta integration schemes to improve temporal integration [8]. On the one side, the entropy viscosity technique dynamically adjusts dissipation based on local flow properties. Specifically, the model spots unstable regions through local entropy-residual analysis, subsequently tackling only those areas where instabilities are detected. Additionally, the formulation also leverages Ducros-type discontinuity sensors [9] to further refine the identification of regions requiring stabilization by distinguishing between compressive (shock-like) and vortical (rotational) structures. The result is an artificial viscosity that enables selective dissipation, preserving fine-scale features in smooth regions while controlling instabilities elsewhere. Furthermore, this artificial viscosity is used to complement gradient-based limiters in TVD schemes, offering a hybrid strategy that enhances stability without significantly compromising accuracy, improving the robustness of SW simulations over long time scales. On the other side, for time integration, we use pseudo-symplectic (Low Storage) Strong Stability Preserving Runge-Kutta methods, which provide a better bounding of the solution while improving energy conservation and long-term stability. These schemes are TVD-implicit and show a significantly larger Region of Absolute Stability (ROS) than commonly used multistep methods such as Adams-Bashforth [8], which makes them particularly suitable for the highly nonlinear flows considered in this study. We validate our approach using a set of benchmark test cases based on [10], demonstrating its effectiveness in controlling numerical dissipation and its improved accuracy as compared to traditional schemes. Finally, we apply our methodology to simulate long-term atmospheric dynamics associated with Jupiter’s Great Red Spot and Saturn’s Great White Spot, illustrating the model’s capability to reproduce long-term, planetary-scale vortices and storms with high fidelity, i.e., preserving reliable dynamics over extended simulation periods with improved stability.   [1] Sánchez-Lavega, A. et al. Jupiter’s Great Red Spot: Strong interactions with incoming anticyclones in 2019, Journal of Geophysical Research: Planets 126 (2021). [2] García-Melendo, E. and Sánchez-Lavega, A. Shallow water simulations of Saturn’s giant storms at different latitudes. Icarus 287, 241-260 (2017). [3] Simon, A., et al. Historical and contemporary trends in the size, drift, and color of Jupiter's Great Red Spot. The Astronomical Journal 155(4), 151 (2018). [4] Sánchez‐Lavega, A. et al. The Origin of Jupiter’s Great Red Spot. Geophysical Research Letters 51, (2024). [5] García-Melendo, E., et al. Atmospheric dynamics of Saturn’s 2010 giant storm. Nature Geoscience 6, 525–529 (2013). [6] Spiga, A., et al. Global climate modeling of Saturn's atmosphere. Part II: Multi-annual high-resolution dynamical simulations. Icarus 335, 113377 (2020). [7] Guermond, J.-L., et al. Entropy viscosity method for nonlinear conservation laws. Journal of Computational Physics 230, 4248–4267 (2011). [8] Higueras, I. and Roldán, T. Efficient SSP low-storage Runge–Kutta methods. Journal of Computational and Applied Mathematics 387, 112500 (2021). [9] Hendrickson, T. R., et al. An Improved Ducros Sensor for the Simulation of Compressible Flows with Shocks. Fluid Dynamics Conference (AIAA, 2018). [10] Williamson, D. L., et al. A standard test set for numerical approximations to the shallow water equations in spherical geometry. Journal of Computational Physics 102(1), 211-224 (1992).