# 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

# Update 7/22

## Charge Exchange

A few more frames. So far, charge exchange is occurring, presumably creating a hot neutral population which is then immediately ionized (on the next line transfer step [not subcycle], by both photoionization and possibly ionization by "reverse" charge exchange). This suggests that the amount of neutral stellar material should be alternating between zero and some small value every step, which means at some point a frame may be output with nonzero amounts of Hhot….

Charge exchange should be creating at most ~10^{-4} Hhot (in CN). The highest-rate locations are in the shadow of the planet. HIIcold reacts with Hhot, and we have large amounts (~10^{-2}-10^{0}) of HIIcold. This reverse charge exchange is already taken into account in the rate calculation, but if it's oscillating we may not have seen the opposite side of the step yet.

All this to say I'm not yet sure where the Hhot is going, and I haven't figured out the appropriate diagnostics to look at to make sure everything is happening as it should be, even if it's not quite the expected result (everything being zero when Hhot is zero, and not having seen a frame with nonzero Hhot yet).

Here are some figures for the latest frame:

rho, 1-X, T

Proportion of species

Clockwise from top left: Ionized cold, ionized hot, neutral hot, neutral cold

Amount of species

"Charge exchange potential"

min(Hcold, HIIhot)/max(Hcold,HIIhot) - so proportion of material that can undergo an exchange, disregarding the reverse interaction.

Charge exchange rate (CU)

Multiply by dt to get delta Hhot

### AMR line transfer ==

Still hasn't started on BlueHive.

## Non-rotating charge exchange

Still waiting on AMR line transfer.

## Paper

Good discussion Thursday. Haven't had much chance to start making figures or new changes yet.

# 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

# Update 7/15

## Charge Exchange

Working on stampede. Here's the first frame:

