MODELS OF VARYING COMPLEXITY: FROM
VOTER NETWORKS TO EXTRASOLAR PLANETS
A Dissertation
Presented to the Faculty of the Graduate School
of Cornell University
in Partial Fulfillment of the Requirements for the Degree of
Doctor of Philosophy
by
Ekaterina Landgren
December 2022
© 2022 Ekaterina Landgren
ALL RIGHTS RESERVED
MODELS OF VARYING COMPLEXITY: FROM VOTER NETWORKS TO
EXTRASOLAR PLANETS
Ekaterina Landgren, Ph.D.
Cornell University 2022
Mathematical modeling and techniques from the field of dynamical systems can
be used to study a variety of phenomena. The applications included in this dis-
sertation span from social to physical sciences. The first chapter presents a con-
ceptual model of voter turnout. Using numerical and analytical techniques, we
explore a variety of network structures and identify the conditions that facilitate
minority factions winning elections. The remainder of the dissertation focuses
on mathematical models in atmospheric science, specifically one-dimensional
and two-dimensional models of exoplanet atmospheres. We compare two one-
dimensional energy balance models with explicit dependence on obliquity to
study the likelihood of different stable ice configurations. We compare the re-
sults of models with different methods of heat transport and different insolation
distributions and show that stable partial ice cover is possible for any obliquity,
provided the insolation distribution is sufficiently accurate. In the final chap-
ters of this dissertation, we present a two-dimensional shallow-water model,
SWAMPE (Shallow-Water Atmospheric Model in Python for Exoplanets). Two-
dimensional shallow-water models fill the gap between minimal-complexity
models which lack longitudinal variation and high-complexity computationally
expensive three-dimensional models. Our model can accurately and rapidly ex-
plore a vast parameter space. The code is flexible, modular, built fully in Python,
and could be easily adapted to model a variety of dissimilar space objects, from
Brown Dwarfs to smaller, terrestrial planets. The code can produce the ther-
mal maps necessary for the generation of phase curves and secondary eclipse
maps, which will help constrain and make predictions for observations. Finally,
we apply SWAMPE to explore the circulation regimes of synchronously rotat-
ing sub-Neptunes as we consider the interactions between planetary rotation
period and radiative forcing from the host star.
BIOGRAPHICAL SKETCH
Ekaterina (Kath) Landgren grew up in Moscow, Russia, before attending
Loomis Chaffee high school in the United States. She obtained a Bachelor of
Science degree from Brown University in Providence, Rhode Island, majoring
in Applied Mathematics and Philosophy. She enrolled in the PhD program at
Cornell’s Center for Applied Mathematics, where she explored her interest in
mathematical modeling and dynamical systems. In the third year of her PhD,
she began to focus on mathematical models of atmospheres—a pursuit which
ultimately led to a fruitful collaboration with exoplanetary scientists. She hopes
to continue working in the interdisciplinary corners of mathematics.
iii
This work is dedicated to anyone searching through this thesis on a quest for
very specific information. I hope you find what you are looking for.
iv
ACKNOWLEDGEMENTS
While this dissertation bears my name only, I am thankful for a supportive
and collaborative environment I have found during the last five years.
Firstly, I would like to express deep gratitude to my mentors. My advisor
Steve Strogatz’s mentorship has enabled me to take an unusual trajectory and
to figure out what I was truly interested in. Like a whole generation of students,
I came to the field dynamical systems by reading his undergraduate textbook,
Nonlinear Dynamics and Chaos, and its charm still has a hold on me. Steve’s
taste in problems and the drive to zero in on most compelling questions have
helped me gain the confidence to pursue and communicate research in a way
that draws people in.
Nikole Lewis was an amazing mentor and role model to me. She took me on
as a student when I knew nothing about exoplanets. Nikole is an exceptional
scientist as well as an exceptional manager. Working with her, I always felt that
my goals are a priority. I am lucky to have enjoyed the environment she created
in the Substellar Explorers Group. Thank you, Nikole, for teaching me about
planetary science, for supporting me in whatever directions I wanted to grow,
for taking a well-thought-out approach to mentorship, and for telling me that it
will be fine (as often as I needed to hear it).
Alice Nadeau will deny this, but I don’t think I would have been able to
graduate without her compassionate mentorship. She provided focus to my
mathematical explorations when I was struggling with the breadth of my in-
terests. Alice guided me through writing my first graduate manuscript, helped
me overcome many mathematical hurdles, and showed incredible persistence
in the face of dysfunctional code. It was so important for me to have a postdoc
mentor who could empathize closely with the challenges of a doctoral program.
v
Alice, thank you for introducing me to mathematics of climate and to PokeLava,
and also for glimpses of your cats on Zoom throughout the pandemic.
Tiffany Kataria has brought both levity and expertise to our weekly meet-
ings, and worked very hard to help me write manuscripts and proposals. Her
sharp insights have been both a resource and an inspiration for much of my sci-
entific writing. She was supportive and enthusiastic even when the project was
at a standstill.
I would like to thank Richard Rand and Toby Ault for serving on my special
committee, and David Bindel, for stepping in.
I am lucky to have had great mentors beyond the projects represented in this
dissertation. I am grateful to my mentors at the applied mathematics depart-
ment at Brown. I would like to thank Bjorn Sandstede, for taking undergradu-
ates seriously, and for his good advice on choosing a graduate program. I thank
John Gemmer for introducing me to the fun of agent-based modeling and for
making a whole summer’s worth of eight a.m. meetings. I would like to thank
my mentors and collaborators from the MCRN summer school, especially John
Gemmer, Mary Silber, and Katie Slyman. I would like to thank Alice Schwarze
for her efforts to create a welcoming and supportive environment for women in
network science.
I have chosen the Center for Applied Math for its community, and I am
glad that my first impressions proved true. I would like to thank Erika Fowler-
Decatur, who is the single most competent person I have ever met, and whose
kindness and genuine interest in CAM and its denizens were a truly sustaining
life force.
Mallory Gaspard, Max Ruth, Gokul Nair, I hope you forgive me for being a
force of distraction in our desk clump. Thank you for all the lunches and dinners
vi
in the CAM space, for weekend company, for rubber-ducking, and for making
CAM a better place. I am deeply grateful to Shriya Nagpal for her compassion,
support, commiseration, great conversations and for reminding me that I am a
dynamicist when I needed a boost. Irena Papst was my CAM mentor in both
formal and informal ways, always led by example, and was there for good con-
versation. Shawn Ong helped me move when we barely knew each other, can
be relied on to make the best jokes, and is a great person to have on a trivia team.
Jonas L. Juul was a wonderful mentor and collaborator. I hope we get to keep
running into each other at complex systems events. I would like to express my
appreciation of Lindsay Mercer, whose presence always brightened CAM, who
is a truly independent thinker, and whose keen eye for justice makes her an asset
to whatever space she is in. I would like to thank a few other CAM students and
postdocs: Max Jenquin for good advice, calming presence, and for all the Satur-
days we’ve spent working on manifolds; Andrew Horning, for elucidating the
depths of Evans for me; Heather Wilber, for being a fantastic role model; Eliz-
abeth Wesson, for their sewing mentorship, cutting wit, and queer community;
Tayler Fernandes Nunez for great workouts, enthusiasm, and commitment to
making math better; Jackson Kulik, for bringing a relaxed air to our space; Car-
los Martinez, for international student solidarity; Kate Meyer, for her advice of
“parking downhill.” Among my astronomy collaborators, I would like to thank
Trevor Foote, for his sober coding advice, Ishan Mishra, for accountability and
commiseration, and Ryan MacDonald, for his unbridled enthusiasm. I had the
privilege of working with brilliant undergraduate researchers throughout my
time at Cornell. I would like to thank Thomas Mitchell, Anushka Naranyan,
and Anna Asch for their collaboration and trust.
I am lucky to have had many friendships that sustained me through the last
vii
five years. Arielle Johnson is always ready to help a friend, is a great conver-
sationalist, and throws the best parties. Aida Feng has brightened many of my
afternoons by sharing her sharp observations and letting me keep in touch with
Germanistik, if only vicariously. Anastasia Bazhenova and Alina Tabernilia
know how to maintain a cross-continental friendship. Debbie Riccioli is a great
hiking buddy for not-very-fast hikers like myself. I’m glad to continue know-
ing Daniela Rakhlina-Powsner as we grow out of some idiosyncrasies and into
others. Sara Hartse inspires me in matters small and large to be a better person.
I would like to thank my parents for sending me away to study across the
ocean, and my brother, who eventually came to terms with it. I thank my
mother-in-law, Karin Landgren, for her excitement and for letting me use her
last name. I would like to express my appreciation of my cat, Fennel, for pro-
viding warmth when I needed it. Finally, and most importantly, I would like to
thank my partner, Benedict – you are my family. I am lucky that we get to walk
through life together.
viii
CONTENTS
Biographical Sketch . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . iii
Dedication . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . iv
Acknowledgements . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . v
Contents . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . ix
List of Tables . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . xi
List of Figures . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . xii
1 Introduction 1
2 How a minority can win: Unrepresentative outcomes in a simple
model of voter turnout 6
2.1 Introduction . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 6
2.2 The model . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 9
2.3 Model Networks . . . . . . . . . . . . . . . . . . . . . . . . . . . . 11
2.3.1 Erdős-Rényi networks . . . . . . . . . . . . . . . . . . . . . 12
2.3.2 Stochastic block networks . . . . . . . . . . . . . . . . . . . 14
2.3.3 Networks with a heavy-tailed degree distribution . . . . . 22
2.3.4 Geometric Random Networks . . . . . . . . . . . . . . . . 24
2.4 Discussion . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 27
3 Comparison of Two Analytic Energy Balance Models 29
3.1 Introduction . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 29
3.2 Governing Equations . . . . . . . . . . . . . . . . . . . . . . . . . . 32
3.2.1 Annual Average Energy Balance . . . . . . . . . . . . . . . 32
3.2.2 Analysis of Temperature Equilibria . . . . . . . . . . . . . 41
3.3 Likelihood of Stable Partial Ice Cover . . . . . . . . . . . . . . . . 46
3.3.1 Defining the Region of Integration . . . . . . . . . . . . . . 46
3.3.2 Relative Likelihood of Stable Partial Ice Cover . . . . . . . 47
3.4 Partial Ice Cover and the Snowball State . . . . . . . . . . . . . . . 53
3.5 Discussion . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 57
3.6 Conclusion . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 61
4 SWAMPE: a Shallow-Water Atmospheric Model in Python for Exo-
planets 64
5 Exploring the interaction of rotation rate and stellar irradiation on syn-
chronously rotating sub-Neptunes 67
5.1 Introduction . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 67
5.2 Methods . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 70
5.2.1 Model . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 70
5.2.2 Model validation . . . . . . . . . . . . . . . . . . . . . . . . 75
5.2.3 Parameter regimes . . . . . . . . . . . . . . . . . . . . . . . 77
5.2.4 Motivation for the choices of Φ . . . . . . . . . . . . . . . . 80
ix
5.3 Results . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 82
5.3.1 Extension to Sub-Neptunes . . . . . . . . . . . . . . . . . . 82
5.3.2 Strong forcing simulations . . . . . . . . . . . . . . . . . . . 85
5.3.3 Medium forcing simulations . . . . . . . . . . . . . . . . . 91
5.3.4 Weak forcing simulations . . . . . . . . . . . . . . . . . . . 95
5.4 Discussion . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 101
5.4.1 Holistic discussion of trends in the regimes . . . . . . . . . 101
5.4.2 Trends in the zonal wind flow . . . . . . . . . . . . . . . . . 106
5.4.3 Idiosyncratic spin-up behavior for long timescales . . . . . 108
5.4.4 Equatorial and Tropical Band Oscillation . . . . . . . . . . 110
5.5 Conclusions . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 111
A Algorithm for generating scale-free networks with homophily 114
B Comparison of the Analytical Equilibrium Solutions 115
C An Open-Source Software Repository for SWAMPE 118
C.1 Documentation . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 118
C.2 Time-stepping . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 118
C.3 Filters . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 123
D Supplementary figures showing geopotential contours and mean-
zonal winds for intermediate scale heights. 124
x
LIST OF TABLES
3.1 Nondimensional parameters for the relaxation to the mean and
diffusion versions of the energy balance model. . . . . . . . . . . 40
5.1 Typical maximal wind speeds, day/night geopotential contrast,
and Rossby number for a representative scale height Φ = 4 × 106
m2/s2 based on temporal averages. . . . . . . . . . . . . . . . . . 83
xi
LIST OF FIGURES
1.1 Hierarchy of atmospheric models. . . . . . . . . . . . . . . . . . . 3
2.1 Illustration of behavioral assumptions for the voter turnout model. 9
2.2 The proportion of unrepresentative outcomes on Erdős-Rényi
random graphs. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 13
2.3 The proportion of unrepresentative outcomes on stochastic block
networks. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 16
2.4 The proportion of unrepresentative outcomes on stochastic block
networks for the special case pin = 1. . . . . . . . . . . . . . . . . . 20
2.5 Figure showing limiting behavior for the proportion of unrepre-
sentative outcomes. . . . . . . . . . . . . . . . . . . . . . . . . . . 21
2.6 Examples of networks with heavy-tailed degree distributions
with varying homophily factors. . . . . . . . . . . . . . . . . . . . 24
2.7 Proportion of minority winning on networks with heavy-tailed
degree distributions. . . . . . . . . . . . . . . . . . . . . . . . . . . 25
2.8 Example of majority and minority node distributions for geo-
metric random networks . . . . . . . . . . . . . . . . . . . . . . . 25
2.9 Proportion of minority victories on a geographic random network. 27
3.1 Illustration of planets exhibiting ice caps and ice belts. . . . . . . 35
3.2 Plots showing the bifurcation diagrams for Snowball transitions. 44
3.3 Plots showing the relative likelihood of stable partial ice cover
compared to partial stable ice cover on Earth. . . . . . . . . . . . 48
3.4 Plots showing stability of the equilibrium ice line locations for
β = 54.47◦ for the diffusion model for ice caps (top row) and ice
belt (bottom row) . . . . . . . . . . . . . . . . . . . . . . . . . . . . 52
3.5 Contour plots showing the minimum value of the nondimen-
sional incoming radiation q for which the Snowball state is not
globally stable. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 54
4.1 Illustration of one time step using the spectral method employed
by SWAMPE . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 65
5.1 Recreation of a hot Jupiter regime in Figure 3 in Perez-Becker and
Showman (2013). Shown in the figure are maps of steady-state
geopotential contours and the wind vector fields for the equili-
brated solutions of the shallow-water model. . . . . . . . . . . . . 76
5.2 Planetary equilibrium temperature and as a function of rotation
period for planets around K0, K5, M0, and M5 stars. . . . . . . . 81
5.3 Recreation of a hot Jupiter regime, with planetary radius turned
down to a = 3a⊕ = 1.3 × 107 m. . . . . . . . . . . . . . . . . . . . . 85
xii
5.4 Averaged maps of geopotential contours and wind vector fields
for the equilibrated solutions of the shallow-water model for the
strong forcing regime (∆Φeq/Φ = 1) at reference geopotential Φ =
4 × 106 m2s−2. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 87
5.5 Mean-zonal winds U for the equilibrated solutions of the
shallow-water model for the strong forcing regime (∆Φeq/Φ = 1)
at reference geopotential Φ = 4 × 106 m2s−2. . . . . . . . . . . . . 88
5.6 Same as Figure 5.4, but at a low scale height corresponding to
Φ = 106 m2s−2. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 90
5.7 Same as Figure 5.5 but at a low scale height corresponding to
Φ = 106 m2s−2. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 91
5.8 Averaged maps of geopotential contours and wind vector fields
for the equilibrated solutions of the shallow-water model for the
medium forcing regime (∆Φeq/Φ = 0.5) at reference geopotential
Φ = 4 × 106 m2s−2. . . . . . . . . . . . . . . . . . . . . . . . . . . . . 93
5.9 Mean-zonal winds U for the equilibrated solutions of the
shallow-water model for the medium forcing regime (∆Φeq/Φ =
0.5) at reference geopotential Φ = 4 × 106 m2s−2. . . . . . . . . . . 94
5.10 Same as Figure 5.8, but at a low scale height corresponding to
Φ = 106 m2s−2. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 96
5.11 Same as Figure 5.9 but at a low scale height corresponding to
Φ = 106 m2s−2. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 97
5.12 Averaged maps of geopotential contours and wind vector fields
for the equilibrated solutions of the shallow-water model for the
weak forcing regime (∆Φeq/Φ = 0.1) at reference geopotentialΦ =
4 × 106 m2s−2. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 98
5.13 Mean-zonal winds U for the equilibrated solutions of the
shallow-water model for the weak forcing regime (∆Φeq/Φ = 0.1)
at reference geopotential Φ = 4 × 106 m2s−2. . . . . . . . . . . . . 99
5.14 Same as Figure 5.12, but at a low scale height corresponding to
Φ = 106 m2s−2. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 100
5.15 Same as Figure 5.13 but at a low scale height corresponding to
Φ = 106 m2s−2. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 101
5.16 Snapshots of spin-up behavior for strong-forcing regime with
Φ = 4 × 106 at rotation period Prot = 10 days and radiative
timescale τrad = 10 days. . . . . . . . . . . . . . . . . . . . . . . . . 109
5.17 Snapshots of oscillatory behavior for a regime corresponding to
a weakly forced planet with a scale height corresponding to Φ =
4 × 106 at rotation period Prot = 1 days and radiative timescale
τrad = 10 days. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 111
xiii
5.18 The estimated contours where τadv = τrad for different day/night
contrasts of the local radiative equilibrium for the strong,
medium, and weak forcing regimes. . . . . . . . . . . . . . . . . . 113
B.1 Comparison of equilibrium temperature profiles for the diffu-
sion model and the relaxation to the mean model. . . . . . . . . . 116
C.1 Screenshot of SWAMPE documentation’s homepage. . . . . . . . 119
D.1 Averaged maps of geopotential contours and wind vector fields
for the equilibrated solutions of the shallow-water model for the
strong forcing regime (∆Φeq/Φ = 1) at reference geopotential Φ =
3 × 106 m2s−2. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 125
D.2 Mean-zonal winds U for the equilibrated solutions of the
shallow-water model for the strong forcing regime (∆Φeq/Φ = 1)
at reference geopotential Φ = 3 × 106 m2s−2. . . . . . . . . . . . . . 126
D.3 Averaged maps of geopotential contours and wind vector fields
for the equilibrated solutions of the shallow-water model for the
strong forcing regime (∆Φeq/Φ = 1) at reference geopotential Φ =
2 × 106 m2s−2. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 127
D.4 Mean-zonal winds U for the equilibrated solutions of the
shallow-water model for the strong forcing regime (∆Φeq/Φ = 1)
at reference geopotential Φ = 2 × 106 m2s−2. . . . . . . . . . . . . 128
D.5 Averaged maps of geopotential contours and wind vector fields
for the equilibrated solutions of the shallow-water model for the
medium forcing regime (∆Φeq/Φ = 0.5) at reference geopotential
Φ = 3 × 106 m2s−2. . . . . . . . . . . . . . . . . . . . . . . . . . . . . 129
D.6 Mean-zonal winds U for the equilibrated solutions of the
shallow-water model for the medium forcing regime (∆Φeq/Φ =
0.5) at reference geopotential Φ = 3 × 106 m2s−2. . . . . . . . . . . 130
D.7 Averaged maps of geopotential contours and wind vector fields
for the equilibrated solutions of the shallow-water model for the
medium forcing regime (∆Φeq/Φ = 0.5) at reference geopotential
Φ = 2 × 106 m2s−2. . . . . . . . . . . . . . . . . . . . . . . . . . . . . 131
D.8 Mean-zonal winds U for the equilibrated solutions of the
shallow-water model for the medium forcing regime (∆Φeq/Φ =
0.5) at reference geopotential Φ = 2 × 106 m2s−2. . . . . . . . . . . 132
D.9 Averaged maps of geopotential contours and wind vector fields
for the equilibrated solutions of the shallow-water model for the
weak forcing regime (∆Φeq/Φ = 0.1) at reference geopotentialΦ =
3 × 106 m2s−2. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 133
xiv
D.10 Mean-zonal winds U for the equilibrated solutions of the
shallow-water model for the weak forcing regime (∆Φeq/Φ = 0.1)
at reference geopotential Φ = 3 × 106 m2s−2. . . . . . . . . . . . . 134
D.11 Averaged maps of geopotential contours and wind vector fields
for the equilibrated solutions of the shallow-water model for the
weak forcing regime (∆Φeq/Φ = 0.1) at reference geopotentialΦ =
2 × 106 m2s−2. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 135
D.12 Mean-zonal winds U for the equilibrated solutions of the
shallow-water model for the weak forcing regime (∆Φeq/Φ = 0.1)
at reference geopotential Φ = 2 × 106 m2s−2. . . . . . . . . . . . . . 136
xv
CHAPTER 1
INTRODUCTION
Applied mathematics is a versatile and diverse field that has enabled me
to explore a variety of disciplines through the lens of mathematical model-
ing. Mathematical modeling is a fundamentally translational work. It trans-
lates real-world phenomena into the language of mathematics, and it requires
back-and-forth communication between mathematicians and scientists. While
the applications included in this thesis range from voter turnout to exoplanet
atmospheres, a dynamical systems perspective is a persistent common thread
throughout this work.
Despite existing on opposite ends of the spectrum from social to physical
science, each part of this thesis is grounded in a dynamical systems approach.
In each of the applications, we first carefully choose how to summarize the full
complexity of the real world in a few relevant mechanisms. Then, using nu-
merical and analytical approaches, we explore the parameter space and analyze
the emergent patterns and phenomena. Mathematical modeling tends to be
most useful where real-world information is limited or incomplete. The work
in this thesis complements the insights we have into the complexities of human
decision-making and the properties of the atmospheres of distant planets.
In Chapter 2, I present a model of voter turnout. Our conceptual, qualitative,
agent-based model examines the effect of voter turnout on election outcomes.
The outcome of an election depends not only on which candidate is more pop-
ular, but also on how many of their voters actually turn out to vote. In our
simple model voters abstain from voting if they think their vote would not mat-
ter. Specifically, they do not vote if they feel sure their preferred candidate will
1
win anyway (a condition we call complacency), or if they feel sure their candi-
date will lose anyway (a condition we call dejectedness). The voters reach these
decisions based on a myopic assessment of their local network, which they take
as a proxy for the entire electorate: voters know which candidate their neigh-
bors prefer and they assume—perhaps incorrectly—that those neighbors will
turn out to vote, so they themselves cast a vote if and only if it would produce
a tie or a win for their preferred candidate in their local neighborhood. We ex-
plore various network structures and distributions of voter preferences and find
that certain structures and parameter regimes favor unrepresentative outcomes
where a minority faction wins, especially when the locally preferred candidate
is not representative of the electorate as a whole.
The rest of this thesis focuses on physical science applications, specifically on
atmospheric models of varying complexity. Figure 1.1 illustrates the complex-
ity hierarchy of atmospheric models and summarizes the advantages of each
complexity level, from one- to three-dimensional models. A key takeaway is
that higher-complexity models are not universally better. Mathematical mod-
eling often involves trade-offs between complexity, specificity, computational
cost, and the ability to diagnose the underlying dynamical mechanisms.
On the lower end of the modeling hierarchy, there are one-dimensional mod-
els. They can capture complex phenomena (e.g. Bell and Cowan (2018)) and
can rapidly explore the parameter space, but they fail to account for longitudi-
nal variation. One distinct advantage of one-dimensional models is their (rel-
ative) mathematical tractability. With just one spatial dimension, these models
lend themselves to explicit mathematical analysis (e.g. Widiasih (2013)), which
brings a robust foundation to the numerical results.
2
Figure 1.1: Hierarchical list of atmospheric models for calculating atmo-
spheric dynamics. Most models used are either one-dimensional and three-
dimensional, and are either low or high in complexity. The work presented
in this thesis spans one-dimensional and 2D models. First, in Chapter 3, we
compare two analytical, mathematically tractable energy-balance models and
analyze their implications for planetary ice cover. Chapter 4 presents a 2D
shallow-water model and the software I have developed to conduct simulations.
Chapter 5 contains application of this software in the context of sub-Neptunes,
a population of planets whose atmospheres will be readily characterized with
current and future telescopes.
In Chapter 3, I present an analysis of two energy-balance models used to
study planetary ice cover. We compare two analytic energy balance models
with explicit dependence on obliquity to study the likelihood of different stable
ice configurations. We compare the results of models with different methods of
heat transport and different insolation distributions. We show that stable par-
tial ice cover is possible for any obliquity, provided the insolation distribution is
sufficiently accurate. Additionally, we quantify the severity of the transition to
the Snowball state as different model parameters are varied. In accordance with
an earlier study, transitions to the Snowball state are more severe for higher val-
ues of the albedo contrast and energy transport across latitudes in both model;
however, we find that the Snowball transition is not equally likely across both
models. This work is general enough to be used to study the likelihood of Snow-
ball transitions on planets within the habitable region of other stars.
3
However, one-dimensional energy balance models have their limitations.
For example, recent observations of giant exoplanets have shown that one-
dimensional models cannot completely describe some of the key processes (e.g.
Feng et al. (2016)) that take place in their atmospheres. More than one physical
dimension is needed to model such systems accurately.
On the opposite end of the complexity hierarchy, there are sophisticated,
three-dimensional models. Such models can capture variation in the physi-
cal space. They are frequently based on primitive equations (e.g. Menou and
Rauscher (2009), Kataria et al. (2016b), Parmentier et al. (2013)) or on the Navier-
Stokes equations (e.g. Cooper and Showman (2006), Dobbs-Dixon and Agol
(2013)) and can be used to understand a variety of radiative, chemical, and dy-
namical processes. However, three-dimensional models tend to be computa-
tionally expensive, sometimes taking months to explore the parameter space.
The difference in capability between one-dimensional and three-dimensional
models leaves a natural gap for two-dimensional shallow-water models, which
can capture the spatial variability as well as run fast enough to rapidly explore
the parameter space and study the underlying dynamical mechanisms.
In Chapter 4, I describe SWAMPE (Shallow-Water Atmospheric Model
in Python for Exoplanets), the software I have developed to conduct
two-dimensional simulations of planetary atmospheres. SWAMPE is two-
dimensional shallow-water equations model. Two-dimensional shallow-water
models fill the gap between minimal-complexity models which lack lon-
gitudinal variation and high-complexity computationally expensive three-
dimensional models. Many previous such models are built in Fortran, which
makes them difficult to adapt to the needs of exoplanet science. The code devel-
4
oped here is flexible, modular, built fully in Python, and could be easily adapted
to model a variety of dissimilar space objects, from Brown Dwarfs to smaller,
terrestrial planets. The Python code can produce the thermal maps necessary
for the generation of phase curves and secondary eclipse maps, which will help
constrain and make predictions for observations. SWAMPE is currently hosted
in an open-source repository and is readily available for use by planetary scien-
tists.
In Chapter 5, I present the application of SWAMPE to a range of possible
synchronously rotating sub-Neptunes. We explore the interaction of planetary
rotation period, radiative forcing, and scale height. The circulation regimes vary
from day-to-night flow to a jet-dominated flow. Planets with lower forcing tend
to exhibit more temporal variability. These simulations will aid in the inter-
pretation of exoplanet observations with the James Webb Space Telescope and
other ground- and space-based facilities.
5
CHAPTER 2
HOW A MINORITY CAN WIN: UNREPRESENTATIVE OUTCOMES IN A
SIMPLE MODEL OF VOTER TURNOUT
Originally published in:
Landgren, E., Juul, J., Strogatz, H. How a minority can win: Unrepresentative
outcomes in a simple model of voter turnout. Physical Review E 104.5 (2021):
054307. https://doi.org/10.1103/PhysRevE.104.054307
2.1 Introduction
Election forecasting is a difficult problem with real-world consequences (Wang
et al., 2015; Rothschild, 2009; Volkening et al., 2020). Part of the difficulty is that
human psychology is murky. How do voters decide which candidate they pre-
fer? What makes them change their minds? And how do they decide whether
to tell pollsters what they really think? More broadly, modeling elections and
voter behavior can shed light on a wide range of puzzling issues about human
decision-making and hot-topic phenomena such as polarization and the forma-
tion of political echo chambers (Baumann et al., 2020; Madsen et al., 2018; Hayat
and Samuel-Azran, 2017; Heider, 1946; Kułakowski et al., 2005; Marvel et al.,
2011; De Marzo et al., 2020; Del Vicario et al., 2017).
There is a rich literature on agent-based opinion dynamics. This literature
includes the “voter model” of probability theory (Holley and Liggett, 1975) and
its many extensions (see Redner (2019) for a review), as well as bounded confi-
dence models (Lorenz, 2007; Hegselmann and Krause, 2002; Hickok et al., 2021).
In such models, agents interact on a network and change their opinions accord-
ing to certain rules. For example, the agents can adopt the opinion of one of
their nearest neighbors chosen at random (Holley and Liggett, 1975), or they
6
can adopt the opinion held by the majority of their neighbors (Krapivsky and
Redner, 2003; Galam, 2008), or they can update their opinion at a nonlinear rate
depending on the opinion distributions of their neighbors (Lambiotte and Red-
ner, 2007, 2008; Slanina et al., 2008). The update rules can also depend on the
state of agents’ opinions (e.g., introducing stubborn (Masuda et al., 2010) or
confident (Volovik and Redner, 2012) voters who do not change their opinions
easily). A key question is when consensus forms among the nodes and what
conditions promote it.
However, opinion dynamics is just one facet of voter behavior. In the real
world, another important factor is voter turnout, defined as the percentage of
eligible voters who cast a ballot in an election. The turnout rate depends on
many socioeconomic, political, and institutional factors, from population size
to campaign expenditures to registration requirements (Geys, 2006; Cancela and
Geys, 2016). The abundance of relevant factors makes predicting voter turnout
difficult.
One factor influencing voter turnout is the closeness of the election (Ender-
sby et al., 2002; Bursztyn et al., 2017). Intuitively, one might expect that close
elections should produce higher turnout, but some scholars dispute that this is
the case (Matsusaka, 1993; Gerber et al., 2020). Here we explore the effect of net-
work structure on individual agents’ perceptions of election closeness and the
consequent impact on turnout and on the election itself.
Certain network structures and opinion distributions can lead to minority
nodes mistakenly believing that they belong to a majority. The phenomenon
whereby local knowledge of the network is not representative of the electorate
as a whole is known as the “majority illusion”; a “minority illusion” is also pos-
7
sible Lerman et al. (2016). We are interested in conditions that allow a minority
to win elections by generating a higher turnout than the majority.
The phenomenon of the minority defeating the majority has been studied
previously in many ways. For example, Iacopini et al. (Iacopini et al., 2021)
examine when a minority can build a critical mass to cause a cascade on hyper-
graphs and become the dominant opinion. In a similar spirit, Touboul (Touboul,
2019) and Juul and Porter (Juul and Porter, 2019) examine how antiestablish-
ment nodes (nodes that prefer to belong to a minority) can spread their influence
and create an antiestablishment majority.
In this chapter we consider a model of voter turnout that allows for majority
and minority illusions. We ask: What network structures enable minority fac-
tions to win? While we do not consider opinion dynamics (our model voters
never change their minds), the mechanisms of voter turnout alone can generate
situations where a small minority can win in a landslide. This counterintuitive
result is one of our main findings. Whether it holds in more realistic models
remains to be seen.
The chapter is laid out as follows. Section 2.2 introduces the model. In Sec-
tion 2.3, we apply the model on a variety of network structures: Erdős-Rényi
networks (2.3.1), stochastic block networks (2.3.2), scale-free networks (2.3.3),
and random geometric networks (2.3.4). Section 2.4 summarizes and discusses
the results.
8
Figure 2.1: The behavioral assumptions of (a) dejectedness, (b) complacency,
and (c) their combination applied to the same ring network with a 5 − 4 split
between orange and purple nodes. (a) The top orange node is surrounded by a
purple node on either side. Thus, in its local (one-hop) neighborhood it is out-
numbered 2 − 1, so its vote cannot tie or win the upcoming election in that local
neighborhood. Making a myopic (and wrong) estimate of the orange opinion’s
chances globally, based solely on its local neighborhood, the orange node be-
lieves its vote cannot affect the upcoming election, so it gets dejected and does
not cast a vote, as indicated by the gray cross. (b) The two bottom nodes are
completely surrounded by other orange nodes. Based on this local information,
they erroneously conclude that the upcoming election is a safe win, become
complacent, and do not vote. (c) Because three orange nodes do not vote, pur-
ple wins the overall election by 4 − 2.
2.2 The model
Our simplified model of voter behavior is intended to spotlight the role of two
social effects: complacency and dejectedness. In the model, voters have fixed opin-
ions and only need to decide whether to participate in an upcoming election.
Whether a node chooses to vote or abstain depends on whether its local neigh-
borhood causes the node to experience complacency, dejectedness, or neither of
these effects. Complacency is the effect where nodes that are surrounded pre-
dominantly by nodes with matching opinions do not bother to vote, because
they are convinced that their preferred candidate is going to win in any case.
Dejectedness is the effect where nodes that are surrounded predominantly by
nodes with opposite opinions tend not to vote, because they are convinced that
9
the situation is hopeless and their preferred candidate is going to lose.
Our model of voter behavior under dejectedness and complacency can be
introduced formally as follows. We assume that N voters live on a network, and
each node has some opinion θ, drawn from a probability distribution f (θ). We
shall assume that only two opinions exist, although studying the more general
case of multiple opinions is a natural direction for future work. In the context of
the model, these opinions can be thought of as preferences for one of two can-
didates in an election, but they could also represent binary referendum options,
or any other binary choice.
Continuing in the spirit of simplicity, we further assume that each node
knows the opinion of all its neighbors. The only question is who will vote.
Whether a node decides to vote or not depends on whether it thinks its vote will
make a difference, which in turn depends on the prevalence of the two opinions
among its neighbors in the network. We assume the following simple-minded
decision rule: A node chooses to cast its ballot if and only if its vote would cause
a tie or a one-vote win in its one-hop network neighborhood (assuming that all
its neighbors choose to vote). More precisely, if a focal node with opinion θ has
kθ neighbors with opinion θ and kϕ neighbors with the opposite opinion ϕ, it will
vote if and only if
0 ≤ kϕ − kθ ≤ 1.
Figure 2.1 illustrates the model. In the example shown, nine nodes live on
a ring graph. Five nodes hold a majority opinion (orange) and four nodes hold
a minority opinion (purple). Figure 2.1(a) illustrates the effect of dejectedness.
When the top orange node decides whether to cast its vote, it sees that both of
its neighbors hold the opposite opinion. Since purple outnumbers orange in the
10
top node’s neighborhood, even if the orange node decides to vote it cannot tie or
win the election locally, so it gets dejected and abstains from voting (as indicated
by the gray cross). Figure 2.1(b) illustrates the effect of complacency. The two
orange nodes at the bottom are completely surrounded by nodes with the same
orange opinion. These two nodes conclude that orange is a local majority, even
without their votes, and thus abstain from voting. Figure 2.1(c) shows the result
of the election: 3 orange nodes abstain from voting, leading to a 4−2 win by the
purple minority.
As this example shows, the election outcome depends on a surprisingly sub-
tle interplay among three factors: the network structure, the proportion of nodes
that hold each of the two opinions, and how the opinions are arranged among
the network nodes. Thus, this result raises several questions: Are some network
structures more likely to result in minority wins than others, at least under our
model? Does homophily (the tendency for neighboring nodes to hold identi-
cal opinions) increase or decrease the likelihood of minority wins? And how
does the minority size affect the likelihood of a minority win? In the remainder
of this chapter, we pursue these questions by simulating our model on various
network topologies, and with different choices for the arrangement of opinions
among the network nodes.
2.3 Model Networks
Our model networks (“electorates”) consist of N nodes (“voters”), of which N+
hold the majority opinion and N− hold the minority opinion. We typically work
with networks of size N = 100, in which case N− can also be interpreted as the
11
minority fraction, defined as the percentage of the electorate that holds the mi-
nority opinion. For each class of networks, we treat N− as a control parameter
and explore how the probability of a minority victory depends on N−. In our
analytical work on stochastic block networks (Section 2.3.2), we also find it con-
venient to express the results as a function of the ratio
N−
α = ≤ 1,
N+
a parameter that quantifies how closely divided the electorate is.
2.3.1 Erdős-Rényi networks
We begin by applying our model to Erdős-Rényi random graphs (Newman,
2018). In these networks, any given pair of nodes is connected by an undirected
edge with probability p. Since the number of nodes, N, and the edge probability,
p, define this family of random graphs, the family is often denoted G(N, p).
Figure 2.2 shows how the average proportion of unrepresentative outcomes
changes as we vary the edge probability p, for fixed network size N = 100 and
three different choices for the minority fraction N−. In Fig. 2.2(a), the majority
nodes outnumber the minority nodes by 80 to 20, a considerable margin. Under
these circumstances it is not easy for the minority to pull off an upset win, but
it is possible, thanks to the complacency of the majority. The probability that
minority wins peaks at around p = 0.25, with a corresponding win probability of
less than 0.2. Figure 2.2(b) shows the corresponding plot when we increase the
fraction of minority nodes to 30 out of 100, and Fig. 2.2(c) does the same for 40
minority nodes. The effects of these changes are mild. The main things to notice
are that as the electorate becomes more nearly evenly split, the peak probability
12
Figure 2.2: Unrepresentative outcomes are rare on Erdős-Rényi random graphs.
The plots show the proportion of minority victories on random graphs drawn
from the family G(100, p) as a function of the edge probability p (x-axis) for three
different scenarios: (a) a minority that is 20% of the entire electorate, (b) 30% mi-
nority, and (c) 40% minority. Each data point in a plot is based on 106 numerical
experiments. For all three scenarios, the peak probability of a minority victory
occurs for an intermediate p. But note that the minority never wins more than
half of the time; the curves lie below the dashed red line at all values of the edge
probability p.
that the minority wins becomes slightly higher and there is a widening of the
range of p-values where minority wins occasionally take place. Still, the main
message of Fig. 2.2 is that unrepresentative outcomes are fairly rare on this class
of random graphs. Indeed, in our simulations of the model on Erdős-Rényi
networks, there is no parameter regime where a minority wins most of the time.
From Figure 2.2, we can make two observations about when a minority can
win: minority victories become more likely for larger minorities and for inter-
mediate values of p. The first observation makes intuitive sense: A minority vic-
tory is less likely when the margin between the number of majority nodes and
the number of minority nodes is wider, because fewer minority nodes means
13
that more majority nodes must abstain from voting in order to ensure a minor-
ity win. Second, to understand why minority wins are most likely for interme-
diate values of p, it is helpful to consider the extreme network structures that
can arise in Erdős-Rényi networks. There are two such extremes. When p = 0,
the network has N components, each consisting of a single node, and no node
has neighbors. In the absence of local information, every node votes, making
unrepresentative outcomes impossible. At the other extreme, when p = 1 the
Erdős-Rényi network becomes a complete graph. On a complete graph, every
node has perfect information about the global state of the network, which leads
to dejectedness for the minority nodes and complacency for the majority nodes
(if the margin is greater than 1). As long as this condition holds true, nobody
votes, and therefore unrepresentative outcomes do not occur in this case either.
2.3.2 Stochastic block networks
In section 2.3.1, we assumed that opinions were distributed uniformly at ran-
dom among the nodes in Erdős–Rényi networks. Distributing the opinions in
this way meant that there was no homophily in the networks. Looking back at
Fig. 2.1, we see that the nodes that are resistant to complacency and dejected-
ness have the majority and minority opinions nearly equally represented in their
local neighborhoods. As such, it is the other nodes, the ones in homophilous
neighborhoods, that tend not to vote and thereby open the door to unrepresen-
tative outcomes. In other words, we expect homophily to play an important role
in enabling the minority to win.
One way to introduce such homophily into randomly generated networks is
14
to create random networks with community structure and assume that nodes
in the same community have the same opinion. We now do exactly this by
simulating our model on “stochastic block networks”(Newman, 2018).
It is helpful to think of stochastic block networks as a generalization of
Erdős-Rényi networks. Whereas in Erdős-Rényi networks, the probability of
forming an edge is the same for any two pairs of nodes, in stochastic block net-
works the node set is partitioned into disjoint subsets. The probability of form-
ing an edge between nodes then depends on the nodes’ respective subsets. If
the nodes are in the same subset, they are part of the same community, and the
probability of them being joined by an edge is high. On the other hand, nodes
in different subsets are assumed to not be part of the same community, and the
probability of an edge between them is low.
Since we are interested in the interactions between majority and minority
nodes, we will use a stochastic block network with two blocks. The probability
of forming an edge between two nodes can be represented as a matrix:
P = p11 p

