Hierarchical Modeling of Solar System Planets with Isca

: We describe the use of Isca for the hierarchical modeling of Solar System planets, with particular attention paid to Earth, Mars, and Jupiter. Isca is a modeling framework for the construction and use of models of planetary atmospheres at varying degrees of complexity, from featureless model planets with an atmosphere forced by a thermal relaxation back to a speciﬁed temperature, through aquaplanets with no continents (or no ocean) with a simple radiation scheme, to near-comprehensive models with a multi-band radiation scheme, a convection scheme, and conﬁgurable continents and topography. By a judicious choice of parameters and parameterization schemes, the model may be conﬁgured for fairly arbitrary planets, with stellar radiation input determined by astronomical parameters, taking into account the planet’s obliquity and eccentricity. In this paper, we describe the construction and use of models at varying levels of complexity for Earth, Mars and Jupiter using the primitive equations and/or the shallow water equations.


Introduction
Planetary atmospheres are extremely diverse, even within our own Solar System. Thus, for example, tidally locked Mercury, in a 3:2 spin resonance [1], has an extremely thin atmosphere mainly of hydrogen and helium with virtually no greenhouse effect and surface temperatures ranging from 100 K to 700 K, whereas Venus has a carbon dioxide atmosphere over 90 times thicker than that of Earth. This, combined with its slow rotation (a sidereal day over 200 Earth days) and its highly reflective upper-tropospheric clouds, give it an almost uniform surface temperature of about 700 K and a troposphere about 60 km thick, as opposed to the 10 km of Earth (e.g., [2,3] and references therein). Mars, perhaps the most Earth-like (or "telluric") of the other planets has a thin carbon dioxide atmosphere with a weak greenhouse effect. Its radius is about half that of Earth and its rotation period is about the same, giving rise to baroclinic instability and large dust storms. The obliquity of Mars is currently about 25 • , similar to that of Earth, although it likely has had much larger values (up to 40 • ) and much smaller values in the distant past [4,5]. The obliquity, combined with the lack of oceans (and thus a small surface heat capacity) gives rise to quite strong seasons, so much so that when a pole is in darkness up to 30% of the atmosphere may be deposited on the surface giving a seasonal varying layer of solid carbon dioxide a meter or so thick (see, e.g., Section 3.4 of [6]).
Moving further from the Sun, we find Jupiter, a gas giant with ten times the radius of Earth and with a rotation period of just 10 h, giving rise to the splendid bands of alternating zonal velocity with speeds around 100 m s −1 with swirling vortices embedded within them, all visible from Earth and, with much more detail, from the recent Juno mission. Beyond Jupiter is Saturn, another gas giant, rather smaller, but surrounded by its beautiful and famous rings. Then, on to the ice giants, Uranus and Neptune, both unique in their own way. Uranus, for example, has an obliquity of 97 • , meaning one of its poles points almost directly at the Sun at the solstice. Its atmosphere is a hydrogen-helium mix, 1. There are so many planets that constructing a comprehensive model for each is simply infeasible. 2. The observational data available for solar system planets are orders of magnitude less than that for Earth in terms of spatial and temporal coverage, even with modern missions to other planets such as the Venus Express or Juno. The data available for exoplanets are orders of magnitude less still. Thus, were we to choose to use a highly complex model to model a given planet, we would be in danger of having to make too many choices about important parameters and processes that are, at the very least, under-constrained by the observations. We would therefore run the risk of over-interpreting matches between model output and observations when the solutions could be highly non-unique, thus running the risk of "over-fitting" the data to the model. 3. Although there is a diversity of planetary atmospheres, they obey the same physical laws. This, and the desire to understand as well as simulate (a desire that also holds for Earth) suggests that we should use models that take advantage of that commonality and build from there, with appropriate complexity for the problem at hand.
There will, however, be no single level of complexity that is appropriate even for a given planet, especially if we are to both understand and compare to observations. Thus, we propose that, as has been implemented to some degree already for Earth (for a discussion, see [10]), a hierarchical approach will be useful. In such an approach, we may use a rather simple class of models that can be applied with minimal changes across a range of planetary atmosphere to understand the basic behavior as various parameters (e.g., rotation rate) are changed, building up in complexity if we wish to study a particular planet in more detail. We choose three Solar System planets to illustrate this approach-Earth, Mars, and Jupiter. The approach for Earth has been discussed elsewhere, particularly in Vallis et al. [11], thus here our description is brief; however, it is a familiar planet and a useful place to begin. Mars is, as noted, a terrestrial planet and a close neighbor, and Jupiter a giant planet. Taken together, they provide a useful introduction to our approach.

