Computational Guided Inquiry for Modeling Earth's Climate (Neshyba, 2023)

Diagnostic Functions

Three types of modeled variable

Here we are laying the groundwork for a revision of Cambio that provides little more organizational structure to our model. To do this, it's useful to distinguish between three different types of variables in climate modeling:

1. Parametric variables are variables that define a given climate model, such as the value of the climate sensitivity parameter and flux constants, that don't change over the course of a climate simulation.
2. Prognostic variables are variables that are integral to, and predicted by, the algorithm that move the model forward in time. In our case, prognostic variables are concentrations of carbon in the carbon reservoirs (atmosphere and oceans), and the fluxes between these reservoirs.
3. Diagnostic variables are values that can be calculated from parametric and prognostic variables, which we might be interested in, but which are not otherwise essential to the model. Examples include the global temperature and ocean pH -- important for the health of the planet, but not needed by the Euler algorithm to calculate concentrations and fluxes as the model progresses from one year to the next.

We should add that it's possible to imagine a diagnostic variable becoming a prognostic one. How? What if we wanted (in a later code revision) to implement a feedback loop in which the global temperature is not only affects th ocean-to-atmosphere flux, but is also influenced by that flux? In that case, we'd say that a feedback loop had converted temperature from diagnostic to prognostic status.

Using dictionaries to structure the code

We'll start this reorganization here by collecting all the parametric variables into a Python dictionary called ClimateParams.

Using Python functions to compute diagnostic variables

Going hand in hand with the structural revisions described above, it will be useful to organize how we compute diagnostic variables. We'll do this here by means of Python functions, one function for each diagnostic variable. As you'll see, some of these will be pretty trivial, like adding a constant to a temperature anomaly to compute an actual temperature. In other cases, you'll see some familiar algorithms repurposed. For example, we're going to use $\sigma_{floor}$ (which we used before to specify the atmosphere-to-land flux) to diagnose Earth's albedo! And in yet other cases, the algorithms will be entirely new.

Computing skills

1. I can translate mathematical formulas into algorithmic form.
2. I'm getting really good at creating Python functions.
3. I can create and use Python dictionaries.

import numpy as np
import pandas as pd
import matplotlib.pyplot as plt
import sys; sys.path.append('/home'); import Conn350Lib as CL

%matplotlib notebook

Describing climate parameters as a dictionary

Executing the cell below creates a dictionary of climate parameters.

# Start with an empty dictionary
ClimateParams = {}

# Preindustrial climate values
ClimateParams['preindust_C_atm'] = 615
ClimateParams['preindust_C_ocean'] = 350
ClimateParams['preindust_albedo'] = 0.3
ClimateParams['preindust_pH'] = 8.2

# Parameter for the basic sensitivity of the climate to increasing CO2
ClimateParams['climate_sensitivity'] = 3/ClimateParams['preindust_C_atm'] # IPCC: 3 degrees for doubled CO2

# Carbon flux constants
ClimateParams['k_la'] = 120
ClimateParams['k_al0'] = 113
ClimateParams['k_al1'] = 0.0114
ClimateParams['k_oa'] = 0.2
ClimateParams['k_ao'] = 0.114

# Parameter for the ocean degassing flux feedback
ClimateParams['DC'] = 0.034 # Pretty well known from physical chemistry

# Parameters for albedo feedback
ClimateParams['albedo_sensitivity'] = -100 # Based on our radiative balance sensitivity analysis
ClimateParams['albedo_transition_temperature'] = 3 # A guess at the T at which significant albedo reduction kicks in
ClimateParams['albedo_transition_interval'] = 1 # Also a guess: Over what temperature range will albedo reduction kick in?
ClimateParams['max_albedo_change_rate'] = 0.0005 # Also a guess, this controls how fast albedo can change in a year
ClimateParams['fractional_albedo_floor'] = 0.9 # Also a guess: .9 would be a 10% maximum reduction in albedo

# Parameters for the atmosphere->land flux feedback
ClimateParams['F_al_transitionT'] = 2 # Pretty much a guess: at what T will photosynthesis be impaired?
ClimateParams['F_al_transitionTinterval'] = 1 # Also a guess
ClimateParams['fractional_F_al_floor'] = 0.9 # Also a guess: .9 would be a 10% maximum reduction in F_al

# Parameter for stochastic processes
ClimateParams['Stochastic_C_atm_std_dev'] = 3