12 , (2.1)
p p 21 22
where pi j is the probability of forming an edge for any pair of nodes from block
i and block j. For the sake of simplicity, we pick one in-block probability (p11 =
p22 = pin) and one inter-block probability (p12 = p21 = pout) to reflect the in-
group/out-group differences. These relative probabilities serve as a homophily
parameter. When pin/pout is high, the network exhibits high homophily, since
nodes are more likely to form edges within their block. When pin/pout is low, the
network exhibits low, or even anti-homophily, since the nodes are more likely to
form edges across blocks. In the special case pin = pout, we obtain Erdős-Rényi
15
Figure 2.3: Introducing community structure to random graphs allows for a
prevalence of unrepresentative outcomes within some parameter regions (the
diagonal yellow regions). The proportion of unrepresentative outcomes is
shown in color as a function of pin and pout on stochastic block networks of
N = 100 nodes with minority blocks of sizes (a) N− = 20, (b) 30, and (c) 40 nodes,
for 106 simulations.
networks.
Numerical experiments on stochastic block networks
Figure 2.3 shows the proportion of unrepresentative outcomes as a function of
pin and pout for networks where majority nodes outnumber minority nodes by
varying amounts. The color represents the proportion of simulations in which
the minority wins. Parameter values leading to unrepresentative outcomes are
conspicuous as the bright yellow regions.
While there are quantitative differences among the three networks, there are
important qualitative similarities. In each of the three panels, most of the pa-
rameter space is colored dark blue, corresponding to the democratic outcomes
16
one would naturally expect. However, there are also yellow diagonal regions
in which the minority wins more than half of the time. The highest probability
of a minority victory occurs close to the midline of the yellow region, where
pout/pin ≈ α. While not visible in the figure, the global maximum occurs on the
right edge of each panel, at the point where when pin = 1 and pout = α. As we
increase the size N− of the minority population, the location of the peak moves
up the pout axis, resulting in an increased slope of the yellow region, while pin
stays pinned at its maximum value, pin = 1.
These results confirm our intuition from the Erdős–Rényi networks: Unrep-
resentative outcomes occur in the intermediate information regime. They do
not thrive on complete networks, nor on fragmented ones with many compo-
nents. Rather, they favor regimes where nodes have an intermediate level of
knowledge about the state of the electorate as a whole.
An intuitive way of understanding Figure 2.3 is to think about the effects of
complacency and dejectedness. In order to avoid these effects, it is necessary to
have both majority and minority opinion nearly equally represented in a node’s
neighborhood. Because there are more majority nodes in the network, at high
pin and intermediate pout settings the minority nodes are most likely to know
almost equal numbers of nodes who agree and disagree with them. However,
in that same setting, the majority nodes are more likely to know more nodes
who agree with them because of the high pin probability, and therefore are more
likely to get complacent. This effect is what allows the minority to win.
17
Analytical results for stochastic block networks with pin = 1: Exact probabil-
ity of a minority victory
For the convenient special case where pin = 1, we can find the probability of a
minority victory exactly. To do so, observe that if the number of majority nodes
exceeds the number of minority nodes by at least two (N+ ≥ N− + 2), then none
of the majority nodes will vote, due to the effects of complacency. Therefore, in
this particular case, the minority will win as long as any minority node votes.
We can compute the probability of that event in a few easy steps as follows.
The first step is to consider the probability that any given minority node
votes. Because pin = 1, the given minority node is certain to be linked to all
the other minority nodes in the electorate and hence is sure to see exactly N−
votes for the minority opinion in its local neighborhood (including its own vote).
Now invoke the decision rule: the given minority node votes if and only if do-
ing so would either cause a tie or a one-vote victory in its local neighborhood.
For those events to happen, the minority node also needs to be connected to
either the same number, N−, of majority nodes, or one less than that number.
Those two events both happen according to binomial probability distributions,
because they involve choosing either N− or N− − 1 majority nodes out of a total
of N+ available. Therefore, the probability that the given minority node votes is
a sum of two binomial terms:
P(a(ny )given minority node votes)
N+ N
= ( p −ou)t (1 − p N+−N−N out) (2.2)−
N+
+ − p
N−−1
out (1 − p N+−(N−−1)N 1 out) .−
The first term expresses the probability that a minority node sees an equal num-
ber of majority and minority nodes (and will vote because it can cause a local
18
tie). The second term represents the probability that the minority node sees
N− − 1 majority nodes (and will vote because it can cause a local minority vic-
tory). All other possibilities are irrelevant: If the minority node sees more than
N− majority nodes, it would become dejected, whereas if it sees fewer than N−−1
majority nodes, it would become complacent.
The next step is to subtract the right hand side of (2.2) from unity, to get the
probability that a given minority node does not vote. Since there are N− such
nodes, and their decisions to vote are all independent, the probability that all of
them do not vote is:
[ P(no minority nodes vo] te) = (2.3)N
1 − −P(any given minority node votes) .
Then, by subtracting this quantity from 1, we obtain the probability that at least
one minority node votes,
P(at least one minority node votes) =
(2.4)
1 − P(no minority nodes vote).
As stated above, this probability is also equal to the probability that the minority
wins. Combining the equations above and replacing N− with αN+ throughout,
we finally arrive at our desired result:
P(mino[ rity(wins))N+
= (1 − 1 − ) pαN+out (1 − pout)N+−αN+αN (2.5)+
− N+ ]− pαN+−1out (1 − pout)N+−(αN+−1) αN+ .αN+ 1
Figure 2.4 shows an excellent match between this analytical prediction and sim-
ulations.
19
Figure 2.4: The proportion of unrepresentative outcomes on stochastic block
networks for the special case pin = 1.. The probability of a minority victory
is plotted as a function of pout, for networks of size N = 10. Results for three
values of N− are shown, corresponding to minority fractions of 20% (red), 30%
(orange), and 40% (yellow). The dotted black lines show the analytical expres-
sion in Eq. (2.5), which agrees with numerical results from 106 simulations (solid
colored lines).
Peak location and probability of a minority victory
In Figure 2.3 we saw that the probability of the minority winning in our simu-
lations on stochastic block networks was reached at high pin and intermediate
values of pout. Continuing to assume fully connected blocks, pin = 1, we can now
calculate at the value of pout that maximizes the probability of a minority victory.
To do so, we differentiate Eq. (2.5) with respect to pout and set the resulting ex-
pression to zero.
After straightforward but extensive algebra, and with the help of Stirling’s
formula, we find that in the limit N+ → ∞ with α held fixed,
pout = α
20
Figure 2.5: The solid curves show the proportion of unrepresentative outcomes
as described by Eq. (2.5) for constant α = 2/3 and N = 5, 10, 20, 100. The
stars indicate the location of the maximum on each curve.The probability pout
that maximizes the proportion of unrepresentative outcomes approaches α as
N increases. The limiting location of the peak, pout = α, is marked by the gray
vertical line.
maximizes the probability of a minority victory.
Figure 2.5 shows how the peak value of pout converges to α as N increases. In
these plots, we fix α = N−/N+ = 2/3 and vary the network size N. Notice that at
the peak, the proportion of unrepresentative outcomes approaches 1 as N goes
to infinity. With further effort, one can show that the peak value of a minority
victory deviates from 1 by an exponentially small term for N ≫ 1:
 √ 2αN 
P(minority wins | pout = α) ∼1 − exp − ( [ π(1 −]α)2)  (2.6)
× 1exp −
(1 − .α)π
Furthermore, the curves in Fig. 2.5 become increasingly sharply peaked as N
increases. To check this, we evaluate Eq. (2.5) in the same way at pout = α + ϵ for
21
ϵ ≪ 1 and find that P(minority wins | pout = α+ ϵ) tends to 0 as N approaches in-
finity. Therefore in the large-N limit, P(minority wins) tends to a discontinuous
function that equals 1 at pout = α and 0 everywhere else.
2.3.3 Networks with a heavy-tailed degree distribution
Erdős–Rényi networks and stochastic block networks are both widely studied.
Their simplicity allowed us to derive analytical results and gain some intuition
for when the minority could win the election in our model. In both models,
however, nodes tend to have very similar numbers of network neighbors. This
homogeneity is different from many real-world networks in which node degrees
can vary a lot (Barabási and Albert, 1999; Newman, 2018; Clauset et al., 2009).
As an example of networks with broad degree distributions we now consider
networks whose degree distributions follow a power law in the limit N → ∞.
Such scale-free networks have been claimed to capture features of many real-
world networks (Albert and Barabási, 2002; Barabási and Albert, 1999). Other
scholars have moderated or even argued against this claim (Clauset et al., 2009;
Broido and Clauset, 2019).
In our investigations of stochastic block networks, we found that the exis-
tence of community structure could in some cases increase the likelihood of
a minority win under our model. To understand the effect of homophily in
more detail, we also incorporate homophily in our simulations of our model
in networks with a heavy-tailed degree distribution. In order to introduce ho-
mophily into the setting of networks with power-law degree distributions, we
introduce a homophily parameter h (0 ≤ h ≤ 1). When h = 0, the node opinions
22
are distributed randomly on the network, whereas when h = 1, the majority
and minority nodes organize into disjoint blocks with no connection between
nodes of different opinions. Our algorithm for generating homophily on net-
works with power-law degree distributions is described in Appendix A. The
algorithm is heavily inspired by algorithms used to create configuration-model
networks (Newman, 2018). In that sense, our networks with power-law degree
distributions can be thought of as a class of configuration-model networks with
homophily.
Figure 2.6 shows examples of the resulting networks. As h increases, the
nodes get a higher preference for connecting to nodes with the same opinion.
When h = 0, the nodes’ local information is most likely to be representative of
the true proportion of opinions across the electorate as a whole. When h = 1, the
nodes’ local information will only reflect the presence of nodes with the same
opinion.
Figure 2.7 shows the proportion of unrepresentative outcomes on a network
with a heavy-tailed degree distribution and size N = 104 with minority fractions
20%, 30%, and 40%. We have chosen a larger network size to avoid undesired
topological correlations (Catanzaro et al., 2005; Boguná et al., 2004). The hori-
zontal axis shows the homophily parameter h. We observe once again that un-
representative outcomes occur most frequently when the homophily parameter
is in the intermediate range. In Figs. 2.7(a) and (b), for homophily parameter
values in range 0.45 ≤ h ≤ 0.95 the minority faction wins more than half of the
time. In Fig. 2.7 (c), the corresponding range is 0.55 ≤ h ≤ 0.9. Surprisingly,
increasing the minority size N− does not yield a larger peak probability of mi-
nority wins for these configuration networks, in contrast to the other network
23
Figure 2.6: Examples of networks with heavy-tailed degree distributions with
homophily factors (a) h = 0 , (b) h = 0.3, (c) h = 0.8, and (d) h = 1. All networks
are of size N = 15 with minority size N− = 5. We use the power law exponent
λ = 2.5 to generate the degree distribution.
structures tested in this chapter.
2.3.4 Geometric Random Networks
In Section 2.3.3, we considered networks with broad degree distributions, a trait
shared by some social networks. A qualitatively different class of networks are
those in which the likelihood of a link between two nodes depends on their
geographical separation. “Geometric random networks” provide some of the
24
Figure 2.7: Proportion of minority winning on networks with heavy-tailed de-
gree distributions and size N = 104, for (a) N− = 2000, (b) 3000, and (c) 4000
nodes as a function of the homophily factor h, for 103 simulations.
Figure 2.8: Example of majority (orange) and minority (purple) node distribu-
tions for geometric random networks with radius (a) r = 0.3, (b) r = 0.5, and (c)
r = 0.8.
25
simplest examples. To generate them, imagine throwing nodes uniformly at
random inside a unit square. We add an edge between any two nodes that lie
within a distance r of each other. A larger value of r results in denser networks,
as illustrated in Fig. 2.8.
In order to incorporate homophily into these sorts of random networks, we
assign minority and majority opinions preferentially to the left and right halves
of the unit square, respectively. With probability equal to the homophily param-
eter h, nodes lie within their preferred half of the square.
We vary the radius of connection r and compute the proportion of unrep-
resentative outcomes. Figure 2.9 shows the results of the simulation. While
the proportion of unrepresentative outcomes peaks in the intermediate radius
range, the peak probability of minority victories moves to the right as ho-
mophily increases. In a low-homophily setting, minority nodes benefit from
low radius to prevent dejectedness (they need to actively avoid knowing major-
ity nodes). In a high-homophily setting, minority nodes benefit from a higher
radius to prevent complacency (they need to ensure they know some majority
nodes). The extreme peak in panel (c) is interesting. It is due to the fact that
in extreme homophily settings, the majority half of the square domain is more
densely populated. Therefore, at low non-zero values of r, majority nodes be-
gin to see other majority nodes and become complacent before minority nodes
begin to see other nodes. This effect results in many disconnected minority
nodes voting. The effect is diminished when the difference between N+ and N−
is lower. While not shown here, our numerical experiments show that the peak
is higher for N− = 20% and lower for N− = 40%.
26
Figure 2.9: Proportion of minority victories on a geographic random network
as a function of the radius of connection r and probability of nodes lying within
their preferred half of the square. (a) h = 0.5 (no homophily), (b) h = 0.75
(moderate homophily), and (c) h = 1 (extreme homophily). The results are for
networks of size N = 100, with minority size N− = 30, for 104 simulations.
2.4 Discussion
In this chapter, we have presented a simple agent-based model of voter turnout.
By simulating the model on a variety of network structures, we found that it
is often possible for a minority faction to win the model election under the ef-
fects of dejectedness and complacency. These unrepresentative outcomes oc-
cur most frequently in the parameter ranges that correspond to intermediate
knowledge of the global state of the electorate, as well as in networks with some
homophily or community structure. We have further shown that unrepresen-
tative outcomes can become more likely in settings where the local distribution
of opinions is not representative of the average global distributions. Intermedi-
ate homophily settings often create regimes in which minority nodes are more
27
likely to overestimate the closeness of an election based on their one-hop net-
work neighborhood, while majority nodes are susceptible to complacency.
In reality, it remains unknown how much complacency and dejectedness
influence whether people cast their vote in elections. It is also unknown to
what extent such complacency and dejectedness would be caused by the im-
mediate social-network neighborhood of the voter; it seems quite possible, for
example, that media reports, forecasting agencies, and other non-local effects
could play even bigger roles in pushing voters to turn out or stay home. All
that one can say with certainty is that voter decision-making is a complex phe-
nomenon with many social, political, and structural factors influencing indi-
vidual choices. Nonetheless, our work suggests that homophily and network
structure can greatly affect vote outcomes in settings where voters choose to ab-
stain or cast their ballots based on the prevalence of opinions in their local social
neighborhood.
There are many extensions of this study that would be intriguing to try
in future work. Some directions could focus on implementing the model in
more general settings such as: Realizing the model on core/periphery networks,
adding more than two opinion states, modifying the decision rule, and perhaps
adding a tension between local and global information in the form of broadcast-
ers or forecasters. Another possibility would be to make the model dynamic.
What do nodes do after having lost an election that they thought was a safe
win? Introducing such dynamics and looking for fixed points, cycles, and other
time-varying states would be interesting.
28
CHAPTER 3
COMPARISON OF TWO ANALYTIC ENERGY BALANCE MODELS
Originally published in:
Landgren, E., Nadeau, A. Comparison of Two Analytic Energy Balance Models
Shows Stable Partial Ice Cover Possible for Any Obliquity. The Planetary Science
Journal 3.4 (2022): 79. https://dx.doi.org/10.3847/PSJ/ac603d
3.1 Introduction
The search for habitable exoplanets, perhaps hosting life, is one of the great en-
deavors of our time. To aid the search, it is important to understand what plan-
etary factors contribute to a planet’s habitability for life as we know it. In this
direction, analytical and computational studies using climate models to investi-
gate different planetary scenarios have been incredibly important for advancing
our scientific understanding of exoplanet climates. Recent work has investi-
gated how a planet’s orbital parameters such as obliquity (e.g. Armstrong et al.
(2014); Rose et al. (2017)), eccentricity (e.g. Kane et al. (2020)), or spin-orbit res-
onance (e.g. tidal locking Checlair et al. (2017, 2019a,b)) may affect the planet’s
climate. Other studies have looked at what role climatic elements such as sea ice
drift (e.g. Yue and Yang (2020)), land albedo (e.g. Rushby et al. (2020)), or ocean
heat transport (e.g. Checlair et al. (2019a)) may play in the long term habitabil-
ity of a planet. This study adds to the rapidly growing body of literature on
exoplanet climate by considering the role of albedo/temperature feedback in a
zonally averaged energy balance model with explicit obliquity dependence. We
compare the model output for two methods of modeling the energy transport
across latitudes and find qualitatively similar results between them.
29
This study is an extension of the work in Rose et al. (2017). In that study, Rose
et al. (2017) analyze an energy balance model with explicit obliquity dependence
and heat transport modeled with a diffusion term. They find an analytical so-
lution for temperature as a function of sine of the latitude and obliquity based
on the methods given in North (1975a). Their work shows that partial stable ice
cover is more likely for ice caps than for ice belts (the model scenario where ice
advances from the equatorial region of the planet rather than the poles, e.g.
as in Figure 3.1) and identifies regions of parameter space where the Snow-
ball catastrophe—where a small change in radiative forcing causes the planet
to quickly become completely ice covered—may be avoided.
Here we compare the model used in Rose et al. (2017) in two ways: (1) to the
same model but with a more accurate approximation to the insolation distribu-
tion and (2) to a similar energy balance model but with heat transport modeled
with relaxation to the global mean average temperature. Comparisons between
energy balance models with diffusive and relaxation terms for heat transport
have been conducted in the past for Earth (e.g. North (1984); Roe and Baker
(2010); Walsh (2017)); however, an analysis agnostic to the specificity of particu-
lar planets has yet to appear in the literature.
In Rose et al. (2017), the incoming stellar radiation is approximated by a
second degree polynomial in sine of latitude and cosine of the obliquity. A
second degree polynomial is sufficient to capture the qualitative behavior of
insolation distributions for planets with very low obliquities (between 0◦ and
45◦) and obliquities close to 90◦ (between 65◦ and 90◦), but does not capture
the qualitative distribution for planets between these ranges.1 To capture the
1Note that annual average insolation distribution for rapidly rotating planets is symmetric
about 90◦ obliquity. The distribution when obliquity is β is the same as when it is 180 − β. We
restrict our attention to obliquities between 0◦ and 90◦.
30
behavior accurately, one needs at least a sixth degree polynomial in sine of the
latitude and cosine of the obliquity (Nadeau and McGehee, 2017; Dobrovolskis,
2021). Here we show that the higher order approximation is a necessary one; we
would not find that stable partial ice cover is possible at any obliquity without
it. In particular, we find that stable partial ice cover is possible for any obliquity
for both diffusive and relaxation to the mean models.
Further, we find that the mode of heat transport does not change the quali-
tative distribution of the likelihood of planets with stable partial ice cover. Here
we show that stable partial ice cover is less likely for high-obliquity planets than
for low-obliquity planets for both diffusive and relaxation to the mean models,
agreeing with the results already shown in Rose et al. (2017). We also show that
low albedo contrast and low efficiency of heat transport favor stable partial ice
cover. The qualitative similarities between the models is encouraging because
although analytical solutions can be found in both cases, the analytical solutions
for the relaxation to the mean model are more tractable as they can be written
using elementary functions.
Following Rose et al. (2017), we consider the stability of equilibrium of the
system as radiative forcing changes in the model. Physically this could be
caused by atmospheric effects, such as changing greenhouse gases. In partic-
ular we consider how different parameters in the model affect the Snowball
catastrophe. We consider the severity of the Snowball catastrophe bifurcation
based on the distance between the bifurcation latitude and the Snowball state.
We identify the regions in parameter space where more severe bifurcations oc-
cur since a severe Snowball catastrophe may have stronger signals in a planet’s
climate record and may be more likely to be observed. We show that the Snow-
31
ball catastrophe is less severe in the relaxation to the mean model compared to
the diffusion model.
The chapter is laid out as follows. In Section 3.2.1, we present the governing
equations and a nondimensionalization of these equations to better quantify the
effects of parameter changes on the behavior of the system. In Section 3.2.2, we
derive the equations that relate the latitude of the saddle node bifurcation to the
corresponding parameter values. In Section 3.3, we calculate the relative likeli-
hood of stable partial ice cover. In Section 3.4, we classify the bifurcation into
the Snowball state based on its severity. A discussion of the results follows in
Section 3.5 and we conclude in Section 5.5. In Appendix B, we briefly analyt-
ically compare the equilibrium solutions of the diffusion and relaxation to the
mean models.
3.2 Governing Equations
3.2.1 Annual Average Energy Balance
We consider a one dimensional energy balance model of the form popularized
by Budyko (1969), Sellers (1969) and North (1975b). These models describe the
time evolution of temperature on a planet depending on the incoming and out-
going radiation and heat transfer across latitudes. We consider the model
∂T
R = Qs(y, β)(1 − α(y, η)) − (A + BT ) + F (T ) (3.1)
∂t
where T = T (y, t) is the annual zonally averaged temperature as a function
32
of sine of latitude y and time t and F (T ) will either be diffusion
F ∂ ∂(T ) = D∇2T = D (1 − y2) T (y, t)
∂y ∂y
or relaxation to the global mean tempera(ture ∫ 1 )
F (T ) = −C(T − T ) = −C T (y, t) − T (ϕ, t)dϕ .
0
As noted in Roe and Baker (2010), both forms of F (T ) are parameterizations of
the divergence of the poleward heat flux. Scientific discussions for the terms
in the relaxation to the mean model can be found in the literature, e.g. Held
and Suarez (1974); Checlair et al. (2017). Readers interested in a mathematical
discussion should see for example Tung (2007); Widiasih (2013); McGehee and
Widiasih (2014); Walsh (2017); Kaper and Engler (2013); and North (1975b). In
the above model, y is the sine of latitude. Hemispheric symmetry is assumed
here so that the latitude y ranges from 0 (the equator) to 1 (the north pole). The
ice line latitude (the boundary between the high and low albedo regions) is de-
noted by η. The mean annual amount of incoming solar radiation (insolation)
is represented by Q. The insolation distribution function s(y, β) depends on the
latitude and on the planetary axial tilt β. The co-albedo function (1 − α(y, η)) de-
termines the proportion of incoming solar radiation absorbed by the planetary
surface at each latitude. Outgoing radiation A + BT is in a linearized form. The
dimensional parameters and their values for the Earth can be found in many ref-
erences such as Tung (2007) or Kaper and Engler (2013). Nondimensionalized
parameter values for Earth are given in Table 3.1.
The albedo function α(y, η) is a piecewise constant function demarcating the
regions of the planet with high and low albedo. Let αp be the albedo polarward
33
of the ice line, and αe be the albedo equatorward of the ice line, then

