Running the GEOCARB model from R

I have provided functions for running the GEOCARB model from R:

To run the model:

run_geocarb(co2_spike, filename, degas_spinup, degas_sim,
            plants_spinup, plants_sim, land_area_spinup, land_area_sim,
            delta_t2x, million_years_ago, mean_latitude_continents)

You need to specify co2_spike (the spike in CO₂ at time zero, measured in billions of tons of carbon).

The other parameters will take default values if you don’t specify them, but you can override those defaults by giving the parameters a value.

The arguments to the function are:

filename

an optional file to save the results of the run to. You can read them back

in using the read_geocarb() function: r run_geocarb(spike = 1000, filename = "test_run.txt") data = read_geocarb("test_run.txt")

degas_spinup and degas_sim

the rates of CO₂ degassing from volcanoes for the spinup and simulation

phases, in trillions of molecules per year.

plants_spinup and plants_sim

TRUE or FALSE values for whether to include the role of plants in

weathering (their roots speed up weathering by making soil more permeable and by releasing CO₂ into the soil), and land_area is the total area of dry land, relative to today.

land_area_spinup and land_area_sim

The amount of land area, compared to today (1.0 means the same amount of

land as today).

delta_t2x

The climate sensitivity (the amount warming for each time CO₂ is

doubled), in degrees Celsius.

million_years_ago

Simulate past climates when the sun was not as bright as today.

The value of this variable is how many million years ago the year zero of the simulation should be. This is not currently working because of a bug in the web version of GEOCARB.

mean_latitude_continents

The mean latitude, in degrees, of the continents.

The default values are: degas = 7.5, plants = TRUE, and land_area = 1 for both the spinup and the simulation.

The default value for delta_t2x is 3.0, million_years_ago is 0, and mean_latitude_continents is 30, which corresponds to today’s world.

mean_latitude_continents and land_area allow you to explore conditions in earth’s past, where the continents had different areas and were located in different parts of the world.

million_years_ago is meant to allow you to explore how the silicate weathering thermostat worked in earth’s past, when the sun was a lot less intense than it is today. However, this part of the model is not working now.

After you run run_geocarb, you would read the data in with read_geocarb(filename). This function will return a data frame with the columns year, co2_total, co2_atmos, alkalinity_ocean, delta_13C_ocean, delta_13C_atmos, carbonate_ocean, carbonate_weathering, silicate_weathering, total_weathering, carbon_burial, degassing_rate, temp_atmos, and temp_ocean.

Refresher on Plotting Several Variables

You may want to go back to the documentation for Lab #2 and refresh your memory about the pivot_longer() function for manipulating data frames and tibbles, and the different ways we can use ggplot to plot several variables on the same plot.

Suppose you have a tibble with columns time, foo and bar, as shown below:

df = tibble(time = seq(100), foo = -1 + 0.1 * time - 0.001 * time^2, 
            bar = sin(time / 10))
kable(head(df), digits = 2)

Now, suppose you want to plot foo and bar on the same graph. You can do

ggplot(df, aes(x = time)) + geom_line(aes(y = foo), color = "darkred") +
  geom_line(aes(y = bar), color = "darkblue")

But it is more elegant to write

df_tidy = pivot_longer(df, cols = -time, names_to = "variable", 
                       values_to = "value")
kable(head(df_tidy))

Now you can plot this:

ggplot(df_tidy, aes(x = time, y = value, color = variable)) + geom_line() +
         scale_color_manual(values = c(foo = "darkred", bar = "darkblue"))

And you can put this all together in a single expression using pipes:

pivot_longer(df, cols = -time, names_to = "variable", 
             values_to = "value") %>%
  ggplot(aes(x = time, y = value, color = variable)) + geom_line() +
  scale_color_manual(values = c(foo = "darkred", bar = "darkblue"))

Modifying the Axes of a Plot