# This displays the dictionary contents
display(ClimateParams)

{'preindust_C_atm': 615,
 'preindust_C_ocean': 350,
 'preindust_albedo': 0.3,
 'preindust_pH': 8.2,
 'climate_sensitivity': 0.004878048780487805,
 'k_la': 120,
 'k_al0': 113,
 'k_al1': 0.0114,
 'k_oa': 0.2,
 'k_ao': 0.114,
 'DC': 0.034,
 'albedo_sensitivity': -100,
 'albedo_transition_temperature': 3,
 'albedo_transition_interval': 1,
 'max_albedo_change_rate': 0.0005,
 'fractional_albedo_floor': 0.9,
 'F_al_transitionT': 2,
 'F_al_transitionTinterval': 1,
 'fractional_F_al_floor': 0.9,
 'Stochastic_C_atm_std_dev': 3}

Your turn

In the cell below, add a new climate parameter called 'ScheduledFlowsFile'. You can give this parameter any string value you like; the idea, however, is that it will be the name of the file that contains your Scheduled Flows (i.e., your $\epsilon(t)$ file).

After adding this parameter, use "display(ClimateParams)" to verify that it was added properly.

# Add the of your scheduled anthropogenic emissions file to the dictionary
### BEGIN SOLUTION
ClimateParams['ScheduledFlowsFile'] = 'Peaks_in_2030.hdf5'
### END SOLUTION

# Display your updated dictionary
### BEGIN SOLUTION
display(ClimateParams)
### END SOLUTION

{'preindust_C_atm': 615,
 'preindust_C_ocean': 350,
 'preindust_albedo': 0.3,
 'preindust_pH': 8.2,
 'climate_sensitivity': 0.004878048780487805,
 'k_la': 120,
 'k_al0': 113,
 'k_al1': 0.0114,
 'k_oa': 0.2,
 'k_ao': 0.114,
 'DC': 0.034,
 'albedo_sensitivity': -100,
 'albedo_transition_temperature': 3,
 'albedo_transition_interval': 1,
 'max_albedo_change_rate': 0.0005,
 'fractional_albedo_floor': 0.9,
 'F_al_transitionT': 2,
 'F_al_transitionTinterval': 1,
 'fractional_F_al_floor': 0.9,
 'Stochastic_C_atm_std_dev': 3,
 'ScheduledFlowsFile': 'Peaks_in_2030.hdf5'}

Diagnostic algorithms in functional form

In this part, you'll be writing Python functions that calculate values of diagnostic variables. For each algorithm, we have describe the algorithm in mathematical terms. Your task is implement the algorithm as a Python function, then run the function through some benchmark tests. The first one below is done for you as a template.

Diagnosing and testing the land-to-atmosphere carbon flux algorithm

The algorithm is

$$ F_{al} = k_{la} $$

The benchmark value is

$$ F_{la} = 120 \ GtC/yr $$

Note that you'll have to extract k_{la} from the ClimateParams dictionary.

def Diagnose_F_la(ClimateParams):
    """Returns the land-to-atmosphere carbon flux from ClimateParam's k_la"""
    
    k_la = ClimateParams['k_la']
    F_la = k_la
    return F_la

# Testing the algorithm
print(Diagnose_F_la(ClimateParams))

120

Diagnosing and testing the ocean surface pH

Algorithm:

$$ pH = -log_{10}(\frac {C_{atm}}{C_{atm,pre-industrial}})+pH_{pre-industrial} $$

Benchmarking: When C_atm is the preindustrial value (615 GtC), the resulting pH should be

$$ pH = 8.2 $$

When C_atm is double the preindustrial value (1230 GtC), the resulting pH should be

$$ pH \approx 7.9 $$

For this function, you'll need to extract some variables from the ClimateParams dictionary:

• preindust_pH
• preindust_C_atm

Notice that we've added the concentration of carbon in the atmosphere to the argument list (as shown below).

Don't forget that numpy's log-base-ten function is np.log10

# Implementing the algorithm as a Python function
def Diagnose_OceanSurfacepH(C_atm,ClimateParams):
    """Returns ocean pH from the carbon amount in the atmosphere"""
    """and ClimateParams' preindust_pH and preindust_C_atm"""

    ### BEGIN SOLUTION
    # Extract the preindustrial pH from ClimateParms
    preindust_pH = ClimateParams['preindust_pH']
    
    # Extract the preindustrial carbon in the atmosphere from ClimateParams
    preindust_C_atm = ClimateParams['preindust_C_atm']
    
    # Calculate the pH according to our algorithm
    pH = -np.log10(C_atm/preindust_C_atm)+preindust_pH

    # Return the diagnosed pH
    return(pH)
    ### END SOLUTION

