# Writing our first climate model

In class, we discussed Earth's global radiation budget, and we plotted the top-of-atmosphere radiative fluxes.  Now we will use our knowledge of energy balance to write a simple climate model!  

The model we are going to create will tell us how much global mean temperature should change for a given change in climate parameters.

First, let's look again at the plot of energy balance we've been studying.  Double-click on the image to see how it was added to the Jupyter notebook as well--this will help you on Assignment 2.



## 1. Finding parameters that fit observations

In class we discussed the Stefan-Boltzmann law of blackbody radiation, and we related the outgoing longwave radiation at the top of the atmosphere to Earth's theoretical surface temperature.  We found that blackbody radiation was not a very good fit to the observed temperature due to the influence of our greenhouse atmosphere.

To account for the atmosphere, we introduced the  **transmissivity**  of the atmosphere, $\tau$, a number less than 1 that represents the greenhouse effect of Earth's atmosphere.  

With the transmissivity included, the **OLR** or Outgoing Longwave Radiation to space is:

$$ \text{OLR} = \tau \sigma T_s^4 $$

where $\sigma$ is the Stefan-Boltzmann constant and $T_s$ is the surface temperature. 

**Exercise: Find a value for $\tau$ that fits the observations of $OLR$ and $T_s$.**


```python
## your python code here
```

From the plot above, the area-averaged incoming solar radiation or **insolation** is 341.3 W m$^{-2}$.

Let's denote this quantity by $Q$.


```python
Q = 341.3  # the insolation
```

#### Planetary albedo

Some of the incoming radiation is not absorbed at all but simply reflected back to space. Let's call this quantity $F_{reflected}$

From observations we have:


```python
Freflected = 101.9  #  reflected shortwave flux in W/m2
```


The **planetary albedo** is the fraction of $Q$ that is reflected.

We will denote the planetary albedo by $\alpha$.

Then we can compute...


```python
## compute alpha and fill in below statement
```

This means that ....% of solar radiation is reflected back to space. 

## 2. Absorbed Shortwave Radiation

The **Absorbed Shortwave Radiation** or ASR is the part of the incoming sunlight that is *not* reflected back to space, i.e. that part that is absorbed somewhere within the Earth system.

Mathematically we write

$$ \text{ASR} = Q - F_{reflected} = (1-\alpha) Q $$

**Exercise: define a new variable (not a function) called ASR_obs.  Compute ASR from the variables you have computed above.**


```python
## your code here
```

What does this value tell you about Earth's global energy balance?  Use what we discussed in class and the information above.

(Your interpretation here)

## 3. Equilibrium temperature

*This is one of the central concepts in climate modeling.*

The Earth system is in **energy balance** when energy in = energy out, i.e. when

$$ \text{ASR} = \text{OLR} $$

Questions we'd like to answer include

- What surface temperature do we need to have this balance?
- By how much would the temperature change in response to other changes in Earth system?
    - Changes in greenhouse gases
    - Changes in cloudiness
    - etc.



The equations above describe a simple climate model.  We can use this model to get an *exact solution* for the **equilibrium temperature**.

**Exercise: Write out the statement of energy balance, replacing ASR and OLR with their formulas we discussed above.**

*(You can write the statement here if you know LaTeX, or you can write it in your notebook)*

**Now: rearrange the statement to solve for the surface temperature, $T_s$.**

*(What expression do you find for $T_s$? You can write it here if you know LaTeX, or write it in your notebook.)*

To produce this expression, we assumed $ASR=OLR$.  Thus, $T_s$ is our equilibrium temperature required for energy balance, $T_{eq}$.  If the Earth had exactly this surface temperature, it would be in radiative balance with the Sun.

**Python implementation: write a reusable function that will compute the equilibrium temperature $T_{eq}$ for any values of albedo, incoming solar flux, and transmissivity.**


```python
## your reusable function here
```

Let's make your function a little bit easier to use.  Add to your function definition so that it has *default values* for each of the parameters.  [Here is a tutorial on default values.](https://www.pythontutorial.net/python-basics/python-default-parameters/)

Choose appropriate default values and assign them below.


```python
## version of your function with default values here
```

Finally, compute the equilibrium temperature with the observed values of each parameter.


```python
## your call to your function
```

How does this compare with the observed global temperature?

*(Your response here)*

## 4. A climate change scenario

Suppose that, due to global warming (changes in atmospheric composition and subsequent changes in cloudiness):

- The longwave transmissitivity decreases to $\tau = 0.57$ 
- The planetary albedo increases to $\alpha = 0.32$

What is the ***new equilibrium temperature***?


```python
Teq_new = ... ## your computation here
#  we'll show an example of formatted print output, limiting to two or one decimal places
print('The new equilibrium temperature is {:.2f} K.'.format(Teq_new))
print('The equilibrium temperature increased by about {:.1f} K.'.format(Teq_new-Teq_observed))
```

