# Fermi Project Update January 2022

The journal article on overleaf is at: https://www.overleaf.com/3756433924fzhskcnhbkmz

**What Else I Need to Do?**

~~Results Section: fix labels on parameter sweep (Fig 2). Use the parameter names on the paper rather than in the code~~- Appendix B: currently has the relevant timescales and ratios written down, what I need to do is figure out which 4/5 parameters or ratios uniquely determine the output. Or influence the output the most.
- Appendix C: currently has density calculations (as function of radius). I need to copy my results showing technology as function of density. Then I need to use both of those to figure out technology as a function of radius.

**What Is Left?**

- send out the paper to collaborators for review
- write introduction, conclusion, and abstract

# Fermi Project Update December 2021

**Goal:** find technology as a function of galactic orbital radius

- Find density as a function of radius

- Get technology as a function of density

2a. Recreate surface plot made with Jonathan (with )

- f is the fraction of systems that are settleable
- is the normalized density of settleable systems within probe range
- is the settlement civilization lifetime
- X is the local fraction of settled systems to total systems

2b. Rerun the same parameter sweep but make

2c. Do some lineouts to get technology and X as function of density

- Interpolate to get density as a function of radius. (not entirely sure how to do this one)

# Fermi Project Update July 2021

## Parameter Sweep: Periodic Box Model

# Fermi Project Update: March 2021

As a reminder, we had previously made the probe ranges and velocities normally distributed by introducing another variable called technology, making this normally distributed through random walks, then making the probe range and velocity relativistic functions of technology.

The current goal we have been working on is introducing a fixed probability of civilization destruction per timestep for the civilizations. We ended up accomplishing this by making the civilization lifetime exponentially distributed, analogous to how we had previously had the abiogenesis time exponentially distributed.

Click here for a summary of this and some parameter sweeps.

# Update March 2021

## Coupled EBM

I submitted the paper to AJ and arXiv.

Our Paper is on arXiv Here

ADS Link

## Fermi Project

# Fermi Project Update 02/18/21

**Goal:** Figure out how to make the probe velocity and range normally distributed within our FORTRAN model.

**Method:** First added a local variable called 'technology' to the program.
Then used random walks to have it be normally distributed around the technological capabilities of the abiogenesis seed (first habited system). Finally, I made the probe range and velocity relativistic functions of this new technology variable. The table below shows the initial values inputted into the model.

Inputs | ||
---|---|---|

v0 | 0.0001c=30 km/s | Initial Probe Velocity |

r0 | 10 lyr | Initial Probe Range |

t0=r0/v0 | 100,000 years | Initial Probe Lifetime |

**Click here to see the PDF which summarizes calculation/implementation of the normally distributed probe ranges/distances**

**Click here for more information about the input variables for the models shown below.**

**Notes**

- The pink boxes surround systems that are uninhabited.
- All simulations begin with 10 systems originally habited.
- Galaxy Model shows time with unit Myr. Total runtime being 1000 Myr.
- In contrast, the periodic box model has a total runtime of 1 Myr. Thus I show frames instead of time, where each frame is approximately 1/1000 Myr.
- x=(number of habited systems)/(total number of systems)

## Galaxy Model

## Periodic Box Model

## Analytic Model

# Fermi Project Update (02/13/2021)

**Goal:** Figure out how to make the probe velocity and range normally distributed within our FORTRAN model.

**Method:** First added a local variable called 'technology' to the program.
Then used random walks to have it be normally distributed around the technological capabilities of the abiogenesis seed (first habited system). Finally, I made the probe range and velocity relativistic functions of this new technology variable. See the pdf below for more details about calculation/implementation.

**Click here to see the PDF which summarizes my progress so far and potential steps.**

### Current Model Output (Colored by Technology, Pink=Unsettled)

Inputs | ||
---|---|---|

v0 | 0.0001c=30 km/s | Initial Probe Velocity |

r0 | 10 lyr | Initial Probe Range |

t0=r0/v0 | 100,000 years | Initial Probe Lifetime |