# Testing the algorithm for the pre-industrial atmospheric carbon amount (615 GtC), and double that (1230 GtC)
### BEGIN SOLUTION
print(Diagnose_OceanSurfacepH(615,ClimateParams))
print(Diagnose_OceanSurfacepH(1230,ClimateParams))
### END SOLUTION

8.2
7.898970004336018

Diagnosing a temperature anomaly

Algorithm:

$$ T_{anomaly} = CS\times(C_{atm}-C_{atm,pre-industrial}) $$

where $CS$ is the climate sensitivity parameter (degrees warming per GtC increase in atmospheric $CO_2$).

Benchmarking: When $C_{atm}=890 \ GtC$ (the current atmospheric carbon amount), we should get $T_{anomaly} \approx +1.3$.

Note that you'll have to extract the following parameters from ClimateParams:

• climate_sensitivity
• preindust_C_atm

# Implementing the algorithm
def Diagnose_T_anomaly(C_atm, ClimateParams):
    """Returns a temperature anomaly from the carbon amount in the atmosphere"""
    """and ClimateParams' climate_sensitivity and preindust_C_atm"""

    ### BEGIN SOLUTION
    # Extract the climate sensitivity and the preindustrial C_atm
    CS = ClimateParams['climate_sensitivity']
    preindust_C_atm = ClimateParams['preindust_C_atm']
    
    # Calculate the temperature anomaly according to our algorithm
    T_anomaly = CS*(C_atm-preindust_C_atm)

    # Return the temperature anomaly
    return(T_anomaly)
    ### END SOLUTION

# Testing the algorithm for 890 GtC)
### BEGIN SOLUTION
print(Diagnose_T_anomaly(890,ClimateParams))
### END SOLUTION

1.3414634146341464

Diagnosing the actual temperature

Algorithm: $T(actual,C) = (T_{anomaly}+14) $. This algorithm doesn't require any climate parameters, so only T_anomaly is in the argument list.

Benchmarking:

• return a temperature of 14 degrees C if the temperature anomaly is zero
• return a temperature of 15 degrees C if the temperature anomaly is +1

# Implementing the algorithm
def Diagnose_actual_temperature(T_anomaly):
    """Computes degrees C from a temperature anomaly"""

    ### BEGIN SOLUTION
    T_C = T_anomaly+14

    # Return our diagnosed temperature in degrees C
    return(T_C)
    ### END SOLUTION

# Testing the algorithm for a temperature anomaly of zero
### BEGIN SOLUTION
print(Diagnose_actual_temperature(0))
### END SOLUTION

# Testing the algorithm for a temperature anomaly of +1 degrees
### BEGIN SOLUTION
print(Diagnose_actual_temperature(1))
### END SOLUTION

14
15

Diagnosing temperature in Fahrenheit from Celsius

Algorithm: $T(actual,F) = T(actual,C) \times 9 /5 + 32$, where $T(actual,C)$ is in degrees Celsius.

Benchmark values:

• a temperature of $32 \ F$ if $T(C,actual)=0 \ C$ (freezing temperature of water)
• a temperature of $212 \ F$ if $T(C,actual)=100 \ C$ (boiling temperature of water)

def Diagnose_degreesF(T_C):
    """Returns an actual temperature in degrees F from an actual temperature in degrees C"""

    ### BEGIN SOLUTION
    # Do the conversion to F
    T_F = T_C*9/5+32

    # Return the diagnosed temperature in F
    return(T_F)
    ### END SOLUTION

# Testing the algorithm for freezing water (T = 0 degrees C)
### BEGIN SOLUTION
print(Diagnose_degreesF(0))
### END SOLUTION

# Testing the algorithm for boiling water (T = 100 degrees C)
### BEGIN SOLUTION
print(Diagnose_degreesF(100))
### END SOLUTION

32.0
212.0

Diagnose_F_ao

Algorithm: $F_{ao} = k_{ao} \times C_{atm}$.

Benchmark values:

• $F_{ao} \approx 70$ if the atmosphere contains 615 GtC
• $F_{ao} \approx 101$ if the atmosphere contains 890 GtC