Most climate models are more complicated than this one, and solving directly for the equilibrium temperature will not be possible.

Instead, we will be able to use the model to calculate the terms in the energy budget (ASR and OLR). To prepare for that style of modelling...

**Python exercise:**
- Write reusable functions that will compute ASR and OLR for any values of albedo, transmissivity, solar flux, and temperature.
- Verify the following:
    - With the new parameter values but the old temperature $T = 288$ K, is ASR greater or lesser than OLR? 
    - Is the Earth gaining or losing energy?
    - How does your answer change if $T = 295$ K (or any other temperature greater than 291 K)?


```python
## your code here
```

## 5. Simulating changes over time

The above exercise shows us that if some properties of the climate system change in such a way that the **equilibrium temperature goes up**, then the Earth system *receives more energy from the sun than it is losing to space*. The system is **no longer in energy balance**.

The temperature must then increase to get back into balance. The increase will not happen all at once! It will take time for energy to accumulate in the climate system. We want to model this **time-dependent adjustment** of the system.

In fact almost all climate models are **time-dependent**, meaning the model calculates **time derivatives** (rates of change) of climate variables.

### An energy balance **equation**

If the Earth system is no longer in energy balance, then there must be some net change in energy within the system over time.  We will write the **total energy budget** of the Earth system as

$$
\begin{align} 
\frac{dE}{dt} &= \text{net energy flux in to system} \\
 &= \text{flux in – flux out} \\
 &= \text{ASR} - \text{OLR}
\end{align}
$$

where $E$ is the **enthalpy** or **heat content** of the total system.

We will express the budget **per unit surface area**, so each term above has units W m$^{-2}$

Note: any internal exchanges of energy between different reservoirs (e.g. between ocean, land, ice, atmosphere) do not appear in this budget – because $E$ is the sum of all reservoirs.

Also note: **This is a generically true statement.** We have just defined some terms, and made the (very good) assumption that the only significant energy sources are radiative exchanges with space.

**This equation is the starting point for _every climate model_.**

But so far, we don’t actually have a MODEL, something that will show us changes in the system over time. We just have a statement of a budget. To use this budget to make a model, we need to relate terms in the budget to state variables of the atmosphere-ocean system.

For now, the state variable we are most interested in is **temperature** – because it is directly connected to the physics of each term above.

### An energy balance **model**

If we now suppose that 

$$ E = C T_s $$

where $E$ is the total heat content as above, $T_s$ is the global mean surface temperature, and $C$ is a constant – the **effective heat capacity** of the atmosphere- ocean column.

then we can rewrite our energy budget equation as...

*(Write down the equation here if you know LaTeX, or write it in your notebook.)*





For our purposes here we are going to use a value of C equivalent to heating 100 meters of water:

$$C = c_w \rho_w H$$

where 

$c_w = 4 \times 10^3$ J kg$^{-1}$ $^\circ$C$^{-1}$ is the specific heat of water,

$\rho_w = 10^3$ kg m$^{-3}$ is the density of water, and

$H$ is an effective depth of water that is heated or cooled.


```python
c_w = 4E3  #  Specific heat of water in J/kg/K
rho_w = 1E3  #  Density of water in kg/m3
H = 100.   #  Depth of water in m
C = c_w * rho_w * H   #  Heat capacity of the model 
print('The effective heat capacity is {:.1e} J/m2/K'.format(C))
```

### Solving the energy balance model

The version of the energy budget equation you rewrote above is a first-order Ordinary Differential Equation (ODE) for $T_s$ as a function of time. It is also **our very first climate model!**

To solve it (i.e. see how $T_s$ evolves from some specified initial condition) we have two choices:

1. Solve it analytically
2. Solve it numerically

Option 1 (analytical) will usually not be possible because the equations will typically be too complex and non-linear. This is why computers are our best friends in the world of climate modeling.

HOWEVER it is often useful and instructive to simplify a model down to something that is analytically solvable when possible. Why? Two reasons:

1. Analysis will often yield a deeper understanding of the behavior of the system
2. Gives us a benchmark against which to test the results of our numerical solutions.

## Sidebar: representing time derivatives on the computer

Recall that the derivative is the **instantaneous rate of change**. It is defined as 

$$ \frac{dT}{dt} = \lim_{\Delta t\rightarrow 0}⁡ \frac{\Delta T}{\Delta t}$$