Sometimes you have a lot of data and you just want to plot a small part of it. Consider the following GEOCARB model run:

geocarb_data = run_geocarb(1000)

## Rows: 140 Columns: 14
## -- Column specification --------------------------------------------------------
## Delimiter: ","
## dbl (14): year, tco2, alk, d13Cocn, pCO2, d13Catm, CO3, WeatC, WeatS, TotW, ...
## 
## i Use `spec()` to retrieve the full column specification for this data.
## i Specify the column types or set `show_col_types = FALSE` to quiet this message.

ggplot(geocarb_data, aes(x = year, y = co2_atmos)) + geom_line() +
  labs(x = "Year", y = "CO2")

This shows us all 2 million years of the model run, but we can’t see the detail of what’s happening near year zero. There are several ways we can zoom our plot in to look only at the region near year zero:

1. Use the xlim and ylim functions to set limits:

   ggplot(geocarb_data, aes(x = year, y = co2_atmos)) + geom_line() + 
     xlim(-500, 1000) + labs(x = "Year", y = "CO2")

   

   If you only want to change one limit, and leave the other at its default, you can put NA for the limit you want to leave alone:

   ggplot(geocarb_data, aes(x = year, y = co2_atmos)) + geom_line() + 
     xlim(NA, 1000) + labs(x = "Year", y = "CO2")

   

2. Use the scale_x_continuous and scale_y_continuous functions

   ggplot(geocarb_data, aes(x = year, y = co2_atmos)) + geom_line() + 
     scale_x_continuous(limits = c(NA, 1E4), labels = label_comma()) +
     labs(x = "Year", y = "CO2")

   

   Using the scale_x_continuous and scale_y_continuous functions lets you also modify the way numbers are formatted on the axis. Here, I used the label_comma() function to insert commas in the thousands and millions places. Other label commands include label_percent. You can read more about these at the web page for the scales package.

3. Another approach to limiting the extent of the plot is to filter the data before you call ggplot

   geocarb_data %>% filter(year >= -500, year <= 1000) %>%
     ggplot(aes(x = year, y = co2_atmos)) + geom_line() + 
     labs(x = "Year", y = "CO2")

   

Exercise 1: Weathering as a Function of CO₂

In the steady state, the rate of weathering must balance the rate of CO₂ degassing from the Earth, from volcanoes and deep-sea vents.

Write this exercise up as a discussion of what happens when the rate of volcanic degassing changes. This rate has changed many times in Earth’s history.

• Discuss how CO₂ and temperature change both in the first thousand years after the rate of degassing changes, and also in the longer term, over the course of one or two million years.
• Discuss what causes the amount of CO₂ in the atmosphere to stabilize after the degassing rate changes.
• Discuss how the size of the change in degassing rage affects the amount of change in CO₂ and temperature between the original climate and where they finally stabilize with the new degassing rate.
• Explain the role of the silicate weathering thermostat in stabilizing the amount of CO₂, and what determines the final stable value of CO₂.

Exercise 2: The long-term fate of fossil fuel CO₂

Write this exercise up as a discussion of what happens if 2000 billion tons of carbon is released into the atmosphere as CO₂. What do you learn from GEOCARB about where that CO₂ ends up and how the earth removes it from the atmosphere. Discuss how long the removal takes and what the implications are for how we should think about CO₂ in comparison to other kinds of pollution.

Exercise 3 (Graduate Students Only): How the Land Plants Changed the Carbon Cycle

The roots of plants accelerate weathering by two processes: First, as they grow, they open up the soil, making it more permeable to air and water. Second, the roots pump CO₂ down into the soil.

Write this exercise up as a report on the effects plants have on atmospheric CO₂ concentrations. If you turn off the plants during the spinup and then turn them on during the simulation, this simulates vascular land plants (plants with roots, stems, etc.) suddenly appearing in a world where they did not previously exist. How would this have affected the global carbon cycle and the composition of the atmosphere?