Greenhouse model#

Task 1: Plug Eq. (7) into Eq. (6) and solve for the radiative equilibrium suface temperature \(T_e\)

# Solve for the radiative equilibrium temperature Te
# Put you own code here
sigma = 5.67e-8
Q = 342
albedo = 0.3

Te = (((1-albedo)*Q)/sigma)**(1/4)

print('Radiative equilibrium temperature: {:.2f}'.format(Te))

Radiative equilibrium temperature: 254.91

Task 2: Where in the atmosphere do we find \(T_e\)?

import numpy as np
import matplotlib.pyplot as plt
import xarray as xr

## NCAR data
url = 'air.mon.ltm.1981-2010.nc'
nc = xr.open_dataset('./files/air.mon.ltm.1981-2010.nc', use_cftime = True)
nc

<xarray.Dataset>
Dimensions:             (level: 17, lat: 73, lon: 144, time: 12, nbnds: 2)
Coordinates:
  * level               (level) float32 1e+03 925.0 850.0 ... 30.0 20.0 10.0
  * lat                 (lat) float32 90.0 87.5 85.0 82.5 ... -85.0 -87.5 -90.0
  * lon                 (lon) float32 0.0 2.5 5.0 7.5 ... 352.5 355.0 357.5
  * time                (time) object 0001-01-01 00:00:00 ... 0001-12-01 00:0...
Dimensions without coordinates: nbnds
Data variables:
    climatology_bounds  (time, nbnds) object 1981-01-01 00:00:00 ... 2010-12-...
    air                 (time, level, lat, lon) float32 ...
    valid_yr_count      (time, level, lat, lon) float32 ...
Attributes:
    description:                     Data from NCEP initialized reanalysis (4...
    platform:                       Model
    Conventions:                    COARDS
    not_missing_threshold_percent:  minimum 3% values input to have non-missi...
    history:                        Created 2011/07/12 by doMonthLTM\nConvert...
    title:                          monthly ltm air from the NCEP Reanalysis
    dataset_title:                  NCEP-NCAR Reanalysis 1
    References:                     http://www.psl.noaa.gov/data/gridded/data...
#  Take global, annual average and convert to Kelvin
nc.air.sel(time='0001-06-01 00:00:00',level=1000).plot(figsize=(20,10))

<matplotlib.collections.QuadMesh at 0x11fa27940>
../_images/nb_greenhouse_6_1.png 
#  a "quick and dirty" visualization of the data
weights = np.cos(np.deg2rad(nc.lat))
air_weighted = nc.air.weighted(weights)

weighted_mean = air_weighted.mean(("lat",'lon','time'))
weighted_mean.plot()

[<matplotlib.lines.Line2D at 0x12c0a7130>]
../_images/nb_greenhouse_7_1.png 
# Create the skew-plot with the metpy module
from metpy.plots import SkewT

fig = plt.figure(figsize=(10,7))

skew = SkewT(fig, rotation=30)
skew.plot(weighted_mean.level, weighted_mean, color='black', linestyle='-')

skew.plot_dry_adiabats(linewidth=0.5)
skew.plot_moist_adiabats(linewidth=0.5)

<matplotlib.collections.LineCollection at 0x28562e620>
../_images/nb_greenhouse_9_1.png

Task 3: What is the surface temperature with the single layer model?

# Solve for the atmospheric surface temperature

# Calc surface temperature Ts
# Put you own code here
Ts = 2**(1/4) * Te

print('Surface temperature: {:.2f}'.format(Ts))

Surface temperature: 303.14

Why does the model overestimate the surface temperature?

Task 5: Write a Python function for \(OLR = U_2 = (1-\epsilon)^2 \sigma T_s^4 + \epsilon(1-\epsilon)\sigma T_0^4 + \epsilon \sigma T_1^4\)

def two_layer_model(Ts, T0, T1, epsilon):
    return ((1-epsilon)**2)*sigma*Ts**4 + epsilon*(1-epsilon)*sigma*T0**4 \
+ epsilon*sigma*T1**4

Task 7: We will tune our model so that it reproduces the observed global mean OLR given observed global mean temperatures. Determine the temperatures for the two-layer model from the following sounding

alt text

Task 8: Find graphically the best fit value of \(\epsilon\)

import numpy as np
import matplotlib.pyplot as plt

import numpy as np
import matplotlib.pyplot as plt

# Assignments
OLR = []        # initialize array
epsilons = []   # initialize array
OLR_obs = 238.5 # observed outgoing long-wave radiation

# Optimize epsilon
# We define a range from 0 to 1 with a 0.01 step and calculate the OLR for each of these epsilon values
for eps in np.arange(0, 1, 0.01):
    OLR.append(OLR_obs - two_layer_model(288, 275, 230, eps))
    # Store the results in the epsilon-array
    epsilons.append(eps)   
    
def find_nearest(array, value):
    """ 
    Auxiliary function to find index of closest value in an array 
    array :: input array
    value :: the value to find in the array
    """
    # This searches the minimum between the values and all values in an array
    # Basically, we enumerate over the array. The enumerator iterates over the array and returns a 
    # tuple (index, value) for each element. We take the value (x[1]) from this tuple and substract the value
    # we are searching for. The element which has the smallest difference is what we are looking for. We finally
    # return the index of this value.
    idx,val = min(enumerate(array), key=lambda x: abs(x[1]-value))
    return idx

# Now, find the closest value to the observed OLR using the previously defined function
idx = find_nearest(OLR, 0)

# Save the optimized epsilon in the epsilons-array
epsilon = epsilons[idx]

# Plot the results
print('The optimized transmissivity is: {:.2f}'.format(epsilons[idx]))
plt.figure(figsize=(15,8))
plt.plot(epsilons,OLR);
plt.scatter(epsilons[idx], OLR[idx], s=50,  color='r')
plt.hlines(0,0,1,linestyles='dotted',color='gray');
    

The optimized transmissivity is: 0.59
../_images/nb_greenhouse_18_1.png 
two_layer_model(288, 275, 230, 0.59)

237.629744608674

Task 9: Write a Python function to calculate each term in the OLR. Plug-in the observed temperatures and the tuned value for epsilon.

def two_layer_terms(Ts, T1, T0, epsilon):
    return ( ((1-epsilon)**2)*sigma*Ts**4, \
            epsilon*(1-epsilon)*sigma*T0**4, \
            epsilon*sigma*T1**4)

# Calculate terms
term1, term2, term3 = two_layer_terms(288, 275, 230, 0.59)

print('Term Ts: {:.2f}'.format(term1))
print('Term T0: {:.2f}'.format(term2))
print('Term T1: {:.2f}'.format(term3))

Term Ts: 65.57
Term T0: 38.38
Term T1: 191.32

Task 10: Changing the level of emission by adding absorbers, e.g. by 10 %. Suppose further that this increase happens abruptly so that there is no time for the temperatures to respond to this change. We hold the temperatures fixed in the column and ask how the radiative fluxes change.

Which terms in the OLR go up and which go down?

# Make simulation here

Task 11: Calculate the radiative forcing for the previous simulation

# Calculate radiative forcing

Task 12: What is the greenhouse effect for an isothermal atmosphere?

# Make simulation here

Task 13: For a more realistic example of radiative forcing due to an increase in greenhouse absorbers, we use our observed temperatures and the tuned value for epsilon. Assume an increase of epsilon by 2 %.

# Make simulation here