Grey radiation modeling with climlab

8. Grey radiation modeling with climlab#

This notebook is part of The Climate Laboratory by Brian E. J. Rose, University at Albany.

1. Introducing `climlab`#

About climlab#

climlab is a specialized python package for process-oriented climate modeling.

It is based on a very general concept of a model as a collection of individual, interacting processes. climlab defines a base class called Process, which can contain an arbitrarily complex tree of sub-processes (each also some sub-class of Process). Every climate process (radiative, dynamical, physical, turbulent, convective, chemical, etc.) can be simulated as a stand-alone process model given appropriate input, or as a sub-process of a more complex model. New classes of model can easily be defined and run interactively by putting together an appropriate collection of sub-processes.

climlab is an open-source community project. The latest code can always be found on github:

climlab/climlab

Installing climlab#

If you’ve followed these instructions from the Climate Laboratory book, then you should be all set – climlab is automatically installed as part of the suite of tools used in The Climate Laboratory.

If you are maintaining your own Python installation (e.g. on a personal laptop), you can always install climlab by doing

conda install -c conda-forge climlab

%matplotlib inline
import numpy as np
import matplotlib.pyplot as plt
import xarray as xr
from numpy import cos, deg2rad, log
import climlab

4. A 30-layer model using the observed temperatures#

#  initialize a grey radiation model with 30 levels
col = climlab.GreyRadiationModel()
print(col)

climlab Process of type <class 'climlab.model.column.GreyRadiationModel'>. 
State variables and domain shapes: 
  Ts: (1,) 
  Tatm: (30,) 
The subprocess tree: 
Untitled: <class 'climlab.model.column.GreyRadiationModel'>
   LW: <class 'climlab.radiation.greygas.GreyGas'>
   SW: <class 'climlab.radiation.greygas.GreyGasSW'>
   insolation: <class 'climlab.radiation.insolation.FixedInsolation'>

col.lev

array([ 16.66666667,  50.        ,  83.33333333, 116.66666667,
       150.        , 183.33333333, 216.66666667, 250.        ,
       283.33333333, 316.66666667, 350.        , 383.33333333,
       416.66666667, 450.        , 483.33333333, 516.66666667,
       550.        , 583.33333333, 616.66666667, 650.        ,
       683.33333333, 716.66666667, 750.        , 783.33333333,
       816.66666667, 850.        , 883.33333333, 916.66666667,
       950.        , 983.33333333])

col.lev_bounds

array([   0.        ,   33.33333333,   66.66666667,  100.        ,
        133.33333333,  166.66666667,  200.        ,  233.33333333,
        266.66666667,  300.        ,  333.33333333,  366.66666667,
        400.        ,  433.33333333,  466.66666667,  500.        ,
        533.33333333,  566.66666667,  600.        ,  633.33333333,
        666.66666667,  700.        ,  733.33333333,  766.66666667,
        800.        ,  833.33333333,  866.66666667,  900.        ,
        933.33333333,  966.66666667, 1000.        ])

# interpolate to 30 evenly spaced pressure levels
lev = col.lev
Tinterp = np.interp(lev, np.flipud(Tglobal.level), np.flipud(Tglobal))
Tinterp
#  Need to 'flipud' because the interpolation routine 
#  needs the pressure data to be in increasing order

array([224.34736633, 211.6633606 , 206.96233453, 208.29142253,
       212.58642578, 217.19396973, 221.78252157, 226.3520813 ,
       231.30422974, 236.08017476, 240.67991638, 245.279658  ,
       249.35979716, 252.92033386, 256.48087056, 259.66789246,
       262.48139954, 265.29490662, 267.81303914, 270.03579712,
       272.25855509, 274.21642049, 275.90939331, 277.60236613,
       279.29533895, 280.98831177, 282.48550415, 283.98269653,
       285.6810201 , 287.4463874 ])

# Initialize model with observed temperatures
col.Ts[:] = Tglobal[0]
col.Tatm[:] = Tinterp

# This should look just like the observations
result_list = [col.state]
name_list = ['Observed, interpolated']
plot_soundings(result_list, name_list);

../_images/228743bdbba42b49016aa682388c5a6d02086555fd11b56fd24ea054d3b337b2.png

Tune absorptivity to get observed OLR#

col.compute_diagnostics()
col.OLR

array([263.15000104])

# Need to tune absorptivity to get OLR = 238.5
epsarray = np.linspace(0.01, 0.1, 100)
OLRarray = np.zeros_like(epsarray)

for i in range(epsarray.size):
    col.subprocess['LW'].absorptivity = epsarray[i]
    col.compute_diagnostics()
    OLRarray[i] = col.OLR

plt.plot(epsarray, OLRarray)
plt.grid()
plt.xlabel('epsilon')
plt.ylabel('OLR')

Text(0, 0.5, 'OLR')