- **On the computer there is no such thing as an instantaneous change.** 
- We are always dealing with *discrete quantities*.
- So we approximate the derivative with $\Delta T/ \Delta t$. 
- So long as we take the time interval $\Delta t$ "small enough", the approximation is valid and useful.
- (The meaning of "small enough" varies widely in practice. Let's not worry about it right now)

So we can write out your model from step 6 as

$$ C  \frac{\Delta T}{\Delta t} \approx \text{ASR} - \text{OLR}$$

where $\Delta T$ is the **change in temperature predicted by our model** over a short time interval $\Delta t$.

We can now use this to make a prediction: 

Given a current temperature $T_1$ at time $t_1$, what is the temperature $T_2$ at a future time $t_2$?

We can write

$$ \Delta T = T_2-T_1 $$
$$ \Delta t = t_2-t_1 $$

and so our model says

$$ C  \frac{T_2-T_1}{\Delta t} = \text{ASR} - \text{OLR} $$

Which we can rearrange to **solve for the future temperature**:

$$ T_2 = T_1 + \frac{\Delta t}{C} \left( \text{ASR} - \text{OLR}(T_1) \right)  $$

We now have a formula with which to make our prediction!

Notice that we have written the OLR as a *function of temperature*. We will use the current temperature $T_1$ to compute the OLR, and use that OLR to determine the future temperature.

## 6. Computational solution of our Energy Balance Model

The quantity $\Delta t$ is called a **timestep**. It is the smallest time interval represented in our model.

Let's use a timestep of 1 year:


```python
dt = 60. * 60. * 24. * 365.   # one year expressed in seconds
```

### Try stepping forward one timestep


```python
# Try a single timestep, using your functions for ASR and OLR.  Input the correct parameter values: alpha=0.32, tau=0.57
T1 = 288.
T2 = T1 + dt / C * ( ASR(...) - OLR(T1...) ) ## edit here to call your functions
print(T2)
```

***Did you get a `NameError` here?***

The code above assumes that we have already defined functions `ASR()` and `OLR()`. If you named them something else in the exercise above, then the code won't work as written.  Try again, making sure to call correct function names to compute ASR and OLR.

**What happened when you stepped forward one year?  Why?**

*(Your answer here)*

**Exercise: step forward one more time step, using the T2 you computed above as your starting value.**


```python
T2 = ...
T3 = ...
print(T3)
```

But we don't want to type this out every time.  That's right: it's time to define another function.

**Exercise: define a function that will step forward one year, given an initial temperature $T_i$**


```python
def step_forward(Ti):
    ## work of the function
    return ... ## write what should be returned
```

The most effective way to write this function is to base it on what we've already done, calling other functions and variables we have already defined.  

**Why do you think we defined default values of the albedo, solar flux, and transmissivity above?  How does that help us in writing step_forward?**

### Automate the timestepping with a loop

Now let's really harness the power of the computer by **making a loop** (and storing values in arrays):


```python
import numpy as np

numsteps = 20
Tsteps = np.zeros(numsteps+1)
Years = np.zeros(numsteps+1)
Tsteps[0] = 288. 
for n in range(numsteps):
    Years[n+1] = n+1
    Tsteps[n+1] = step_forward( Tsteps[n] )
print(Tsteps)
```

What did we just do?

- Created an array of zeros
- set the initial temperature to 288 K
- repeated our time step 20 times. 
- Stored the results of each time step into the array.

#### Python fact: the `for` statement executes a statement (or series of statements) a specified number of times (a loop!)

## 7. Make a plot

Let's present our results in a more visually interpretable way, using matplotlib to create a plot.

We can easily make simple graphs with the function `plt.plot(x,y)`, where `x` and `y` are arrays of the same size. But we must import it first. 

This uses a graphics library called [matplotlib](https://matplotlib.org). This is the workhorse of scientific plotting in Python, and we will be using it all the time!

When we plan on plotting, our notebooks will usually start with this:
```
%matplotlib inline
import numpy as np
import matplotlib.pyplot as plt
```
...unless we are going to use xarray plotting only.


```python
# a special instruction for the Jupyter notebook
#   Display all plots inline in the notebook
%matplotlib inline  
#  import the plotting package
import matplotlib.pyplot as plt
```


```python
plt.plot(Years, Tsteps)
plt.xlabel('Years')
plt.ylabel('Global mean temperature (K)');
```

**Exercise: Use the matplotlib documentation, or other online resources, to change something about your plot.  Maybe you want to use a different color line, or different line style, or line width?**


```python
## Your revised version of the plot here
```

**Finally: Describe what happens in your simulation over time.**

*(Your answer here)*

This lab draws heavily from Climlab resources by Brian Rose (not linked at top to avoid ruining the surprise!)


```python

```

f9ba13dc4219396ac1de4da91f9f424ec58dbf02
ipynb
Jupyter Notebook
03-energy_balance_model.ipynb
atakoudes/climdyn-labs
81b83f6109107e9957f40f5a5090a6ce3b1c3477
MIT
03-energy_balance_model.ipynb
atakoudes/climdyn-labs
81b83f6109107e9957f40f5a5090a6ce3b1c3477
MIT
03-energy_balance_model.ipynb
atakoudes/climdyn-labs
81b83f6109107e9957f40f5a5090a6ce3b1c3477
MIT
2022-02-24T18:42:08.000Z
2022-03-17T19:04:24.000Z
Qwen/Qwen-72B
1. YES
2. YES