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.

\frac{dN}{dt}=\alpha (T,E)N(1-\frac{N}{N_{max}})
\alpha (T,E) = \alpha_{birth} (T,E) - \alpha_{death} (T,E)=E(B_{0}e^{-(\frac{T-T_{max}}{\Delta t})^{2}}-D_{0}e^{(\frac{T-T_{max}}{\Delta t})^{2}})

• \alpha (T,E) = relative growth rate
  • T = temperature
  • E = technological efficiency (ie: measurement of how much technology a civilization has)
• N = current population
• N_{max} = maximum population \approx 20\ billion\ ppl
• T_{max} = global temperature resulting in maximum growth rates \approx 75^{o}F \approx 297K
• B_{0} = initial growth rate \approx 1/20
• D_{0} = initial death rate \approx 1/70
  • assuming a 70 year avg lifespan
• \Delta t = width of gaussian (or exponential ramping distance) \approx 10

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.

C\frac{\partial T}{\partial t}= S(1-a)-IR+\frac{1}{cos(L)}\frac{\partial}{\partial L}(Dcos(L)\frac{\partial T}{\partial L})

• 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 = 2.1*10^{8} \frac{J}{m^{2}K}
• D = diffusive parameter describing efficiency of energy transport = 0.38 \frac{W}{m^{2}K}
• 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

\frac{dP_{co2}}{dt}=C_{o}EN-T_{r}

• C_{o} = annual average pco2 contributions from 1 million people \approx 3.459*10^{-4} \frac{ppm\ from\ million\ ppl}{yr}
• T_{r} = constant reduction of CO2 by trees = 1 \frac{ppm}{yr}

Results (E=1)

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

1. a higher peak growth rate
   • the more technology you have, the lower your infantile death rates becomes, which raises the relative growth rate
2. 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

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⁰) 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.

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

• Minimum Temperature: 70^{o}F=21^{o}C
• Maximum Temperature: 140^{o}F=60^{o}C

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

\qquad \qquad \qquad \qquad \qquad \qquad \quad \qquad \frac{dN}{dt}=\alpha (T,E)N(1-\frac{N}{N_{max}})

• N = current population
• N_{max} = maximum population
• \alpha (T,E) = relative growth rate
  • T = temperature
  • E = technological efficiency (ie: measurement of how much technology a civilization has)

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:

\quad \alpha_{net} (T,E) = \alpha_{birth} (T,E) - \alpha_{death} (T,E)=E(B_{0}e^{-(\frac{T-T_{max}}{\Delta t})^{2}}-D_{0}e^{(\frac{T-T_{max}}{\Delta t})^{2}})

• T_{max} = global temperature resulting in maximum growth rates \approx 98^{o}F
• B_{0} = initial growth rate \approx 1/40
• D_{0} = initial death rate \approx 1/70
  • assuming a 70 year avg lifespan
• \Delta t = width of gaussian (or exponential ramping distance)
  • \approx 33\ for\ T<T_{max}
  • \approx 50\ for\ T\geq T_{max}

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

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

C\frac{\partial T(\ell ,t)}{\partial t}= \bar{S}(\ell, L)(1-\alpha (T))-IR+\frac{1}{cos(\ell)}\frac{\partial}{\partial \ell}(Dcos(\ell)\frac{\partial T(\ell ,t)}{\partial \ell})

• 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

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

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:

\qquad \qquad \qquad \qquad \qquad \qquad S_{sun}=\frac{L}{4\pi d^{2}}=Solar\ Constant\ \Longrightarrow S_{rel} \varpropto \frac{1}{d^{2}}

• The Habitable ("Goldilocks") Zone is defined to be any orbital distance that can result in liquid water:
  • Water Freezes: 0^{o}C=273.15K (blue dotted line)
  • Water Boils: 100^{o}C=373.15K (red dotted line)
• I then run our uncoupled EBM, with default, deterministic, Earth-like inputs, for three years, and then plot the final global temperatures, discretized by latitude
  • 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

1. Add complexity to the population logistic growth equation
2. Then couple the updated EBM to the updated population model
3. Explore Results!

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:

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?

• 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.

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 .

​Movies and frames

​density-edgeon-throughparticles

​density-faceon-box

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.

​Movies and frames

​density-edgeon-throughparticles

​density-faceon-box

​Temp-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.

​Movies and franes

​density-edgeon-throughparticles_frame16

​density-faceon-box

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.

​Movies and frames

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.

​Movies and frames

Same effect observed in the orbiting case: ​density-faceon

a zoomed-in view: ​density-faceon-zoom

Parameters

• R_p = 1000 km
• v_w = 500 km/s
• \beta = 1.85
• M = 6.119

Define critical resistivity

\eta_c = v_w R_p = 5 \times 10^5 km^2/s = 10^11 m^2/s = 10^5 \Omega m

Table of Runs

• Resistivity = 0 and 10^6 (1/6 of critical)
• Magnetic field direction = By, Bz, none
• Planet absorption coefficient = -1, 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.