I believe the neutral planetary material top-left is from the higher density (higher recombination rate, possibly also from higher temperature) in the stellar-planetary wind shock. Note that no charge exchange has occurred yet (but there's little contact between neutral planetary and ionized stellar material so far).

Trying to make a VisIt expression for charge exchange potential, but the following gives a "cannot divide by zero" error:

if(gt(Hcold,1e-2), HIIhot/Hcold, 0)

Note that the division should never happen if there could be a zero in the denominator, so I'm not sure why VisIt is complaining here…

## AMR line transfer

Fixed the missing condition on the displaced child mask update. Still waiting on newest test on BlueHive.

## Non-rotating charge exchange timing tests

Have the steady states ready. Need to calculate the stellar winds - use original distance (vs. larger orbital distance e.g. O-type)?

## Radiation Pressure Paper

Got referee response. Sent around preliminary thoughts, will try to start addressing them this week. To recap my email:

- Collisional nature of the wind (not hard-body collisional but electrostatically collisional)
- Clearer distinction between our focus on the cometary tail vs. the observed velocities
- Discussing more thoroughly the absorption of momentum by ionized hydrogen

# 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!

# Update 7/8

## Charge Exchange

As usual, Stampede is continually frustrating (related: Jonathan, you mentioned a plan to migrate from gitlab.circ to git.circ—any ETA on that?). Still making progress, though.

## AMR timing

Baseline test complete (single frame of HD209458b steady state). What we believe to be the most efficient version is queued on BlueHive. Here are the relevant stats from the baseline run:

Total Runtime | 43101 s |

Total AMR Time | 40543 s |

Relative Time Across All Levels: LineTransfer | 56.10% |

Relative Time Within MaxLevel: LineTransfer | 61.06% |

Relative Times of Each Level: 4 | 91.87% |

## Other tests

- Charge exchange timing for non-rotating parameter space planets
- AMR ray tracking

## Misc

Posted my suggested progress through Toro on a blog post. Feel like there should be a better place to put this - perhaps an "Introduction for New Students"-type page would be useful?

# Suggestions for Working Through Toro

- Read Chs 1-3 on the Euler equations and solutions to hyperbolic equations in general. Some of this will likely be familiar.
- Read Ch 4 and write the exact Euler solver described.
- Read Ch 5 and write the Godunov upwind scheme, for the linear systems described.
- Read Ch 6 and update the Godunov upwind scheme to use the exact Euler solver.
- Skip 7&8
- Read Ch 9 and write the approximate-state solvers (good place to find out if you're following good programming practices).
- Read/skim Chs 10-12 and write at least the HLLC solver (which is the one used in AstroBEAR). Note that there have been some significant improvements to the HLLC chapter in the third edition. Write the rest if you're interested, though I wouldn't recommend the Osher solver (it's interesting mathematically, but doesn't seem to offer much over the other solvers).
- Read Ch 16 and write a 2D/3D solver (test problem in Ch 17). Another place where you'll really notice if you've been following good practice.
- For second-order solvers: Read Ch 14 (skim 13 for an introduction to the same concepts for scalar equations) and write the WAF and MUSCL schemes (esp. MUSCL-Hancock, used in AstroBEAR).

Note that the ends of the code-related chapters have example code. I'd strongly suggest looking at some modern Fortran style guides rather than following the old style given in these, not least because it's much easier to follow your own code if it's written with meaning.

# Update 07/01

# CEJet

## From last week: 3 test runs planned

4 levels of AMR with 64 base grid resolution. Each run is planed for 10 days in simulation time, and takes about 2 days to complete on bluehive.

### Run 45 - low ambient, small separation

Secondary starts right against the primary envelope: a0 = 49 R_sun.

Use lower ambient density and pressure profile: rho_amb = 1e-10 g/cc, p_amb = 1961 dyn/cm2 .

density-edgeon-throughparticles

### Run 46 - high ambient, small separation

Secondary starts right against the primary envelope: a0 = 49 R_sun.

Use ambient profile as the CE fiducial run (run 143): rho amb = 6.67e-9 g/cc, p_amb = 1e5 dyn/cm2

This and the first run will compare the effect of ambient material over the jet. Also a comparison with the CE fiducial run.

density-edgeon-throughparticles

### Run 47 - low ambient, large separation

Start with the secondary farther out: a0 = 98 R_sun (meant to do 72.5 R_sun, but I messed up), and use lower ambient profile.

density-edgeon-throughparticles_frame16

## Issue: unstable primary, 2 tests without jet

Secondary starts right against the primary envelope: a0 = 49 R_sun. Use lower ambient density and pressure profile: rho_amb = 1e-10 g/cc, p_amb = 1961 dyn/cm2 .

### test 1 - low ambient, small separation

no mass in the secondary (M2=0),

but the two particles are given the same initial velocity as the orbiting case.

See the primary blows up in the opposite direction of its motion: density-faceon

### test 2 - low ambient, small separation

with secondary (M2=0.978 Msun),

this is like Run 45 without the jet.

Same effect observed in the orbiting case: density-faceon

a zoomed-in view: density-faceon-zoom

# Update 7/1

## Charge Exchange

Modified to track stellar and planetary ionization rates separately. I believe this is the most efficient solution to dividing the ionization and recombination properly among the species, since a single rate can't account for the difference in proportion between the recombination and ionization rates.

rates(i,j,k,iPlanetIoniz) = f*(1d0-emtau)/dx*(q(iHcold)/q(iH)) rates(i,j,k,iStellarIoniz) = f*(1d0-emtau)/dx*(q(iHhot)/q(iH)) rates(i,j,k,iPlanetIoniz) = rates(i,j,k,iPlanetIoniz) - RecombinationRate(T, ne*(1-ambFrac), q(iHIIcold)*(1-ambFrac)) rates(i,j,k,iStellarIoniz) = rates(i,j,k,iStellarIoniz) - RecombinationRate(T, ne*(1-ambFrac), q(iHIIhot)*(1-ambFrac))

### Non-rotating timing on BlueHive

Need to get the data from Baowei. Then will do this after AMR tests are done.

## AMR line transfer

The first timing test (projected-grid, for comparison) is queued. I'm using a single frame (220→221) of the steady state of HD209458b, no radiation pressure.

### Tracking rays

Lowest testing priority right now, but ready once other tests are complete. Contemplating how best to implement plotter for output in Mathematica/Matlab/python.

## Misc

Working through pre-req to fall course. My website is up.