(y ) 
αe 0 < y < η
α , η = 
αe+αp
2 y = η
αp η < y < 1
In the case of ice caps, the polar region has a higher albedo than the equatorial
region with αp > αe. For ice belts, the situation is reversed and αp < αe (Figure
3.1). In the following, we will denote high albedo with αh (ice covered regions)
and low albedo with αl (non-ice covered regions). In some energy balance mod-
els of Earth, the albedo values are set to αl = .32 and αh = .62 with αe = αl
and αp = αh, (e.g. in Widiasih (2013)). Ice-line-dependent albedo is used in
McGehee and Widiasih (2014); Widiasih (2013); Walsh (2017); Barry et al. (2017).
This is in contrast with the temperature-dependent albedo used in Rose et al.
(2017); North (1975b). Here we use the ice-dependent version because while the
two forms are equivalent in the diffusion model (Cahalan and North, 1979), the
temperature dependent version in the relaxation to the mean model can result
in spurious temperature solutions (Walsh and Rackauckas, 2015).
We model the annual average changes in the temperature profile by using
zonally averaged mean annual insolation in equation (3.1). The mean annual
insolation distribution s(y, β) depends only on obliquity β and latitude y. The
insolation distribution is given by the integral (e.g. Ward (1974); McGehee and
Widiasih (2014)) ∫ √
2 2π √
s(y, β) = 1 − ( 1 − y22 sin β cos γ − y cos β)2dγπ 0
where γ is the planet’s longitude and which is not amenable to analytical calcu-
lations in the models we consider. Instead, Nadeau and McGehee (2021) show
34
Figure 3.1: Planets at low obliquity (left) tend to exhibit ice caps, while planets
at high obliquity (right) tend to exhibit ice belts. The ice line latitude is marked
η.
that we may write s(y, β) as a Legendre series
∑∞
s(y, β) = a2n p2n(cos β)p2n(y)
n=0
where p2n is the Legendre polynomial of(degr)e(e 2n and)( )
(−1)n(4n n+ 1) ∑ 2n 2n + 2k 1/2
a2n = 22n−1 n − k 2k k + 1
k=0
using the standard notatio(n )
x x(x − 1) · · · (x − k + 1)
= .
k k!
Truncating the series for some fixed N gives a degree 2N polynomial approx-
imation σ2N to the true insolation distribution, where increasing N decreases
the root mean square error of the approximation (Nadeau and McGehee, 2021).
Nadeau and McGehee (2017, 2021) show that the sixth degree approximation
35
given by ∑3
σ6(y, β) = a2n p2n(cos β)p2n(y)
n=0
(3.2)
= 1 + a2 p2(cos β)p2(y) + a4 p4(cos β)p4(y)
+ a6 p6(cos β)p6(y),
where a2 = −5/8, a4 = −9/64, and a6 = −65/1024, and
p2(y) = (3y2 − 1)/2,
p4(y) = (35y4 − 30y2 + 3)/8,
p6(y) = (231y6 − 315y4 + 105y2 − 5)/16
is the smallest degree approximation which captures the qualitative shape of the
distribution for all obliquities and approximates the true distribution to within
1.6% for all obliquities. In particular, it is the smallest degree approximation
which captures the characteristic ‘W’ shape of insolation distributions for plan-
ets with obliquity between approximately 45◦ and 65◦ (Nadeau and McGehee,
2017). Dobrovolskis (2021) shows that extending the approximation to degree 8
or 10 improves the approximation significantly at the poles for obliquities close
to zero but only slightly at other latitudes or for higher obliquities. In Rose et al.
(2017), the degree two approximation is used (a4 = 0 and a6 = 0). Notice that
√
when cos β = 3/3 (i.e. when β ≈ 54.74◦), the degree two approximation is iden-
tically equal to 1 for all values of y. Here we use a generic truncation σ2N in our
analysis. In Section 3.3 we present results for the cases where N = 1 and N = 3.
The position of the ice lines, denoted by η, depends on the mean annual tem-
perature of the ice line. Ice–albedo feedback is incorporated by the dynamic ice
line equation that is coupled with equation (3.1). We use the following dynamic
ice line equation, first formulated in Widiasih (2013) for ice caps on Earth as
dη
= ρ(T (η, t) − Tc). (3.3)dt
36
For ice belts, the right-hand side should be multiplied by negative one. The
physical boundaries at the pole and the equator are built into the model, i.e. η
cannot be greater than 1 or less than 0. A mathematical treatment of the result-
ing nonsmooth system using a projection rule and a Filippov framework can
be found in Barry et al. (2017), where the invariance of the physically possible
region is shown.
An intuitive interpretation of Equation (3.3) notes that the mean annual tem-
perature at the ice line is denoted by T (η, t). The critical temperature Tc is the
highest temperature at which multiyear ice can be present. If the ice line tem-
perature is above Tc, then the ice cover shrinks. If the temperature at the ice line
is below Tc, the ice cover grows. The response constant ρ controls the speed of
the ice line response to a change in temperature. We are interested in the equi-
librium position of the ice line, which is obtained when the ice line temperature
is exactly Tc.
We nondimensionalize the system using transformations analogous to those
in Rose et al. (2017), namely
2πt A + BT (y)
τ = ωt = , T ∗ = ,
tyear A + BTc
where ω = 2π/tyear rescales time to be dimensionless, such that τ = 1 corresponds
to one year. The nondimensionalized temperature T ∗ is proportional to temper-
ature and outgoing longwave radiation. At the ice line, T ∗ is always equal to 1,
i.e. T ∗(η) = 1.
The nondimensionalized parameters are summarized in Table 3.1. The pa-
rameter transformations (which are the same for both diffusion and relaxation
37
to the mean) are
Rω (1 − α
q l
)Q
γ = , = ,
B A + BTc
1 − α
α = 1 − h− , ζ = cos(β), and1 αl
ρ(A + BTc)
λ = .
Bω
These parameters have the following physical interpretations:
• γ: Seasonal heat capacity of the system relative to the outgoing radiation
over one year.
• q: A measure of radiative forcing balance. It is directly proportional to the
annual average incoming solar radiation and inversely proportional to the
outgoing radiation at critical temperature Tc.
• α: A measure of albedo contrast that changes the ice—albedo feedback.
The minimum value of α = 0 means that the high albedo regions, αh, and
the low albedo regions, αl, have the same albedo. The maximum value
of α = 1 is not necessarily the maximal albedo contrast, but is instead the
range of states where αh = 1 (the regions with high albedo are infinitely
reflective and absorb no energy) and αl ≠ 1 (the regions with low albedo
absorb some energy).
• λ: A measure of the speed of ice line response to the changes in tempera-
ture.
The nondimensional parameter for heat transport in the diffusion model is
D
δ =
B
38
and it is a measure of the efficiency of heat transport across latitudes (Rose et al.,
2017; Stone, 1978). The nondimensional parameter for heat transport in the re-
laxation to the mean model has the same form
C
µ =
B
and similar interpretation. Note that despite similar forms, µ does not directly
correspond to δ. The discrepancy is caused by the fact that the diffusion coeffi-
cient in the diffusion model is not a linear scaling of the horizontal heat transfer
coefficient C in relaxation to the mean model. North (1975a) comments that the
divergence operator and the relaxation to the mean operator have the same ef-
fect on degree two polynomials when δ = µ/6. When N = 3, the relationship
between these parameters is more complicated. We consider the relationship
between equilibrium solutions to the two models in more detail the Appendix.
Throughout this work, we use δ = µ/6 when directly comparing the behavior of
the two models for fixed values of heat transport parameter.
Rewriting the albedo function with the nondimensionalization for ice caps
yields 
∗ 

(y ) 
1 0 < y < η
α , η =  2−α 2 y = η
.
1 − α η < y < 1
For belts, the positions of 1 and 1− α are swapped with α∗(y, η) = 1 for y > η and
α∗(y, η) = 1 − α for y < η.
The nondimensionalized version of the temperature model (3.1) when heat
transport is modeled as relaxation to the annual average temperature is given
by
∂T ∗
γ = qσ ∗2N(y, ζ)α (y, η) − T ∗ − µ(T ∗ − T ∗) (3.4)
∂τ
39
Table 3.1: Nondimensional Parameters for (3.4) and (3.5).2
Parameter Definition Brief descrip- Value
tion for
Earth
γ RωB Seasonal heat 6.13
capacity
q (1−αl)QA BT Radiative 1.27+ c
forcing
α 1 − 1−αh1− Albedo con- 0.44αl
trast
ζ cos(β) Cosine of 0.92
obliquity
λ ρ(A+BTc)B Ice line re- variesω
sponse
N 2N is the 1
degree of the
insolation ap-
proximation
µ CB Relaxation 1.6
to mean effi-
ciency of heat
transport
δ DB Diffusion effi- 0.31
ciency of heat
transport
and when heat transport is modeled as diffusion it is given by
∂T ∗
γ = qσ2N(y, ζ)α∗(y, η) − T ∗ + δ∇2T ∗. (3.5)
∂τ
Note that T ∗ is a function of y, η and τ but those dependencies are being sup-
pressed in the above equations. In the relaxation to the mean model, T ∗ depends
on η and τ (see Section 3.2.2 for more details). The nondimensionalized ice-line
equation for ice caps is
∂η
= λ(T ∗(η, τ) − 1), (3.6)
∂τ
where T ∗(η, τ) is the temperature evaluated at the ice line y = η. Note that as
with the dimensional equation (3.3), the right-hand side must be multiplied by
40
−1 for ice belts. In the following we do not explicitly consider dynamics of the
ice line and so λ plays no role in the following analysis.
3.2.2 Analysis of Temperature Equilibria
We expect to observe planets that have existed for a long time, and therefore we
expect in the simplest case to find the corresponding planets at temperature—
ice equilibrium. We focus on the equilibria of the ice line η as they undergo
bifurcations in radiative forcing q and the associated hysteresis loops in our
non-smooth system. Previously, bifurcations in A and Q have been considered
for Earth’s range of parameter values (Widiasih (2013)) and bifurcations in q
have been studied in Rose et al. (2017). We follow the framework developed
in Rose et al. (2017) and study bifurcations in q by considering a illustrative
sample of combinations of the physical parameters. In Section3.3 we will also
conduct a parameter sweep in the obliquity, albedo contrast, and efficiency of
heat transport to determine likelihood of stable partial ice cover.
In order to compute the parameter values that correspond to the bifurca-
tions, we first need to find an expression for the equilibrium temperature pro-
file T ∗(y, τ). Below we compute the equilibrium temperature profile for the re-
laxation to the mean model. This derivation can be found in many places (e.g.
McGehee and Widiasih (2014)), so we only summarize it here for convenience.
Following these calculations, we summarize the method to compute the equilib-
rium for the diffusion model which is shown in detail in North (1975a) and also
summarized in Rose et al. (2017). We briefly compare these analytical results in
the Appendix B.
41
Relaxation to the Mean Equilibrium Temperature
∗
To find T ∗eq(y, η), the equilibrium temperature at each latitude, we set
∂T
t = 0 and∂
first find T ∗eq(η), the global average∫equilibrium temperature, by averaging over1
the latitudinal range 0 to 1. Since T ∗(y, η, τ)dy = T ∗(η, τ), the last term in (3.4)
0
cancels, and we can explicitly so∫lve for T ∗eq(η), namely1
T ∗eq(η) = qσ2N(y, ζ)(α∗(y, η))dy.
0
Note that Teq(η) depends on η since α∗(y, η) depends on η, while the dependence
on y is integrated out.
Since α∗(y, η) differs for ice belts and ice caps, T ∗eq(η) also differs. Both solu-
tions are polynomials depending on obliquity and ice edge latitude. For ease of
notation, we let ∫ η ∑N
Σ2N(η) = σ2N(y, ζ)dy = a2n p2n(ζ)P2n(η)
∫ 0 n=0
where Pi(y) = pi(y)dy is the antiderivative of the Legendre polynomial pi(y).
Then the global average equilibrium temperature is
 ( )q((1 − α) + αΣ)2N(η) , ice capsT ∗eq(η) =  . (3.7)q 1 − αΣ2N(η) , ice belts
Note that T ∗eq(η) is proportional to the nondimensionalized radiative forcing q.
For fixed η, the more radiative forcing the planet receives, the warmer its mean
equilibrium temperature.
Once T ∗eq(η) has been found, we may solve for T ∗eq(y, η) in Equation (3.4) with
∂T ∗
t = 0. Indeed, for y ≠ η∂
∗ 1 ( )Teq(y, η) = qσ2N(y, ζ)α∗(y, η) + µT ∗eq(η)1 + µ
42
Due to the discontinuity in α∗(y, η), the temperature profile is discontinuous at
the ice line. We define the value at the ice line to be the average of the left and
right limits of T ∗eq as y approaches η, namely
lim T ∗ ∗∗ y→η+ eq(y, η) + limy→η− Teq(y, η)Teq(η, η) = .2
For ice caps the left and right hand limits are
∗ qσ2N(η, ζ) + µT
∗
eq(η)
lim Teq(y, η) = , (3.8)y→η− 1 + µ
∗
lim T ∗
qσ2N(η, ζ)(1 − α) + µTeq(η)
eq(y, η) = (3.9)y→η+ 1 + µ
and for ice belts limy→η− and limy→η+ are swapped. In both cases, the equilibrium
temperature at the ice line is given by
∗ qσ2N(η, ζ)(2 − α) + 2µT
∗
eq(η)
Teq(η, η) = ;2(1 + µ)
however, the function for T ∗eq(η) differs for ice caps and ice belts (see Equation
(3.7)). Note that this definition of the temperature profile at the ice line is equiv-
alent to using our definition of α∗(η, η).
Ice line equilibria occur when T ∗eq(η, η) = 1. We consider the response of
equilibria to changes in radiative forcing q. The relaxation to the mean model
allows us to solve exactly for the unique value of radiative forcing q as a function
of the other parameters, namely
2(1 + µ)
qη(ζ, α, µ) = − , (3.10)σ2N(η, ζ)(2 α) + 2µTx(η, ζ)
where 
T 
(1 − α) + αΣ2N(η, ζ) ice caps
x(η, ζ) =  .1 − αΣ2N(η, ζ) ice belts
43
Relaxation to the Mean Diffusion (Rose 2017)
1.0 1.0
qfree q ●free SN bifurcation
0.8 0.8
0.6 ● SN bifurcation 0.6 ● SN bifurcation
0.4 0.4
0.2 0.2
qsnow qsnow
0.0 0.0
1.0 1.2 1.4 1.6 1.8 1.0 1.2 1.4 1.6 1.8
q q
Figure 3.2: Plots showing the bifurcation diagrams demonstrating qfree, qsnow,
and the saddle node point ∂qη = 0 for Earth’s parameter values from Table 3.1
∂η
using the relaxation to the mean model used in this work (left) and the diffusion
model used in Rose et al. (2017) (right). Solid curves indicate stable ice line
locations. Dashed curves indicate unstable ice line locations.
Note that Tx(η, ζ) = T ∗eq(η)/q.
For ice caps, a planet is ice free when η = 1, and in a Snowball state when
η = 0. Following the conventions from Rose et al. (2017), for ice caps, we let qfree
denote the lowest value of q for which q1 is stable. Stability of the ice free state
is inferred from where qη intersects the line η = 1. When q > q1, the nondimen-
sional temperature T ∗eq at the pole is greater than 1. Similarly we let qsnow denote
the location where qη intersects the line η = 0. We can think of this as the highest
value of q for which q0 is stable. Note that for ice belts, q1 is qsnow, and q0 is qfree.
See Figure 3.2 for a bifurcation diagram demonstrating qfree and qsnow. It should
be noted that the definitions of qfree and qsnow used in this work take the average
of the two sides of the discontinuity at the ice line. Alternative definitions can
be used. We elaborate on the implications of this choice in Section 3.5.
44
η
η
Diffusion Equilibrium Temperature
For the remainder, we will use the tilde symbol (˜) over a function or variable to
denote that it is associated with the diffusion version of the model.
North (1975a) analytically solves the diffusive energy balance equation (3.5)
with ∂T/∂t = 0 for Earth with a degree two approximation to the insolation
distribution. He finds that solutions to the diffusive equation can be expressed
in terms of hypergeometric functions. This analytical solution is generalized to
the exoplanet case with degree two approximation to the insolation distribution
in Rose et al. (2017).
For an arbitrary approximation to the insolation approximation, one may
use the same methods described in North (1975a) and Rose et al. (2017). For
ease of notation, let
∑N a2n p2n(ζ)H2N(x) = p2n(x).1 + 2n(2n + 1)δ
n=0
As noted in North (1975a), qH2N(x) is the particular solution to the diffusion
equation on the interval [0, η) and q(1 − α)H2N(x) is the particular solution to
the diffusion equation on the interval (η, 1]. To find the general solution, one
must find the solution to the homogeneous equation. This derivation is given in
North (1975a), so we do not give it here. Instead we report the results of those
computations which yield the equilibrium temperature at the ice line
q (H2N(η) − αF2N(η)) , ice capsT̃ ∗eq(η) = q ((1 − α)H2N(η) + αF2N(η)) , ice belts
45
where
P′ν(η)H
′
2N(η) − H2N(η)Pν(η)F2N(η) = f (η),P′ νν(η) fν(η) − f ′ν (η)Pν(η)
fν(x) = 2F1(−ν/2, (1 + ν)/2, 1/2, x2),
Pν(x) = 2F1((1 + ν)/2,−ν/2, 1, 1 − x2),
1/2
ν = −1 (1 − 4/δ)+ ,
2 2
and 2F1 is the hypergeometric function. Although Pν(x) and fν(x) are complex-
valued for δ < 4, taking the real parts of these functions yields linearly indepen-
dent solutions as noted in North (1975a). Solving for q gives
− 
H2N(η) − αF2N(η), ice caps
q̃η(ζ, α, δ) 1 = (1 − α)H2N(η) + αF2N(η), ice belts
When N = 1, the results from Rose et al. (2017) are recovered.
3.3 Likelihood of Stable Partial Ice Cover
3.3.1 Defining the Region of Integration
Planets with stable partial ice cover are potential candidates for a Snowball
catastrophe bifurcation. We focus on quantifying the likelihood of stable par-
tial ice cover depending on planetary obliquity. We follow Rose et al. (2017) and
compute an estimate of the likelihood of stable partial ice cover based on the
size of the region in parameter space that admits these types of solutions. Loca-
tions where qη or q̃η have a zero derivative demarcate regions of stable partial ice
cover because the same location in the transformed plot with η on the vertical
axis is a saddle node bifurcation point (Figure 3.2).
46
Taking the derivative of q with respect to η yields
∂qη −2(1 + µ)
= ×
∂η [(2(− α)σ2N(η, ζ) + 2µTx(η)]2 ) (3.11)
(2 − ∂α) σ2N(η, ζ) + 2µ(±ασ2N(η, ζ))
∂η
for caps and belts, where ±ασ2N(η, ζ) is the derivative of Tx(η, ζ) for caps and
belts, respectively. Since Tx(η) ≥ 0, this definition is well defined. Parameter
values that result in ∂qη = 0 give the location of the saddle node. Setting Equa-
∂η
tion (3.11) to zero and solving for α yields ∂
αcrit(ζ, µ, η) = 
2 ∂ησ2N (η,ζ)
∂ for caps,
∂ησ2N (η,ζ)−2µσ2N (η,ζ) (3.12)
2 ∂∂ησ2N (η,ζ)
∂ for belts
∂ησ2N (η,ζ)+2µσ2N (η,ζ)
which is the critical value of the albedo contrast at the saddle node bifurcation
latitude. Stable partial ice cover is possible whenever α < αcrit(ζ, µ, η). At α =
αcrit(ζ, µ, η), the saddle node bifurcation occurs. Note that since µ is theoretically
unbounded, the function αcrit can become arbitrarily small. In the next section,
we use αcrit to find the relative likelihood of stable partial ice cover.
Following the method in Rose et al. (2017) but with the general degree 2N
insolation approximation, we find that in the diffusion model.H′2N (η)F′2N ( ) for caps,ηα̃crit(ζ, µ, η) =  (3.13)H′2N (η)H′ ( )−F′ for belts.2N η 2N (η)
Recall that H2N and F2N both also depend on ζ and δ.
3.3.2 Relative Likelihood of Stable Partial Ice Cover
In this section, we integrate over the parameter region in order to estimate the
relative likelihood of stable partial ice cover using the same method as was de-
47
Diffusion, 2nd degree Diffusion, 6th degree
1.4 1.4
N = 1 N = 3
1.2 1.2
1.0 ● 1.0 ●
PDF0
0.8 0.8
PDF1
0.6 PDF2 0.6
0.4 0.4
0.2 0.2
0.0 0.0
0 20 40 60 80 0 20 40 60 80
RelaxationObtloiquthityeβM, deeagnre, e2snd degree RelaxatioOnbtloiquthityeβM, deeagnre,e6sth degree
1.4 1.4
N = 1 N = 3
1.2 1.2
1.0 ● 1.0 ●
0.8 0.8
PDF0
0.6 PDF1 0.6
0.4 PDF2 0.4
0.2 PDF3 0.2
PDF4
0.0 0.0
0 20 40 β 60 80 0 20 40 β 60 80Obliquity , degrees Obliquity , degrees
Figure 3.3: Plots showing the reFliagtuirvee4l:ikliekleilhihoooodd of stable partial ice cover com-
pared to partial stable ice cover on Earth. Plot show how the likelihood changes
for different degree insolation approximations for Top row: the diffusion model
and Bottom row: the relaxation to the mean model. The top left plot is quali-
tatively similar to the likelihood plot in Rose et al. (2017), but includes the like-
lihood contribution from inaccessible edges. Likelihoods are given for PDF0
(blue), PDF1 (orange), PDF2 (green) for both models and of PDF3 (red) and
PDF4 (purple) for the relaxation to the mean model. The Earth obliquity is
marked with a blue circle on each plot.
scribed in Rose et al. (2017). We summarize it here for convenience.
We assume that there exists a true d4istribution for each parameter and that
the parameters are independent of each other. However, since the true probabil-
ity distributions are not known, we consider a handful of candidate probability
distributions. We integrate the composite probability density function over the
region of the domain where stable edges are present to obtain the likelihood of
stable partial ice cover for a given value of obliquity ζ. We normalize our results
by the likelihood value for the obliquity of the Earth, ζ = cos(23.5◦). From the in-
48
Relative Likelihood Relative Likelihood
Relative Likelihood Relative Likelihood
dependence assumption, we can write the overall probability density function
hplanet as follows:


hq(q)hµ(µ)hα(α), relaxation to mean,
hplanet = hq(q)hδ(δ)hα(α), diffusion.
We follow Rose et al. (2017) in our choice of candidate probability functions
to test in order to facilitate the comparison with their results and include two
additional probability functions to account for difference in δ and µ. Since q,
µ, and δ are nonnegative and unbounded, log-normal distributions are used to
incorporate the possibility of a long tail. The general form of the probability
density function of log-normal distribution is
 ( ( ))− ln x−
2
θ
1 