Isca, the Modeling Framework
Isca [11] is a framework for building atmospheric models at various degrees of complexity that utilizes the the Flexible Modeling System (FMS, https://www.gfdl.noaa.gov/fms/) software infrastructure developed at the Geophysical Fluid Dynamics Laboratory (GFDL), USA (see also [12]). Most of the model is written in Fortran, but it is configured, and the output generally analyzed, with Python. Isca itself is not a model and, although there is a certain modeling philosophy associated with its use, Isca does not impose strict modeling choices on the user. Roughly speaking, the following options are currently available: The above options are all publicly available on the main GitHub site (although the SOCRATES radiation code must be downloaded separately from the UK Met Office). Various other options are under development. In particular, the next release of Isca will have a fully dynamical ocean and a choice of dynamical cores, including a scheme on the cubed-sphere. We now describe how these options can be used for Earth, Mars, and Jupiter.

Earth
Various models of Earth's atmosphere have been constructed using Isca and are described in more detail elsewhere (e.g., [11,[28][29][30][31]). Here, we give three examples covering the range of possibilities. For each of these, and for those listed for other planets, we provide a parameter list in Table 1.

A Simple Thermal-Relaxation Model
For our simplest model of Earth's atmosphere, we use a primitive equation dynamical core with Newtonian cooling of temperature and Rayleigh drag in the boundary layer. Traditionally, models with Newtonian cooling of temperature tend to have a constant-in-time relaxation temperature, T eq , e.g., the classical Earth example of Held and Suarez [14]; however, it is useful to incorporate Earth's seasonal cycle into such a model. We therefore construct a seasonally-varying equilibrium temperature field using the method described in Section 3 of Vallis et al. [11]. Table 1. Parameters for each of the runs shown in the present work . E-S, E-I, and E-C are the simple, intermediate, and complex models of Earth described in Section 3. M-S, M-I, and M-C are the Mars models described in Section 4. J-S and J-I are the simple and intermediate models of Jupiter from Section 5. Other short-hand used are Prim., meaning primitive-equation dynamics; SW, meaning shallow-water dynamics; SBM, meaning the simplified Betts-Miller convection scheme; τ Conv , being the timescale for convective relaxation; and t 30 , being the wall-time in minutes taken to run 30-Earth-days in this configuration using N procs processors. None of the model configurations displayed in this paper use Q-fluxes. All of the parameters below may be easily changed as desired. This method uses Isca's orbital-calculation routines to calculate the incoming top-of-atmosphere shortwave radiation, S(θ, φ, t), where θ and φ are longitude and latitude, respectively. This is then averaged over the diurnal cycle to generate a daily-mean insolation,Ŝ(φ, t), which is then used to calculate the temperature at the tropopause

Planet
which is assumed to be in radiative equilibrium with the net incoming shortwave radiation, where S is the incoming top-of-atmosphere shortwave, σ is the Stefan-Boltzmann constant, and α is the planet's surface albedo. From this temperature, the height of the tropopause is defined, following Equation (A11) of Vallis et al. [32], where Γ is a prescribed lapse rate, C ≈ log 4 ≈ 1.4, τ s is the surface longwave optical depth, and H a is the scale height of the main infrared absorber. This expression comes from assuming grey radiative transfer in the infrared, an infrared optical depth of τ = τ s e −z/H a that varies exponentially with geometric height above the surface (z) with a scale height of H a , hydrostatic balance, and a constant tropospheric lapse rate, Γ (for further details see, Appendix of [32]). The same constant lapse rate used in the calculation of the tropopause height, Γ, is then used to extend temperatures from the tropopause down to the surface, defining a near-surface air temperature, These near-surface air temperatures are then used to calculate a ground temperature, T g , where where C g is the heat capacity of the surface and A is a constant, for example 4σT 3 0 , with T 0 a representative surface value. (We actually use the first formulation in this paper, although the second is equally or perhaps more appropriate.) This calculation allows the surface to have a finite heat capacity, enabling the surface temperatures to lag the insolation and giving a more realistic seasonality. This surface energy balance is simplified compared with that, e.g., for a mixed-layer ocean, as it does not directly account for shortwave absorption by the ground, or explicitly sensible or latent heat fluxes from the ground into the atmosphere, or the full longwave downward flux at the surface. All of these processes are implicitly accounted for by the form of Equation (4), since in steady state it gives T s = T g , with no temperature discontinuity at the surface, which strict radiative equilibrium has.
The ground temperatures T g are then extrapolated upward to give the atmospheric radiativeconvective-equilibrium temperature distributions, T eq (θ, φ, z, t) using Thus, the radiative-convective equilibrium temperatures have a seasonal dependence, one that may lag the distribution of incoming solar radiation due to the finite heat capacity of the surface.
A decision that one must make with this setup is how to deal with temperatures above the tropopause, where the lapse rate Γ is no longer appropriate. For the simulation presented in Figure 1a, we impose a minimum value on T eq (θ, φ, z, t) of 200 K, which gives a constant-temperature stratosphere, and is the same stratospheric condition as that imposed by Held and Suarez [14]. In such a configuration, the lack of shortwave absorption by the atmosphere means that the only way the atmosphere knows about the seasonal cycle is through the atmosphere's absorption of longwave emission by the ground (as seen, for example, in Equation (5)). However, the ground temperatures lags the insolation by an amount that depends on the ground's heat capacity, as discussed by Schneider [33], with the implications of this lag for simplified models discussed by Jucker [34]. The size of this lag on Earth is roughly 40 • , and also roughly 40 days (see, e.g., Figure 1 of [35]). A simple scaling argument (given in Section 3.2) predicts this to occur with a mixed-layer depth of about 5 m. However, such shallow mixed-layer depths allow the ITCZ to travel significantly further from the equator than on the real Earth [35]. This in turn suppresses the summer Hadley cell and promotes equatorward motion aloft in the winter hemisphere, leading to an unrealistic reverse-direction subtropical jet. Using a depth of 20 m gives a better climatology in the winds and temperatures but still allows a too-large migration of the ITCZ across the equator and subsequent reverse-direction subtropical jet, and also results in seasonal lag of approximately 90 days, meaning that the model's boreal Spring looks more like boreal Winter. We therefore show zonal-mean zonal winds averaged over March, April, and May in Figure 1a in the case with 20 m depth for comparison with the real-world December, January, and February zonal-mean zonal wind from the Japanese Reanalysis product, JRA-55 [36], in Figure 1d. Aside from the time lag and the aforementioned reverse-direction subtropical jet, there are many similarities between these fields including the strength and location of the midlatitude jet streams and the presence of the trade winds at low latitudes. Evidently an Earth-like circulation with a seasonal cycle can be produced with a very simple model, although there are several factors to consider regarding the representation of seasonality. The seasonality issues are not a particular consequence of using a simple radiative forcing; rather, they arise because a slab or mixed-layer lower boundary condition cannot properly account for ocean heat transport and its seasonal variability. This also holds in more complex models, as we see below.

An Intermediate Complexity Model
One step up in the hierarchy is a so-called "intermediate-complexity" model, namely a grey-radiation aquaplanet with a simple inclusion of moisture. The setup is essentially that used by Frierson et al. [15], but here we add a seasonal cycle to the radiation. The model has a number of useful features, specifically that it produces Earth-like zonal-winds, temperatures, and precipitation patterns at a relatively low computational cost.
Associated with its simplicity the model has a number of shortcomings-its representation of the stratosphere is poor, as can be seen by the lack of polar vortex in Figure 1a as compared with the wind peak seen at low pressures in the polar regions in Figure 1c,d. However, and not surprisingly, with a simple mixed layer the problem of seasonality discussed above remains, as we now describe.
The evolution equation for our mixed-layer ocean's temperature, T g is where C g is the heat capacity per unit area, SW is the net shortwave at the surface, LW is the net longwave cooling of the surface, and SH and LH are the sensible and latent heat fluxes out of the surface. A simple analysis of this equation can be done by ignoring the SH and LH terms, and assuming that the ground temperature T g is close to the overlying atmospheric temperature T 0 . This means that we can expand the longwave cooling of the surface (LW ∼ σT 4 g − σT 4 0 ) around T g ∼ T 0 , meaning that LW ≈ 4σT 3 0 T g . If SW = S 0 exp iωt , then the lag between the insolation and the surface temperatures is approximately φ ≈ arctan (C g ω/4σT 3 0 ). This scaling predicts that the realistic lag of about ∼40 days will occur with a mixed-layer depth of 5 m, which is clearly much shallower than the true ocean's mixed-layer and as noted in the previous section leads to unreasonable zonal-mean winds in the tropics. Thus, as with the simple model, we run our intermediate complexity model with a mixed-layer depth of 20 m, and again we show the zonal-mean zonal winds averaged over March, April, and May in Figure 1b for comparison with the real-world December, January, and February average in panel Figure 1d. The migration of the ITCZ is smaller in this setup compared with that in the simple model, meaning we do not see a reversed subtropical jet, and aside from the obvious time difference the simulation compares fairly well with observations To obtain a still more realistic simulation, we include continents, a more realistic radiation scheme, and impose the sea-surface temperature, as we now describe.

A More Comprehensive Model
In Figure 1c, we show results from a configuration at the more complex end of those available within Isca. It has realistic Earth topography, and a prescribed seasonally varying but annually repeating sea-surface temperatures taken from the AMIP experiment [37]. These additions give the model the ocean-temperature signatures of ocean currents, such as the Gulf stream and the Kurushio, and keep the seasonal cycle of ocean temperatures realistic without worrying about mixed-layer depths and the associated seasonal lag, with the result that the models seasonality compares well with reanalysis products, as shown in Figure 2. The land in this simulation does not have its temperatures prescribed; rather, it is allowed to evolve freely with a heat capacity corresponding to 2 m of water. Compared to the ocean, it has a higher albedo and reduced moisture availability, implemented in the form of an evaporative resistance. We use the comprehensive Socrates radiation code, with the same number of shortwave and longwave bands as are used in the UK Met Office's Global-Atmosphere 7 configuration [19]. This scheme was also used in Isca by Thomson and Vallis [38]. The result of all this is to significantly improve the zonal-mean zonal winds in winter, as seen in the comparison with the JRA-55 reanalysis shown in Figure 1. There is significant improvement in the stratosphere, where shortwave absorption by ozone improves the representation of stratospheric temperature gradients and therefore the zonal winds. Despite its high level of complexity compared with other models within Isca, the model is still simplified, and consequently faster to run, than a comprehensive atmosphere-only climate model. We are missing processes related to clouds, aerosols, sea ice, and gravity wave drag, although the inclusion of the first and last of these processes is currently underway. Despite these simplifications, however, the climate of this simulation has many realistic features in the mean, including the meridional overturning circulation shown in Figure 3, and in the atmospheric variability, with the model displaying a reasonable North-Atlantic Oscillation (see [28], for a full description).
An alternative formulation of similar complexity that we do not discuss here is to use a mixed-layer ocean with Q-fluxes, namely seasonally-dependent source and sink terms in the surface mixed-layer temperature equation that account for the effects of horizontal transport. The Q-fluxes allow the seasonal cycle to be properly simulated while avoiding the disadvantages of fixed surface temperatures (see [28,30], for further details). An important conclusion emerging from all of the above formulations is that simulations of the seasonal cycle require a careful treatment of the surface boundary condition, whether it is a fixed temperature, a mixed layer or a full ocean. Figure 2. Zonal winds at 250 hPa in our most realistic Earth model (a,b) and JRA-55 (c,d). The left column is averaged over December, January, and February, and the right column is averaged over June, July, and August.  Figure 3. Meridional overturning streamfunction in our most realistic Earth model (a,b) and from JRA-55 (c,d). The left column is averaged over December, January, and February, and the right column is averaged over June, July, and August.

Mars
Mars is the best-observed planet within the solar system aside from Earth, so much so that there are now several Martian reanalysis products available that combine temperature observations from satellites with state-of-the art Martian GCMs to create a "best-guess" version of Martian atmospheric conditions, e.g., the ensemble Kalman filter approach of [39] and the "Mars Analysis Correction Data Assimilation" (MACDA [40]). Nevertheless, several unanswered questions remain. For example, Mitchell et al. [41] found that Mars' polar vortex has an annular shape in potential vorticity (PV), suggesting that it should be barotropically unstable, but its observed steadiness suggests that it somehow remains stable. Recent work has found that Mars' short radiative timescale (a consequence of its thin atmosphere) could be responsible for its stabilization [42], but that work was done using a single-layer shallow water model. Experiments with a more comprehensive model are required to understand the influence of such things as radiative transfer, topography, Mars' thermal tides, and global dust storms, which were found to influence the vortex in Mitchell et al. [41].
Our approach to this problem is, again, hierarchical. We begin with a rather basic Martian general circulation model and add features gradually and systematically until a picture is built up as to which features are necessary to produce a circulation that compares well with observations. The details of this work will be reported elsewhere but here we present an overview of the range of models we are using, all constructed using Isca.
One aspect that is common to all our models constructed for Mars is that they measure the time of year using "Areocentric longitude", L S , which is a measure of the angle that Mars is through its orbit around the Sun (see e.g., figure A.1 of [6]), and is a simple offset of the orbit's true anomaly so that perihelion happens at L S = 251 • , and aphelion happens at L S = 71 • . All plots in the following sections therefore have their times referred to using L S .

A Simple Model
As with our Earth model hierarchy, our simplest Mars model uses Newtonian relaxation of temperatures with a seasonally-varying equilibrium temperature. Here, instead of calculating the equilibrium temperatures using radiative-convective equilibrium, we use pure radiative equilibrium. On Mars, the radiative relaxation timescale is quite short and convection has less influence on the time-mean, meaning that radiative equilibrium is a reasonable first step. To find the equilibrium temperatures, we first calculate the incoming shortwave radiation received by Mars using Isca's routines for orbital calculations, which account for both the eccentricity of Mars' orbit and its orbital inclination (obliquity). We then make the simplifying assumption that Mars' atmosphere is transparent in the shortwave, meaning that the outgoing longwave radiation (OLR) required for the planet to be in radiative balance at each location is where α is the surface albedo, θ is longitude, φ is latitude, and t is time. We assume that the atmosphere is in radiative equilibrium, is grey in the infrared, and we make the two-stream approximation, meaning that radiation only travels vertically, and is not scattered. To calculate the radiative equilibrium temperatures, we use the Schwartzchild equations [43,44], where D is the downward flux of infrared radiation, U is the upward flux of infrared radiation, τ is the infrared optical depth, which increases from zero at the top of the atmosphere to τ s = 0.2 at the surface, and B = σT 4 (τ) is the thermal emission from the atmosphere at the vertical point with optical depth τ. The value of τ s = 0.2 is chosen as a rough approximation to the infra-red absorption in relatively dust-free atmospheric conditions, and is similar to values found in Table 1 of Martin et al. [45]. We suppose that τ varies with height as τ = τ s e −z/H , where z is the height above the surface, and H is the atmospheric scale height, thus neglecting pressure broadening for the sake of simplicity.
In radiative equilibrium, the net radiative flux between the layers is zero, meaning that meaning that the surface temperatures at z = 0, T s (θ, φ, t) are As in the Earth model case, we then introduce a ground temperature, T g , which satisfies where C g is the heat capacity of the surface, and again we use the first formulation here. To then account for the effects of this ground temperature on the atmosphere, we take the lowest-level atmospheric temperatures to be T g and use Equation (11) to recalculate the OLR to be This means that Equation (7) is no longer satisfied at any given time, with radiative imbalance possible due to the storage of energy in the ground. However, this balance will be satisfied in a time-mean, as required by conservation of energy. We then recalculate the atmospheric radiative equilibrium temperatures using Equation (10) with the new expression for OLR. These atmospheric temperatures then form the basis for the Newtonian cooling.
We make two additional simplifications. The first is that we do not include the diurnal cycle in the insolation, and account only for the seasonal variations; this is quite significant for Mars, where thermal tides caused by the diurnal cycle can play a significant role [6]. Second, we impose a lower bound on T g of 149 K, following the simplified Mars model described by Haberle et al. [46]. This lower bound is designed to replicate the effects of CO 2 condensation, which at Martian pressures occurs at 149 K. This lower bound represents the heating produced by this condensation, which would prevent the temperatures from decreasing too much below that value. The lower bound is particularly important in a radiative equilibrium temperature calculation for it stops the equilibrium temperature dropping to absolute zero during the polar night.
Despite these substantial simplifications, the formulation reproduces structure and magnitude of the Mars winter-time zonal-mean zonal wind quite well-compare Figure 4a,d. There are some significant differences, particularly in the eastward zonal jet in the southern hemisphere in Figure 4a that is not present in Figure 4d. However, much of the broad zonal-jet structure is well captured in this simple framework. In addition, it is noted that this model well captures Mars' seasonal cycle in temperatures without a seasonal lag compared with the observations. This is because Mars is a rocky planet without an ocean (at least today), meaning that a realistic seasonal cycle can be well approximated by using a small surface heat capacity.

An Intermediate Complexity Model
The next level up in complexity in our Martian hierarchy is to move to an interactive surface that can exchange heat and momentum with the atmosphere, and to use a gray radiative transfer instead of Newtonian relaxation. In the longwave, we choose a surface longwave optical depth of τ lw0 = 0.2, consistent with the simple model, with the vertical variation given by As an approximation for relatively dust-free atmospheric conditions, we choose a shortwave surface optical depth of τ sw0 = 0.2 (see, e.g., Section 5.

of [6]), with a vertical variation of
meaning the scale height of this absorption is four times smaller than the scale height for infra-red optical depth. This is chosen to concentrate shortwave absorption at lower levels, where dust levels are highest in non-dust-storm conditions (e.g., Figure 2.2b of [6]). This has the effect of warming the lower atmosphere compared with the results using only well-mixed gases. We do not impose a lower bound on the temperature, preferring to keep the system energetically closed. The zonal-mean zonal wind is now a reasonable match to the MACDA reanalysis (see Figure 4b). The southern-hemisphere eastward jet is somewhat stronger than in the reanalysis, but the much stronger westward jet in the upper atmosphere in the tropics now compares well with observations, suggesting temperature gradients are more accurate. This model also accounts for Mars' diurnal cycle, meaning that thermal tides due to the strong diurnal forcing of surface temperatures are present. We do not account for any dust-related shortwave absorption in this model, and no topography is present, meaning that the driving of tides and semi-diurnal tides by shortwave absorption by dust, and through interaction with topography, are absent (see discussion in Section 5.3 of [6]).

A More Comprehensive Model
The final step up our hierarchy involves use of the comprehensive Socrates radiative transfer code, but otherwise keeping the model the same as the grey radiation one. We take the mass-mixing ratio of CO 2 on Mars to be 0.949, and also account for 0.0168 of N 2 . The Socrates spectral files we use for Mars were created for the ROCKE-3D model [47], and include CO 2 − CO 2 collision-induced absorption [48][49][50], CO 2 sub-Lorentzian line wings [50,51], and CO 2 self-broadening.
We also include Martian topography, with height data taken from measurements made by the Mars Orbiter Laser Altimeter (MOLA) aboard NASA's Mars Global Surveyor satellite [52]. The zonal winds for this setup are shown in Figure 4c. As in the simpler model we do not, in the version described here, account for the radiative effects of atmospheric dust, nor the heating effects of CO 2 condensation. Despite these simplifications, the zonal-mean atmospheric temperatures compare reasonably well with those in MACDA, as shown in Figure 5. We capture many of the salient features of Mars atmosphere, including the temperature asymmetry between northern-and southern-hemisphere summers. The absolute north-south temperature gradients are a little out from MACDA, and this may well be due to our lack of representation of dust, and our lack of latent heat release from CO 2 condensation in the winter polar regions. The comprehensive model includes both new radiation and topography compared with the intermediate model. The combination of these two effects changes the time-mean zonal winds by weakening the winter hemisphere jet slightly, and also by weakening the southern-hemisphere eastward jet, with the former looking less like observations, and the latter looking more like observations. The weakening of the winter-hemisphere jet is mostly due to the changes in radiation scheme, as flattening the topography in the comprehensive model still results in a weaker jet compared with the intermediate model (not shown). The summer-hemisphere changes are mostly due to the topography (not shown). The weakening of the winter-hemisphere jet must be due to unrealistic horizontal temperature gradients compared with the intermediate model and MACDA. This may well be to do with Socrates not accounting for the effects of dust, where the intermediate model's optical depths attempt to do this, but further work is required to investigate this further.
In the comprehensive model, we do capture some significant features of Mars' surface pressure variations with season, as shown in Figure 6, even without accounting for the well-known seasonal cycle in total atmospheric mass on Mars due to CO 2 condensation. This omission may well be part of the reason for the large pressure differences between Isca and MACDA, particularly in the northern polar regions.^Ž

Comparison of Models and Observations
The three simulations shown in Figure 4 are an example of the type of Martian hierarchy that can be constructed using Isca. The setups here do not form a "continuous" hierarchy, as more than one feature has been changed between each of them, but constructing such a hierarchy would be relatively straightforward. Isca has the ability to interchange each of the features described; for example, the lower bound on temperature motivated by CO 2 condensation could be implemented in other simulations, and topography could be added to the simplest simulations. Isca may be a powerful tool but the actual construction of a scientifically meaningful hierarchy is a difficult scientific task; it is made easier, but not automatic, by Isca. The use of Isca to better understand Mars' polar vortex is underway and will be presented for publication elsewhere.

Jupiter
Jupiter's multiple jet streams have long fascinated atmospheric dynamicists of all persuasions, with a myriad of studies looking at various possibilities for how these jets form (e.g., [53][54][55][56][57][58]), how they remain stable despite appearing to be barotropically unstable, and how the jets seen at cloud level are related to flows deeper inside the planet, with some studies seeking to combine answers to these questions [59][60][61].
It is interesting to note that many different model types have produced jets that are somewhat Jupiter-like. These models range from single-layer quasi-geostrophic models (e.g., [53,61]), through shallow-water models (e.g., [56]), to those including varying degrees of three-dimensionality, including the primitive-equation models of Schneider and Liu [62] through to anelastic models [63] and finally to deep-convection models such as those of Heimpel et al. [55] and Heimpel et al. [60]. It is not always clear, however, how each of these rungs in the Jupiter hierarchy connect together, and so a joined-up approach of looking at the same question across a range of models would certainly be of value in this area.
An interesting open question about Jupiter's atmosphere is the nature of the interaction between Jupiter's outer "weather-layer" and the deeper layers. This question is particularly timely in the light of new observations of Jupiter's gravitational field by NASA's Juno spacecraft. These observations suggest that the jet streams visible at cloud-level continue below the base of the weather layer (O(10)bars) deep into the planet's interior, to a depth on the order of 3000 km [64], although the uniqueness of this conclusion has been questioned [65].
The fluid beneath the weather layer has a much higher density than that in the weather layer, so that any deep flows have a significantly higher inertia, and therefore almost certainly longer timescales. Given this separation of timescales, one modeling approach is to prescribe the deep flows as a time-constant boundary condition on weather-layer only models. The simplest models of this type are so-called "1.5-layer" models, with the half layer being the prescribed deep flow, with an active layer evolving above, with examples including those of Dowling and Ingersoll [59] and Thomson and McIntyre [61]. This type of approach is less common in more comprehensive models, with the notable exception of the EPIC model of Dowling [66]. However, conducting such experiments through a hierarchy of models seems likely to yield interesting insights, not only into the interaction of Jupiter's layers, but also into the consistency, or not, of the array of Jupiter models in use. With this in mind, the Isca framework has a range of Jupiter models from the simple shallow-water models through to 3D primitive-equation dynamics, all with the option of prescribed deep flow, as we now describe.

A Simple Model
The simplest model in our Jovian hierarchy within Isca is a global 1.5-layer shallow-water model, as described by Thomson [67]. This model's parameters include the equilibrium layer depth (which is a control on the Rossby deformation scale, L D ), the strength of radiative relaxation, the deep jet profile and the properties of the forcing. As an example of the type of results produced by this model, Figure 7a shows the zonal-mean zonal wind profiles in the northern hemisphere of a run with deep jets, small deformation scale, and intermediate-strength radiative damping. This integration, as do most in a Jupiter-like regime, shows significant matching between the weather-layer jets and the prescribed deep jets (see [67] for more detailed analysis). This includes the equatorial regions, where the upper-layer jets super-rotate, reflecting the super-rotating deep jets that were prescribed. The integration also produces very straight weather-layer jets, i.e. they run around latitude circles, as shown in Figure 8a, although this is partly due to the prescribed deep jets also having zonal symmetry. Still, the straightness matches well with the observed straightness of Jupiter's weather-layer jets, as is clear from watching movies of Jupiter's cloud motions (e.g., https://voyager.jpl.nasa.gov/mission/science/jupiter/).
The fact that the weather-layer and deep jets match so closely in this model suggests that it would be difficult for significant vertical shear to exist between these layers on the real Jupiter. Clearly, there are many complicating factors not present in this simple model including possibilities for baroclinic instability, a proper treatment of radiative transfer, etc., all of which might have a significant impact on the results.

An Intermediate Complexity Model
Our intermediate-complexity Jupiter model is a significant jump from the simplest model. It is a primitive equation model with 60 vertical levels, a grey radiative transfer similar to that described by Schneider and Liu [62], a vertical diffusion scheme based on that described by Young et al. [68], and a uniform bottom-boundary heat flux of 5.7 Wm −2 . The model differs from other conventional Jupiter GCMs in that it has prescribable jets in the bottom model level, which is at 14.5 bars in the present configuration. These jets are created using a Rayleigh drag back towards a prescribed jet profile on a short timescale of 10 Earth days. Enforcing such jets is not necessary to achieve barotropic Jupiter-like flows within the weather layer [69]. However, we include prescribable deep jets as an option within Isca to investigate how the weather-layer's jet flows respond to the presence of deep jets, much like in the simple model described previously.
Zonal-mean zonal-wind profiles from an example run with this model are shown in Figure 7b. Despite the inclusion of several important processes compared with the simple model, the results are actually very similar, with the jets within the weather-layer matching the prescribed deep jets rather closely, including the superrotating equatorial jet. The barotropic structure of the jet streams in this model is shown in more detail in Figure 7c. The jets in this model are also very straight, as shown in Figure 8b. Although consistent with the simple model, the straightness is more surprising given the increase in the number of processes that could cause longitudinal variations, particularly baroclinic instability. Further work is underway looking at the factors affecting the vertical coupling of Jupiter's weather-layer and deep jets using both the intermediate and simple models, and results are being prepared for publication.

Toward a More Comprehensive Model
Although, compared to Earth models, our grey-radiation Jupiter model is of "intermediate complexity", it is actually one of the more comprehensive models of that planet, attesting to the difficulty of building models of Jupiter that are both comprehensive and meaningful. Compared to our model, the models of Young et al. [68] and Palotai and Dowling [70] include more sophistical treatments of moist processes and representations of clouds, and the ambitious deep atmosphere model of Heimpel et al. [60] simulates both the weather layer and the deep jets, albeit while making various assumptions about the density structure and using simple physics schemes.
Clearly, there is scope for more comprehensive models of Jupiter. The role of moisture, for example, is far from understood. Observational studies have argued that moist convection due to water could help carry Jupiter's internal heat across the weather layer and out to space [71,72]. On the other hand, Jupiter's moist adiabatic lapse rate is very close to the dry adiabatic lapse rate due to the small molecular weight of Jupiter's hydrogen atmosphere, suggesting rather weak moist effects; suffice it to say that the question of the importance of moisture is an open one. Properly incorporating all the effects of moisture is quite a challenge because its physics is extremely complex, potentially involving computationally-expensive microphysical calculations and small spatial scales. This, combined with the geophysical-fluid-dynamics problem of understanding deep and shallow jets, might well cause the rational investigator to despair. A hierarchical approach may be the only way to proceed, and effecting that is a long-term goal of the authors.

Discussion and Conclusions
In this paper we describe the use of Isca to construct models of varying complexity for three Solar System planets, Earth, Mars, and Jupiter. Hierarchical modeling of Earth's atmosphere and climate is a well established field, with notable contributions and endorsements by Schneider and Dickinson [73] and Hoskins [74]; reinforced by Held [75]; discussed by Vallis [76]; and reviewed by Maher et al. [10]. Hierarchical modeling is, in our view, even more important if we wish to understand and simulate (not always the same thing) planets other than our own, for the reasons discussed in the Introduction. There is a huge diversity of planetary atmospheres and data are sufficiently sparse (relative to that for Earth) that comprehensive modeling is difficult because of this lack of observational constraints. However, data are becoming available, both from Extra Large Telescopes on Earth, orbiting telescopes (Hubble and hopefully James Webb) (as reviewed in, e.g., [9], and references therein), and from missions to Solar System planets, such as the Exo-Mars Trace Gas Orbiter [77], Juno [78], and Dragonfly [79], the proposed exploration of Titan's atmosphere with a drone. Thus, quantitative comparison with data is becoming possible and therefore is, in conjunction with seeking a basic understanding, a goal of our modeling activities.
The twin goals of understanding and of quantitative simulation cannot, in our view, be achieved with a single level of complexity. Our goal in constructing Isca is to make a tool that enables the scientist to construct models of appropriate complexity for the task at hand. Sometimes that task will necessitate construction of a suite of models of a single planet, and sometimes that suite, taken as a whole, may essentially constitute our theory of the planetary atmosphere, or some aspect of it.
In this paper, we focus on Earth, Mars, and Jupiter because the modeling of their atmospheres is quite a severe test for any modeling system. Venus and Titan are other candidate atmospheres of different kinds that we expect to address in the future. Exoplanets are another obvious target for hierarchical modeling, and, although we used Isca in some preliminary studies (e.g., [80]), much more can be done. We concentrate mainly on the use of the primitive equations, although for Jupiter we include a shallow water model. The various models then differ mainly in the sophistication of their "physics" and their parameterizations of subgrid processes, varying from Newtonian relaxation to multi-band radiation schemes, and from aquaplanets and flat rocky planets, to planets with continents and topography. We do not claim that the choices we have made here are necessarily the "best" ones; rather, our intent is to give a flavor of what can be done and how it may be carried further. More detailed descriptions of some of the models will appear in the specialist literature.
In addition to its capabilities for simulating planets, Isca is also an easy-to-use framework, with runscripts written in Python, an automated code testing suite provided, and a selection of test cases that sample the range of possible planets configurable with Isca. We would therefore like to encourage other users to build and/or use models of any and all planets using Isca, including those described here if desired-the framework is sufficiently flexible that quite different models of the same planet could be constructed. The framework is fully open source and available online, with versions configured to run on many different architectures, including native versions for Mac computers, several example setups on Linux machines, including supercomputers, and the ability to run within a Docker environment. The main requirements for model compilation are a Fortran compiler (e.g., GFortran), the Message-Passing Interface (MPI), and a netCDF library for output. We hope that the above features, alongside Isca's scientific capabilities, will make it a useful community tool for modeling of planetary atmospheres.