../_images/41686d8dd1771c1331e72766bcafc3c3ccd6bdd976d368e50abd7291a0496160.png

The necessary value seems to lie near 0.055 or so.

We can be more precise with a numerical root-finder.

def OLRanom(eps):
    col.subprocess['LW'].absorptivity = eps
    col.compute_diagnostics()
    return col.OLR - 238.5

# Use numerical root-finding to get the equilibria
from scipy.optimize import brentq
# brentq is a root-finding function
#  Need to give it a function and two end-points
#  It will look for a zero of the function between those end-points
eps = brentq(OLRanom, 0.01, 0.1)
print( eps)

0.05369073166506245

col.subprocess.LW.absorptivity = eps
col.subprocess.LW.absorptivity

Field([0.05369073, 0.05369073, 0.05369073, 0.05369073, 0.05369073,
       0.05369073, 0.05369073, 0.05369073, 0.05369073, 0.05369073,
       0.05369073, 0.05369073, 0.05369073, 0.05369073, 0.05369073,
       0.05369073, 0.05369073, 0.05369073, 0.05369073, 0.05369073,
       0.05369073, 0.05369073, 0.05369073, 0.05369073, 0.05369073,
       0.05369073, 0.05369073, 0.05369073, 0.05369073, 0.05369073])

col.compute_diagnostics()
col.OLR

array([238.5])

5. Radiative forcing in the 30-layer model#

Let’s compute radiative forcing for a 2% increase in absorptivity.

#  clone our model using a built-in climlab function
col2 = climlab.process_like(col)
print(col2)

climlab Process of type <class 'climlab.model.column.GreyRadiationModel'>. 
State variables and domain shapes: 
  Ts: (1,) 
  Tatm: (30,) 
The subprocess tree: 
Untitled: <class 'climlab.model.column.GreyRadiationModel'>
   LW: <class 'climlab.radiation.greygas.GreyGas'>
   SW: <class 'climlab.radiation.greygas.GreyGasSW'>
   insolation: <class 'climlab.radiation.insolation.FixedInsolation'>

col2.subprocess['LW'].absorptivity *= 1.02
col2.subprocess['LW'].absorptivity

Field([0.05476455, 0.05476455, 0.05476455, 0.05476455, 0.05476455,
       0.05476455, 0.05476455, 0.05476455, 0.05476455, 0.05476455,
       0.05476455, 0.05476455, 0.05476455, 0.05476455, 0.05476455,
       0.05476455, 0.05476455, 0.05476455, 0.05476455, 0.05476455,
       0.05476455, 0.05476455, 0.05476455, 0.05476455, 0.05476455,
       0.05476455, 0.05476455, 0.05476455, 0.05476455, 0.05476455])

#  Radiative forcing by definition is the change in TOA radiative flux,
# HOLDING THE TEMPERATURES FIXED.
col2.Ts - col.Ts

Field([0.])

col2.Tatm - col.Tatm

Field([0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0.,
       0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0.])

col2.compute_diagnostics()
col2.OLR

array([236.65384137])

The OLR decreased after we added the extra absorbers, as we expect. Now we can calculate the Radiative Forcing:

RF = -(col2.OLR - col.OLR)
print( 'The radiative forcing is %.2f W/m2.' %RF)

The radiative forcing is 1.85 W/m2.

6. Radiative equilibrium in the 30-layer model#

re = climlab.process_like(col)

#  To get to equilibrium, we just time-step the model forward long enough
re.integrate_years(1.)

Integrating for 365 steps, 365.2422 days, or 1.0 years.
Total elapsed time is 0.9993368783782377 years.

#  Check for energy balance
print( 'The net downward radiative flux at TOA is %.4f W/m2.' %(re.ASR - re.OLR))

The net downward radiative flux at TOA is -0.0015 W/m2.

result_list.append(re.state)
name_list.append('Radiative equilibrium (grey gas)')
plot_soundings(result_list, name_list)

<AxesSubplot:xlabel='Temperature (K)', ylabel='Pressure (hPa)'>

../_images/ff1ae3f8ac0e047354b8518e721be5f9186db90ae2426457ce2917a6bc0ad26a.png

Some properties of the radiative equilibrium temperature profile:

The surface is warmer than observed.
The lower troposphere is colder than observed.
Very cold air is sitting immediately above the warm surface.
There is no tropopause, no stratosphere.

7. Radiative-Convective Equilibrium in the 30-layer model#

We recognize that the large drop in temperature just above the surface is unphysical. Parcels of air in direct contact with the ground will be warmed by mechansisms other than radiative transfer.

These warm air parcels will then become buoyant, and will convect upward, mixing their heat content with the environment.

We parameterize the statistical effects of this mixing through a convective adjustment.

At each timestep, our model checks for any locations at which the lapse rate exceeds some threshold. Unstable layers are removed through an energy-conserving mixing formula.

