Shallow water simulations are used to present a unified study of three major storms on Saturn (nicknamed as Great White Spots, GWS) at different latitudes, polar (1960), equatorial (1990), and mid-latitude (2010) (Sánchez-Lavega, 2004; Sánchez-Lavega et al., 2011). In our model, the three GWS are initiated by introducing a Gaussian function pulse at the latitude of the observed phenomena with controlled horizontal size and amplitude. This function represents the convective source that has been observed to trigger the storm. A growing disturbance forms when the pulse reacts to ambient winds, expanding zonally along the latitude band of the considered domain. We then compare the modeled potential vorticity with the cloud field, adjusting the model parameters to visually get the closest aspect between simulations and observations. Simulations of the 2010 GWS (planetographic latitude ∼+40º, zonal velocity of the source ∼−30m s−1) indicate that the Coriolis forces and the wind profile structure shape the disturbance generating, as observed, a long region to the east of the convective source with a high speed peripheral anticyclonic circulation, and a long-lived anticyclonic compact vortex accompanied by strong zonal advection on the southern part of the storm forming a turbulent region. Simulations of the equatorial 1990 GWS (planetographic latitude +12º–+5º, zonal velocity of the source 365–400m s−1) show a different behavior because of the intense eastward jet, meridional shear at the equatorial region, and low latitude dynamics. A round shaped source forms as observed, with the rapid growth of a Kelvin-Helmholtz instability on the north side of the source due to advection and to the strong meridional wind shear, whereas at the storm latitude the disturbance grows and propagates eastward. The storm nucleus is the manifestation of a Rossby wave, while the eastward propagating planetary-scale disturbance is a gravity-Rossby wave trapped around the equator. The simulated 1960 GWS disturbance (planetographic latitude +56º, zonal velocity 4m s−1) formed a chain of periodic oval spots that mimic the few available observations of the phenomenon. For the mid and high latitude storms, simulations predict a strong injection of negative relative vorticity due to divergence of the upwelling storm material, which may produce large anticyclones on the anticyclonic side of the zonal profile, and a quick turbulent expansion on the background cyclonic regions. In general, simulations indicate that negative relative vorticity injected by storms determines the natural reaction to zonal winds at latitudes where Coriolis forces are dominant.