h(x) = √ exp m2 ,
σ 2π(x − θ)  2σ 
where σ is referred to as the shape parameter, m is the scale parameter, and θ
is the location parameter. In our two new PDFs we use gamma distributions for
µ. As with the log-normal distribution, the gamma distribution is a maximum
entropy probability distribution and so minimizes the amount of prior informa-
tion included in the distribution. The use of the gamma distribution allows us
to create probability density functions that make Earth’s value of µ = 1.6. While
the gamma distribution has similar properties to those of the log-normal distri-
bution, they are different enough to make them good candidates for sensitivity
analysis (Wiens, 1999; Iaci, 2000). The general form of the probability density
function for the gamma distribution is
baxa−1e−bx
h(x) =
Γ(a)
49
where a is the shape parameter, b is the rate parameter, and Γ is the gamma
function. Since α ∈ [0, 1], uniform and beta distributions are used. The beta
distribution has the probability density function
hα(α) = 6α(1 − α)
and favors values of α close to the value for Earth (α = 0.44) compared to ex-
tremes of α = 0 or 1. The probability density function parameters are chosen so
that Earth’s parameter values are not unlikely.
For PDF0, hα is uniform on [0, 1]; hµ and hδ are log-normal on [0,∞] with
shape parameter 1.0, scale parameter 1.0, and location parameter 0; hq is log-
normal on [0,∞] with shape parameter 0.5, scale parameter 1.0, and location
parameter 0.
PDF1 is the same as PDF0 except hµ and hδ are log-normal on [0,∞] with
shape parameter 2.0, scale parameter e, and location parameter 0. Compared to
PDF0 and PDF2, PDF1 makes very small values and large values of µ or δ more
likely.
PDF2 is the same as PDF0 except hα is parabolic beta-distribution on [0, 1],
with mode at 0.5. Compared to PDF0 and PDF1, PDF2 favors intermediate val-
ues of albedo contrast compared to extreme values.
PDF3 is the same as PDF0 except hµ and hδ have a gamma distribution on
[0,∞] with shape parameter 4.0, rate parameter 2.5. Compared to PDF0, PDF3
makes the µ value for Earth (1.6) more likely.
PDF4 is the same as PDF1 except hµ and hδ have a gamma distribution on
[0,∞] with shape parameter 4.0, rate parameter 1.85. Compared to PDF3, PDF4
makes larger values of µ more likely.
50
The (non-normalized) lik∫el∫ihoo∫d is given by
∫1 ∞ αcrit h (q , µ, α)dα dµ dηP (β) 0 0 0ice = ∞ ∫ ∞ ∫ planet η1 . (3.14)
hplanet(q, µ, α)dα dµ dq0 0 0
Note that the denominator in equation (3.14) is equal to 1 since hplanet is a prob-
ability density function. For higher values of µ and δ, the corresponding val-
ues of αcrit and αcrit become vanishingly small for all latitudes η. Therefore,
even though the region of integration is unbounded, numerical integration con-
verges. The results of the integration are summarized in Figure 3.3.
The qualitative results are similar across all models. Planets at low obliquity
have a higher likelihood of stable partial ice cover than planets at high obliquity,
and planets with moderate obliquity have the lowest likelihood of stable partial
ice cover. The apparent discontinuities in the relative likelihood plots for the
relaxation to the mean model in Figure 3.3 are due to cut-offs for the obliquities
where we use a model with caps or belts. As shown by Dobrovolskis (2021), the
minima of the insolation distribution occur between the poles and the equator
when the obliquity is between approximately 45◦ and 65.355◦. This behavior
can be thought of as a transitional regime between ice caps (where the insola-
tion minima occur at the poles) and ice belts (where the insolation minimum
occurs at the equator). While we do not specifically consider planets with vary-
ing obliquity in this work, we use both the caps and belts models for this range
of obliquities.
In agreement with the degree two diffusion model in Rose et al. (2017), the
likelihood of stable partial ice cover goes to zero for the degree 2 relaxation to
the mean model. It is easy to see the reason for this by noting that the derivative
of the insolation function appears in the numerator of αcrit and recalling that
√
the insolation function σ2 is constant for ζ = 3/3 which is when the obliquity
51
10.0030 10.0030 10.0030
↵ = 0.2 ↵ = 0.44 ↵ = 0.7
3.0025 3.0025 3.0025
1.0020 1.0020 1.0020
0.3015 0.3015 0.3015
0.1010 0.1010 0.1010
0.035 0.035 0.035
0.01 0.01 0.01
0 20 40 60 80 0 20 40 60 80 0 20 40 60 80
0 20 40 60 80 0 20 40 60 80 0 20 40 60 80
10.0030 10.0030 10.0030
↵ = 0.27 ↵ = 0.44
3.0025 3.0025 3.0025
1.0020 1.0020 1.0020
0.3015 0.3015 0.3015
0.1010 0.1010 0.1010
0.035 0.035 0.035
0.01 0.01 0.01
0 20 40 60 80 0 20 40 60 80 0 20 40 60 80
0 20 40 60 80 0 20 40 60 80 0 20 40 60 80
Obliquity  Obliquity  Obliquity 
Figure 1: Example tikz
Diffusion Caps 0.05,0.04 (Rose 2017) Diffusion Caps 0.44,0.04 (Rose 2017)
11.0 1.0
0.80.8 0.8
0.60.6 0.6
0.40.4 0.4
↵ = 0.05 ↵ = 0.44
 = 0.04  = 0.04
0.20.2 0.2
0. 0.0 0.0
0.96 0.98 1.00 1.02 1.04 1.06 1.0 1.2 1.4 1.6 1.8
0.96 1.00 1.04
q 1.0 1.4 1.8q
Diffusion Belt 0.05,0.04 (Rose 2017) Diffusion Belt 0.44,0.04 (Rose 2017)
11.0 1.0
0.80.8 0.8
0.60.6 0.6
0.40.4 0.4
↵ = 0.05 ↵ = 0.44
 = 0.04  = 0.04
0.20.2 0.2
0. 0.0 0.0
0.96 0.98 1.00 1.02 1.04 1.06 1.0 1.2 1.4 1.6 1.8
0.96 1.00 1.04
q 1.0 1.4 1.8q q q
Figure 3.4: Plots showing stability of the equilibrium ice line locations for β =
Figure 2: 5F4.4ig7◦ufroer t4hemdiaffunsuiosncrmipodtelvfeorrsiicoencap-sd(tio↵pursoiwo)nancdapicse b(etlot p(b)otatonmdrobwe)l.t (bottom)
The model takes N = 3, the minimum approximation needed to capture the
qualitative characteristics of the insolation distribution at this obliquity. Solid
curves indicate stable ice line locations. Dashed curves indicate unstable ice
line locations. Note that the horizontal scales are different between the left and
right columns.
1
52
heat transport eciency µ/6 heat transport eciency 
⌘ ⌘
η η
η η
is 54.74◦. For the models with higher degree insolation approximations, the
likelihood is never identically zero but nevertheless attains its minimum at mid-
obliquities for all tested probability distribution functions.
For example, in Figure 3.4 we plot bifurcation diagrams for the diffusion
model at β = 54.47◦. For N = 3 shown in the figure, small values of α and δ
admit stable partial ice edges for both ice caps and ice belts, although it is only
for the smallest values of these parameters that the edges are accessible in a hys-
teresis loop. Once α ≈ 0.3, there are no longer any stable partial ice edges in the
diffusion model. In contrast, no value of albedo contrast nor efficiency of heat
transport will cause the diagram for N = 1 to have a saddle node bifurcation.
3.4 Partial Ice Cover and the Snowball State
In addition to assessing the likelihood of partial ice cover, we quantify the sever-
ity of the Snowball state. Typically, one is interested in determining the inner
and outer edges of the habitable zone based on the system parameters; however,
as Rose et al. (2017) notes in Section 5.3 (and references therein), the models are
too simplistic to give good estimates for this range. It is possible to quantify
whether a bifurcation from no or partial ice cover to the Snowball state occurs
for the system parameters.
The most severe Snowball bifurcation occurs when qfree < qsnow (similarly
q̃free < q̃snow), and for all 0 < η < 1, qfree ≤ qη ≤ qsnow. Although there may be
stable ice line equilibria between 0 and 1, they would inaccessible by varying
q through a hysteresis loop in this situation. A less severe bifurcation occurs
when the ice line continuously transitions from the ice free state to small stable
53
10.0030 10.0030 10.0030
↵ = 0.2 ↵ = 0.44 ↵ = 0.7
3.0025 3.0025 3.0025
1.0020 1.0020 1.0020
0.3015 0.3015 0.3015
0.1010 0.1010 0.1010
0.035 0.035 0.035
0.01 0.01 0.01
0 20 40 60 80 0 20 40 60 80 0 20 40 60 80
0 20 40 60 80 0 20 40 60 80 0 20 40 60 80
10.0030 10.0030 10.0030
↵ = 0.27 ↵ = 0.44
3.0025 3.0025 3.0025
1.0020 1.0020 1.0020
0.3015 0.3015 0.3015
0.1010 0.1010 0.1010
0.035 0.035 0.035
0.01 0.01 0.01
0 20 40 60 80 0 20 40 60 80 0 20 40 60 80
0 20 40 60 80 0 20 40 60 80 0 20 40 60 80
Obliquity  Obliquity  Obliquity 
0.90 0.96 F1.0ig2ur1e.081:1.E14xa1m.2p0le1t.2ik6z1.32 1.38 1.42 1.46
Figure 3.5: Contour plots showing the minimum value of the nondimensional
incoming radiation qDiffuosiorn Cwapsh0i.0c5,h0.04t(hRoeseS20n17)owbDifafusliolnsCatpas t0.e44,i0.s04n(Roset2g017l)obally stable. Regions
11.0 1.0
above the dashed curves are where the Snowball catastrophe occurs directly
from the ice free s0t.a80.8te. First row: The d0i.8ffusion model with N = 3. Second row:
The relaxation to t0h.60.e6 mean model with 0N.6 = 3.
0.40.4 0.4
↵ = 0.05 ↵ = 0.44
 = 0.04  = 0.04
ice caps (or ice bel0t.2)0.2 and then the ice lin0.2e drops to the Snowball state. Behavior
0. 0.0 0.0
of solutions is the ty0.960.9p6 ic0.9a8 l p1.001.0ass1.a02 ge1.04thr1.o06 ugh1.0 a s1a.2 dd1l.4 1.6 1.80 1.04q 1.0 1.4e node1.b8q ifurcation (Strogatz,
Diffusion Belt 0.05,0.04 (Rose 2017) Diffusion Belt 0.44,0.04 (Rose 2017)
11.0 1.0
2018). The least severe scenario is where no saddle node bifurcation occurs
0.80.8 0.8
and stable ice line equilibrium depends continuously on the parameter q. This
0.60.6 0.6
occurs when qsnow0.<40.4 qfree, and for all 0 <0.4 η < 1, q↵ = 0.05 ↵ = 0.44 snow ≤ qη ≤ qfree.
 = 0.04  = 0.04
0.20.2 0.2
In Figure 3.5, w0.e0.0 plo 0.00.96 0.98t co1.00nto1.0u2 rs1.04for1.06the smallest value of the nondimensional1.0 1.2 1.4 1.6 1.8
0.96 1.00 1.04
q 1.0 1.4 1.8q q q
incoming radiation q for which the Snowball state is not globally attracting (i.e.
Figure 2: Figure 4 manuscript version - di↵usion caps (top) and belt (bottom)
the only stable solution). In Rose et al. (2017), this value was referred to as qhab
and computed by
1
qhab = min{qfree, qstab}
where qstab is the minimum q value for which an ice line equilibrium exists be-
tween 0 and 1. Frequently this corresponds to the saddle node bifurcation which
54
heat transport eciency µ/6 heat transport eciency 
⌘ ⌘
η η
η η
causes the Snowball catastrophe (if present, e.g. in Figure 3.4 only the lower
left plot exhibits this bifurcation). If the model has stable partial ice edges be-
tween 0 and 1 but no saddle node bifurcation causing the Snowball catastrophe,
qstab = qsnow. When qfree > qstab, the hysteresis loop generated by varying q causes
a drop from the ice free state to a partially ice covered state instead of directly
to the Snowball state.
In Figure 3.5, the regions above the dashed black curves are where the Snow-
ball catastrophe occurs directly from the ice free state. In an extension of the
work in Rose et al. (2017), here we consider the models with the degree six
approximation to the insolation distribution instead of the degree two approx-
imation for insolation in the diffusion model. For ease of comparison with the
results in Rose et al. (2017), we restrict the ice cap models to obliquities between
0◦ and 55◦ and the ice belt models to obliquities between 55◦ and 90◦.
It is straightforward to compare the top row of Figure 3.5 to Figure 8 in Rose
et al. (2017). Analytically, the only difference is the degree of the insolation dis-
tribution used. On the whole, the figures are very similar, having nearly the
same contours. They differ only slightly in the maximum and minimum values:
the sixth degree approximation decreases the extrema achieved for qhab for the
model with ice caps and increases the extrema for the model with an ice belt rel-
ative to the second degree approximation. This is due to the increased accuracy
of the insolation approximation at the poles (minima for ice caps, maxima for an
ice belt). The largest difference is that for the smallest obliquities in the degree
six diffusion model, the Snowball catastrophe from the ice free state occurs for
larger values of the heat transport parameter. This means the Snowball catas-
trophe is less severe for low obliquity planets than in the degree two diffusion
55
model. The transition from ice free directly to the Snowball state is minimally
affected for obliquities above approximately 30◦.
In the bottom row of Figure 3.5 we plot the contours for the degree six relax-
ation to the mean model. The relaxation to the mean model has stark contrasts
with the diffusion model in the row above it. For low to moderate albedo con-
trast the minima and maxima of both models are relatively similar. However
for large albedo contrast, the maximum values achieved for the relaxation to the
mean model are significantly greater. Analytically this is caused by the depen-
dence of qfree on the albedo contrast. In the diffusion model qfree is independent
of α, but in the relaxation to the mean model it is not. As α increases so does
qfree for the relaxation to the mean model (this can be seen in equation (3.10)).
In both models, increasing the albedo contrast increases the severity of the
Snowball catastrophe for all obliquities (the dashed black curves decrease), but
the increase in severity in the relaxation to the mean model is minimal. Al-
though for all obliquities in either model lower heat transport efficiency always
guards against the most severe transition (from ice free directly to Snowball).
In the diffusion model for the range of heat transport plotted this most severe
transition is always present for the mid-obliquities and increasingly present for
higher obliquities as the albedo contrast increases. To avoid the ice free directly
to Snowball transition requires decreasing the heat transport efficiency at least
another order of magnitude than what is plotted (i.e. to at least 0.001 for β = 90◦
and α = 0.7). Conversely, in the relaxation to the mean model there is only a
small range of obliquities where this most severe transition is present for the
range of heat transport considered.
56
3.5 Discussion
In this chapter we have analyzed a one dimensional energy balance model with
heat transport modeled by relaxation to the global mean temperature and dif-
fusion. The relaxation to the mean and diffusion versions of the Budyko-Sellers
model provide two similar ways to model the multitude of processes involved
in energy transfer between latitudes. Both methods ensure that energy is trans-
ported from latitudes that are “hot” to ones that are “cold.” The diffusive heat
transport is a local process that necessitates special treatment at the poles, while
relaxation to the mean global temperature is a global process that does not re-
quire special boundary conditions (Widiasih, 2013). This work may be inter-
preted for any rapidly rotating rocky planet with some physical mechanism of
heat transport. This work applies to planets where the temperature affects the
albedo, in particular, we assume that higher temperatures decrease the albedo
as they do for ice/water on Earth.
We have also considered the effects of different approximations to the an-
nual average insolation distribution on the results of the models. The second
degree approximation of the insolation distribution used in Rose et al. (2017)
does not capture the qualitative distribution of mid-obliquity planets. Planets
with obliquities between approximately 45◦ and 65◦ have a characteristic ‘W’
shape that requires a degree six (or higher) polynomial approximation (Nadeau
and McGehee, 2017).
A main result from Rose et al. (2017) is that the likelihood of stable partial ice
cover goes to zero at 55◦ obliquity. This result is due to the insolation approxi-
mation used in their study, which is constant for 55◦ obliquity. In the definition
57
of α̃crit (equation (3.13)), if H2N were constant then α̃crit = 0 for all values of the
arguments. When N = 1 as in Rose et al. (2017), the function H2N is constant
exactly when the insolation is constant, i.e. when the obliquity is 54.47◦. This
means that the integral in the numerator of the likelihood calculation (equation
(3.14)) is zero. We see a similar problem with the relaxation to the mean model
when N = 1. Here it is clear to see that when the insolation is constant then
αcrit = 0 (see Equation (3.12)). Taking a higher degree approximation in either
model, as we do here, avoids these problems.
In Figure 3.3 we saw that low obliquity planets are more likely to have stable
partial ice cover than those with high obliquity, and planets at middle obliquities
are least likely to have stable partial ice cover, which is qualitatively similar
to the likelihood computations in Rose et al. (2017). As noted above, we find
that stable partial ice cover is possible at all obliquities and that, in particular,
the relative likelihood of finding a planet with partial stable ice cover ranges
between 7% and 32%.
Comparing the relaxation to the mean model to the diffusion model, we note
that the latter predicts lower likelihood of stable partial ice cover at lower obliq-
uities. This effect is due to the fact that the diffusion model can exhibit a second
saddle node bifurcation at high values of η, close to the poles (called the small
ice cap instability) for Earth’s obliquity. As the obliquity decreases to zero the
small ice cap instability shrinks in size and can disappear completely for small
to moderate values of the albedo contrast α, increasing the range of η where sta-
ble ice edges are possible. In contrast, the relaxation to the mean model does not
exhibit the small ice cap instability at Earth’s obliquity for degree two or degree
six insolation approximation. This means that the likelihood remains relatively
58
flat as the obliquity decreases.
The relaxation to the mean model also exhibits pronounced differences be-
tween PDF0, PDF1 and PDF2 at high obliquities. For high values of the obliq-
uity β, the likelihood of stable partial ice cover is lower for PDF2. The difference
between PDF2 and other tested probability density functions is due to the dif-
ferences in αcrit between the relaxation to the mean model and the diffusion
model. Since PDF2 changes the distribution of α from a uniform distribution to
a parabolic beta distribution, the shape of αcrit results in a more pronounced dif-
ference for the relaxation to the mean model than for the diffusion model. The
resulting gap between the likelihood curves conveys decreased certainty about
the likelihood of stable partial ice cover on high obliquity planets. Note also that
for PDF3 and PDF4, the relaxation to the mean model exhibits behavior similar
to that for PDF1.
We quantify the effects of albedo contrast and efficiency of heat transport on
the presence of hysteresis loops in radiative forcing. We find that the severity of
the Snowball bifurcation increases as the albedo contrast α and the efficiency of
heat transport (δ or µ) increase.
The above behavior can be explained by the effect of albedo contrast and ef-
ficiency of heat transport on the planetary climate mechanism. When the albedo
contrast α is low, ice is not much more reflective than the non-frozen regions (ei-
ther ground or water), resulting in suppressed ice–albedo feedback. When δ or
µ are low, the near-absence of heat transport across latitudes limits the interac-
tion between the ice regions and the water regions of the planet, thus reducing
the likelihood of the Snowball catastrophe. When α is high, the reflectivity of ice
is much higher than that of the non-frozen regions, expediting the ice—albedo
59
processes. When µ is high, the heat transport across latitudes makes it difficult
to maintain a difference in temperatures between ice regions and water regions,
leading to an ice free or a Snowball planet.
In the range of obliquities where both ice caps and ice belts may be stable
it is not possible to transition continuously between the ice free and Snowball
states as q is varied, i.e. there will always be a hysteresis loop when varying q.
The hysteresis will either contain a saddle node bifurcation or the most severe
Snowball bifurcation from ice free to completely ice covered. A future study
might explore whether the lack of a region without hysteresis in the parameter
space is a contributing factor for the decrease in the likelihood of stable partial
ice cover for these obliquities in Figure 3.3.
It should be noted that dealing with the discontinuity at the ice line affects
the definitions of qfree and qsnow. For ice caps, physically, the definition qfree = q1
can be interpreted as vanishingly small ice caps, and qsnow = q0 as vanishingly
small equatorial strip of water. This choice is consistent with the literature on
the relaxation to the mean model (Widiasih, 2013; McGehee and Widiasih, 2014).
This definition ensures that the bifurcation diagram for the relaxation to the
mean model is continuous at η = 0 and η = 1. However, one could consider an
alternative definition of qfree and qsnow: instead of taking the average of the two
sides of the discontinuity at the ice line, one could use only the warm branch
(Equation (3.8)) of the temperature equilibrium to define qfree and only the cold
branch (Equation (3.9)) to define qsnow. This choice yields q
(1+µ)
free = σ2N (η,ζ) and+µ
q (1+µ)snow = (1− )( ( ) ) . Such a definition might be more physically intuitive, andα σ2N η,ζ +µ
removing the dependence of qfree on α might lead to a higher estimate for the
likelihood of ice-free solutions. However, analyzing the implications of these
60
alternatives would require extensive mathematical analysis of bistability of ice
free and partial ice states that we defer to future work.
The robustness of the Snowball catastrophe and the parameter regimes
where an energy balance model might be applicable has been debated. In the
GCM simulations conducted by Ferreira et al. (2014), the Snowball catastrophe
occurs only for particular ocean regimes. Wagner and Eisenman (2015) show
that meridional heat transfer may increase ice cover stability. Rose and Marshall
(2009) have extended the energy balance models to include ocean heat transport
and meridional structure and have found that Snowball catastrophe is possible
in those models.
If this work were to be applied to an observed planet, the obliquity ζ and
the albedo contrast α could perhaps be measured, and albedo signatures could
be compared to the predictions of our model. The parameter q will be more
difficult to estimate. While the mean annual insolation Q can be derived from
information about the star and orbit of the planet, the dependence of q on the
atmospheric parameters A and B make determining q challenging. The param-
eters δ and µ, efficiency of heat transport, would also be difficult to measure
directly because of our limited knowledge of rates of heat transport on different
planets.
3.6 Conclusion
In this chapter we have analyzed a one dimensional energy balance model with
two different methods of modeling the heat transport. The models have explicit
dependence on the planet’s obliquity through the annual average insolation dis-
61
tribution which we approximate with degree two and degree six polynomials.
We pay particular attention to the planet’s obliquity, radiative forcing, albedo
contrast, and efficiency of heat transport and find:
1. With an improved approximation to the insolation distribution function,
planets at all values of obliquity exhibit nonzero likelihood of stable par-
tial ice cover regardless of mode of heat transport. Minimum likelihood
ranges from 7% to 32% relative to Earth’s likelihood of stable partial ice
cover, depending on the model and probability density functions for the
parameters. Maximum likelihood ranges from 110% to 140% relative to
Earth.
2. Several results from the earlier study by Rose et al. (2017) are seen here in
both the diffusion and relaxation to the mean models:
(a) Models with high obliquity are less likely to have stable partial ice
cover than ones with low obliquity but are more likely than models
with moderate obliquity.
(b) Low albedo contrast and low efficiency of heat transport favor stable
partial ice cover.
(c) High albedo contrast and high efficiency of heat transport favor se-
vere Snowball catastrophe in a hysteresis loop caused by changes in
radiative forcing.
3. In the relaxation to the mean model, a larger range of heat transport ef-
ficiency supports partial stable ice cover at high obliquities even for high
albedo contrast. Whereas in the diffusion model stable partial ice cover
a high obliquities and high albedo contrast is possible only for very low
heat transport efficiency.
62
Both models discussed in this work are highly simplified, and so direct im-
plications for real physical systems are tenuous when taken individually. We
may interpret the results with higher confidence in regions of parameter space
where their solutions give similar results. Places where they differ indicate need
for further investigations with more complex models.
For example, due to the disagreement of the models shown in Figure 3.5, a
future study may explore the transition to the Snowball state for high obliquity
planets in a general circulation model. Other extensions of this work would be
to look at climactic histories of planets by introducing obliquity or eccentricity
variations in time into the energy balance model used in this study. Variations
in a planet’s orbital parameters change the behavior of planetary ice cover over
geological time and which are easiest to model over long time periods with
simple energy balance models like the ones presented here.
63
CHAPTER 4
SWAMPE: A SHALLOW-WATER ATMOSPHERIC MODEL IN PYTHON
FOR EXOPLANETS
Current efforts to model exoplanet atmospheres primarily focus on minimal-
complexity one-dimensional and high-complexity three-dimensional models.
In the previous chapter, we have explored one-dimensional energy-balance
models. While they are useful for characterizing rapidly rotating planets with
little longitudinal variation, such as terrestrial planets, they can fail to account
for important atmospheric processes (Feng et al., 2016). On the other hand,
three-dimensional models have the capacity to explore the entirely of the phys-
ical space (typically in latitude, longitude, and pressure). They are frequently
based on primitive equations (e.g. Menou and Rauscher (2009), Kataria et al.
(2016b), Parmentier et al. (2013)) or on the Navier-Stokes equations (e.g. Cooper
and Showman (2006), Dobbs-Dixon and Agol (2013)) and can be used to un-
derstand a variety of radiative, chemical, and dynamical processes. Three-
dimensional models such as ROCKE-3D Way et al. (2017) can be tuned to a
variety of exoplanets. However, three-dimensional models tend to be computa-
tionally expensive, sometimes taking months to explore the parameter space.
The difference in capability between one-dimensional and three-dimensional
models leaves a natural gap for two-dimensional shallow-water models, which
can capture the spatial variability as well as run fast enough to rapidly ex-
plore the parameter space and study the dynamical mechanisms. In particular,
shallow-water models have been used to study solar system planets (including
Earth), e.g. Ferrari and Ferreira (2011), Brueshaber et al. (2019). Outside of the
solar system, shallow-water models have been used to understand a variety of
64
Figure 4.1: Illustration of the sequence of steps performed at each time step
using the spectral method employed by SWAMPE based on Hack and Jakob
(1992).
atmospheric phenomena of hot Jupiters, such as atmospheric variability Menou
et al. (2003) and superrotation Showman and Polvani (2011). They have also
been used to make observational predictions for hot Jupiters (e.g. Langton and
Laughlin (2008), Perez-Becker and Showman (2013)). However, many of these
models are written in Fortran, which makes them difficult to adapt for the var-
ied needs of exoplanetary science.
SWAMPE offers a fully Python, open-source implementation of the 2D
shallow-water system, based on the approach outlined in Hack and Jakob
(1992). The code employs a spectral method, which takes advantage of the
convenient properties of spherical harmonics to speed up the computation. A
high-level overview of the spectral method is illustrated in Figure 4.1. In the
65
physical space, the variables are resolved on a latitude-longitude grid. These
values are then transformed into the wavenumber space (the space of coeffi-
cients for spherical harmonic functions). The time-stepping takes place in the
wavenumber space, but the variables are transformed into the physical space so
that nonlinear terms can be updated.
This package does not require multiple cores (in fact, SWAMPE simulations
can be run on a personal laptop), and is flexible and modular. SWAMPE is
designed to be easily modified to model dissimilar space objects, from Brown
Dwarfs to terrestrial, potentially habitable exoplanets. SWAMPE provides the
capability to conduct wide parameter sweeps and to produce maps of the ther-
mal and wind properties of the planets in latitude and longitude, which can
be used to help constrain and make predictions for observations of their atmo-
spheres.
Chapter 5 presents applications of SWAMPE to the atmospheres of syn-
chronously rotating sub-Neptunes. In particular, Section 5.2.1 provides the
governing equations implemented in the software package. In Appendix C, I
present an open source software repository for SWAMPE that makes this tool
available to the planetary science community and provide the time-stepping al-
gorithm used for SWAMPE.
66
CHAPTER 5
EXPLORING THE INTERACTION OF ROTATION RATE AND STELLAR
IRRADIATION ON SYNCHRONOUSLY ROTATING SUB-NEPTUNES
5.1 Introduction
Despite being some of the most common planet types in our galaxy (Deeg
and Belmonte, 2018), sub-Neptunes, extrasolar planets with sizes smaller than
Neptune, have no analogs in the solar system. They can be readily character-
ized with current facilities like the Hubble and James Webb Space Telescopes.
Rocky super-Earths and gaseous mini-Neptunes (planets 1.6-4 times the radius
of Earth), collectively referred to as sub-Neptunes, are prime candidates for
modeling that can reveal a variety of atmospheric phenomena. In order to an-
swer questions about potential habitability of exoplanets, it is important to de-
velop a robust understanding of a variety of dynamic processes that can take
place in exoplanetary atmospheres. These types of planets exist in a wide phase
space in terms of temperature (most of them falling in the 400-1200 K range
(Crossfield and Kreidberg, 2017)), chemical composition, and aerosol composi-
tion. Investigating the atmospheric circulation of planets across the entire pa-
rameter space and providing observational predictions can be time-consuming
with high-complexity three-dimensional (3D) general circulation models, espe-
cially those that try to treat the full range of chemical and physical processes
taking place in their atmospheres. On the other end of model complexity, in
many cases 1D models can be insufficient in capturing the inherently 3D pro-
cesses taking place in exoplanetary atmospheres (e.g. Feng et al. (2016), Kataria
et al. (2016a),Line and Parmentier (2016)).
67
For these reasons, shallow-water models are useful tools that can capture
key dynamical processes despite being more simplified than 3D models. Such
models represent the dynamics of a thin layer of a fluid of constant density and
variable thickness. The shallow-water equations are a coupled system of hor-
izontal momentum and mass conservation, which govern the evolution of the
zonal and meridional velocity components as well as the layer thickness. These
equations are typically solved numerically as a function of longitude, latitude,
and time. For a description of the shallow-water model focused on application
to exoplanetary atmospheric circulation, see, e.g. Showman et al. (2010).
There is a rich heritage of shallow-water models that have been used to study
solar system planets. For example, shallow-water 2D GCMs have been used
to study topography variation on Earth (?), (Kraucunas and Hartmann, 2007),
planetary-scale waves on Venus Iga and Matsuda (2005), latitudinal transport
and superrotation on Titan (Luz and Hourdin, 2003), and polar vortex dynamics
on Jupiter, Saturn, Uranus, and Neptune (Brueshaber et al., 2019) as well as
on Mars (Rostami and Zeitlin, 2017). Shallow-water models have been used
to model the atmospheric phenomena of objects beyond the solar system as
well. Perez-Becker and Showman (2013) explore the interaction of radiative
and drag timescales using a 2D model, Zhang and Showman (2014) explore
the circulation of Brown Dwarfs, and Penn and Vallis (2018) apply the shallow-
water framework to terrestrial planets, Hammond and Pierrehumbert (2018)
use it to consider tidally locked exoplanets, and Ohno and Zhang (2019) model
nonsynchronized, eccentric-tilted exoplanets.
The atmospheric properties of sub-Neptunes have been probed in several
studies using 3D GCMs. Zhang and Showman (2017) explore the effects of
68
bulk composition on the atmospheres of sub-Neptunes, and find that the effect
of molecular weight dominates in determining day-night temperature contrast
and the width of the superrotating jet. Wang and Wordsworth (2020); Christie
et al. (2022) study the atmospheric circulation of the sub-Neptune GJ 1214 b us-
ing 3D GCMs, and find that a long convergence time is needed to simulate long
radiative timescale, and that high metallicity and clouds with large vertical ex-
tents are needed to match observation. Dos Santos et al. (2022) estimate the
atmospheric escape rate and outflow temperature on the warm Neptune HAT-
P-11b by fitting retrieval models to the observed He transmission spectrum.
In this study, we use an open-source shallow water model, SWAMPE, to
investigate the atmospheric dynamics of sub-Neptunes. We contribute to the
growing literature on sub-Neptunes by systematically investigating the inter-
action of planetary rotation rate and radiative forcing through their effect on
planetary atmospheric dynamics. Our goal is to explore the atmospheric cir-
culation of sub-Neptunes and to identify the key differences in the dominant
atmospheric properties of this type of planets compared to hotter Jovian atmo-
spheres. While idealized, the shallow-water framework is well-suited for the
analysis of the dynamics of large-scale atmospheric flow, and we employ it to
explore the interaction of planetary rotation period, radiative timescale, and ir-
radiation. This simulation ensemble can guide further studies with more so-
phisticated models.
The chapter is organized as follows. Section 5.2.1 introduces our dynamical
model. In Section 5.2.3, we present the parameter regimes used for our simula-
tions. We present our simulation results in Section 5.3. We discuss the results in
Section 5.4.
69
5.2 Methods
5.2.1 Model
We study atmospheric circulation on sub-Neptunes using an idealized shallow-
water model we developed, SWAMPE (Shallow-Water Atmospheric Model in
Python for Exoplanets). The atmosphere is assumed to be a fluid that is incom-
pressible and hydrostatically balanced. For a derivation of the shallow water
equations from the continuity equation and the equation of motion, see, e.g.
Kaper and Engler (2013). As the name suggests, the shallow-water simplifi-
cation of the Navier-Stokes equations assumes that the fluid (in our case, the
planetary atmosphere) is spread in a thin layer (i.e. the horizontal scale is much
larger than the vertical scale). This idealized system can nevertheless capture
large-scale phenomena on a sphere in response to gravitational and rotational
accelerations. Within this framework, the wind patterns exhibited by the model
can be explained by considering the interaction of standing, planetary-scale
waves with the mean flow (e.g. Showman and Polvani, 2011; Showman et al.,
2013).
The idealized shallow-water system models the atmosphere as two layers:
an buoyant upper layer at constant density ρupper and an infinitely deep con-
vective layer at a higher constant density ρlower. Sub-Neptunes are likely to have
deep atmospheres (> 103 km, (see e.g. Kite et al., 2020)) atmospheres), and so the
shallow-water framework is appropriate for modeling these types of exoplan-
ets. We focus on modeling horizontal velocities in the upper layer, but vertical
velocities are captured in the model as mass transfer between layers.
70
The equations governing the top layer are:
dV
+ f k × V + ∇Φ = FV, (5.1)dt
dΦ
+ Φ∇ · V = FΦ, (5.2)dt
where V = ui + vj is the horizontal velocity vector and i and j are the unit
eastward and northward directions, respectively. The quantities FV and FΦ re-
fer to the forcing terms in the momentum equation and mass continuity equa-
tion, respectively. In this work, these terms are tied to stellar insolation, but in
general they represent any source or sink terms. The free surface geopotential
is given by Φ ≡ gh where g is the gravitational acceleration and h is geopotential
height. The geopotential height represents the height of the pressure surface of
the modeled upper layer. The geopotential can be separated into the constant
component Φ = gH (where H is the scale height) and the time-dependent devi-
ation Φ. We follow this convention here (following, e.g. Hack and Jakob, 1992;
Perez-Becker and Showman, 2013).
The geopotential forcing is represented by linear relaxation to the local ra-
diative equilibrium state Φeq, given by:
Φ + ∆Φeq cos λ cos ϕ on the dayside,Φeq =  (5.3)0 on the nightside,
where ∆Φeq is the maximal day-night contrast assuming the thickness of the up-
per layer in radiative equilibrium. The ratio ∆Φeq/Φ serves as a proxy for the
equilibrium day-night temperature contrast. The radiative equilibrium thick-
ness is highest at the substellar point, which occurs at latitude 0◦, longitude
71
0◦. This equation represents the radiative heat the planet receives from its star,
hence the forcing is only present on the dayside. Here we assume that sub-
Neptune planets are orbiting close to their host stars and synchronously rotat-
ing with the rotational period equal to the orbital period (1:1). The bulk of the
known sub-Neptune population has orbital periods of less than 30 days (Bean
et al., 2021) that place them generally within the tidal locking distance of their
host stars (e.g Barnes, 2017; Leconte et al., 2015). Hence the substellar point
and the radiative equilibrium profile are constant in time in our model. Linear
relaxation to the local radiative equilibrium takes the following form:
Φeq − Φ
FΦ = ≡ Q, (5.4)
τrad
where τrad is the radiative time scale. In this study, the timescale τrad is a free
parameter. We give an overview of the values of τrad we sample in Section 5.2.3.
If the geopotential is below the local radiative equilibrium, the forcing will
act to increase it, and vice versa. In the shallow-water framework, regions of
higher geopotential correspond to regions of heating, since mass is transferred
into the upper layer in the presence of radiative forcing.
In the momentum equation, the mass transfer between the two layers in the
model is governed by the following expression:

F 

V =  QV
− Q > 0,
Φ (5.5)
0 Q ≤ 0.
This term captures the effect of momentum advection from the lower layer on
the upper layer. The term R captures the mass transferred from the lower layer
to the upper layer through heating (Q > 0). The mass transfer from the up-
per layer to the lower layer due to cooling, on the contrary, does not affect the
72
specific momentum of the upper layer, and so R is zero wherever Q < 0. The
expression for R used here is also used in Perez-Becker and Showman (2013);
Showman and Polvani (2011); Showman et al. (2012); Shell and Held (2004).
We solve Equations 5.1 and 5.2 using a method based on Hack and Jakob
(1992) implemented in Python. We use a spectral method with a global, spher-
ical geometry, which allows us to perform differentiation using a truncated se-
quence of spherical harmonics, resulting in relatively fast computation time.
The non-linear terms are evaluated in physical space, and are then transformed
into the wavenumber (spectral) space using a spectral transform, which is a
composition of the Fast Fourier Transform and Gaussian quadrature. We use
the SciPy (Virtanen et al., 2020) implementation of the fast Fourier transform
and our implementation of the Gaussian quadrature due to frequent manipu-
lation of Fourier terms prior to the computation of Gaussian quadrature. In
wavenumber space, linear terms and derivatives are computed. Time-stepping
is performed and prognostic variables are evaluated in spectral space. The prog-
nostic variables are then transformed back into physical space, where the diag-
nostic and non-linear terms are evaluated.
Following Hack and Jakob (1992), we use Gaussian latitudes (µ = sin ϕ), uni-
formly spaced longitude, and U = u cos ϕ, V = v cos ϕ for zonal and merid-
ional winds, since Robert (1966) showed that the variables U and V are better
suited to computations in spectral space. While the quantities of interest are the
geopotentialΦ and the horizontal wind vector field V, we rely on the equivalent
vorticity/divergence form of the shallow-water equations, where the prognos-
tic variables are geopotential Φ, absolute vorticity η, and divergence δ. In this
form, the shallow-water equations have terms that have simple representations
73
in spectral space. Somewhat counterintuitively, U and V become diagnostic
variables, albeit readily determined from η and δ at any time step. In physi-
cal space, the variables can be represented as functions of longitude, latitude,
and time: prognostic variables Φ(λ, ϕ, t), η(λ, ϕ, t), δ(λ, ϕ, t), and diagnostic vari-
ables U(λ, ϕ, t), V(λ, ϕ, t). Note that there is no explicit dependence on pressure,
since the shallow-water system is a barotropic formulation. Readers interested
in the derivation of the vorticity/divergence form and the diagnostic relation-
ship should consult, e.g. Hack and Jakob (1992).
We use the spectral truncation T42, corresponding to a spatial resolution of
128 × 64 in longitude and latitude, with 2.81◦ cell size. Given that we wish to
study global-scale flows in sub-Neptune atmospheres, we opt to use the resolu-
tion which provides the best balance between computational efficiency and spa-
tial resolution. Previous studies using GCMs and shallow-water models agree
that the dynamics on tidally-locked exoplanets are global in scale (e.g. Show-
man and Polvani, 2011; Showman et al., 2015; Heng et al., 2011; Perez-Becker
and Showman, 2013; Zhang and Showman, 2017). We provide scale analysis for
the simulations in our ensemble in Section 5.4. We use modified Euler’s method
time-stepping scheme, rather than a more common semi-implicit scheme, due
to stability considerations (following Langton, 2008a). We apply two filters to
further ensure numerical stability: the modal splitting filter (Hack and Jakob,
1992), and a ∇6 hyperviscosity filter (based on Gelb and Gleeson, 2001). Models
are integrated from an initially flat layer (no deviation from Φ) and zero winds.
We select a time step that meets the CFL condition (Boyd, 2001). In general, sim-
ulations of highly irradiated planets require a shorter time step due to higher
wind speeds. The time steps used in our simulations range from 30 s to 180 s.
Each model is integrated for 1000 Earth days. We monitor the RMS windspeeds
74
of each simulation to ensure that they reach a steady state, which typically oc-
curs between 50 and 100 simulated days, but may take longer for simulations of
weakly forced planets.
Our model is consistent with the shallow-water framework as it has been
previously applied to hot Jupiters. Specifically, the momentum advection term
FV matches the term R used in Shell and Held (2004); Showman and Polvani
(2011); Showman et al. (2013). The equations (5.1) and (5.2) are consistent with
the form given in Vallis (2017). These equations also reduce to the form given
in Perez-Becker and Showman (2013) for a regime where the atmospheric drag
timescale τdrag = ∞, equivalent to the absence of atmospheric drag. Perez-Becker
and Showman (2013) have explored how the addition of atmospheric drag inter-
acts with the radiative timescale. Since the focus of this work is on the interplay
of radiative timescale and rotation rate, which has not yet been explored, we
have chosen to omit the drag term. We assume that the modeled layer is suf-
ficiently deep that there is no interaction with the planetary surface. Note also
that while Perez-Becker and Showman (2013) use the geopotential height as the
prognostic variable, the form we use in our model is equivalent since the geopo-
tential height and the geopotential differ only by a factor of reduced gravity g.
5.2.2 Model validation
In order to validate our model, we replicate three hot Jupiter regimes from
Perez-Becker and Showman (2013), shown in Figure 5.1, which matches the
three middle panels in the top row of Figure 3 in Perez-Becker and Showman
(2013). While the model we are replicating includes the effects of atmospheric
75
rad = 0.1 rad = 1 rad = 10
50
0
50
100 0 100 100 0 100 100 0 100
1.2e+04 4.4e+05 8.6e+05 1.3e+06 1.7e+06 2.1e+06 2.6e+06 3.0e+06 3.4e+06 3.8e+06
Figure 5.1: Recreation of a hot Jupiter regime in Figure 3 in Perez-Becker and
Showman (2013). Shown in the figure are maps of steady-state geopoten-
tial contours and the wind vector fields for the equilibrated solutions of the
shallow-water model. The planetary parameters are as follows: planetary ra-
dius a = 8.2 × 107 m, planetary rotation rate Ω = 3.2 × 10−5 rad/s, no drag,
reference geopotential Φ = 4 × 106, for strong radiative forcing (∆Φeq/Φ = 1).
Three radiative time scales are depicted: τrad=0.1 days (left panel), τrad=1 day (cen-
ter panel), and τrad=10 days (right panel). We have subtracted the constant value
of Φ from each panel. Panels share the same color scale for geopotential, but
winds are normalized independently in each panel. The substellar point is lo-
cated at (0◦, 0◦) for each panel.
drag, in this regime, at τdrag = ∞, which corresponds to the absence of drag.
We’ve chosen the color bar to match that of Figure 3 in Perez-Becker and Show-
man (2013). The other parameters here are ∆Φeq/Φ = 1 (high contrast regime),
geopotential at scale height Φ = gH = 4 × 106 m2 s−2. The planetary radius
a = 8.2 × 107 m and the rotation frequency Ω = 3.2 × 10−5 s−1 (equivalent to
Prot = 2.27 days) were selected in Perez-Becker and Showman (2013) based on
those of HD 189733b. Our simulations match the expected behavior, replicating
the transition from shorter τrad, where the planetary behavior resembles the ra-
diative forcing profile accompanied by a day-to-night flow, to longer τrad, where
the longitudinal gradients are small and the structure exhibits a hot equatorial
band.
76
5.2.3 Parameter regimes
The modeling phase space of sub-Neptunes is vast. For example, the population
of sub-Neptunes detected by Kepler exhibits orbital periods ranging from 0.4 to
300 Earth days Bean et al. (2021). These planets exhibit a variety of possible
atmospheric compositions and orbit a range of host stars. Therefore we select
to vary the radiative timescale to capture the breadth of possible atmospheric
compositions, and vary the intensity of incoming stellar radiation to model the
possible range of stellar host types and orbital distances. Under the constraint of
synchronous rotation, which necessitates planets to be orbiting within the tidal
locking distance, we focus on the interplay of rotation period with the radiative
timescale and the strength of instellation.
While sub-Neptune planets range in size from 1.6 to 4 times the radius of
Earth, here we select a representative sub-Neptune planetary radius equal to
three Earth radii: a = 3a⊕ = 1.91 × 107 m (Bean et al., 2021). The reduction
in radius compared to a Jupiter-sized planet RJ = 8.2 × 107 leads to significant
changes in the atmospheric circulation by affecting the advective timescale and
the Coriolis force. The advective timescale τadv ≈ a/U decreases as the plan-
etary radius decreases, in turn increasing the efficiency of heat transport and
leading to lower overall temperature contrasts. While the Coriolis parameter
f = 2Ω sin ϕ does not depend on the planetary radius, its meridional gradient
β = 2Ω cos ϕ/a increases as the planetary radius decreases, affecting the forma-
tion and velocity of planetary-scale Rossby waves. The Rossby waves form as
a response to potential vorticity gradient, which is caused by differential rota-
tion (different latitudes moving at different velocities). The phase speed ω of
the Rossby waves in the zonal direction is always westward and directly pro-
77
portional to the medional gradient of the Coriolis parameter β (see, e.g. Vallis,
2006). Therefore we expect that decreasing the radius would lead to the rise of
stronger westward Rossby waves compared to the hot Jupiter regime.
The reference geopotential Φ comes from separating the geopotential into
the time invariant spatial mean Φ and the time-dependent deviation from the
mean Φ (for a derivation from the continuity and momentum equations, see,
e.g. Hack and Jakob (1992)). The reference geopotential Φ can additionally be
thought of as geopotential at scale height H, Φ = gH, as in Perez-Becker and
Showman (2013). We nominally select g = 9.8 m/s2, but it should be noted
that g is included as a factor in geopotential and does not appear separately in
the model. Within the context of a shallow-water model, different values of Φ
can be interpreted as probing atmospheric layers with different scale heights,
which in turn depend on the mean molecular weight of the atmosphere and
the average temperature of the layer. The sample of reference geopotential Φ,
varying 4 × 106 to 106 m2s−2, is calculated based on a representative sample of
equilibrium temperatures ranging from 400 K to 1200 K and the mean molecu-
lar weights ranging from 2.3 to 5.5. This range of mean molecular masses cor-
responds to atmospheric metallicities ranging from 1 to 200 times solar values
(e.g. Goyal et al., 2019), which is consistent with the expected range of values
for exoplanets with masses between three and thirty Earth masses (Wakeford
and Dalba, 2020; Welbanks et al., 2019). We give the details of this calculation in
Section 5.2.4. Our model is agnostic of atmospheric composition, but we have
sampled a broad range of values of Φ, which incorporate the dependence on at-
mospheric composition and temperature via the scale height. Further, we have
modeled those layers under a variety of forcing conditions, which are related
to photospheric opacities and stellar insolation that set photospheric pressures
78
and temperatures.
We vary the planetary rotation period Prot to be 1, 5, and 10 Earth days to
reflect a variety of plausible rotation periods of synchronously rotating sub-
Neptunes (Yang et al., 2014). We are assuming Prot equal to the orbital period
Porb. Figure 5.2 illustrates the equilibrium temperatures in the sampled region
as a function of planetary rotation periods. All the planets in our sample lie
within the tidal locking distance (varying from 0.1 to 0.8 AU). Their short orbital
periods imply that our simulations could be matched to plausible targets for at-
mospheric characterization observations with facilities like JWST. Our samples
of reference geopotentials and rotation periods fall within a plausible window
for a variety of potential insolations.
Gravity is held constant at g = 9.8 m/s. While this value is based on Earth,
there is reason to believe that this gravity value might be representative of the
known exoplanet population (Ballesteros and Luque, 2016). This value is the
reduced gravity, which measures the restoring force at the interface between the
two layers. We follow Perez-Becker and Showman (2013) in their interpretation
of g in the shallow-water context.
The radiative timescale τrad is the time for gas parcels to reach local radiative
equilibrium. We follow Showman and Guillot (2002), where
∆p cp
τrad ∼ , (5.6)g 4σT 3
where ∆p is the thickness of the pressure layer in Pascals, cp is the specific heat
capacity at pressure p, σ is the Stefan-Boltzmann constant, and T is the temper-
ature of the layer. The radiative timescale represents how quickly the parcels
adjust to radiative changes. While τrad is a free parameter in this study, varying
the radiative timescale between 0.1, 1, and 10 Earth days corresponds to probing
79
a variety of pressure level and temperature combinations. While the radiative
time constant τrad decreases sharply with increasing temperature, note that the
estimate breaks down in deep, optically thick atmospheres, e.g. deep inside hot
Jupiter atmospheres, τrad is long. Our samples span the range from a few hours,
typical for hot Jupiters (Showman et al., 2010) to a few days, approximately
corresponding to Earth (Showman et al., 2008b). Showman and Guillot (2002)
predicted that on hot Jupiters, large day/night fractional differences will persist
when τrad ≪ τadv and small day/night fractional differences will appear when
τrad ≫ τadv.
We investigate three different relative radiative forcing amplitudes, which
prescribe the radiative equilibrium temperature differences between the day-
side and the nightside. We let ∆Φeq/Φ = 1, 0.5, and 0.1, corresponding to hot,
warm, and cool sub-Neptune planets. In conjunction with our choices of ref-
erence geopotential Φ, this choice allows us to span the vast parameter space
of potential stellar forcings and atmospheric scale heights with our simplified
model. In the context of hot Jupiters, the high-contrast regime has been ex-
plored by Perez-Becker and Showman (2013), and the medium-contrast regime
has been explored by Liu and Showman (2013). Both of the above studies also
include a low-contrast regime, with ∆Φeq/Φ varying from 0.001 to 0.1.
5.2.4 Motivation for the choices of Φ
We have chosen the mean geopotential Φ by calculating the scale heights for
a variety of equilibrium temperatures and planetary atmospheric compositions.
While the shallow-water model is agnostic about atmospheric compositions, we
80
2000
K0
K5
1600 M0M5
1200
sampled region
800
400
1 5 10
Planetary rotation/orbital period Prot [days]
Figure 5.2: Planetary equilibrium temperature and as a function of rotation pe-
riod for planets around K0 (purple), K5 (blue), M0 (yellow), and M5 (red) stars.
The chosen grid of rotation periods of 1, 5, and 10 days and of planetary equi-
librium temperatures of 400, 800, and 1200 K generates a realistic sample. The
sample region is highlighted in grey. The stellar parameters were taken from
Zombeck (1990).
have selected the sampled scale heights and forcing amplitudes to encompass a
variety of plausible scenarios. We select the mean geopotential to correspond to
scale height H, which we compute as follows:
kTeq
Φ = gH = , (5.7)
m
where k is the Boltzmann constant, Teq is the equilibrium temperature, and
m (g/mol) is the mean molar mass of the atmospheric constituents. The range
of scale heights was computed using temperatures in the 400 − 1200 K range
and mean molecular weights of 2.3-5.5 times the mass of a hydrogen atom, cor-
responding to 1-200 multiples of solar metallicity. The resulting values for Φ
were relatively evenly spread out between 6 × 105 m2s−2 and 4.3 × 106 m2s−2,
so the samples of 106 m2s−2, 2 × 106m2s−2, 3 × 106 m2s−2, and 4 × 106 m2s−2 ef-
81
Equilibrium temperature Teq [K]
ficiently sample the parameter space, despite combining two parameters into
one. The models at high Φ can be interpreted as corresponding to planets with
a hydrogen-dominated atmosphere, similar to that of Jupiter. A decrease in Φ
can be thought of as a transition to a more metal-rich atmosphere. While the
strength of stellar forcing is not independent of reference geopotential in our
model, our simulations span three forcing regimes for each scale height.
5.3 Results
5.3.1 Extension to Sub-Neptunes
After benchmarking our model in the hot Jupiter regime, we transition the
model toward a regime appropriate for sub-Neptunes by first decreasing the
planetary radius a. We decrease the planetary radius from that of Jupiter to
equal three times Earth radius (a = 1.9× 107 m) while keeping the rest of the pa-
rameters identical to our validation regime: a hot Jupiter based on HD 189733b,
explored in Section 5.2.2. Compared to our replication of a hot Jupiter regime,
reducing the planetary radius still preserves the transition from day-to-night
flow for short radiative timescales (0.1 days) to a jet-dominated flow for long
radiative timescales (10 days). Compared to Figure 5.1, the sub-Neptune sim-
ulations shown in Figure 5.3 exhibit lower geopotential height contrast both in
the zonal direction (between the dayside and the nightside) and in the merid-
ional direction (between the equator and the pole). We also observe that for
short and medium values of τrad (0.1 and 1 days, left and center panels, respec-
tively), the nightside cyclonic vortices are more confined in latitude. For the
82
High forcing
τrad = 0.1 days τrad = 1 day τrad = 10 days
Vmax ≈ 1700 m/s Vmax ≈ 300 m/s Vmax ≈ 160 m/s
Prot = 1 day ∆Φ ≈ 1.1 × 106 m2/s2 ∆Φ ≈ 5.1 × 103 m2/s2 ∆Φ ≈ 3.4 × 102 m2/s2
Ro=0.62 Ro=0.11 Ro=0.06
Vmax ≈ 900 m/s Vmax ≈ 570 m/s Vmax ≈ 80 m/s
P = 5 days ∆Φ ≈ 5.7 × 105 m2/s2 ∆Φ ≈ 2.4 × 105 m2/s2 ∆Φ ≈ 8.4 × 102rot m2/s2
Ro=1.59 Ro=1.03 Ro=0.15
Vmax ≈ 600 m/s Vmax ≈ 330 m/s Vmax ≈ 145 m/sU
P 10 days max
≈ 520 m/s
= 5 2rot ≈ × 3 m2/s2 ∆Φ ≈ 1.1 × 10 m /s
2 ∆Φ ≈ 7.1 × 103 m2/s2
∆Φ 7.2 10
Ro=1.87 Ro=1.17 Ro=0.53
Medium Forcing
τrad = 0.1 days τrad = 1 day τrad = 10 days
Vmax ≈ 1200 m/s Vmax ≈ 160 m/s Vmax ≈ 120 m/s
Prot = 1 day ∆Φ ≈ 6.4 × 105 m2/s2 ∆Φ ≈ 7.1 × 103 m2/s2 ∆Φ ≈ 3.4 × 102 m2/s2
Ro=0.44 Ro=0.06 Ro=0.04
Vmax ≈ 600 m/s V ≈ 490 m/s V ≈ 80 m/s
Umax ≈ 675 m/s max maxP 4 2 2 2 2 2rot = 5 days ∆Φ ≈ 2.4 × 10 m /s ∆Φ ≈ 2.0 × 10 m /s∆Φ ≈ 7.2 × 104 m2/s2
Ro=1.42 Ro=0.87 Ro=0.15
Vmax ≈ 430 m/s Vmax ≈ 280 m/s Vmax ≈ 200 m/s
Prot = 10 days ∆Φ ≈ 1.3 × 104 m2/s2 ∆Φ ≈ 8.7 × 103 m2/s2 ∆Φ ≈ 1.5 × 104 m2/s2
Ro=1.53 Ro=0.99 Ro=0.71
Low Forcing
τrad = 0.1 days τrad = 1 day τrad = 10 days
Vmax ≈ 330 m/s Vmax ≈ 50 m/s Vmax ≈ 60 m/s
P = 1 day ∆Φ ≈ 4.1 × 104 m2/s2 ∆Φ ≈ 1.4 × 103 m2/s2 ∆Φ ≈ 6 × 10 m2 2rot /s
Ro=0.10 Ro=0.02 Ro=0.02
Vmax ≈ 440 m/s Vmax ≈ 80 m/s Vmax ≈ 100 m/s
Prot = 5 days ∆Φ ≈ 5.3 × 104 m2/s2 ∆Φ ≈ −5 × 10 m2/s2 ∆Φ ≈ 2.4 × 102 m2/s2
Ro=0.79 Ro=0.16 Ro=0.20
Vmax ≈ 310 m/s Vmax ≈ 160 m/s Vmax ≈ 100 m/s
Prot = 10 days ∆Φ ≈ 1.4 × 104 m2/s2 ∆Φ ≈ 8.0 × 103 m2/s2 ∆Φ ≈ 6 × 10 m2/s2
Ro=1.13 Ro=0.82 Ro=0.37
Table 5.1: Typical maximal wind speeds, day/night geopotential contrast, and
Rossby number for a representative scale height Φ = 4 × 106 m2/s2 based on
temporal averages. The geopotential contrast is computed by subtracting the
nightside average from the dayside average. When maximal wind speeds Vmax
values differ from maximal values of the zonal component Umax, the zonal max-
imum Umax is also included.
83
short timescale τrad = 0.1 days, the vortices extend from ≈ 10◦ to ≈ 50◦ latitude,
with wind speeds |V|≈ 103 m/s and the geopotential contrast of ≈ 106 within
the vortices. For the medium radiative timescale τrad = 1 day, the vortices extend
from ≈ 15◦ to ≈ 70◦ latitude, with wind speeds |V|≈ 8×102 m/s and geopotential
contrast ≈ 6 × 105 within the vortices. While the rotation rate Ω (rad/s) was not
changed from that one used in Perez-Becker and Showman (2013), the change
in planetary radius leads to an inversely proportional change in the beta param-
eter β ∂ f= y 2Ω cos ϕ0/a. The beta parameter is in turn proportional to the speed∂
of Rossby waves. In the left panel, for τrad = 0.1 days, the vortices are centered
at ∼ 30◦ latitude. In the center panel, where τrad = 1 day, the vortices are larger
and father from the equator, at ∼ 40◦ latitude, but exhibit lower geopotential
contrast. The center panel simulation also exhibits a secondary, warmer pair of
cyclonic vortices. Reducing the planetary radius by 75% has led to the emer-
gence of significant changes in global scale atmospheric flow patterns and wind
speeds, showing distinct differences in global circulation patterns between hot
Jupiters and sub-Neptunes.
We now proceed to model sub-Neptunes by performing a parameter sweep
of rotation periods Prot and radiative timescales τrad in three forcing regimes:
strong forcing (∆Φeq/Φ = 1, Section 5.3.2), medium forcing (∆Φeq/Φ = 0.5, Sec-
tion 5.3.3), and weak forcing (∆Φeq/Φ = 0.1, Section 5.3.4). In Table 5.1, we sum-
marize the maximal wind speeds Vmax, and the geopotential contrast between
the dayside and the nightside. The maximal equatorial zonal speeds Umax are
similar to the maximal overall wind speeds. While our parameter sweep in-
cludes four scale height values, a representative scale height corresponding to
Φ = 4 × 106 was chosen to be represented in the table.
84
rad = 0.1 rad = 1 rad = 10
50
0
50
100 0 100 100 0 100 100 0 100
0.15 0.30 0.45 0.60 0.75 0.90 1.05 1.20 1.35
1e6
Figure 5.3: Recreation of a hot Jupiter regime in Figure 3 in Perez-Becker and
Showman (2013), with planetary radius turned down to a = 3a = 1.3 × 107⊕
m. Shown are the geopotential contour maps and the wind vector fields.
The remaining planetary parameters are as follows: planetary rotation rate
Ω = 3.2 × 10−5 rad/s, no drag, reference geopotential Φ = 4 × 106, for strong
radiative forcing (∆Φeq/Φ = 1). Three radiative time scales are depicted: τrad=0.1
days (left panel), τrad=1 day (center panel), and τrad=10 days (right panel). We have
subtracted the constant value of Φ from each panel. Panels share the same color
scale for geopotential, but winds are normalized independently in each panel.
The substellar point is located at (0◦, 0◦) for each panel.
5.3.2 Strong forcing simulations
In this section, we explore what we are calling “strong forcing” simulations in
the context of our ensemble. For these simulations the prescribed local radiative
equilibrium day-night forcing contrast is equal to the the nightside geopotential
at radiative equilibrium (Φ = gH, where H is the scale height). This implies that
the equilibrium day-night contrast is of the same magnitude as the equilibrium
nightside geopotential. Note that changing the scale height results in changing
the amplitude of the local radiative equilibrium, and so depending on the value
of Φ, the stellar forcing varies.
First we consider Φ = 4 × 106 m2s−2, corresponding to a large scale height.
While our simplified shallow-water model is agnostic to the exact specification
85
of the planetary composition, a large scale height could model, for example,
a predominantly H/He atmosphere with ≈ 1× solar metallicity. The resulting
geopotential profile for this regime is shown in Figure 5.4. The corresponding
mean zonal wind profiles are shown in Figure 5.5. In the top left panel of Figure
5.4, for the short-timescales regime, where radiative timescale τrad = 0.1 days
and rotation period Prot = 1 day, strong forcing is combined with strong Cori-
olis force. This regime exhibits high day-night contrast combined with a high
equator-to-pole contrast. We observe atmospheric superrotation with a strong
eastward equatorial jet (the mean zonal wind speeds reach ∼ 0.6 km s−1 at the
equator) that shifts the dayside hotspot to the east. A combination of eastward
Kelvin waves and westward Rossby waves results in equatorial superrotation
that gives rise to a hotspot shape and offset from the substellar longitude char-
acteristic for hot Jupiters (Showman et al., 2020). Indeed, this strongly forced,
short-timescale regime is closest to a hot Jupiter regime among the ones we ex-
plore. On the night side, the simulation exhibits standing polar vortices. In
contrast, in the right column of Figure 5.4, long radiative timescale τrad = 10
days results in the formation of a jet-dominated pattern with little longitudinal
variation, regardless of the rotation period. The equator-to-pole variation is sig-
nificantly lower than that of the prescribed radiative equilibrium. No vortices
are present. In this regime, the radiative timescale is a few times longer than
the advective timescale, and since the dayside radiative forcing is the primary
driver of longitudinal variation, this ratio results in longitudinally near-uniform
geopotential distribution. squashing any longitudinal variation.
In the top center panel of Figure 5.4, the simulation exhibits a transition be-
tween the one dominated by the day-night flow, which occurs for shorter radia-
tive timescales, and the one dominated by the jet flow, which occurs for longer
86
eq/ = 1, = 4x10^6
rad [days]
0.1 1 10
50
1.77 × 106
0
1.57 × 106
50
1.38 × 106
50 1.18 × 10
6
0 9.83 × 105
50 7.86 × 105
5.90 × 105
50 3.93 × 105
0 1.97 × 105
50
0.00 × 100
100 0 100 100 0 100 100 0 100
Figure 5.4: Averaged maps of geopotential contours and wind vector fields for
the equilibrated solutions of the shallow-water model for the strong forcing
regime (∆Φeq/Φ = 1) at reference geopotential Φ = 4 × 106 m2s−2. The con-
stant value of Φ has been subtracted from each panel. Each panel in the grid
corresponds to a unique combination of planetary rotation period Prot and ra-
diative timescale τrad. Panels share the same color scale for geopotential, but
winds are normalized independently in each panel. The substellar point is lo-
cated at (0◦, 0◦) for each panel.
radiative timescales. There is some longitudinal variation, but the hottest parts
of the planet are no longer constrained to the dayside. There is high equator-to-
pole variation. Polar vortices are present, extending to ≈ 60◦ latitude, centered
at the terminator 90◦ east of the substellar point. The equator-to-pole contrast is
significantly lower than that of the prescribed radiative equilibrium. The aver-
age zonal wind speeds exhibit peaks of ≈ 150 m/s at ∼ 40◦ − 50◦ latitude, and
these wind speeds are much lower than the ≈ 600 m/s eastward equatorial jet
exhibited by the τrad = 0.1 days regime.
87
Prot [days]
10 5 1
eq/ = 1, = 4x10^6
rad = 0.1 rad = 1 rad = 10
50
0
50
50
0
50
50
0
50
0 200 400 600 0 200 400 600 0 200 400 600
Figure 5.5: Mean-zonal winds U for the equilibrated solutions of the shallow-
water model for the strong forcing regime (∆Φeq/Φ = 1) at reference geopotential
Φ = 4 × 106 m2s−2. Each panel in the grid corresponds to a unique combination
of planetary rotation period Prot and radiative timescale τrad.
As Prot increases from 1 to 5 days, the eastward equatorial jet becomes
weaker, with lower associated mean zonal wind speeds. The models shown
in the middle row of Figure 5.4 exhibits a robust day-to-night flow pattern.
For τrad = 0.1 days (left column), compared to faster rotation rates, the night-
side cyclonic vortices are smaller and closer to the equator (centered at around
∼ 20◦ − 30◦ latitude). The cyclonic vortices correspond to the westward mean
zonal winds around those latitudes. The equatorial zonal winds are eastward,
and the amplitude is lower than that for shorter rotation period values. When
Prot = 10 days, the simulation exhibits a cool region west of the hot spot, ap-
proximately at the terminator. When Prot = 5 days and τrot = 1 days, we observe
cyclonic nightside vortices centered at ∼ 30◦ latitude just west of the termina-
tor. We also observe a secondary pair of warmer but still cyclonic vortices just
88
Prot = 10 Prot = 5 Prot = 1
east of the terminator (∼ 100◦ longitude). Apart from the cyclonic vortices, the
equator-to-pole contrast is relatively low. Further increasing the rotation period
Prot to 10 days, in the bottom center panel, this regime exhibits nightside cyclonic
vortices, and the vortices are shifted west compared to the middle panel with
the same radiative timescale and a shorter rotation period. This behavior cre-
ates a geopotential nightside asymmetry, where the west half of the nightside
exhibits larger geopotential than the east half. Overall, the shift from Prot = 1
day to Prot = 5 days creates a more dramatic change in the atmospheric circula-
tion patterns than the shift from Prot = 5 days to Prot = 10 days. The shift from
short radiative timescales to longer ones corresponds to the shift from a day-to-
night flow with nightside vortices to a jet-dominated flow, as well as from a low
time-variability to a higher time-variability.
As expected, the wind speeds are strongest (≈ 600 m/s) when both radia-
tive and rotational timescales are short, and weakest (≈ 10 − 20 m/s) when
both timescales are longest. Most of our simulations in the high-forcing region
in the parameter space exhibit an equatorial jet pattern, but when the radiative
timescale is long, but the rotation period is short (top center and top right in Fig-
ure 5.5), the average zonal wind speeds exhibit two tropical peaks. We discuss
this transition in more detail in Section 5.4.2.
As we reduce the scale height (and hence Φ, the relative forcing also de-
creases, since our forcing term is defined in terms of the fractional difference.
Qualitatively, the patterns in atmospheric dynamics similar to the one shown in
Figure 5.4.1 The mean zonal winds for Φ = 3 × 106 m2s−2 are similar in speed
and amplitude to those occurring at higher scale height (Φ = 4 × 106, Figure
1We show the geopotential and mean-zonal wind grids for intermediate scale heights corre-
sponding to Φ = 3 × 106 m2s−2 and Φ = 2 × 106 m2s−2 in Appendix D.
89
eq/ = 1, = 1x10^6
rad [days]
0.1 1 10
50
5.97 × 105
0
5.30 × 105
50
4.64 × 105
50 3.98 × 105
0 3.31 × 105
50 2.65 × 105
1.99 × 105
50 1.33 × 105
0 6.63 × 104
50
0.00 × 100
100 0 100 100 0 100 100 0 100
Figure 5.6: Same as Figure 5.4, but at a low scale height corresponding toΦ = 106
m2s−2. Note the color-scale has been adjusted. Wind speeds are normalized
individually for each panel.
5.5), ranging from ≈ 10 − 20 m/s for the longest timescales to ≈ 600 m/s for the
shortest timescales. As we continue to decrease the scale height to Φ = 2 × 106
m2s−2, the qualitative behavior remains the same while the geopotential ampli-
tude and the mean zonal wind speed decrease, ranging from ≈ 10 − 20 m/s for
the longest timescales to ≈ 500 m/s for the shortest timescales.
We show the high-forcing simulations corresponding to the smallest scale
height in our ensemble (with Φ = 1 × 106 m2s−2) in Figure 5.6. Note that the col-
orbar limits for the geopotential figures are the same across panels but differ for
each figure in order to better show the structure of each scale height regime. It is
interesting to note that as the scale height decreases andΦ = 1×106 m2s−2, the hot
spot in the short time scales regime (where τrad = 0.1 days, Prot = 1 day) changes
90
Prot [days]
10 5 1
eq/ = 1, = 1x10^6
rad = 0.1 rad = 1 rad = 10
50
0
50
50
0
50
50
0
50
100 0 100 200 300 100 0 100 200 300 100 0 100 200 300
Figure 5.7: Same as Figure 5.5 but at a low scale height corresponding to Φ =
106 m2s−2. Note the limits on the horizontal axis have been adjusted between
figures.
its shape from one matching hot Jupiters (Figure 5.4) to a more elongated one in
Figure 5.6. This change is driven by a decrease in Rossby deformation radius,
which is proportional to the scale height.
5.3.3 Medium forcing simulations
In this section, we explore the medium-forcing regime, where the prescribed lo-
cal radiative equilibrium day-night contrast is equal to half of the mean geopo-
tential (∆Φeq/Φ = 0.5). While the strong forcing regime (∆Φeq/Φ = 1) has been
widely explored in literature in the context of hot Jupiters (e.g. Perez-Becker and
Showman, 2013; Showman and Polvani, 2011), the medium forcing regime has
not been investigated as thoroughly (e.g. Liu and Showman, 2013). This regime
91
Prot = 10 Prot = 5 Prot = 1
is relevant to the modeling of exoplanet atmospheres since many sub-Neptunes
orbit closely around cooler stars, such as M-dwarfs, and are synchronous rota-
tors (Bean et al., 2021). Such sub-Neptunes are also advantageous targets for
atmospheric characterization with observational facilities due to the favorable
ratio between the planet and star radii (Triaud, 2021).
The trends in qualitative behavior are the same as for high contrast (strong
forcing) regime, with only a few differences. The amplitude of the resulting
geopotential profile appears to scale linearly with the day-night contrast of the
local radiative equilibrium.
For large scale height (Φ = 4 × 106 m2s−2), the geopotential profile is shown
in Figure 5.8, and the corresponding mean zonal wind profiles are shown in
Figure 5.9. We observe a transition from a high day-night contrast for shorter
radiative and rotational timescales (all three panels in the left column (τrad = 0.1
days, as well as center panel and the bottom panel in the middle row (τrad = 1
day and Prot = 5 and 10 days)) to a jet-dominated flow (especially pronounced
in the top center panel (τrad = 1 day, Prot = 5 days) and in all three panels in the
right column (τrad = 10 days) of Figure 5.8). In the left column (τrad = 0.1 days),
the simulations exhibit a hot spot. As the rotation period increases, the super-
rotation increases. In the center and bottom center panels (τrad = 1 day,Prot = 5
and 10 days), the models exhibit nightside cyclonic vortices. In the center left
panel (τrad = 0.1 days, Prot = 5 days, we observe a long-timescale oscillation with
momentum advection across the equator, which generates a cyclonic vortex in
the southern hemisphere in the averaged geopotential profile shown in Figure
5.8 and the mean-zonal wind profile in Figure 5.9.
While the qualitative trends in response to changing radiative timescale and
92
eq/ = 0.5, = 4x10^6
rad [days]
0.1 1 10
50
9.75 × 105
0
8.67 × 105
50
7.59 × 105
50 6.50 × 105
0 5.42 × 105
50 4.34 × 105
3.25 × 105
50 2.17 × 105
0 1.08 × 105
50
0.00 × 100
100 0 100 100 0 100 100 0 100
Figure 5.8: Averaged maps of geopotential contours and wind vector fields for
the equilibrated solutions of the shallow-water model for the medium forcing
regime (∆Φeq/Φ = 0.5) at reference geopotential Φ = 4 × 106 m2s−2. The constant
value of Φ has been subtracted from each panel. Each panel in the grid corre-
sponds to a unique combination of planetary rotation period Prot and radiative
timescale τrad. Panels share the same color scale for geopotential, but winds
are normalized independently in each panel. The substellar point is located at
(0◦, 0◦) for each panel.
rotation rate are similar to the corresponding scale height with strong forcing,
there are several quantitative differences. The overall geopotential contrast
tends to be smaller. The greatest geopotential contrast is ≈ 6.4 × 105 m2s−2
(for τrad = 0.1 days, Prot = 1 day). Compared to the strong-forcing regime,
the hotspot for the medium-forcing regime extends farther west at high lati-
tudes. The lowest geopotential contrast is ≈ 2 × 102m2s−2 (for τrad = 10 days,
Prot = 5 days). The qualitative wind speed trends for medium-forcing models
follow those for strong forcing, but the speeds are significantly lower (Figure 5.9
shows the mean zonal wind speeds ranging from ≈ 500 m/s to ≈ 50 m/s). In the
93
Prot [days]
10 5 1
eq/ = 0.5, = 4x10^6
rad = 0.1 rad = 1 rad = 10
50
0
50
50
0
50
50
0
50
0 200 400 0 200 400 0 200 400
Figure 5.9: Mean-zonal winds U for the equilibrated solutions of the shallow-
water model for the medium forcing regime (∆Φeq/Φ = 0.5) at reference geopo-
tential Φ = 4× 106 m2s−2. Each panel in the grid corresponds to a unique combi-
nation of planetary rotation period Prot and radiative timescale τrad.
Prot = 1 day, τrad = 1 day regime, the nightside vortices are still present but the jet
structure is more defined compared to the strong forcing regime. Overall, the
atmospheric dynamics of the medium forcing regime exhibit lower day-night
contrasts, weaker equatorial winds speeds, and different location of vortices,
compared to the strong forcing regime. However, the transitions between day-
to-night flow and jet-dominated flow occur at the same points in the τrad—Prot
parameter space as for the strong forcing regime.
As scale height decreases, the transitions from day-night flow to jet-
dominated flow occur at the same values of τrad and Prot as for the larger scale
height. The geopotential contrast appears to decrease linearly in response to a
decrease in forcing (recall that forcing is tied to scale height in our model). The
decrease in geopotential contrast and wind speeds is also pronounced when
94
Prot = 10 Prot = 5 Prot = 1
compared to the strong forcing regime at the corresponding height. For shortest
timescales (for τrad = 0.1 days, Prot = 1 day), the hotspot is elongated in the east-
west direction, and secondary eastward jets emerge in addition to the equatorial
jet. The elongation is due√to the decrease in the Rossby radius of deformation,
which is proportional to Φ. This behavior is illustrated in the top left panel
of Figure 5.10, which shows medium-forcing simulations for our lowest simu-
lated scale height, corresponding to Φ = 106 m2s−2. Compared to both larger
scale height regimes at medium forcing and same scale height regime at strong
forcing, both geopotential contrast and wind speeds are lower, while the jet
structure and transitions from day-night to jet-dominated flow are preserved.
The right columnn, corresponding the the long radiative timescale τrad = 10
days, exhibits jet-dominated flow. In the transitional regime (τrad = 0.1 days,
Prot = 1 day), the anticyclonic vortices are farther from the equator (centered at
≈ 45◦ latitude) compared to medium-forcing regimes at larger scale heights. In
comparison to the strong forcing regime at the same scale height, this regime
exhibits lower geopotential contrast (≈ 50% decrease) and lower wind speeds.
5.3.4 Weak forcing simulations
For these simulations, the local radiative equilibrium amplitude is equal to 0.1
of the mean geopotential (∆Φeq/Φ = 0.1), approaching the linear regime com-
monly assumed for Earth. Low-contrast regimes have been previously explored
in the context of hot Jupiters (e.g. Perez-Becker and Showman, 2013; Shell and
Held, 2004) due to their analytical tractability. Similarly to the medium-forcing
regime, this regime models sub-Neptunes orbiting closely around cooler stars.
95
eq/ = 0.5, = 1x10^6
rad [days]
0.1 1 10
50
2.98 × 105
0
2.64 × 105
50
2.31 × 105
50 1.98 × 105
0 1.65 × 105
50 1.32 × 105
9.92 × 104
50 6.61 × 104
0 3.31 × 104
50
0.00 × 100
100 0 100 100 0 100 100 0 100
Figure 5.10: Same as Figure 5.8, but at a low scale height corresponding to Φ =
106 m2s−2. Note the color-scale has been adjusted. Wind speeds are normalized
individually for each panel.
For large scale height (Φ = 4 × 106 m2s−2), the geopotential profile is shown
in Figure 5.12, and the corresponding mean zonal wind profiles are shown in
Figure 5.13. The maximal geopotential contrast is ≈ 1.75 × 105 m2s−2, and mean
zonal the wind speeds are in the ≈ 30 − 100 m/s range. When τrad = 0.1 days, in
the left column of Figure 5.12, we observe day-to-night flow, similar to strong
and medium forcing regimes. The short-timescales regime (τrad = 0.1 days,
Prot = 1 day, upper left panel) is qualitatively similar to the regime where τrad = 1
day, Prot = 1 day for strong and medium forcing. This regime is a transition
between day-to-night and jet-dominated flow. The similarity occurs because
increasing τrad has the same effect on Equation 5.3 as decreasing ∆Φeq/Φ. This
regime exhibits visually identifiable cyclonic and anticyclonic Rossby gyres. As
rotation period increases, the cyclonic vortices remain, while the anticyclonic
96
Prot [days]
10 5 1
eq/ = 0.5, = 1x10^6
rad = 0.1 rad = 1 rad = 10
50
0
50
50
0
50
50
0
50
50 0 50 100 150 50 0 50 100 150 50 0 50 100 150
Figure 5.11: Same as Figure 5.9 but at a low scale height corresponding to Φ =
106 m2s−2. Note the limits on the horizontal axis have been adjusted between
figures.
vortices break down. As a result, the regimes where τrad = 0.1 days, Prot = 5
and Prot = 10 days (center left and bottom left of Figure 5.12) exhibit superro-
tated hot spots. As the radiative timescale increases to τrad = 1 day, for Prot = 1
day and Prot = 5 days, we observe regimes that are in transition between day-
to-night flow and jet-dominated flow. There is some longitudinal variation, but
the geopotential contrast is lower than for shorter radiative timescales. The cen-
ter panel exhibits hot band around the equator, but still maintains longitudinal
variation in the form of cyclonic polar vortices at the terminator. For the longest
radiative timescale, τrad = 10 days (right column of Figure 5.12), the simulations
exhibit jet-dominated flow for all rotation periods. Planets at faster rotation
rates exhibit more equator-to-pole geopotential contrast. The combination of
long radiative timescale (τrad = 10 days) and short rotation period (Prot = 1 day)
produces a regime with time variation, where a hot equatorial band alternates
97
Prot = 10 Prot = 5 Prot = 1
eq/ = 0.1, = 4x10^6
rad [days]
0.1 1 10
50
1.69 × 105
0
1.50 × 105
50
1.32 × 105
50 1.13 × 105
0 9.40 × 104
50 7.52 × 104
5.64 × 104
50 3.76 × 104
0 1.88 × 104
50
0.00 × 100
100 0 100 100 0 100 100 0 100
Figure 5.12: Averaged maps of geopotential contours and wind vector fields
for the equilibrated solutions of the shallow-water model for the weak forcing
regime (∆Φeq/Φ = 0.1) at reference geopotential Φ = 4 × 106 m2s−2. The constant
value of Φ has been subtracted from each panel. Each panel in the grid corre-
sponds to a unique combination of planetary rotation period Prot and radiative
timescale τrad. Panels share the same color scale for geopotential, but winds
are normalized independently in each panel. The substellar point is located at
(0◦, 0◦) for each panel.
with hot tropical bands with a periodicity of ∼ 50 hours. When both timescales
are long, we observe some temporal variability as the hotspot is advected west
by westward-propagating Rossby waves.
As we decrease the scale height, the qualitative behavior stays similar, with
transitions occurring at the same points in the parameter space. For Φ = 3 ×
106 m2s−2, the highest geopotential contrast decreases to ≈ 1.5 × 105, and the
mean zonal wind speeds decrease slightly, resulting in ≈ 25 − 100 m/s range.
For this and smaller scale heights, the radiative timescale (τrad = 10 days) and
98
Prot [days]
10 5 1
eq/ = 0.1, = 4x10^6
rad = 0.1 rad = 1 rad = 10
50
0
50
50
0
50
50
0
50
0 50 100 0 50 100 0 50 100
Figure 5.13: Mean-zonal winds U for the equilibrated solutions of the shallow-
water model for the weak forcing regime (∆Φeq/Φ = 0.1) at reference geopoten-
tial Φ = 4 × 106 m2s−2. Each panel in the grid corresponds to a unique combina-
tion of planetary rotation period Prot and radiative timescale τrad.
short rotation period (Prot = 1 day) regime no longer exhibits a separation of the
equatorial band into tropical bands, although the model still exhibits some time
variation. For Φ = 2× 106 m2s−2, the maximal geopotential contrast is ≈ 1.1× 105
m2s−2. The mean zonal winds decrease to ≈ 10 m/s in the shortest timescales
regime (τrad = 0.1 days and Prot = 1 day) and ≈ 80 m/s in the longest timescales
regime (τrad = 0.1 days and Prot = 1 day).
Finally, at the lowest scale height in our ensemble, Φ = 106 m2s−2, we once
again observe similar atmospheric patterns. The geopotential profile is shown
in Figure 5.14, and the corresponding mean zonal wind profiles are shown in
Figure 5.15. The maximal geopotential contrast is ≈ 7 × 104 m2s−2. The mean
zonal winds are in a range similar to Φ = 2 × 106 m2s−2, in the ≈ 10 − 80
m/s range. In the middle column of Figure 5.14, the transitional regimes dif-
99
Prot = 10 Prot = 5 Prot = 1
eq/ = 0.1, = 1x10^6
rad [days]
0.1 1 10
50
6.31 × 104
0
5.61 × 104
50
4.91 × 104
50 4.21 × 104
0 3.51 × 104
50 2.81 × 104
2.10 × 104
50 1.40 × 104
0 7.01 × 103
50
0.00 × 100
100 0 100 100 0 100 100 0 100
Figure 5.14: Same as Figure 5.12, but at a low scale height corresponding to Φ =
106 m2s−2. Note the color-scale has been adjusted. Wind speeds are normalized
individually for each panel.
fer slightly in structure, while still exhibiting a transtion from day-to-night flow
to a jet-dominated flow. Specifically, in the τrad = 1 day and Prot = 1 regime, the
anticyclonic vortices move poleward, centered at ≈ 40◦ latitude (compared to
≈ 25◦ latitude for Φ = 2 × 106 m2s−2).
Overall, compared to strong and medium forcing simulations, the weak forc-
ing exhibits lower geopotential contrast, lower wind speeds, more east-west
elongated structures, and more (periodic) time variability. Note also that for
strong and medium forcing, the highest mean zonal wind speeds occur when
both radiative timescale and rotation period are short, and the lowest wind
speeds occur when both timescales are long. This pattern is reversed in the
weak forcing regime. We discuss this transition further in 5.4.2. Weak forcing
100
Prot [days]
10 5 1
eq/ = 0.1, = 1x10^6
rad = 0.1 rad = 1 rad = 10
50
0
50
50
0
50
50
0
50
20 0 20 40 60 80 20 0 20 40 60 80 20 0 20 40 60 80
Figure 5.15: Same as Figure 5.13 but at a low scale height corresponding to Φ =
106 m2s−2. Note the limits on the horizontal axis have been adjusted between
figures.
induces weak winds, which result in a smaller advective term, especially at long
radiative timescale τrad=10 days.
5.4 Discussion
5.4.1 Holistic discussion of trends in the regimes
In this section, we discuss the overall trends in the atmospheric patterns in our
simulations, specifically horizontal winds and geopotential. Since atmospheric
heating will generally increase the geopotential Φ and atmospheric cooling will
generally reduce Φ, the qualitative changes in geopotential can be thought of as
changes in temperature, with higher thickness of the upper layer corresponding
101
Prot = 10 Prot = 5 Prot = 1
to higher temperatures. Using the ideal gas relationship Φ ≈ RT , where R is
the specific gas constant, we can approximate changes in temperature. For our
range of mean molecular weights, the approximate range for values of R is 750—
1800 J kg−1K−1. We can estimate that, for example, a change of 105 m2s−2 in
the geopotential can correspond to a temperature change in the ≈ 50 − 130 K
range, depending on assumed mean molecular weight. The geopotential and
horizontal winds are driven by the interacting processes of irradiation, rotation,
and advection.
The strength of stellar irradiation in our simulation can be analyzed through
two parameters: the radiative timescale τrad and the day/night temperature
contrast of the local radiative equilibrium. The radiative timescale is one of
the principal determiners of the resulting atmospheric patterns. Simulations
at the shortest radiative timescale tend to exhibit a hot spot on the day side
and cyclonic vortices on the nightside. The hot spot is shifted eastward from
the substellar longitude, similar to the equatorial superrotation on hot Jupiters.
High day-night heating contrasts result in standing, planetary-scale Rossby and
Kelvin waves as shown in Showman and Polvani (2011). As the timescale τrad
increases, the hot spot curves westward in the mid-latitudes under the influence
of the Coriolis force, breaking up into cyclonic vortices that equilibriate west of
the substellar point. As the radiative timescale τrad increases further to 10 days,
the behavior transitions from a day-to-night flow to a jet-dominated flow. Vor-
tices tend to be present for short and medium values of τrad but not for τrad = 10
days. Temporal variability is also primarily determined by τrad. Almost no vari-
ability is present for the shortest radiative timescales, where the radiative forc-
ing dominates all the other terms of the shallow-water system. Large temporal
variability occurs at longest radiative timescales, typically in the equatorial and
102
tropical latitudes.
The other parameter controlling the strength of irradiation is the prescribed
equilibrium day/night temperature contrast. The amplitude of the resulting
day-night thickness contrast appears to scale linearly with the fractional day-
night equilibrium difference ∆Φeq/Φ. This behavior makes sense as a response
to linear Newtonian relaxation. The fractional day-night difference controls the
amplitude of the dayside cosine bell of the prescribed local radiative equilib-
rium. Similarly, as we lower the scale height, lower mean geopotential Φ cor-
responds to lower resulting day-night thickness contrasts. The equator-to-pole
contrast decreases as the radiative timescale or the rotation period decrease.
Overall the response in our simulations to changes in radiative forcing is
consistent with the argument given in Showman et al. (2013). When radiative
forcing is weak (e.g. in the weak forcing regime or due to longer τrad), the heating
of air parcels at they cross from the day side to the night side will be too small
to induce significant day-night contrast, and the dominant driver of the flow
will be the meridional (latitudinal) gradient in the zonal-mean radiative heating.
As expected, we generally observe zonal symmetry with primary temperature
variations occuring between the poles and the equator.
The horizontal wind speeds in our simulations are controlled by radiative
forcing. In the strong and medium forcing regimes, the largest mean zonal wind
speeds tend to occur when both the radiative time scale and the rotation period
are short. In this case, the mean zonal wind speeds form a global-scale eastward
equatorial jet that moves with the speed of a few hundred ms−1. For longer
timescales τrad and Prot, the mean zonal wind speeds are U ≈ 100 ms−1. In the
weak forcing regime, the mean zonal wind speeds are generally lower and can
103
drop to U ≈ 10 ms−1. In the weak forcing regime, the mean zonal winds tend to
be lower when the planetary rotation period is short. We discuss this behavior
in more detail in Section 5.4.2.
We now turn to consider the balance of advective and radiative forces using
scale analysis. Here we compare the radiative timescale τrad and the advective
timescale τadv = L/U. For regimes where τrad = 0.1 days and τrad = 1 day, τrad <
τadv, and so the day-night differences tend to be more pronounced. However,
when τrad = 10, we tend to observe τrad > τadv, where the advective behavior
dominates exhibiting planetary-scale westward Rossby waves. These waves
are especially pronounced for regimes where Prot = τrad = 10 days. For the
intermediate scale τrad = 1 day, the strong forcing regime tends to have τrad >
τadv, while the weak forcing regime tends to have τrad < τadv, and the medium
forcing regime is transitional. This behavior is summarized in Figure 5.18. The
figure shows estimated contours of τrad = τadv based on the values for each of
the hot, warm, and cool regimes. The values were averaged over all reference
values of Φ.
The importance of rotation varies across our simulations. Overall, the strong
and medium forcing regimes tend to be more driven by insolation, while in the
weak forcing regime, rotation has a more pronounced effect. We quantify the
dominance of rotation using the Rossby number Ro = U/2ΩL, which gives the
ratio of advective to Coriolis forces in the horizontal momentum equation. U is
the characteristic windspeed and L is the characteristic horizontal length scale.
We use L ∼ a for a global-scale flow. Rapidly rotating planets with Prot = 1
day and planets with τrad = 10 days exhibit Ro ≪ 1 regardless of irradiation.
The highest Rossby numbers (∼ 2) are exhibited by slowly rotating planets with
104
short radiative timescales. When rotation period is 5 or 10 days and the radia-
tive timescale τrad = 1 day, irradiation plays more of a role in influencing the
Rossby number. Strongly forced planets tend to exhibit Ro ∼ 1, weakly forced
planets exhibit Ro ≪ 1, and planets with medium forcing exhibit a transition
between the two. As expected, the effects of rotation are more pronounced in
the atmospheres of rapidly rotating planets and planets with lower insolation.
Overall, the strong and medium forcing regimes are similar, while the weak
forcing regime exhibits lower Rossby numbers. Changes in scale height do not
affect the Rossby number significantly. Rossby numbers for a representative
scale height corresponding to Φ = 4 × 106 m2s−2 are represented in Table 5.1.
Due to low insolation, the dynamics of the weak forcing regime tend to be
qualitatively different from the dynamics of the strong and medium forcing
regimes at the same timescales. Low stellar irradiation results in lower wind
speeds and lower contrasts, and the effects of rotation lead to the formation of
more vortices. Instead of eastward superrotation, the weak forcing regimes tend
to exhibit westward shifts driven by Rossby waves.
Rossby radius of deformation is greater than the planetary radius for Prot = 5
days and Prot = 10 days. This predicts the formation of planetary-scale pat-
terns in the corresponding regimes. This prediction is confirmed by our sim-
ulations, which exhibit planetary-scale behavior in both geopotential and hori-
zontal wind patterns. For Prot = 1 day, the Rossby radius of deformation ranges
from 6.8 × 106 for Φ = 106 to 1.3 × 107 for Φ = 4 × 106. This corresponds to
the shorter widths of the equatorial jets, which extend up to 20◦ − 45◦ latitude.
Compared with our simulations, this estimate tends to hold for τrad = 0.1 days,
Prot = 1 day, but as the radiative timescale increases, the equatorial jet becomes
105
suppressed. We discuss this behavior in further detail in Section 5.4.2.
The Rhines length π (Ua/2Ω)1/2 is longer than the planetary radius is the
majority of simulations in our ensemble and shorter than the planetary radius
(∼ a/2) when longer radiative timescale τrad combines with shorter rotation pe-
riod Prot. Correspondingly, the simulations in this parameter regime tend to
exhibit banded jet-stream structures rather than vortices. This mechanism is
described in Showman et al. (2010).
5.4.2 Trends in the zonal wind flow
In our simulations, we have identified several important trends in the zonal
wind flow, which we discuss here. In the high forcing, high scale height regime,
the shape of the mean zonal winds depends strongly on the interaction of the
rotation rate and radiative timescale. For example, in the top row of Figure 5.5,
the shortest timescales regime (Prot = 1 day, τrad = 0.1 days) exhibits a strong
equatorial jet. However, as the radiative timescale increases to τrad = 1 day
and τrad = 10 days, the equatorial jet becomes suppressed and replaced by two
symmetrical jets with peaks at ≈ 40◦ latitude. In this regime, as τrad increases,
the forcing weakens. This in turn weakens the equatorial jet, which typically
emerges as a response to longitudinal variation in radiative forcing (e.g. Show-
man and Polvani, 2011). In these conditions, strong Coriolis effect results in the
emergence of symmetric tropical waves, similar to what is observed in idealized
models of Earth troposphere (e.g. Kraucunas and Hartmann, 2007).
The transition from a single equatorial jet to symmetric tropical jets appears
only in the presence of a high rotation rate. Note that in the middle row of Figure
106
5.5 (where the rotation period Prot = 5 days), the symmetric tropical jet pattern
is weaker and only present for the longest radiative timescale (τrad = 10 days).
In the bottom row (at rotation period Prot = 10 days), the equatorial jet weakens
but does not disappear as τrad increases. Our simulations suggest that in the sub-
Neptune regime, the rotaiton period must be sufficiently short and the radiative
timescale must be sufficiently long in order to see symmetric tropical jets.
This qualitative behavior persists across scale height variations in the high-
forcing and medium-forcing regimes. As scale height decreases, the overall
mean zonal wind speed decreases, but the transitions in qualitative behavior oc-
cur at the same values of τrad and Prot. The medium-forcing regime exhibits the
same qualitative behavior at all the simulated scale heights, except for τrad = 10
days and Prot = 5 days, where we observe a pattern which is neither a strong
equatorial jet nor symmetric tropical jets. Here a single westward jet has a flat
peak that extends across the equator, exhibiting an intermediate, transitional
regime. High- and medium-forcing regimes exhibit more similarities to each
other than to the weak-forcing regime.
Another important trend occurs at shortest timescales (τrad = 0.1 days,
Prot = 1 day). In the strong-forcing regime, when the scale height is large (at
Φ = 4 × 106 m2/s2), the simulation exhibits westward jets north and south of
the eastward equatorial jet. The peaks of the westward jets correspond to the
minima of mean zonal winds U, which occur at ≈ 50◦ latitude. However, as
scale height decreases, the minima move toward the equator. This behavior is a
response to the narrowing of the equatorial jet, since decreasing the s√cale height√
also decreases the Rossby radius of deformation (λR = gH/ f = Φ/ f ). For
Φ = 106 m2/s2, the minima occur at ≈ 35◦ latitude. Curiously, in this regime,
107
secondary westward jets emerge polarward of the eastward jets. The westward
jets are similar in speed to the eastward jets (≈ 35 m/s). In the medium-forcing
regime, the secondary jets are present at both Φ = 106 (≈ 35 m/s) and Φ = 2×106
(with speeds of ≈ 20 m/s) and are much slower than the eastward equatorial jet.
In contrast, in the weak-forcing regime, at shortest timescales, the secondary jets
are of the same order of magnitude as the eastward equatorial jet. Scale analy-
sis predicts that the lower (≈ 30 m/s) wind speeds in this regime (compared to
the short-timescales regime for strong and medium forcing) would lead to the
√
emergence of three jets (Njet = 2Ωa/U), which is consistent with our simula-
tions. The combination of fast rotation rate and weak forcing produces a zonal
wind pattern similar to that exhibited by the models of hot Jupiter WASP-43b in
presence of high drag (Kataria et al., 2015). Further, Zhang and Showman (2017)
conclude that sub-Neptunes with higher molecular weight atmospheres tend to
exhibit narrower equatorial jets. This is consistent with our simulations, where
high molecular weight correspond to lower scale heights and lower values for
the Rossby radius of deformation.
5.4.3 Idiosyncratic spin-up behavior for long timescales
When both τrad and Prot are long, the model exhibits peculiar spin-up behavior
driven by the westward Rossby wave. During spin-up, the hot spot can be
advected around the planet several times before the oscillations subside. This
behavior is illustrated in Figure 5.16 for Φ = 4 × 106 in a strong forcing regime,
but occurs for long radiative timescales for a variety of Φ in the strong and
medium forcing regimes. The left column of Figure 5.16 captures the Rossby
wave advecting the hotspot westward all the way around the planet with a
108
3000 h 5000 h
50
0
50
3150 h 5150 h
1.36 × 106
50 1.21 × 106
0 1.06 × 106
50 9.08 × 105
3300 h 5300 h 7.56 × 105
5
50 6.05 × 10
5
0 4.54 × 10
50 3.03 × 10
5
1.51 × 105
3450 h 5450 h 0.00 × 100
50
0
50
100 0 100 100 0 100
Figure 5.16: Snapshots of spin-up behavior for strong-forcing regime with Φ =
4 × 106 at rotation period Prot = 10 days and radiative timescale τrad = 10 days.
The hotspot is advected fully around the planet during spin-up (left column
shows cyclic behavior with a period of ≈ 450 hours). The right colum shows the
equilibrated behavior. While some oscillations are still present, the equilibrated
behavior shows much less time variability.
period of ∼ 450 hours. However, this oscillation is transient, and eventually
the hot spot is not exhibiting subrotation but is no longer advecting around the
plant. This behavior highlights the need for longer simulations, especially at
long timescales where time varibility is more likely.
109
5.4.4 Equatorial and Tropical Band Oscillation
For simulations with a short rotation period (Prot = 1 day) and long radiative
timescale (τrad = 10 days), we observe a persistent oscillation in the weak forcing
regime with a scale height corresponding toΦ = 4×104 m2s−2. Figure 5.17 shows
the snapshots from this regime: a hot equatorial band (top panel) alternates with
two hot tropical bands (bottom panel). The meridional component of the wind
vector field oscillates between north and south. The period of these oscillations
are ≈ 55 hours, or between two and three days. Note that Figure 5.17 shows
snapshots corresponding to the average geopotential profile shown in the top
right panel of Figure 5.12. The oscillation arises due to gravity waves forming
and propagating in the north-south direction, driving the mass of the buoyant
layer into alternating equatorial and tropical bands.
The behavior of a single hot equatorial band separating into two hot tropical
bands is unique among our simulations. When Prot = 1 days and τrad = 10 days,
simulations in the same weak-forcing regime at different scale heights as well as
planets with strong and medium forcing all show similar qualitative behavior,
but much less north-to-south variability in the meridional component of the
winds, and retain the hot equatorial band throughout.
A unique combination of fast rotation, weak radiative forcing, and large
scale height gives the opportunity for gravity waves to emerge. Firstly, a long
radiative timescale τrad = 10 days, combined with the weak-forcing regime leads
to the lowest radiative forcing among our simulations.√The shallow-water grav-√
ity wave speed is proportional to scale height gH = Φ, and so for Φ = 4×104
m2s−2, is the highest in our ensemble. For this regime, the high rotation rate
(Prot = 1 day) drives the momentum advection by mean-meridional circulation,
110
50 1.36 × 106
0 1.21 × 106
1.06 × 106
50 9.08 × 105
7.56 × 105
50 6.05 × 10
5
4.54 × 105
0 3.03 × 105
5
50 1.51 × 10
0.00 × 100
100 0 100
Figure 5.17: Snapshots of oscillatory behavior for a regime corresponding to a
weakly forced planet with a scale height corresponding to Φ = 4×106 at rotation
period Prot = 1 days and radiative timescale τrad = 10 days. Over the course of
the oscillation, a hot equatorial band (top panel) alternates with two hot tropical
bands (bottom panel) with a period of ≈ 55 hours.
[ ]
or v∗ f − ∂uy (see, e.g. Showman and Polvani, 2011). The interaction of the above∂
forces results in a persistent oscillation.
5.5 Conclusions
Decreasing the planetary radius from a hot Jupiter to a sub-Neptune leads to
more moderate geopotential contrasts, even at equivalent radiative forcing. This
phenomenon can partly be explained by the following mechanism. Reduced
planetary radius leads to shorter advective timescales: a parcel of gas moving at
the same zonal speed will move around the planet and pass through the dayside
more frequently on a smaller planet rather than a larger one. Shorter advective
timescales result in more uniform temperatures.
111
The strong and medium forcing regimes tend to exhibit similar behavior:
equatorial superrotation for short radiative timescales and jet-dominated flow
with little longitudinal variation for long radiative timescales. In the weak forc-
ing regime, due to weaker insolation, the resulting patterns are more driven by
the Coriolis force and are more likely to exhibit a westward shift of the hot spot
as well as a jet-dominated structure.
The transition from day-to-night flow to jet-dominated flow is tied predom-
inantly to the radiative timescale τrad, irrespective of rotation period Prot. Short
radiative timescales lead to a day-to-night flow, while long radiative timescales
lead to jet-dominated flow. The medium value of τrad = 1 day is transitional. For
short rotation periods, the intermediate value τrad = 1 day leads to jet-dominated
flow, but as rotation period increases, the regimes tend to exhibit more longitu-
dinal variation.
The geopotential contrast tends to scale linearly with the amplitude of the
prescribed local radiative equilibrium. The contrast is highest when both the
radiative timescale and the rotation period are short, and is lowest when both
are long.
Mean zonal wind speeds range from ∼ 600 ms−1 for the hot planet regime to
U < 100 ms−1 for the cool regime. In the low forcing regime, when the rotation
period is short, the mean zonal wind speeds can be as low ∼ 10 ms−1. Using
estimates for the wind speeds from our model, we compute the Rossby number
and Rhines scale. Both estimates indicate that we should expect global-scale
phenomena. This behavior is confirmed by our simulations.
Temporal variability is most prevalent when τrad and Prot are both long. When
112
Contours where rad = adv
101
eq/ = 0.1
100
eq/ = 0.5
rad > adv
eq/ = 1
rad < adv
10 1
100 101
Rotation period Prot [days]
Figure 5.18: The estimated contours where τadv = τrad for different day/night
contrasts of the local radiative equilibrium for the strong (solid line), medium
(dashed line), and weak (dot-dashed line) forcing regimes. The blue crosses
mark the simulated regimes. In order to compute the contours, the values for
τrad and τadv were averaged across the four scale heights in our ensemble.
τrad is short, we see almost no variation, even with a long rotation period. With
an eye toward observations, τrad might be difficult to determine a priori, how-
ever, planets at lower rotation periods are more likely to have short radiative
timescales less temporal variability and can therefore be more favorable candi-
dates for observation using current facilities.
113
Radiative timescale rad [days]
APPENDIX A
ALGORITHM FOR GENERATING SCALE-FREE NETWORKS WITH
HOMOPHILY
We generate a scale-free network with homophily using the following algo-
rithm:
1. Fix n nodes.
2. Draw degrees from a power law distribution.
3. Generate a vector of length n assigning a binary opinion: 0 to majority
nodes and 1 to minority nodes.
4. Initialize two empty stacks: the majority stack and the minority stack.
5. For each node:
If a node is a minority node, add its index to the minority stack the num-
ber of times corresponding to its degree.
If the node is a majority node, add its index to the majority stack the num-
ber of times corresponding to its degree.
6. Shuffle the majority and minority stacks.
7. While the minority stack is non-empty:
pop node1 from the top of the minority stack. generate a random number
between 0 and 1. If the random number is less than the homophily factor
h, draw an edge between node1 and the first node in the minority stack
(node2). If the random number is greater, draw an edge between node1
and the top node in the majority stack.
8. If the majority stack is nonempty by the time the minority stack is empty,
connect the remaining majority nodes in pairs.
114
APPENDIX B
COMPARISON OF THE ANALYTICAL EQUILIBRIUM SOLUTIONS
In this section, we discuss a brief comparison of the analytical equilibrium
solutions. As mentioned in Section 3.2.1, North (1975a) comments that the di-
vergence operator and the relaxation to the mean operator have the same effect
on degree two polynomials when δ = µ/6. These operators are the same on the
function subspace of degree two polynomials; however, solutions to the models
live in larger spaces of functions, sometimes called solution spaces. In particu-
lar, the equilibrium temperature solutions must be in the solution space so the
relaxation model solution space contains piecewise polynomials of degree 2N
and the diffusion model solution space contains bounded differentiable hyper-
geometric functions on the interval [0,1]. Degree two polynomials on [0,1] are a
subset of both spaces. Even though solutions to the models are not restricted to
the function subspace of degree two polynomials, the relation between δ and µ
of δ = µ/6 remains a remarkably good approximation for comparing solutions
between the models. We find that in the small and large limits of the heat trans-
port efficiency parameter with other parameters fixed, the equilibrium temper-
ature solutions from both models approach the same functions.
As µ and δ approach zero, solutions of the relaxation to the mean model
approach solutions to the diffusion model. This can be seen directly by con-
sidering Equations (3.4) and (3.5). As the parameters µ and δ approach zero,
heat transport becomes negligible and the models become an energy balance
equation where the latitude is simply a parameter. Convergence of the solu-
tions across all latitudes as µ and δ approach 0 appears to be rather slow due
to the continuous diffusive solution approaching the jump discontinuity in the
115
RTM Caps 0.05,0.24 RTM Caps 0.44,0.24
11.0 1.0
0.80.8 0.8
0.60.6 0.6
0.40.4 0.4
↵ = 0.05 ↵ = 0.44
µ/6 = 0.04 µ/6 = 0.04
0.20.2 0.2
0. 0.0 0.0
0.96 0.98 1.00 1.02 1.04 1.06 1.0 1.2 1.4 1.6 1.8
0.96 1.00q 1.04 1.0 1.4 1.8q
RTM Belt 0.05,0.24 RTM Belt 0.44,0.24
11.0 1.0
0.80.8 0.8
0.60.6 0.6
0.40.4 0.4
↵ = 0.05 ↵ = 0.44
µ/6 = 0.04 µ/6 = 0.04
0.20.2 0.2
0. 0.0 0.0
0.96 0.98 1.00 1.02 1.04 1.06 1.0 1.2 1.4 1.6 1.8
0.96 1.00
qq
1.04 1.0 1.4 1.8
q q
Figure 3: Figure 4 supplementary version - relaxation to the mean caps (top) and belt (bottom)
1.4  = µ/6 = 0.0001 1.4  = µ/6 = 0.01
1.2 1.2
1.0 1.0
0.8 0.8
0.6 0.6
0.4 0.4
0.2 0.2
0.0 0.0
0.0 0.2 0.4 0.6 0.8 1.0 0.0 0.2 0.4 0.6 0.8 1.0
1.4  = µ/6 = 0.04 1.4  = µ/6 = 0.31
1.2 1.2
1.0 1.0
0.8 0.8
0.6 0.6
0.4 0.4
0.2 0.2
0.0 0.0
0.0 0.2 0.4 0.6 0.8 1.0 0.0 0.2 0.4 0.6 0.8 1.0
1.4  = µ/6 = 1 1.4  = µ/6 = 10
1.2 1.2
1.0 1.0
0.8 0.8
0.6 0.6
0.4 0.4
0.2 0.2
0.0 0.0
0.0 0.2 0.4 0.6 0.8 1.0 0.0 0.2 0.4 0.6 0.8 1.0
⌘ ⌘
FiguFreig4u: re B.1: Comparison of equilibrium temperaturComparison of equilibrium temperature profiles T ⇤e profiles T
∗ for the diffusion
for the di↵ueqsion model (black dashed
model (black dashed line) and the relaxation to theeqmean model (gray solid line).
line) and the relaxation to the mean model (gray solid line). The plots are for Earth’s albedo
The plots are for Earth’s albedo contrast and obliquity, N = 3, ice line latitude
contrast and obliquity, N = 3, ice line latitude ⌘ = 0.1.
η = 0.1.
relaxation to the mean solution at the ic2e line.
This behavior is illustrated in Figure B.1, wherein we compare the equilib-
rium temperature profile solutions of the two models with N = 3 and Earth
values for all parameters except δ and µ.
As µ and δ become large, solutions of the relaxation to the mean model again
approach solutions to the diffusion model. Intuitively this occurs because as
the transport coefficients become large, heat is redistributed almost instanta-
neously, and the solutions approach the global mean temperature (see Equation
(3.7)). Numerically we see that the solutions approach each other faster than
they approach the global mean. Convergence for large µ and δ appears to be
much faster compared to when the parameters approach zero. For instance, for
116
⌘ ⌘
T ⇤ T ⇤ T ⇤ η η
η η
the parameters used in Figure B.1, the solutions are within 3% of each other for
all latitudes when δ = µ/6 = 1 and within 1% of each other for all latitudes when
δ = µ/6 = 10. Further mathematical investigation of these properties is a natural
direction for future work.
117
APPENDIX C
AN OPEN-SOURCE SOFTWARE REPOSITORY FOR SWAMPE
I have created a repository for SWAMPE that is publicly available on
Github at https://github.com/kathlandgren/SWAMPE and documented
at https://swampe.readthedocs.io. The Python package has been de-
veloped in accordance with the standards for open-source research software, as
defined by the Journal of Open-Source Software1. SWAMPE has also been up-
loaded to the Python Package Index (PyPI) and is available for download and
installation.
C.1 Documentation
The documentation contains installation instructions, tutorial notebooks, and
API reference, i.e., all the modules and their functions, along with the input pa-
rameters they require and the output they generate. All functions in SWAMPE
have been documented extensively following the latest Python docstring con-
vention.2 The homepage of the documentation website is shown in Figure C.1.
C.2 Time-stepping
While Hack and Jakob (1992) describe a semi-implicit time-stepping scheme,
we follow Langton (2008b) and implement a modified Euler’s method scheme.
In practice, this method proves more stable, especially for strongly irradiated
planets. In this section, we derive the modified Euler formulation.
1https://joss.theoj.org/
2https://peps.python.org/pep-0257/
118
119
Figure C.1: Screenshot of SWAMPE documentation’s homepage.
The modified Euler’s method for differential equations has the form:
yn+1 = yn + (∆t/2)[ f (tn, yn) + f (tn+1, yn + ∆t f (tn, yn)]. (C.1)
Equivalently, we can write:
K1 = ∆t f (tn, yn) (C.2)
K2 = ∆t f (tn+1, yn + K1) (C.3)
K1 + Ky y 2n+1 = n + . (C.4)2
We will now derive the coefficients for our timestepping scheme.
Firstly, we can write the unforced shallow-water system to split it into the
linear and nonlinear components as follows:
       ηm  
n  0 0 0d  
 
 m   n(n+1)
 
 η
m
 
n
 
  E (t) 
dt δn  = 0 0 a2  δ
m   n
Φm 0 −Φ 0 Φm
+ D(t)  (C.5)
n n P(t)

 
E (t)
 

where 

D(t)  represents the nonlinear time-dependent components. We willP(t)
evaluate the first component of the right hand side implicitly, while evaluating
the second component explicitly.
The (unforced) nonlinear components can be expressed as follows:
1 ∂A 1 ∂B
E (t) = − − , (C.6)
a(1 − µ2) ∂λ a ∂µ
1 ∂B 1 ∂A
D(t) = 2− 2 − − ∇ E, (C.7)a(1 µ ) ∂λ a ∂µ
120
− 1 ∂C 1 ∂DP(t) = − 2 − . (C.8)a(1 µ ) ∂λ a ∂µ
Let FΦ be the geopotential forcing (for SWAMPE, due to stellar irradiation,
but more general in theory). Let FU = Fu cos ϕ and FV = Fv cos ϕ be momentum
forcing. Then the (forced) nonlinear components are as follows:
− 1 ∂ − − 1 ∂E (t) = − 2 (A Fa(1 ) V) (B + F ), (C.9)µ ∂λ a U∂µ
1 ∂ − 1 ∂D(t) = (B + F ) (A − F ) − ∇2− 2 U V E, (C.10)a(1 µ ) ∂λ a ∂µ
1 ∂C 1 ∂D
P(t) = − − 2 − + FΦ. (C.11)a(1 µ ) ∂λ a ∂µ
Following the notation of the modified Euler’s method, we write K1 =
∆t f (t, yt):
K1η = ∆t(E (t)), (C.12)
( )
1 n(n + 1)K = ∆t Φm(t)δ 2 n +D(t) , (C.13)a
( )
K1Φ = ∆t −Φδm(t)n +P(t) . (C.14)
Then we can write the K2 = ∆t( f (t + 1, yt + K1)) coefficients.
K2η = ∆t(E (t + 1)), (C.15)
121
( )
2 n(n + 1)Kδ = ∆t D(t + 1) + 2 (Φ
m
n + K
1
a Φ
) , (C.16)
( )
K2Φ = ∆t P(t + 1) − Φ(δmn + K1δ ) . (C.17)
Expanding( the equations for K2 and K2δ , we obtain:Φ )
2 n(n + 1) n(n + 1) m − n(n + 1)K mδ = ∆t D(t + 1) + 2 (P(t)) +a a2 Φn Φ a2 δn , (C.18)
( )
n(n + 1)
K2Φ = ∆t P(t + 1) − Φ(D(t)) − Φδmn − Φ ma2 Φn . (C.19)
We evaluate the time-dependentterm s explicitly, assuming   

 E (t)  
  E (t + 1) 
D(t) 
 = 
D(t + 1) 

P(t) P(t + 1)
(C.20)
to first order. This is what is done in the semi-implicit method in Hack and
Jakob (1992). An alternative variant would be to approximate η, δ, Φ, U, and V
by a different method, such as forward Euler’s method or a semi-implicit one.
This would result in a higher computational cost and hopefully higher accuracy
as well, while maintaining the stability properties of modified Euler’s method.
Note that in the current implementation, η time-stepping is equivalent to for-
ward Euler’s method, since η does not depend linearly on other state variables,
only nonlinearly in the E (t) term. Writing (K1 + K2)/2 in order to evaluate the
modified Euler schem( e as in (C.4), we can s(implify:
K1
))
δ + K
2
δ n(n + 1) m 1 n(n + 1)= ∆t 2 Φn +D(t) + 2 (P(t) − Φδ
m , (C.21)
2 a 2 a n
122
and
K1 + K2 ( ) ( )
Φ Φ ∆t n(n + 1)
= ∆t −Φδmn +P(t) − Φ D(t) + . (C.22)2 2 a2
C.3 Filters
To ensure numerical stability, SWAMPE applies two filters at each time step.
The first is a modal-splitting filter as described in Hack and Jakob (1992). The
second is a sixth-degree hyperviscosity filter. We use the formulation based on
Gelb and Gleeson (2001). Either filter can be turned off, and its coefficients can
be changed by the user.
123
APPENDIX D
SUPPLEMENTARY FIGURES SHOWING GEOPOTENTIAL CONTOURS
AND MEAN-ZONAL WINDS FOR INTERMEDIATE SCALE HEIGHTS.
In this appendix, we present the supplementary figures showing geopotential
contours and mean-zonal winds for intermediate scale heights corresponding
to reference geopotential Φ = 3 × 106 m2s−2 and Φ = 2 × 106 m2s−2.
124
eq/ = 1, = 3x10^6
rad [days]
0.1 1 10
50
1.43 × 106
0
1.27 × 106
50
1.11 × 106
50 9.53 × 105
0 7.94 × 105
50 6.35 × 105
4.77 × 105
50 3.18 × 105
0 1.59 × 105
50
0.00 × 100
100 0 100 100 0 100 100 0 100
Figure D.1: Averaged maps of geopotential contours and wind vector fields
for the equilibrated solutions of the shallow-water model for the strong forcing
regime (∆Φeq/Φ = 1) at reference geopotential Φ = 3 × 106 m2s−2. The constant
value of Φ has been subtracted from each panel. Each panel in the grid corre-
sponds to a unique combination of planetary rotation period Prot and radiative
timescale τrad. Panels share the same color scale for geopotential, but winds
are normalized independently in each panel. The substellar point is located at
(0◦, 0◦) for each panel.
125
Prot [days]
10 5 1
eq/ = 1, = 3x10^6
rad = 0.1 rad = 1 rad = 10
50
0
50
50
0
50
50
0
50
0 200 400 600 0 200 400 600 0 200 400 600
Figure D.2: Mean-zonal winds U for the equilibrated solutions of the shallow-
water model for the strong forcing regime (∆Φeq/Φ = 1) at reference geopotential
Φ = 3 × 106 m2s−2. Each panel in the grid corresponds to a unique combination
of planetary rotation period Prot and radiative timescale τrad.
126
Prot = 10 Prot = 5 Prot = 1
eq/ = 1, = 2x10^6
rad [days]
0.1 1 10
50
1.05 × 106
0
9.33 × 105
50
8.16 × 105
50 7.00 × 105
0 5.83 × 105
50 4.66 × 105
3.50 × 105
50 2.33 × 105
0 1.17 × 105
50
0.00 × 100
100 0 100 100 0 100 100 0 100
Figure D.3: Averaged maps of geopotential contours and wind vector fields
for the equilibrated solutions of the shallow-water model for the strong forcing
regime (∆Φeq/Φ = 1) at reference geopotential Φ = 2 × 106 m2s−2. The constant
value of Φ has been subtracted from each panel. Each panel in the grid corre-
sponds to a unique combination of planetary rotation period Prot and radiative
timescale τrad. Panels share the same color scale for geopotential, but winds
are normalized independently in each panel. The substellar point is located at
(0◦, 0◦) for each panel.
127
Prot [days]
10 5 1
eq/ = 1, = 2x10^6
rad = 0.1 rad = 1 rad = 10
50
0
50
50
0
50
50
0
50
0 200 400 0 200 400 0 200 400
Figure D.4: Mean-zonal winds U for the equilibrated solutions of the shallow-
water model for the strong forcing regime (∆Φeq/Φ = 1) at reference geopotential
Φ = 2 × 106 m2s−2. Each panel in the grid corresponds to a unique combination
of planetary rotation period Prot and radiative timescale τrad.
128
Prot = 10 Prot = 5 Prot = 1
eq/ = 0.5, = 3x10^6
rad [days]
0.1 1 10
50
7.61 × 105
0
6.77 × 105
50
5.92 × 105
50 5.08 × 105
0 4.23 × 105
50 3.38 × 105
2.54 × 105
50 1.69 × 105
0 8.46 × 104
50
0.00 × 100
100 0 100 100 0 100 100 0 100
Figure D.5: Averaged maps of geopotential contours and wind vector fields for
the equilibrated solutions of the shallow-water model for the medium forcing
regime (∆Φeq/Φ = 0.5) at reference geopotential Φ = 3 × 106 m2s−2. The constant
value of Φ has been subtracted from each panel. Each panel in the grid corre-
sponds to a unique combination of planetary rotation period Prot and radiative
timescale τrad. Panels share the same color scale for geopotential, but winds
are normalized independently in each panel. The substellar point is located at
(0◦, 0◦) for each panel.
129
Prot [days]
10 5 1
eq/ = 0.5, = 3x10^6
rad = 0.1 rad = 1 rad = 10
50
0
50
50
0
50
50
0
50
100 0 100 200 300 400 100 0 100 200 300 400 100 0 100 200 300 400
Figure D.6: Mean-zonal winds U for the equilibrated solutions of the shallow-
water model for the medium forcing regime (∆Φeq/Φ = 0.5) at reference geopo-
tential Φ = 3× 106 m2s−2. Each panel in the grid corresponds to a unique combi-
nation of planetary rotation period Prot and radiative timescale τrad.
130
Prot = 10 Prot = 5 Prot = 1
eq/ = 0.5, = 2x10^6
rad [days]
0.1 1 10
50
5.24 × 105
0
4.66 × 105
50
4.07 × 105
50 3.49 × 105
0 2.91 × 105
50 2.33 × 105
1.75 × 105
50 1.16 × 105
0 5.82 × 104
50
0.00 × 100
100 0 100 100 0 100 100 0 100
Figure D.7: Averaged maps of geopotential contours and wind vector fields for
the equilibrated solutions of the shallow-water model for the medium forcing
regime (∆Φeq/Φ = 0.5) at reference geopotential Φ = 2 × 106 m2s−2. The constant
value of Φ has been subtracted from each panel. Each panel in the grid corre-
sponds to a unique combination of planetary rotation period Prot and radiative
timescale τrad. Panels share the same color scale for geopotential, but winds
are normalized independently in each panel. The substellar point is located at
(0◦, 0◦) for each panel.
131
Prot [days]
10 5 1
eq/ = 0.5, = 2x10^6
rad = 0.1 rad = 1 rad = 10
50
0
50
50
0
50
50
0
50
0 100 200 300 0 100 200 300 0 100 200 300
Figure D.8: Mean-zonal winds U for the equilibrated solutions of the shallow-
water model for the medium forcing regime (∆Φeq/Φ = 0.5) at reference geopo-
tential Φ = 2× 106 m2s−2. Each panel in the grid corresponds to a unique combi-
nation of planetary rotation period Prot and radiative timescale τrad
132
Prot = 10 Prot = 5 Prot = 1
eq/ = 0.1, = 3x10^6
rad [days]
0.1 1 10
50
1.39 × 105
0
1.24 × 105
50
1.08 × 105
50 9.28 × 104
0 7.73 × 104
50 6.19 × 104
4.64 × 104
50 3.09 × 104
0 1.55 × 104
50
0.00 × 100
100 0 100 100 0 100 100 0 100
Figure D.9: Averaged maps of geopotential contours and wind vector fields
for the equilibrated solutions of the shallow-water model for the weak forcing
regime (∆Φeq/Φ = 0.1) at reference geopotential Φ = 3 × 106 m2s−2. The constant
value of Φ has been subtracted from each panel. Each panel in the grid corre-
sponds to a unique combination of planetary rotation period Prot and radiative
timescale τrad. Panels share the same color scale for geopotential, but winds
are normalized independently in each panel. The substellar point is located at
(0◦, 0◦) for each panel.
133
Prot [days]
10 5 1
eq/ = 0.1, = 3x10^6
rad = 0.1 rad = 1 rad = 10
50
0
50
50
0
50
50
0
50
25 0 25 50 75 100 25 0 25 50 75 100 25 0 25 50 75 100
Figure D.10: Mean-zonal winds U for the equilibrated solutions of the shallow-
water model for the weak forcing regime (∆Φeq/Φ = 0.1) at reference geopoten-
tial Φ = 3 × 106 m2s−2. Each panel in the grid corresponds to a unique combina-
tion of planetary rotation period Prot and radiative timescale τrad
134
Prot = 10 Prot = 5 Prot = 1
eq/ = 0.1, = 2x10^6
rad [days]
0.1 1 10
50
1.05 × 105
0
9.35 × 104
50
8.18 × 104
50 7.02 × 104
0 5.85 × 104
50 4.68 × 104
3.51 × 104
50 2.34 × 104
0 1.17 × 104
50
0.00 × 100
100 0 100 100 0 100 100 0 100
Figure D.11: Averaged maps of geopotential contours and wind vector fields
for the equilibrated solutions of the shallow-water model for the weak forcing
regime (∆Φeq/Φ = 0.1) at reference geopotential Φ = 2 × 106 m2s−2. The constant
value of Φ has been subtracted from each panel. Each panel in the grid corre-
sponds to a unique combination of planetary rotation period Prot and radiative
timescale τrad. Panels share the same color scale for geopotential, but winds
are normalized independently in each panel. The substellar point is located at
(0◦, 0◦) for each panel.
135
Prot [days]
10 5 1
eq/ = 0.1, = 2x10^6
rad = 0.1 rad = 1 rad = 10
50
0
50
50
0
50
50
0
50
25 0 25 50 75 25 0 25 50 75 25 0 25 50 75
Figure D.12: Mean-zonal winds U for the equilibrated solutions of the shallow-
water model for the weak forcing regime (∆Φeq/Φ = 0.1) at reference geopoten-
tial Φ = 2 × 106 m2s−2. Each panel in the grid corresponds to a unique combina-
tion of planetary rotation period Prot and radiative timescale τrad.
136
Prot = 10 Prot = 5 Prot = 1
BIBLIOGRAPHY
Réka Albert and Albert-László Barabási. Statistical mechanics of complex net-
works. Reviews of Modern Physics, 74(1):47, 2002.
JC Armstrong, R Barnes, S Domagal-Goldman, J Breiner, TR Quinn, and
VS Meadows. Effects of extreme obliquity variations on the habitability of
exoplanets. Astrobiology, 14(4):277–291, 2014.
Fernando J Ballesteros and Bartolo Luque. Walking on exoplanets: Is star wars
right?, 2016.
Albert-László Barabási and Réka Albert. Emergence of scaling in random net-
works. Science, 286(5439):509–512, 1999.
Rory Barnes. Tidal locking of habitable exoplanets. Celestial Mechanics
and Dynamical Astronomy, 129(4):509–536, December 2017. doi: 10.1007/
s10569-017-9783-7.
Anna M Barry, Esther Widiasih, and Richard McGehee. Nonsmooth frame-
works for an extended budyko model. Discrete and Continuous Dynamical Sys-
tems Series B, 22(6):2447–2463, 2017. doi: 10.3934/dcdsb.2017125.
Fabian Baumann, Philipp Lorenz-Spreen, Igor M Sokolov, and Michele Starnini.
Modeling echo chambers and polarization dynamics in social networks. Phys-
ical Review Letters, 124(4):048301, 2020.
Jacob L Bean, Sean N Raymond, and James E Owen. The nature and origins
of sub-neptune size planets. Journal of Geophysical Research: Planets, 126(1):
e2020JE006639, 2021.
137
Taylor J Bell and Nicolas B Cowan. Increased heat transport in ultra-hot jupiter
atmospheres through h2 dissociation and recombination. The Astrophysical
Journal Letters, 857(2):L20, 2018.
Marián Boguná, Romualdo Pastor-Satorras, and Alessandro Vespignani. Cut-
offs and finite size effects in scale-free networks. The European Physical Journal
B, 38(2):205–209, 2004.
John P Boyd. Chebyshev and Fourier spectral methods. Courier Corporation, 2001.
Anna D Broido and Aaron Clauset. Scale-free networks are rare. Nature Com-
munications, 10(1):1–10, 2019.
Shawn R Brueshaber, Kunio M Sayanagi, and Timothy E Dowling. Dynamical
regimes of giant planet polar vortices. Icarus, 323:46–61, 2019.
Mikhail I Budyko. The effect of solar radiation variations on the climate of the
earth. Tellus, 21(5):611–619, 1969.
Leonardo Bursztyn, Davide Cantoni, Patricia Funk, and Noam Yuchtman. Polls,
the press, and political participation: The effects of anticipated election close-
ness on voter turnout. 2017.
Robert F Cahalan and Gerald R North. A stability theorem for energy-balance
climate models. Journal of the Atmospheric Sciences, 36(7):1178–1188, 1979.
João Cancela and Benny Geys. Explaining voter turnout: A meta-analysis of
national and subnational elections. Electoral Studies, 42:264–275, 2016.
Michele Catanzaro, Marián Boguná, and Romualdo Pastor-Satorras. Generation
of uncorrelated random scale-free networks. Physical Review E, 71(2):027103,
2005.
138
Jade Checlair, Kristen Menou, and Dorian S Abbot. No snowball on habitable
tidally locked planets. The Astrophysical Journal, 845(2):132, 2017.
Jade H Checlair, Stephanie L Olson, Malte F Jansen, and Dorian S Abbot. No
snowball on habitable tidally locked planets with a dynamic ocean. The As-
trophysical Journal Letters, 884(2):L46, 2019a.
Jade H Checlair, Andrea M Salazar, Adiv Paradise, Kristen Menou, and Do-
rian S Abbot. No snowball cycles at the outer edge of the habitable zone for
habitable tidally locked planets. The Astrophysical Journal Letters, 887(1):L3,
2019b.
D A Christie, N J Mayne, R M Gillard, J Manners, E Hébrard, S Lines, and K Ko-
hary. The impact of phase equilibrium cloud models on GCM simulations of
GJ 1214b. Monthly Notices of the Royal Astronomical Society, 09 2022. ISSN 0035-
8711. doi: 10.1093/mnras/stac2763. URL https://doi.org/10.1093/
mnras/stac2763. stac2763.
Aaron Clauset, Cosma Rohilla Shalizi, and Mark E J Newman. Power-law dis-
tributions in empirical data. SIAM Review, 51(4):661–703, 2009.
Curtis S Cooper and Adam P Showman. Dynamics and disequilibrium carbon
chemistry in hot jupiter atmospheres, with application to hd 209458b. The
Astrophysical Journal, 649(2):1048, 2006.
Ian JM Crossfield and Laura Kreidberg. Trends in atmospheric properties of
neptune-size exoplanets. The Astronomical Journal, 154(6):261, 2017.
Giordano De Marzo, Andrea Zaccaria, and Claudio Castellano. Emergence of
polarization in a voter model with personalized information. Physical Review
Research, 2(4):043117, 2020.
139
Hans J. Deeg and Juan Antonio Belmonte. Handbook of Exoplanets. 2018. doi:
10.1007/978-3-319-55333-7.
Michela Del Vicario, Antonio Scala, Guido Caldarelli, H Eugene Stanley, and
Walter Quattrociocchi. Modeling confirmation bias and polarization. Scientific
Reports, 7(1):1–9, 2017.
Ian Dobbs-Dixon and Eric Agol. Three-dimensional radiative-hydrodynamical
simulations of the highly irradiated short-period exoplanet hd 189733b.
Monthly Notices of the Royal Astronomical Society, 435(4):3159–3168, 2013.
Anthony R Dobrovolskis. On a planet’s annual mean insolation. Icarus, 363:
114297, 2021.
Leonardo A. Dos Santos, Aline A. Vidotto, Shreyas Vissapragada, Munazza K.
Alam, Romain Allart, Vincent Bourrier, James Kirk, Julia V. Seidel, and
David Ehrenreich. p-winds: An open-source Python code to model plane-
tary outflows and upper atmospheres. , 659:A62, March 2022. doi: 10.1051/
0004-6361/202142038.
James W Endersby, Steven E Galatas, and Chapman B Rackaway. Closeness
counts in canada: Voter participation in the 1993 and 1997 federal elections.
The Journal of Politics, 64(2):610–631, 2002.
Y Katherina Feng, Michael R Line, Jonathan J Fortney, Kevin B Stevenson, Jacob
Bean, Laura Kreidberg, and Vivien Parmentier. The impact of non-uniform
thermal structure on the interpretation of exoplanet emission spectra. The
Astrophysical Journal, 829(1):52, 2016.
Raffaele Ferrari and David Ferreira. What processes drive the ocean heat trans-
port? Ocean Modelling, 38(3-4):171–186, 2011.
140
David Ferreira, John Marshall, Paul A O’Gorman, and Sara Seager. Climate at
high-obliquity. Icarus, 243:236–248, 2014.
Serge Galam. Sociophysics: A review of galam models. International Journal of
Modern Physics C, 19(03):409–440, 2008.
Anne Gelb and James P Gleeson. Spectral viscosity for shallow water equations
in spherical geometry. Monthly Weather Review, 129(9):2346–2360, 2001.
Alan Gerber, Mitchell Hoffman, John Morgan, and Collin Raymond. One in a
million: Field experiments on perceived closeness of the election and voter
turnout. American Economic Journal: Applied Economics, 12(3):287–325, 2020.
Benny Geys. Explaining voter turnout: A review of aggregate-level research.
Electoral Studies, 25(4):637–663, 2006.
Jayesh M. Goyal, Hannah R. Wakeford, Nathan J. Mayne, Nikole K. Lewis, Ben-
jamin Drummond, and David K. Sing. Fully scalable forward model grid
of exoplanet transmission spectra. , 482(4):4503–4513, February 2019. doi:
10.1093/mnras/sty3001.
James J Hack and Ruediger Jakob. Description of a global shallow water model
based on the spectral transform method. NCAR, 1992.
Mark Hammond and Raymond T. Pierrehumbert. Wave-mean Flow Interac-
tions in the Atmospheric Circulation of Tidally Locked Planets. , 869(1):65,
December 2018. doi: 10.3847/1538-4357/aaec03.
Tsahi Hayat and Tal Samuel-Azran. “you too, second screeners?” second screen-
ers’ echo chambers during the 2016 us elections primaries. Journal of Broad-
casting & Electronic Media, 61(2):291–308, 2017.
141
Rainer Hegselmann and Ulrich Krause. Opinion dynamics and bounded confi-
dence models, analysis, and simulation. Journal of Artificial Societies and Social
Simulation, 5(3), 2002.
Fritz Heider. Attitudes and cognitive organization. The Journal of Psychology, 21
(1):107–112, 1946.
Isaac M Held and Max J Suarez. Simple albedo feedback models of the icecaps.
Tellus, 26(6):613–629, 1974.
Kevin Heng, Dargan M. W. Frierson, and Peter J. Phillipps. Atmospheric circula-
tion of tidally locked exoplanets: II. Dual-band radiative transfer and convec-
tive adjustment. , 418(4):2669–2696, December 2011. doi: 10.1111/j.1365-2966.
2011.19658.x.
Abigail Hickok, Yacoub Kureh, Heather Z Brooks, Michelle Feng, and Mason A
Porter. A bounded-confidence model of opinion dynamics on hypergraphs.
arXiv preprint arXiv:2102.06825, 2021.
Richard A Holley and Thomas M Liggett. Ergodic theorems for weakly inter-
acting infinite systems and the voter model. The Annals of Probability, pages
643–663, 1975.
Ross Joseph Iaci. The gamma distribution as an alternative to the lognormal distribu-
tion in environmental applications. University of Nevada, Las Vegas, 2000.
Iacopo Iacopini, Giovanni Petri, Andrea Baronchelli, and Alain Barrat. Vanish-
ing size of critical mass for tipping points in social convention. arXiv preprint
arXiv:2103.10411, 2021.
Shin-ichi Iga and Yoshihisa Matsuda. Shear instability in a shallow water model
142
with implications for the venus atmosphere. Journal of the atmospheric sciences,
62(7):2514–2527, 2005.
Jonas S Juul and Mason A Porter. Hipsters on networks: How a minority group
of individuals can lead to an antiestablishment majority. Physical Review E, 99
(2):022313, 2019.
Stephen R Kane, Zhexing Li, Eric T Wolf, Colby Ostberg, and Michelle L Hill.
Eccentricity driven climate effects in the kepler-1649 system. The Astronomical
Journal, 161(1):31, 2020.
Hans Kaper and Hans Engler. Mathematics and climate. SIAM, 2013.
Tiffany Kataria, Adam P Showman, Jonathan J Fortney, Kevin B Stevenson,
Michael R Line, Laura Kreidberg, Jacob L Bean, and Jean-Michel Désert.
The atmospheric circulation of the hot jupiter wasp-43b: Comparing three-
dimensional models to spectrophotometric data. The Astrophysical Journal, 801
(2):86, 2015.
Tiffany Kataria, David K. Sing, Nikole K. Lewis, Channon Visscher, Adam P.
Showman, Jonathan J. Fortney, and Mark S. Marley. THE ATMOSPHERIC
CIRCULATION OF a NINE-HOT-JUPITER SAMPLE: PROBING CIRCULA-
TION AND CHEMISTRY OVER a WIDE PHASE SPACE. The Astrophys-
ical Journal, 821(1):9, apr 2016a. doi: 10.3847/0004-637x/821/1/9. URL
https://doi.org/10.3847/0004-637x/821/1/9.
Tiffany Kataria, David K Sing, Nikole K Lewis, Channon Visscher, Adam P
Showman, Jonathan J Fortney, and Mark S Marley. The atmospheric circu-
lation of a nine-hot-jupiter sample: Probing circulation and chemistry over a
wide phase space. The Astrophysical Journal, 821(1):9, 2016b.
143
Edwin S. Kite, Jr. Fegley, Bruce, Laura Schaefer, and Eric B. Ford. Atmosphere
Origins for Exoplanet Sub-Neptunes. , 891(2):111, March 2020. doi: 10.3847/
1538-4357/ab6ffb.
Paul L Krapivsky and Sidney Redner. Dynamics of majority rule in two-state
interacting spin systems. Physical Review Letters, 90(23):238701, 2003.
Ian Kraucunas and Dennis L Hartmann. Tropical stationary waves in a nonlin-
ear shallow-water model with realistic basic states. Journal of the atmospheric
sciences, 64(7):2540–2557, 2007.
Krzysztof Kułakowski, Przemysław Gawroński, and Piotr Gronek. The heider
balance: A continuous approach. International Journal of Modern Physics C, 16
(05):707–716, 2005.
Renaud Lambiotte and Sidney Redner. Dynamics of vacillating voters. Journal
of Statistical Mechanics: Theory and Experiment, 2007(10):L10001, 2007.
Renaud Lambiotte and Sidney Redner. Dynamics of non-conservative voters.
EPL (Europhysics Letters), 82(1):18007, 2008.
Jonathan Langton. Atmospheric dynamics on strongly irradiated Jovian planets. Uni-
versity of California, Santa Cruz, 2008a.
Jonathan Langton. Atmospheric dynamics on strongly ir-
radiated Jovian planets. PhD thesis, 2008b. URL
https://www.proquest.com/dissertations-theses/
atmospheric-dynamics-on-strongly-irradiated/docview/
304661389/se-2.
Jonathan Langton and Gregory Laughlin. Hydrodynamic simulations of un-
evenly irradiated jovian planets. The Astrophysical Journal, 674(2):1106, 2008.
144
Jérémy Leconte, Hanbo Wu, Kristen Menou, and Norman Murray. Asyn-
chronous rotation of Earth-mass planets in the habitable zone of lower-
mass stars. Science, 347(6222):632–635, February 2015. doi: 10.1126/science.
1258686.
Kristina Lerman, Xiaoran Yan, and Xin-Zeng Wu. The “majority illusion” in
social networks. PloS One, 11(2):e0147617, 2016.
Michael R. Line and Vivien Parmentier. THE INFLUENCE OF NONUNIFORM
CLOUD COVER ON TRANSIT TRANSMISSION SPECTRA. The Astrophys-
ical Journal, 820(1):78, mar 2016. doi: 10.3847/0004-637x/820/1/78. URL
https://doi.org/10.3847/0004-637x/820/1/78.
Beibei Liu and Adam P. Showman. Atmospheric Circulation of Hot Jupiters:
Insensitivity to Initial Conditions. , 770(1):42, June 2013. doi: 10.1088/
0004-637X/770/1/42.
Jan Lorenz. Continuous opinion dynamics under bounded confidence: A sur-
vey. International Journal of Modern Physics C, 18(12):1819–1838, 2007.
D Luz and F Hourdin. Latitudinal transport by barotropic waves in titan’s
stratosphere.: I. general properties from a horizontal shallow-water model.
Icarus, 166(2):328–342, 2003.
Jens Koed Madsen, Richard M Bailey, and Toby D Pilditch. Large networks
of rational agents form persistent echo chambers. Scientific Reports, 8(1):1–8,
2018.
Seth A Marvel, Jon Kleinberg, Robert D Kleinberg, and Steven H Strogatz.
Continuous-time model of structural balance. Proceedings of the National
Academy of Sciences, 108(5):1771–1776, 2011.
145
Naoki Masuda, N. Gibert, and S. Redner. Heterogeneous voter models. Phys.
Rev. E, 82:010103, Jul 2010. doi: 10.1103/PhysRevE.82.010103. URL https:
//link.aps.org/doi/10.1103/PhysRevE.82.010103.
John G Matsusaka. Election closeness and voter turnout: Evidence from califor-
nia ballot propositions. Public Choice, 76(4):313–334, 1993.
Richard McGehee and Esther Widiasih. A quadratic approximation to budyko’s
ice-albedo feedback model with ice line dynamics. SIAM Journal on Applied
Dynamical Systems, 13(1):518–536, 2014.
Kristen Menou and Emily Rauscher. Atmospheric circulation of hot jupiters: a
shallow three-dimensional model. The Astrophysical Journal, 700(1):887, 2009.
Kristen Menou, James YK Cho, Sara Seager, and Bradley MS Hansen. “weather”
variability of close-in extrasolar giant planets. The Astrophysical Journal Letters,
587(2):L113, 2003.
Alice Nadeau and Richard McGehee. A simple formula for a planet’s mean
annual insolation by latitude. Icarus, 291:46–50, 2017.
Alice Nadeau and Richard McGehee. Approximating annual mean incoming
solar radiation. Journal of Mathematical Analysis and Applications, To appear.:
XX–XX, 2021.
M. E. J. Newman. Networks. Oxford University Press, Oxford, UK, 2nd edition,
2018. ISBN 9780198805090.
Gerald R North. Analytical solution to a simple climate model with diffusive
heat transport. Journal of the Atmospheric Sciences, 32(7):1301–1307, 1975a.
146
Gerald R North. Theory of energy-balance climate models. Journal of the Atmo-
spheric Sciences, 32(11):2033–2043, 1975b.
Gerald R North. The small ice cap instability in diffusive climate models. Journal
of the atmospheric sciences, 41(23):3390–3395, 1984.
Kazumasa Ohno and Xi Zhang. Atmospheres on Nonsynchronized Eccentric-
tilted Exoplanets. I. Dynamical Regimes. , 874(1):1, March 2019. doi: 10.3847/
1538-4357/ab06cc.
Vivien Parmentier, Adam P Showman, and Yuan Lian. 3d mixing in hot jupiters
atmospheres-i. application to the day/night cold trap in hd 209458b. Astron-
omy & Astrophysics, 558:A91, 2013.
James Penn and Geoffrey K. Vallis. Atmospheric Circulation and Thermal
Phase-curve Offset of Tidally and Nontidally Locked Terrestrial Exoplanets. ,
868(2):147, December 2018. doi: 10.3847/1538-4357/aaeb20.
Daniel Perez-Becker and Adam P. Showman. Atmospheric Heat Redistribution
on Hot Jupiters. , 776(2):134, October 2013. doi: 10.1088/0004-637X/776/2/
134.
Daniel Perez-Becker and Adam P Showman. Atmospheric heat redistribution
on hot jupiters. The Astrophysical Journal, 776(2):134, 2013.
Sidney Redner. Reality-inspired voter models: A mini-review. Comptes Rendus
Physique, 20(4):275–292, 2019.
Andre J Robert. The integration of a low order spectral form of the primitive
meteorological equations. Journal of the Meteorological Society of Japan. Ser. II,
44(5):237–245, 1966.
147
Gerard H Roe and Marcia B Baker. Notes on a catastrophe: A feedback analysis
of Snowball Earth. Journal of climate, 23(17):4694–4703, 2010.
Brian EJ Rose and John Marshall. Ocean heat transport, sea ice, and multiple
climate states: Insights from energy balance models. Journal of the atmospheric
sciences, 66(9):2828–2843, 2009.
Brian EJ Rose, Timothy W Cronin, and Cecilia M Bitz. Ice caps and ice belts: The
effects of obliquity on ice- albedo feedback. The Astrophysical Journal, 846(1):
28, 2017.
M Rostami and V Zeitlin. Understanding dynamics of martian winter polar
vortex with “improved” moist-convective shallow water model. In Journal of
Physics: Conference Series, volume 936, page 012005. IOP Publishing, 2017.
David Rothschild. Forecasting elections: Comparing prediction markets, polls,
and their biases. Public Opinion Quarterly, 73(5):895–916, 2009.
Andrew J Rushby, Aomawa L Shields, Eric T Wolf, Marysa Laguë, and Adam
Burgasser. The effect of land albedo on the climate of land-dominated planets
in the trappist-1 system. The Astrophysical Journal, 904(2):124, 2020.
William D Sellers. A global climatic model based on the energy balance of the
earth-atmosphere system. Journal of Applied Meteorology, 8(3):392–400, 1969.
Karen M. Shell and Isaac M. Held. Abrupt Transition to Strong Superrotation
in an Axisymmetric Model of the Upper Troposphere. Journal of Atmospheric
Sciences, 61(23):2928–2935, December 2004. doi: 10.1175/JAS-3312.1.
A. P. Showman and T. Guillot. Atmospheric circulation and tides of “51 Pegasus
b-like” planets. , 385:166–180, April 2002. doi: 10.1051/0004-6361:20020101.
148
Adam Showman, Kristen Menou, and J. Cho. Atmospheric circulation of hot
jupiters: A review of current understanding. Astronomical Society of the Pacific
Conference Series, 398:419–441, 11 2008b.
Adam P Showman and Lorenzo M Polvani. Equatorial superrotation on tidally
locked exoplanets. The Astrophysical Journal, 738(1):71, 2011.
Adam P Showman, James YK Cho, and Kristen Menou. Atmospheric circulation
of exoplanets. Exoplanets, 526:471–516, 2010.
Adam P. Showman, Jonathan J. Fortney, Nikole K. Lewis, and Megan Shabram.
Doppler signatures of the atmospheric circulation on hot jupiters. The Astro-
physical Journal, 762(1):24, dec 2012. doi: 10.1088/0004-637x/762/1/24. URL
https://doi.org/10.1088/0004-637x/762/1/24.
Adam P Showman, Jonathan J Fortney, Nikole K Lewis, and Megan Shabram.
Doppler signatures of the atmospheric circulation on hot jupiters. The Astro-
physical Journal, 762(1):24, 2013.
Adam P. Showman, Nikole K. Lewis, and Jonathan J. Fortney. Three-
dimensional atmospheric circulation of warm and hot jupiters: effects of or-
bital distance, rotation period, and nonsynchronous rotation. The Astrophys-
ical Journal, 801(2):95, mar 2015. doi: 10.1088/0004-637x/801/2/95. URL
https://doi.org/10.1088/0004-637x/801/2/95.
Adam P Showman, Xianyu Tan, and Vivien Parmentier. Atmospheric dynamics
of hot giant planets and brown dwarfs. Space Science Reviews, 216(8):1–83,
2020.
František Slanina, Katarzyna Sznajd-Weron, and Piotr Przybyła. Some new re-
149
sults on one-dimensional outflow dynamics. EPL (Europhysics Letters), 82(1):
18006, 2008.
Peter H Stone. Constraints on dynamical transports of energy on a spherical
planet. Dynamics of atmospheres and oceans, 2(2):123–139, 1978.
Steven H Strogatz. Nonlinear dynamics and chaos with student solutions manual:
With applications to physics, biology, chemistry, and engineering. CRC press, 2018.
Jonathan D Touboul. The hipster effect: When anti-conformists all look the
same. Discrete & Continuous Dynamical Systems-B, 24(8):4379, 2019.
Amaury H. M. J. Triaud. Small Star Opportunities. In Nikku Madhusudhan,
editor, ExoFrontiers; Big Questions in Exoplanetary Science, pages 6–1. 2021. doi:
10.1088/2514-3433/abfa8fch6.
Ka-Kit Tung. Topics in mathematical modeling. Princeton University Press Prince-
ton, NJ, 2007.
Geoffrey Vallis. Atmospheric and Oceanic Fluid Dynamics. 01 2006. doi: 10.2277/
0521849691.
Geoffrey K Vallis. Atmospheric and oceanic fluid dynamics. Cambridge University
Press, 2017.
Pauli Virtanen, Ralf Gommers, Travis E Oliphant, Matt Haberland, Tyler Reddy,
David Cournapeau, Evgeni Burovski, Pearu Peterson, Warren Weckesser,
Jonathan Bright, et al. Scipy 1.0: fundamental algorithms for scientific com-
puting in python. Nature methods, 17(3):261–272, 2020.
Alexandria Volkening, Daniel F Linder, Mason A Porter, and Grzegorz A Rem-
150
pala. Forecasting elections using compartmental models of infection. SIAM
Review, 62(4):837–865, 2020.
Daniel Volovik and Sidney Redner. Dynamics of confident voting. Journal of
Statistical Mechanics: Theory and Experiment, 2012(04):P04003, 2012.
Till JW Wagner and Ian Eisenman. How climate model complexity influences
sea ice stability. Journal of Climate, 28(10):3998–4014, 2015.
HR Wakeford and PA Dalba. The exoplanet perspective on future ice giant ex-
ploration. Philosophical Transactions of the Royal Society A, 378(2187):20200054,
2020.
James Walsh. Diffusive heat transport in budyko’s energy balance climate
model with a dynamic ice line. Discrete and Continuous Dynamical Systems
B, 22(7):2687–2715, 2017.
James Walsh and Christopher Rackauckas. On the budyko-sellers energy bal-
ance climate model with ice line coupling. Discrete & Continuous Dynamical
Systems-B, 20(7):2187, 2015.
Huize Wang and Robin Wordsworth. Extremely long convergence times in
a 3d GCM simulation of the sub-neptune gliese 1214b. The Astrophysical
Journal, 891(1):7, feb 2020. doi: 10.3847/1538-4357/ab6dcc. URL https:
//doi.org/10.3847/1538-4357/ab6dcc.
Wei Wang, David Rothschild, Sharad Goel, and Andrew Gelman. Forecasting
elections with non-representative polls. International Journal of Forecasting, 31
(3):980–991, 2015.
William R Ward. Climatic variations on mars: 1. astronomical theory of insola-
tion. Journal of Geophysical Research, 79(24):3375–3386, 1974.
151
Michael J Way, Igor Aleinov, David S Amundsen, MA Chandler, TL Clune, An-
thony D Del Genio, Yuka Fujii, Maxwell Kelley, Nancy Y Kiang, Linda Sohl,
et al. Resolving orbital and climate keys of earth and extraterrestrial environ-
ments with dynamics (rocke-3d) 1.0: a general circulation model for simulat-
ing the climates of rocky planets. The Astrophysical Journal Supplement Series,
231(1):12, 2017.
Luis Welbanks, Nikku Madhusudhan, Nicole F. Allard, Ivan Hubeny, Fernand
Spiegelman, and Thierry Leininger. Mass–metallicity trends in transiting ex-
oplanets from atmospheric abundances of hsub2/subo, na, and k. The Astro-
physical Journal Letters, 887(1):L20, dec 2019. doi: 10.3847/2041-8213/ab5a89.
URL https://doi.org/10.3847/2041-8213/ab5a89.
Esther R Widiasih. Dynamics of the budyko energy balance model. SIAM Jour-
nal on Applied Dynamical Systems, 12(4):2068–2092, 2013.
Brian L Wiens. When log-normal and gamma models give different results: a
case study. The American Statistician, 53(2):89–93, 1999.
Jun Yang, Gwenaël Boué, Daniel C Fabrycky, and Dorian S Abbot. Strong de-
pendence of the inner edge of the habitable zone on planetary rotation rate.
The Astrophysical Journal Letters, 787(1):L2, 2014.
Wenshuo Yue and Jun Yang. Effect of sea-ice drift on the onset of snowball
climate on rapidly rotating aqua-planets. The Astrophysical Journal Letters, 898
(1):L19, 2020.
Xi Zhang and Adam P. Showman. Atmospheric Circulation of Brown Dwarfs:
Jets, Vortices, and Time Variability. , 788(1):L6, June 2014. doi: 10.1088/
2041-8205/788/1/L6.
152
Xi Zhang and Adam P Showman. Effects of bulk composition on the atmo-
spheric dynamics on close-in exoplanets. The Astrophysical Journal, 836(1):73,
2017.
MV Zombeck. Book-review-handbook of space astronomy and astrophysics.
Astronomy Quarterly, 7:256, 1990.
153