Note that you'll have to extract k_ao from the ClimateParams dictionary.

def Diagnose_F_ao(C_atm, ClimateParams):
    """Returns the atmosphere-to-ocean carbon flux from the carbon amount in the atmosphere"""
    """and ClimateParams' k_ao"""

    ### BEGIN SOLUTION
    # Calculate the F_ao based on k_ao and the amount of carbon in the atmosphere
    k_ao = ClimateParams['k_ao']
    F_ao = k_ao*C_atm
    
    # Return the diagnosed flux
    return F_ao
    ### END SOLUTION

# Testing the algorithm for 615 GtC
### BEGIN SOLUTION
print(Diagnose_F_ao(615, ClimateParams))
### END SOLUTION

# Testing the algorithm for 890 GtC
### BEGIN SOLUTION
print(Diagnose_F_ao(890, ClimateParams))
### END SOLUTION

70.11
101.46000000000001

Diagnosing the ocean-to-atmosphere flux

Algorithm: $F_{oa} = k_{oa} (1+T_{anomaly} \times DC) \times C_{ocean}$ where $DC$ is the degassing coefficient for $CO_2$ in water.

Benchmark values:

• $F_{oa} = 70 \ GtC/yr$ if the temperature anomaly is zero and the ocean carbon amount is 350 GtC
• $F_{oa} \approx 93 \ GtC/yr$ if the temperature anomaly is 1 degree and the ocean carbon amount is 450 GtC

Note that you'll have to extract the following from the ClimateParams dictionary:

• k_oa
• DC

def Diagnose_F_oa(C_ocean, T_anomaly, ClimateParams):
    """Returns the ocean-to atmosphere carbon flux from the carbon amount in the ocean"""
    """and ClimateParams' k_oa and DC"""

    ### BEGIN SOLUTION
    DC = ClimateParams['DC']
    k_oa = ClimateParams['k_oa']
    F_oa = k_oa*(1+DC*T_anomaly)*C_ocean
    
    # Return the diagnosed flux
    return F_oa
    ### END SOLUTION

# When the temperature anomaly is zero and the ocean carbon amount is 350 GtC ...
### BEGIN SOLUTION
print(Diagnose_F_oa(350, 0, ClimateParams))
### END SOLUTION

# When the temperature anomaly is 1 degree and the ocean carbon amount is 450 GtC ...
### BEGIN SOLUTION
print(Diagnose_F_oa(450, 1, ClimateParams))
### END SOLUTION

70.0
93.06

Diagnosing the atmosphere-to-land flux

Algorithm: $F_{al} = k_{al0} + k_{al1} \times \sigma_{floor} \times C_{atm}$. Here, $\sigma_{floor}$ is a function that goes from 1 to a smaller value (not zero), so it controls how low the atmosphere-to-land flux can become. You invoke it using syntax like

CL.sigmafloor(T_anomaly, F_al_transitionT, F_al_transitionTinterval, floor)

Benchmark values:

• $F_{al} \approx 120 \ GtC/yr $ if the concentration of carbon in the atmosphere is 615, and the temperature anomaly is zero.
• $F_{al} \approx 123 \ GtC/yr $ if the concentration of carbon in the atmosphere is 890 GtC, and the temperature anomaly is +1 degrees.

You'll see in the cell below that three arguments are supplied as input to this diagnostic function. In addition, note that you'll have to extract several parameters from ClimateParams:

• k_al0
• k_al1
• F_al_transitionT
• F_al_transitionTinterval
• fractional_F_al_floor

def Diagnose_F_al(T_anomaly, C_atm, ClimateParams):
    """Returns the atmosphere-to-land carbon flux from the temperature anomaly, carbon amount in the atmosphere,"""
    """and ClimateParams' k_al0, k_al1, F_al_transitionT, F_al_transitionTinterval, and fractional_F_al_floor"""
    
    # Extract parameters we need from ClimateParameters, and calculate a new flux
    ### BEGIN SOLUTION
    k_al0 = ClimateParams['k_al0']
    k_al1 = ClimateParams['k_al1']
    F_al_transitionT = ClimateParams['F_al_transitionT']
    F_al_transitionTinterval = ClimateParams['F_al_transitionTinterval']
    floor = ClimateParams['fractional_F_al_floor'] 
    F_al = k_al0 + k_al1*CL.sigmafloor(T_anomaly,F_al_transitionT,F_al_transitionTinterval,floor)*C_atm
    
    # Return the diagnosed flux
    return F_al
    ### END SOLUTION