**Figure:** The above gif shows the temporal evolution of a model galaxy. Beginning with 10 habited solar systems (ie: abiogensis seeds), these systems send probes out to nearby systems, thus making those systems inhabited and repeating the process until the entire galaxy is filled with life. The uninhabited systems are shown as the pink boxes. Note that a selection effect results in the systems on the outer edge of the galaxy gaining high technological abilities before systems near the center.

# Coupled EBM Update 01/09

## Updates

Senior Thesis has Finally Been Finished

## ToDo for Paper

Preparing for submission into the astrophysical journal. Other than basic formatting and spellchecks, it seems like our paper needs three appendices. I'll make these this week, most of the information is in my senior thesis so I just need to rewrite/reformat it.

- Gamma/Theta/Beta Derivation
- Derive timescales (t_C,t_G,t_T)
- This is important for when we derive the dimensionless timescales

- Use these to derive gamma/theta/beta, and N_A
- This is important for when we plot N_A in the contour plots

- Talk about the difference between gamma and gammaEff (gamma scales with initial climate sensitivity (dTdP|T=T_0) while gammaEff scales with climate sensitivity when temperature has increased by the temperature tolerance (dTdP|T=T0+dT)
- This is important for plotting gammaEff on the contourplot and on the scatter plot vs the decline time. It helps reduce the scatter greatly for the latter.

- Derive timescales (t_C,t_G,t_T)
- Derive dimensionless timescale for high gamma (tau_coll=1/max(sqrt(2),theta)
- tauColl=1/sqrt(2) for low theta
- tauColl=1/theta for high theta

- Explain scatter in gamma vs time to decline plot
- Show plot of gamma vs decline time compared against analytical predictions for both gamma and gammaEff, to show how the effective gamma reduces the scatter
- Show plot of gammaEff vs decline time colored by both T0/a/P0, to show how the points that diverge from our analytical predictions are those with low initial temperature (T0) and large orbital distance (a), yet P0 is still approximately constant for constant gamma

# Coupled EBM Project Update 09/26

### Outline

- Introduction
- Numerical Model
- Modeling Earth

- Analytic Model
- Dimensionless Model

- Numerical Experiment Setup
- 5 Trajectories

- Full Suite of Results
- Parameter Sweeps
- KDE plots of decline time for different dT's

- Discussion

### Questions

- Formatting: for emulateapj LaTeX document class, does anyone know how to have a figure on the top/bottom of the page that spans both columns?
- When talking about CO2, should I say P(CO_2) or pCO_2?
- For our timescales, call them…
- t_c, t_g, t_t?
- Or t_C, t_G, t_T?
- Or t_clim, t_grow, t_tech?

### ToDo

- Finish the dT sweeps for the new version of the model, currently running.
- Finish filling in the details of section 1-4.
- Finish writing up section 5.
- Start/Finish writing up section 6.

## Current Version of the Paper is Here

# Update 08/01

# Update 07/27

# Movie Night! (07/29, 7:30PM)

Vote here if you haven't voted yet.

As of now we have a two way tie for first place.

- Hitchhikers Guide to the Galaxy
- Inception (tied for first)
- Forbidden Planet

# Fermi Project

**Goal:** Figure out how to make the probe velocity and range normally distributed within our FORTRAN model.

**Click here to see the PDF which summarizes my progress so far and my next steps.**

### Current Model Output (Colored by Technology, Pink=Unsettled)

Inputs | ||
---|---|---|

v0 | 0.0001c=30 km/s | Initial Probe Velocity |

r0 | 10 lyr | Initial Probe Range |

t0=r0/v0 | 100,000 years | Initial Probe Lifetime |

# Coupled EBM

## Parameter Sweeps (20 different temperatures and distances)

### Contour Plots (Earth-Like Values)

### KDE Plots

# Update 07/16

# Fermi Project

**Goal:** Figure out how to make the probe velocity and range normally distributed within our FORTRAN model.

**Click here to see the PDF which summarizes my progress so far and my next steps.**

### Current Model Output (Colored by Technology, Pink=Unsettled)

Inputs | ||
---|---|---|

v0 | 0.0001c=30 km/s | Initial Probe Velocity |

r0 | 10 lyr | Initial Probe Range |

t0=r0/v0 | 100,000 years | Initial Probe Lifetime |

# Coupled EBM

## Comparing Experiments

**dPdT:** Amount of pCO2 that must be added to the atmosphere in order to change the climate by 1 degree Kelvin.

### Experiment #1: Constant Composition (P0=284 ppm)

### Experiment #2: Constant Temperature (T0=287.09 K)

## Earth-Like Inputs

- t0=1820
- We choose this because our data starts here.

- P0=284 ppm
- This was the global CO2 concentrations during 1820.

- N0=1.129 billion people
- This was the global population during 1820.

- Nmax=10 Billion people
- This one is based on current trends and global projections

- dT=5K=Temperature Sensitivity
- We chose this because it seemed reasonable in regards to what kind of temperature range humans can survive in.

- dP=200 ppm=Carbon Dioxide Sensitivity
- Similar to dT, we choose this because it seemed reasonable

- A0=0.04 1/yr=Initial Per-Capita Birth Rate
- We chose this based on the assumption that the average time for a generation is 25 years, so that A0=1/25=0.04 1/yr

- B0=0.036 1/yr=Initial Per-Capita Death Rate
- We chose this value in order to tune our model to the data.

- C=0.000275 ppm/(10
^{6}ppl*yr)=Per-Capita Carbon Footprint- We chose this based on present-day values

### Are these numbers reasonable?

The relevant quantity for how population grows is the difference between A0 and B0, this quantity is called the per-capita (net) growth rate. The relevant quantity for how CO2 concentrations change as a function of population growth is C, the per-capita carbon footprint. We can compare our models values against global values to determine how accurate our model is.

# Fermi Project Update 07/02

**Goal:** Figure out how to make the probe velocity and range normally distributed within our FORTRAN model.

**Summary of Progress**:
Click here to see the PDF which summarizes my progress so far and my next steps.

** Video Of Current Model Output (Colored by Technology, Pink=Unsettled): **
http://www.pas.rochester.edu/~esavitch/fermiProject/results2.ogv

# Fermi Project Update 06/28

**Goal:** Figure out how to make the constant parameters that are currently in the Fermi FORTRAN program instead normally distributed.

**Summary of Progress**:
Click here to see the PDF which summarizes my progress so far and my next steps.

** Video Of Current Model Output (Colored by Technology): **
http://www.pas.rochester.edu/~esavitch/fermiProject/results.ogv

# Fermi Project Update 06/14

**Goal:** Figure out how to make the constant parameters that are currently in the Fermi FORTRAN program instead normally distributed.

**Summary of Progress**:
Click here to see the PDF which summarizes my progress so far and my next steps.

# Coupled EBM Project Update 02/24

Explanation of Numerical and Analytic Models, with examples using earth-like inputs

### Results of Linear Fits of pCO2 vs Temp for Both Experiments

http://www.pas.rochester.edu/~esavitch/planet_Plots/both_exp_dTdP.png

### 4x5 Arrays for Experiment One: "Constant Composition"

Array of Plots

Array of Phase Plots

### 4x5 Arrays for Experiment Two: "Constant Temperature"

Question: for these, should I have the initial pco2 be the value that makes the initial temperature be approximately like earth's? Also, this initial pco2 should be constant across maxPops and variable across distances?

# Planet-Civ Evolution Project Update 2/16

### Updates: New Method for Linear Regressions

# Planet-Civ Evolution Project Update 2/9

## Current State of the Equations and Example Outputs

Current Sate of the Equations and Example Outputs Using Earth-Like Inputs

### 4x5 Grid of Plots with Earth-Like dT

- dT parameter is a measure of the temperature range that a civilization can survive in
- Horizontal plots vary by orbital distance
- 0.94AU, 0.97AU, 1.00AU, 1.03AU, 1.06AU

- Vertical plots vary by carrying capacity (maximum population)
- 10 Billion, 40 Billion, 70 Billion, 100 Billion

Plots Showing the Temporal Evolution of Global Temperature and Population (colored by pCO2)

Phase Plots of Global Population vs Global Temperature (colored by pCO2)

### 4x5 Grid of Plots with Earth-Like dT/2 (High Fragility)

(0.5dT) Plots Showing the Temporal Evolution of Global Temperature and Population (colored by pCO2)

(0.5dT) Phase Plots of Global Population vs Global Temperature (colored by pCO2)

### 2x4 Grid of Plots with Earth-Like dT*100 (Low Fragility)

# Planet-Civ Evolution Project Update 1/27/20

Current state of the Equations and example plots using Earth-like inputs

4x4 Grid of Temperature/Population vs time: Horizontal Plots vary by Orbital Distance, Vertical Plots vary by the Maximum Allowed Population (Carrying Capacity)

4x4 Grid of Population vs Temperature: Horizontal Plots vary by Orbital Distance, Vertical Plots vary by the Maximum Allowed Population (Carrying Capacity)

# Planet-Civ Evolution Project Update 11/6

## Energy Balance Model (EBM)

The EBM reaches equilibrium at Teq. Below are plots of this for varying distances and partial CO2 pressures (pCO2):

3D Plot (top view)

3D Plot (side view)

Contour Plot

## Numerical Model

Code Location on GitHub
(**Initial Conditions**: t=1880, pCO2=284ppm, N=1.13 Billion People)

### General Simulation Steps

- Run program with Earth like inputs as they were before the Anthropocene (in 1880), until an equilibrium temperature is reached.

- Set the populations peak growth rate to this temperature as shown here.

- Then turn on population growth; exponential rise in population results in exponentially increasing partial CO2 pressures, ultimately triggering a global climate change. Once the global temperatures pass a certain point, the growth rate of the population becomes negative, which results in mass extinction.

### Goal #1

Repeat steps 1-3 for **four different values of initial partial CO2 pressure**, ranging from a tenth of ours to one hundred times ours. At each new value of pCO2, determine the habitable zone for that point, defined here to be where water can be liquid (273-373K).

Then repeat steps 1-3 for **4 different values of orbital distance**, spanning the entirety of that habitable zone.

Results Shown Here (as 4x4 grid)

### Goal #2

Choose **3 different orbital distances**. For each one, increase pCO2 until the equilibrium temperature reaches 287.09K (the temperature of earth at 1AU).

Results Shown Here (Earth included for reference)

### Next Steps

- Model 1: Introduce constant carrying capacity (maximum population)
- Model 2: Introduce a temperature-dependence to the carrying capacity, to emulate reduction of habitable regions due to global climate change.
- Model 3: Introduce a time-dependence to the carrying capacity, to emulate environmental sensitivity to climate change. So environments with low sensitivity could exceed their carrying capacities for a longer time than ones with higher sensitivity before climate change begins.

## Analytical Model

Current state of the Equations and Nondimensionalization Process

Phase Diagrams of Analytic Model

# Planet-Civ Evolution Project Update 10/21

Current state of the Equations

Link to Phase Diagrams of Analytic Model

# Plots of our Habitable Zone (273K-373K)

# Coupled EBM Project Update 10/13

# Goal

Make plots for the equilibrium temperatures as a function of pCO2 and orbital distance.

(*Equilibrium Temperature* = temperature at which model balances incoming solar radiation with outgoing long-wave terrestrial radiation, without coupling)

## Visualizing How our Model Behaves at Different Orbital Radii

# Coupled EBM Project Update 10/4

## Equations As of Now

- T = Temperature
- P = pCO2
- = per-capita carbon footprint

- N = Population
- = per-capita growth rate, dependent on temperature
- =temperature resulting in optimum growth rate
- = exponential ramping distance
- = proxy for technology, kept constant here

## Methodology

- Run energy balance model until equilibrium is reached, with an average global temperature of
- Center the optimum growth rate for the civilization to that temperature, ie: set
- Turn on coupling functionality, which allows population to rise, causing more CO2 to be in the atmosphere, which inevitably raises global temperatures past a habitable point

## Example outputs for Varying Orbital Distances

## notes

- 1.524AU is the orbital distance of mars
- 273.15K is the freezing temperature of water
- 373.15K is the boiling temperature of water

# Coupled EBM Project Update 9/21

# Model Output vs Global (True) Values

**Population and pCO2 Data From:** Frank, Adam, and Woodruff Sullivan. "Sustainability and the astrobiological perspective: Framing human futures in a planetary context." Anthropocene 5 (2014): 32-41.

** Temperature Data From: ** https://data.giss.nasa.gov/gistemp/graphs/graph_data/Global_Mean_Estimates_based_on_Land_and_Ocean_Data/graph.txt

How I Converted:

- Got temperature anomaly data in Celsius from NASA Goddard's website (both smoothed and raw)
- Converted the smoothed anomaly data into smoothed global temperature dat (in celsius) using
**Tony Del Genio's advice**:*People only report global temperature anomalies because surface weather stations only cover a fraction of the globe, and thus various assumptions are made by different groups to fill in all the places in the world with no data, so the actual temperature of Earth is only known to within a degree or so. People usually define the anomaly relative to the 1951-1980 reference period, whose mean anomaly should be close to zero. If you*(will post his talk at the bottom of this post)**add a number like 14.5 deg. C to that**, that should be as close to a real temperature as you can get. - Then converted from celsius to kelvin (+273.15 deg. C), to compare that data to our model's output

# Visualization of the Anthropocene With Actual Global Data

The above plot has been made using the data from Frank and Sullivan, 2014. From this plot, and with the assumption that the Anthropocene started at around 1945, it seems safe to define the Anthropocene physically as the point at which the relationship between human population and their carbon footprint switched from linear to exponential.

## Tony's Talk: Albedos, Temperatures, and Habitability of Rocky Exoplanets

# Coupled EBM Project Update 9/9

# Code Location on GitHub

https://github.com/ethansavitch/coupled_ebm

## Jacob's Answers To My Questions

**Purpose of the Stochastic Parameter?**

The stochastic parameter in the EBM equation was used in my glacial-interglacial cycles paper (cited on your poster), which was done to explore the concept of "stochastic resonance" in a climate system. The idea is that the Earth system has multiple stable climate states (ice free, large ice cap, small ice cap, ice covered), and that variability in the Earth system itself might contribute to periodic transitions between states (which for small random variability would be between the large and small ice cap solutions). The stochastic parameter thus is a hand-wavy way of approximating the effect of sub-model-scale processes that would exert a random noise upon equilibrium climate. The net effect of this is that the hysteresis loop decreases in width, since it becomes slightly easier to access different states. **There's no need to include the stochastic parameter for accuracy**. I included a longer explanation here in case it is useful at a later stage in the project.

**What's the purpose of the soladj=solar forcing adjustment parameter?**

The solar forcing adjustment is used to time step the model through a prescribed increase in stellar brightness as the star ages through the main sequence. If you look near line 265 of the code, you will see that when soladj is true, the model reads in the eccentricity, precession, obliquity, and solar constant from the file data/insolaout.dat. This data is from Laskar's astronomical solutions for Earth's past and future (http://vo.imcce.fr/insola/earth/online/earth/earth.html). If soladj is false, then these values are all treated as fixed with time. **You only need the solar forcing adjustment feature if you are interested in long time integrations with realistic representation of solar brightening and orbital evolution.**

**How best to incorporate weathering?**

If you want CO2 weathering as described by Williams & Kasting (Icarus, 1997), then you could uncomment line 637, where I have added a comment "removed co2 adjusting.". However, the timescales we are talking about are on the order of thousands to tens of thousands of years, where shorter-term variation in atmospheric CO2 will dominate. The **carbonate-silicate cycle operates with a residence time of about half a million years**, and only matters for Earth climate when we consider timescales into the hundreds of thousands of years (see, e.g., Berner, Lasaga, & Garrels, Am. J. Sci., 1983). So I would say leave the code as-is, and **don't worry about carbonate-silicate weathering**. If you do want to experiment with the effects of weathering, then the solution in 3b of your blog post sounds like a good approach to me. (constant annual per-capita reduction in CO2, parameterized with )

## ToDo Next

- Clean/Fix Jupyter Notebook. All parts are currently working except for the namelist functionality.
- Compare output of coupled EBM with global average temperatures.
- Increase accuracy of the growth rates temperature dependencies, which will hopefully make our output more in line with the global averages.

# Questions for Jacob

## Shorthand

In this blogpost, in order to specify between the original EBM and the second version of the EBM, I will denote:

**D1**= original, faster driver, with the simplified radiative transfer- more accurate for small changes in radiation

**D2**= second, larger, slower driver, which allows explicit pCO2 values- by including a polynomial parameterization of outgoing longwave radiation and planetary albedo

# Questions

- in D2 there is no stochastic functionality/variables (as there is in D1)
- Will the addition of stochastic calculations in D2 help or hurt the accuracy of our model?

- in D1 there is a variable called soladj=solar forcing adjustment…
- What is its purpose?
- It is not in D2, will that considerably affect the accuracy of the results?

- in the version of D2 you sent me, is there any built-in, uncommented, functionality which dictates how pCO2 changes with time?
- I've added an ODE for pCO2 which allows me to specify as an input into D2 a civilizations annual, per-capita addition of atmospheric CO2:
- Should I leave this be, and turn on the programs built-in reduction of pCO2 by weathering?
- if yes, how?

- Should I reduce a term from my existing pCO2 ODE (
- allowing me to have an input in D2, ( =annual reduction of pCO2 by Earths natural weathering processes)

) and turn off the programs built-in functionality?

- Should I leave this be, and turn on the programs built-in reduction of pCO2 by weathering?

- I've added an ODE for pCO2 which allows me to specify as an input into D2 a civilizations annual, per-capita addition of atmospheric CO2:

# Coupled EBM Meeting Update 08/05

## Summary of Current State of Equations

Link to REU Paper describing Summary of Work Done This Summer (so far)

**Progress**: Calibrated our model with population and CO2 data from the past two centuries, summed up in the table below.

(yrs) 0 140 190 198 t=time (yrs) 1820 1960 2010 2018 N=Population ( ppl)1.042 3.032 6.834 7.594 P=pCO2 (ppm) 284 316 390 413

## Poster Showing Result of Modelling Earth

## Next Steps

- increase the accuracy of modeling earth
- find true temperature dependence on the growth rate
- find more accurate value for (=reduction of CO2 by earth's natural weathering processes)

- determine the climate sensitivity of various energy resources
- can let us determine the inevitability of global warming

- put a dTdt dependence on the deathrate
- to quantify acclimation

- introduce parameter E as a proxy for technology
- will raise peak of relative growth rate (r)
- will increase the average carbon footprint ( )

# Coupled EBM Project Update 7/29

# Overview of Current Model

## Change in Population

To model rise in population, we begin with a basic logistic growth equation, which models rise in any population with a carrying capacity, by having exponential rise slow as the population size approaches its maximum.

- T = temperature
- E = technological efficiency (ie: measurement of how much technology a civilization has)

= relative growth rate
- = current population
- = maximum population
- = global temperature resulting in maximum growth rates
- = initial growth rate
- assuming a 70 year avg lifespan

= initial death rate
- = width of gaussian (or exponential ramping distance)

## Change in Temperature

In order to model rise in temperature, we use an energy balance model, which balances the incoming solar radiation with the outgoing longwave radiation, in order to give an approximate planetary temperature.

- S(1-a) = amount of solar flux being absorbed by the planet
- S = diurnally averaged solar flux
- a = albedo of the planet (1 for perfect reflector)

- IR = outgoing infrared flux, modeled with a higher-order polynomial parameterization of outgoing longwave radiation
- C = effective heat capacity of surface and atmosphere =
- D = diffusive parameter describing efficiency of energy transport =
- L = latitude (18 bands total)

## Change in CO2

We model changes in CO2 using a constant, annual, per-capita CO2 increase, assuming constant technology (ie: constant E). We can eventually progress this model by having E progress with time, as it would for any realistic technological civilization. If I were to leave the model here, we would inevitable end up with extinction. Thus, the extra term that is subtracted represents the constant reduction of CO2 by trees

- = annual average pco2 contributions from 1 million people
- = constant reduction of CO2 by trees

## Results (E=1)

**E** acts to quantify a civilizations technological capacities, and affects our model in two ways.
Higher E results in…

- a higher peak growth rate
- the more technology you have, the lower your infantile death rates becomes, which raises the relative growth rate

- higher per capita CO2 contributions
- the more technologically advanced a civilization, the more evenly distributed the technology becomes. (ex: cars were invented in the late 1700's, yet weren't commonplace until the 1900's)

## Next Steps

- Put a dTdt dependence on the death rate as a way to quantify acclimation
- Put a time dependence on E as a way to quantify the technological advancements of energy-harvesting civilizations
- Include more complete models for the environmental impacts of various resources, so that given a set of planetary/climate parameters, we could calculate climate sensitivity for the use of various resources
- Put more accurate temperature dependence onto the relative growth rate

# Coupled EBM Project Update 7/22

The above graph details what range of temperatures an unprotected human can survive in…

- Minimum Temperature:
- Maximum Temperature:

# Coupling Population Growth Equation to EBM

**Goal:** Couple the logistic growth equation to the energy balence model (EBM) by putting a temperature dependence on the relative growth rate

- = current population
- = maximum population
- T = temperature
- E = technological efficiency (ie: measurement of how much technology a civilization has)

= relative growth rate

**Method:** the behavior we desire can be modelled by the difference of two gaussian functions, corresponding to the average birth rate and death rate of a civilization:

- = global temperature resulting in maximum growth rates
- = initial growth rate
- assuming a 70 year avg lifespan

= initial death rate

# Coupled EBM Project Update 7/15

**EBM** (Energy Balance Model) = climate model that balances the incoming
solar radiation with outgoing infrared radiation (IR) to approximate
the temporal progression of planetary temperature

**Simplifications**(from Williams & Kasting (1997))- ebm averages heat capacity, surface albedo, and temperature over 10 deg wide latitudinal areas…
- causes us to overestimate the amplitude of the seasonal cycle over the continents
- also in doing this we ignore meridional land-sea temperature gradients which affect dynamic heat transport and weather
- becomes more of an issue with higher obliquity, as these types of planets have very large latitudinal temperature gradient

- ebm models latitudinal heat transport as diffusion, where in reality, wind patterns and energy transport are much more complicated

- ebm averages heat capacity, surface albedo, and temperature over 10 deg wide latitudinal areas…

**Initial IR calculation:** linear relationship for outgoing IR, based off of North et al. 1981 (A+BT)

**New IR calculation:** higher-order polynomial parameterization of outgoing longwave radiation and planetary albedo, from Williams & Kasting (1997) paper, obtained from fits to Jim Kasting's radiative-convective model earthtem_ir

**Benefits:**- lets us specify atmospheric co2 levels with variable=pco2 (unit bar)

**Limitations:**- error is 4.56 W*m
^{-2 } - this version of ir can be applied for…
- 190 K < Temperatures < 370 K
- .00001 bars < pco2 < 10 bars

- error is 4.56 W*m

# Exploring the EBM

**Goal:** Use our Energy Balance Model (EBM) to approximate the width of the habitable zone of our solar system

### Method

- Assuming all other parameters of the Earth are held constant, we change the effective orbital distance (d) by altering the input corresponding to the relative solar constant
- the solar constant (S) is normally a function of stellar Luminosity and orbital distance, but with Luminosity being constant, it reduces to an inverse square law:

- The Habitable ("Goldilocks") Zone is defined to be any orbital distance that can result in liquid water:
- Water Freezes: (blue dotted line)
- Water Boils: (red dotted line)

- I then run our
*uncoupled*EBM, with*default, deterministic, Earth-like inputs, for*, and then plot the final global temperatures, discretized by latitude**three years**- black dotted line is the final globally averaged temperature

The table below summarizes the findings from the plots above…

Orbital Distance (AU) Relative Solar Constant (relsolcon) New Solar Constant (Wm ^{-2) }Final Globally Averaged Temperature (Kelvin) .61 2.69 3,658.4 371.94 Too Hot (left plot) 1 1 1,360 287 Just Right (middle plot) 1.2 0.69 938.4 274.05 Too Cold (right plot)

**Conclusion:** This data lets us give an approximate width for our habitable zone of 0.59 AU

*(Williams and Kastings (1997) estimated a width of our habitable zone of approximately 0.4 AU)*

## Next Steps

- Add complexity to the population logistic growth equation
- Then couple the updated EBM to the updated population model
- Explore Results!

# EBM Project Update

## Project Description

Figure out how to run Jacob's Fortran Energy Balance Model through Mathematica. Put in inputs within Mathematica notebook, then receive/extract all data through the same notebook.

## Idea 1: Use MathLink (status: not working)

- MathLink
- allows external programs both to call Mathematica, and to be called by Mathematica. By using MathLink, you can, for example, treat Mathematica essentially like a subroutine embedded inside an external program. You can also use MathLink to let Mathematica call individual functions inside an external program. (MathLink (as well as LibraryLink) are designed to be used in conjunction with programs written in C/C++. )

Where I'm at here is I have outlined the 4 steps/goals that if accomplished would allow me to successfully use MathLink. I have successfully accomplished the first two steps, yet the other two remain unworking…

- Successfully test MathLink using the example file they have built into Mathematica by compiling the example C++ program and template file into an executable file which is then run through Mathematica.
*SUCCESS*

- Figure out how to call a Fortran 90 subroutine from a C++ program.
*SUCCESS*

- Successfully test MathLink's ability to call an external FORTRAN subroutine from Mathematica by moving the meat of the example C++ program to a FORTRAN test subroutine, and then compile the example C++ program, the test FORTRAN subroutine, and the example template file into an executable file which is then run through Mathematica.
*UNSUCCESSFUL*(I keep on running into various bugs)

- Do the same as step 3, yet instead of the example C++ program, create a wrapper which will call the EBM driver.f90 instead of the test FORTRAN subroutine.
*UNSUCCESSFUL*(can only work if I figure out step 3)

## Idea 2: Use LibraryLink (status: semi-working)

- LibrayLink
- provides a powerful way to connect external code to the Wolfram Language, enabling high-speed and memory-efficient execution. It does this by allowing dynamic libraries to be directly loaded into the Wolfram Language kernel so that functions in the libraries can be immediately called from the Wolfram Language.

I think I finally have LibraryLink working at least well enough for us to use it in conjunction with Jacob's EBM model. When I run it through Mathematica, I can now…

- Successfully extract:
- globally averaged temperature (double)
- minimum global temperatures by lat band (array size 18)
- maximum global temperatures by lat band (array size 18)

- Unsuccessfully extract:
- average global temperatures by lat band (array size 18)

I think this may not be an issue for our purposes, as it seems like we only need the globally averaged temperature. In addition, if we want an approximated average temp array, I can calculate that just by averaging the min and max temp arrays, as I did in the Mathematica notebook.

### Gif Demonstrating that LibraryLink is Working by Showing how the Temperature Arrays change with Obliquity*…

*(the avg temp array used here was calculated by averaging the min and max temp arrays, it was NOT calculated directly from the EBM program)