This process is assumed to be fast relative to radiative heating. In the model, it is instantaneous.

Add the convective adjustment as an additional subprocess#

#  Here is the existing model
print(re)

climlab Process of type <class 'climlab.model.column.GreyRadiationModel'>. 
State variables and domain shapes: 
  Ts: (1,) 
  Tatm: (30,) 
The subprocess tree: 
Untitled: <class 'climlab.model.column.GreyRadiationModel'>
   LW: <class 'climlab.radiation.greygas.GreyGas'>
   SW: <class 'climlab.radiation.greygas.GreyGasSW'>
   insolation: <class 'climlab.radiation.insolation.FixedInsolation'>

#  First we make a new clone
rce = climlab.process_like(re)
#  Then create a new ConvectiveAdjustment process
conv = climlab.convection.ConvectiveAdjustment(state=rce.state, 
                                               adj_lapse_rate=6.)
#  And add it to our model
rce.add_subprocess('Convective Adjustment', conv)
print(rce)

climlab Process of type <class 'climlab.model.column.GreyRadiationModel'>. 
State variables and domain shapes: 
  Ts: (1,) 
  Tatm: (30,) 
The subprocess tree: 
Untitled: <class 'climlab.model.column.GreyRadiationModel'>
   LW: <class 'climlab.radiation.greygas.GreyGas'>
   SW: <class 'climlab.radiation.greygas.GreyGasSW'>
   insolation: <class 'climlab.radiation.insolation.FixedInsolation'>
   Convective Adjustment: <class 'climlab.convection.convadj.ConvectiveAdjustment'>

This model is exactly like our previous models, except for one additional subprocess called Convective Adjustment.

We passed a parameter adj_lapse_rate (in K / km) that sets the neutrally stable lapse rate – in this case, 6 K / km.

This number is chosed to very loosely represent the net effect of moist convection.

#  Run out to equilibrium
rce.integrate_years(1.)

Integrating for 365 steps, 365.2422 days, or 1.0 years.
Total elapsed time is 1.9986737567564754 years.

#  Check for energy balance
rce.ASR - rce.OLR

Field([0.0007796])

result_list.append(rce.state)
name_list.append('Radiatve-Convective equilibrium (grey gas)')

plot_soundings(result_list, name_list)

<AxesSubplot:xlabel='Temperature (K)', ylabel='Pressure (hPa)'>

../_images/1e9135b03d38533f0f3b94ad46997337dada5188c3e64afbdc61755baf8a5fc7.png

Introducing convective adjustment into the model cools the surface quite a bit (compared to Radiative Equilibrium, in green here) – and warms the lower troposphere. It gives us a MUCH better fit to observations.

But of course we still have no stratosphere.

Credits#

This notebook is part of The Climate Laboratory, an open-source textbook developed and maintained by Brian E. J. Rose, University at Albany.

It is licensed for free and open consumption under the Creative Commons Attribution 4.0 International (CC BY 4.0) license.

Development of these notes and the climlab software is partially supported by the National Science Foundation under award AGS-1455071 to Brian Rose. Any opinions, findings, conclusions or recommendations expressed here are mine and do not necessarily reflect the views of the National Science Foundation.

Grey radiation modeling with climlab

Contents

8. Grey radiation modeling with climlab#

1. Introducing `climlab`#

About climlab#

Installing climlab#

2. Using `climlab` to implement the two-layer leaky greenhouse model#

Validation#

Initialize a model in `climlab`#

3. The observed annual, global mean temperature profile#

4. A 30-layer model using the observed temperatures#

Tune absorptivity to get observed OLR#

5. Radiative forcing in the 30-layer model#

6. Radiative equilibrium in the 30-layer model#

7. Radiative-Convective Equilibrium in the 30-layer model#

Add the convective adjustment as an additional subprocess#

8. Putting stratospheric ozone in the grey-gas model#

Take the global average of the ozone climatology, and plot it as a function of pressure (or height)#

The take home message#

Credits#

Grey radiation modeling with climlab

Contents

8. Grey radiation modeling with climlab#

1. Introducing climlab#

About climlab#

Installing climlab#

2. Using climlab to implement the two-layer leaky greenhouse model#

Validation#

Initialize a model in climlab#

3. The observed annual, global mean temperature profile#

4. A 30-layer model using the observed temperatures#

Tune absorptivity to get observed OLR#

5. Radiative forcing in the 30-layer model#

6. Radiative equilibrium in the 30-layer model#

7. Radiative-Convective Equilibrium in the 30-layer model#

Add the convective adjustment as an additional subprocess#

8. Putting stratospheric ozone in the grey-gas model#

Take the global average of the ozone climatology, and plot it as a function of pressure (or height)#

The take home message#

Credits#

1. Introducing `climlab`#

2. Using `climlab` to implement the two-layer leaky greenhouse model#

Initialize a model in `climlab`#