# Testing the algorithm at 615 and 890 GtC
### BEGIN SOLUTION
print(Diagnose_F_al(0, 615, ClimateParams))
print(Diagnose_F_al(1, 890, ClimateParams))


# Not part of the solution, but interesting to look at:
T_anomaly = np.linspace(0,5)
F_al = Diagnose_F_al(T_anomaly, 890, ClimateParams)
plt.figure()
plt.plot(T_anomaly,F_al)
plt.grid(True)
plt.xlabel('Temperature anomaly')
plt.ylabel('Atmosphere-to-land flux when C_atm=890 GtC')
### END SOLUTION

120.00926644390488
123.09788170907404

Text(0, 0.5, 'Atmosphere-to-land flux when C_atm=890 GtC')

Diagnose_albedo

Algorithm: $\alpha = \alpha_{pre-industrial} \times \sigma_{floor}$. Here, $\sigma_{floor}$ is the same function used by Diagnose_F_al; in this context, it controls how low the albedo can get.

Benchmark values:

• $\alpha = 0.30 $ if the temperature anomaly is zero
• $\alpha \approx 0.27 $ if the temperature anomaly is +4 degrees

Note that you'll have to extract several parameters from ClimateParams:

• albedo_transition_temperature
• albedo_transition_interval
• preindust_albedo
• fractional_albedo_floor

def Diagnose_albedo(T_anomaly, ClimateParams):
    """Returns the albedo from a temperature anomaly"""
    """and ClimateParams' albedo_transition_temperature, albedo_transition_interval, """
    """preindust_albedo, and fractional_albedo_floor"""
        
    ### BEGIN SOLUTION
    # Extract parameters we need from ClimateParameters, and calculate a new albedo
    transitionT = ClimateParams['albedo_transition_temperature'] 
    transitionTinterval = ClimateParams['albedo_transition_interval'] 
    preindust_albedo = ClimateParams['preindust_albedo']
    floor = ClimateParams['fractional_albedo_floor']
    albedo = CL.sigmafloor(T_anomaly,transitionT,transitionTinterval,floor)*preindust_albedo
    
    # Return the albedo
    return albedo
    ### END SOLUTION

# Testing the algorithm for temperature anomaly of zero
### BEGIN SOLUTION
print(Diagnose_albedo(0, ClimateParams))
### END SOLUTION

# Testing the algorithm for temperature anomaly of +4 degrees
### BEGIN SOLUTION
print(Diagnose_albedo(4, ClimateParams))
### END SOLUTION

0.2999962981627204
0.271422776195327

Diagnosing Delta_T from albedo

Algorithm: $ \Delta T_{from \ albedo} = AS \times (\alpha - \alpha_{preindust})$ where $AS$ is the albedo sensitivity parameter.

Benchmark values:

• $\Delta T_{from \ albedo} = 0$ if the albedo is the pre-industrial value, 0.3.
• $\Delta T_{from \ albedo} = +1 \ C $ if the albedo is 0.29.

Parameters you'll have to extract from ClimateParams include albedo_sensitivity and preindust_albedo.

def Diagnose_Delta_T_from_albedo(albedo,ClimateParams):
    """Returns the planetary temperature increase from the albedo"""
    """and ClimateParams' albedo_sensitivity and preindust_albedo"""
    
    ### BEGIN SOLUTION
    # Extract parameters we need and make the diagnosis
    AS = ClimateParams['albedo_sensitivity']    
    preindust_albedo = ClimateParams['preindust_albedo']
    Delta_T_from_albedo = (albedo-preindust_albedo)*AS
    
    return Delta_T_from_albedo
    ### END SOLUTION

# Testing the algorithm when the albedo = 0.3
### BEGIN SOLUTION
print(Diagnose_Delta_T_from_albedo(.3,ClimateParams))
### END SOLUTION

# Testing the algorithm when the albedo = 0.29
### BEGIN SOLUTION
print(Diagnose_Delta_T_from_albedo(.29,ClimateParams))
### END SOLUTION

-0.0
1.0000000000000009

Refresh/save/validate

Double-check everything is OK, and press the "Validate" button (as usual).

Close/submit/logout

Close, submit, and log out.