Tutorial 6: Radiative-Convective Equilibrium

Week 1, Day 5, Climate Modeling

Content creators: Jenna Pearson

Content reviewers: Yunlong Xu, Will Gregory, Peter Ohue, Derick Temfack, Zahra Khodakaramimaghsoud, Peizhen Yang, Younkap Nina Duplex, Ohad Zivan, Chi Zhang

Content editors: Brodie Pearson, Abigail Bodner, Ohad Zivan, Chi Zhang

Production editors: Wesley Banfield, Jenna Pearson, Chi Zhang, Ohad Zivan

Our 2023 Sponsors: NASA TOPS and Google DeepMind

Tutorial Objectives

Building on the understanding of a one-dimensional radiative balance model from previous tutorial, in this tutorial students will learn about radiative-convective-equilibrium. Much of the code shown here was taken from The Climate Laboratory by Brian Rose. Students are encouraged to visit this website for more tutorials and background on these models.

By the end of this tutorial students will be able to:

• Implement a radiative-convective equilibrium model using the python package climlab.

• Understand how this model builds off the one-dimensional radiative balance model used in the previous tutorials.

Setup

# installations ( uncomment and run this cell ONLY when using google colab or kaggle )

# note the conda install takes quite a while, but conda is REQUIRED to properly download the 
# dependencies (that are not just python packages)

# !pip install condacolab &> /dev/null           
# import condacolab
# condacolab.install()

# !mamba install -c anaconda cftime xarray numpy &> /dev/null    # for decoding time variables when opening datasets
# !mamba install -c conda-forge metpy climlab &> /dev/null

# imports

# google colab users, if you get an xarray error, please run this cell again

import xarray as xr  # used to manipulate data and open datasets
import numpy as np  # used for algebra/arrays
import urllib.request  # used to download data from the internet
import climlab  # one of the models we are using
import matplotlib.pyplot as plt  # used for plotting
import metpy  # used to make Skew T Plots of temperature and pressure
from metpy.plots import SkewT  # plotting function used widely in climate science
import pooch
import os
import tempfile
from IPython.display import HTML
from matplotlib import animation

/home/wesley/miniconda3/envs/climatematch/lib/python3.10/site-packages/climlab/convection/akmaev_adjustment.py:142: NumbaDeprecationWarning: The 'nopython' keyword argument was not supplied to the 'numba.jit' decorator. The implicit default value for this argument is currently False, but it will be changed to True in Numba 0.59.0. See https://numba.readthedocs.io/en/stable/reference/deprecation.html#deprecation-of-object-mode-fall-back-behaviour-when-using-jit for details.
  Akmaev_adjustment = jit(signature_or_function=Akmaev_adjustment)

# @title Figure settings
import ipywidgets as widgets  # interactive display

%config InlineBackend.figure_format = 'retina'
plt.style.use(
    "https://raw.githubusercontent.com/ClimateMatchAcademy/course-content/main/cma.mplstyle"
)

# @title Helper functions

def pooch_load(filelocation=None, filename=None, processor=None):
    shared_location = "/home/jovyan/shared/Data/tutorials/W1D4_ClimateModeling"  # this is different for each day
    user_temp_cache = tempfile.gettempdir()

    if os.path.exists(os.path.join(shared_location, filename)):
        file = os.path.join(shared_location, filename)
    else:
        file = pooch.retrieve(
            filelocation,
            known_hash=None,
            fname=os.path.join(user_temp_cache, filename),
            processor=processor,
        )

    return file

# @title Plotting Functions

#  make the videos at the end of the tutorial
plt.rcParams["animation.html"] = "jshtml"

# these three functions are used to make videos at the end of the tutorial


def initial_figure(model):
    with plt.ioff():  # will hide the inital figure which will plot separate from the video otherwise
        fig = plt.figure(figsize=(6, 6))
        lines = []

        skew = SkewT(fig, rotation=30)
        #  plot the observations
        skew.plot(
            Tglobal.level,
            Tglobal,
            color="black",
            linestyle="-",
            linewidth=2,
            label="Observations",
        )
        lines.append(
            skew.plot(
                model.lev,
                model.Tatm - climlab.constants.tempCtoK,
                linestyle="-",
                linewidth=2,
                color="C0",
                label="RC model (all gases)",
            )[0]
        )
        skew.ax.legend()
        skew.ax.set_ylim(1050, 10)
        skew.ax.set_xlim(-60, 75)
        # Add the relevant special lines
        skew.plot_dry_adiabats(linewidth=1.5, label="dry adiabats")
        # skew.plot_moist_adiabats(linewidth=1.5, label = 'moist adiabats')
        skew.ax.set_xlabel("Temperature ($^\circ$C)", fontsize=14)
        skew.ax.set_ylabel("Pressure (hPa)", fontsize=14)
        lines.append(
            skew.plot(
                1000,
                model.Ts - climlab.constants.tempCtoK,
                "o",
                markersize=8,
                color="C0",
            )[0]
        )

    return fig, lines


def animate(day, model, lines):
    lines[0].set_xdata(np.array(model.Tatm) - climlab.constants.tempCtoK)
    lines[1].set_xdata(np.array(model.Ts) - climlab.constants.tempCtoK)
    # lines[2].set_xdata(np.array(model.q)*1E3)
    # lines[-1].set_text('Day {}'.format(int(model.time['days_elapsed'])))
    # This is kind of a hack, but without it the initial frame doesn't appear
    if day != 0:
        model.step_forward()
    return lines


# to setup the skewT and plot observations
def make_basic_skewT():
    fig = plt.figure(figsize=(9, 9))
    skew = SkewT(fig, rotation=30)
    skew.plot(
        Tglobal.level,
        Tglobal,
        color="black",
        linestyle="-",
        linewidth=2,
        label="Observations",
    )
    skew.ax.set_ylim(1050, 10)
    skew.ax.set_xlim(-90, 45)
    # Add the relevant special lines
    # skew.plot_dry_adiabats(linewidth=1.5, label = 'dry adiabats')
    # skew.plot_moist_adiabats(linewidth=1.5, label = 'moist adiabats')
    # skew.plot_mixing_lines()
    skew.ax.legend()
    skew.ax.set_xlabel("Temperature (degC)", fontsize=14)
    skew.ax.set_ylabel("Pressure (hPa)", fontsize=14)
    return skew


# to setup the skewT and plot observations
def make_skewT():
    fig = plt.figure(figsize=(9, 9))
    skew = SkewT(fig, rotation=30)
    skew.plot(
        Tglobal.level,
        Tglobal,
        color="black",
        linestyle="-",
        linewidth=2,
        label="Observations",
    )
    skew.ax.set_ylim(1050, 10)
    skew.ax.set_xlim(-90, 45)
    # Add the relevant special lines
    skew.plot_dry_adiabats(linewidth=1.5, label="dry adiabats")
    # skew.plot_moist_adiabats(linewidth=1.5, label = 'moist adiabats')
    # skew.plot_mixing_lines()
    skew.ax.legend()
    skew.ax.set_xlabel("Temperature (degC)", fontsize=14)
    skew.ax.set_ylabel("Pressure (hPa)", fontsize=14)
    return skew


# to add a model derived profile to the skewT figure
def add_profile(skew, model, linestyle="-", color=None):
    line = skew.plot(
        model.lev,
        model.Tatm - climlab.constants.tempCtoK,
        label=model.name,
        linewidth=2,
    )[0]
    skew.plot(
        1000,
        model.Ts - climlab.constants.tempCtoK,
        "o",
        markersize=8,
        color=line.get_color(),
    )
    skew.ax.legend()

# @title Video 1: Radiative-Convective Equilibrium

from ipywidgets import widgets
from IPython.display import YouTubeVideo
from IPython.display import IFrame
from IPython.display import display


class PlayVideo(IFrame):
  def __init__(self, id, source, page=1, width=400, height=300, **kwargs):
    self.id = id
    if source == 'Bilibili':
      src = f'https://player.bilibili.com/player.html?bvid={id}&page={page}'
    elif source == 'Osf':
      src = f'https://mfr.ca-1.osf.io/render?url=https://osf.io/download/{id}/?direct%26mode=render'
    super(PlayVideo, self).__init__(src, width, height, **kwargs)


def display_videos(video_ids, W=400, H=300, fs=1):
  tab_contents = []
  for i, video_id in enumerate(video_ids):
    out = widgets.Output()
    with out:
      if video_ids[i][0] == 'Youtube':
        video = YouTubeVideo(id=video_ids[i][1], width=W,
                             height=H, fs=fs, rel=0)
        print(f'Video available at https://youtube.com/watch?v={video.id}')
      else:
        video = PlayVideo(id=video_ids[i][1], source=video_ids[i][0], width=W,
                          height=H, fs=fs, autoplay=False)
        if video_ids[i][0] == 'Bilibili':
          print(f'Video available at https://www.bilibili.com/video/{video.id}')
        elif video_ids[i][0] == 'Osf':
          print(f'Video available at https://osf.io/{video.id}')
      display(video)
    tab_contents.append(out)
  return tab_contents


video_ids = [('Youtube', 'KV198mGq1GM'), ('Bilibili', 'BV1gP411Q7f3')]
tab_contents = display_videos(video_ids, W=730, H=410)
tabs = widgets.Tab()
tabs.children = tab_contents
for i in range(len(tab_contents)):
  tabs.set_title(i, video_ids[i][0])
display(tabs)

Section 1: Reproducing Data from the Last Tutorial’s One-dimensional Radiative Equilibrium Model Using Climlab

filename_sq = "cpl_1850_f19-Q-gw-only.cam.h0.nc"
url_sq = "https://osf.io/c6q4j/download/"
ds = xr.open_dataset(
    pooch_load(filelocation=url_sq, filename=filename_sq)
)  # ds = dataset

filename_ncep_air = "air.mon.1981-2010.ltm.nc"
url_ncep_air = "https://osf.io/w6cd5/download/"
ncep_air = xr.open_dataset(
    pooch_load(filelocation=url_ncep_air, filename=filename_ncep_air)
)  # ds = dataset

# take global, annual average
weight_factor = ds.gw / ds.gw.mean(dim="lat")
Qglobal = (ds.Q * weight_factor).mean(dim=("lat", "lon", "time"))

# use 'lev=Qglobal.lev' to create an identical vertical grid to water vapor data
mystate = climlab.column_state(lev=Qglobal.lev, water_depth=2.5)

radmodel = climlab.radiation.RRTMG(
    name="Radiation (all gases)",  # give our model a name!
    state=mystate,  # give our model an initial condition!
    specific_humidity=Qglobal.values,  # tell the model how much water vapor there is
    albedo=0.25,  # this the SURFACE shortwave albedo
    timestep=climlab.constants.seconds_per_day,  # set the timestep to one day (measured in seconds)
)

# need to take the average over space and time
# the grid cells are not the same size moving towards the poles, so we weight by the cosine of latitude to compensate for this
coslat = np.cos(np.deg2rad(ncep_air.lat))
weight = coslat / coslat.mean(dim="lat")

Tglobal = (ncep_air.air * weight).mean(dim=("lat", "lon", "time"))

# take a single step forward to the diagnostics are updated and there is some energy imbalance
radmodel.step_forward()

# run the model to equilibrium (the difference between ASR and OLR is a very small number)
while np.abs(radmodel.ASR - radmodel.OLR) > 0.001:
    radmodel.step_forward()

skew = make_basic_skewT()
add_profile(skew, radmodel)
skew.ax.set_title("Pure Radiative Equilibrium", fontsize=18);

Downloading data from 'https://osf.io/c6q4j/download/' to file '/tmp/cpl_1850_f19-Q-gw-only.cam.h0.nc'.

---------------------------------------------------------------------------
KeyboardInterrupt                         Traceback (most recent call last)
Cell In[8], line 4
      1 filename_sq = "cpl_1850_f19-Q-gw-only.cam.h0.nc"
      2 url_sq = "https://osf.io/c6q4j/download/"
      3 ds = xr.open_dataset(
----> 4     pooch_load(filelocation=url_sq, filename=filename_sq)
      5 )  # ds = dataset
      7 filename_ncep_air = "air.mon.1981-2010.ltm.nc"
      8 url_ncep_air = "https://osf.io/w6cd5/download/"

Cell In[4], line 10, in pooch_load(filelocation, filename, processor)
      8     file = os.path.join(shared_location, filename)
      9 else:
---> 10     file = pooch.retrieve(
     11         filelocation,
     12         known_hash=None,
     13         fname=os.path.join(user_temp_cache, filename),
     14         processor=processor,
     15     )
     17 return file

File ~/miniconda3/envs/climatematch/lib/python3.10/site-packages/pooch/core.py:239, in retrieve(url, known_hash, fname, path, processor, downloader, progressbar)
    236 if downloader is None:
    237     downloader = choose_downloader(url, progressbar=progressbar)
--> 239 stream_download(url, full_path, known_hash, downloader, pooch=None)
    241 if known_hash is None:
    242     get_logger().info(
    243         "SHA256 hash of downloaded file: %s\n"
    244         "Use this value as the 'known_hash' argument of 'pooch.retrieve'"
   (...)
    247         file_hash(str(full_path)),
    248     )

File ~/miniconda3/envs/climatematch/lib/python3.10/site-packages/pooch/core.py:803, in stream_download(url, fname, known_hash, downloader, pooch, retry_if_failed)
    799 try:
    800     # Stream the file to a temporary so that we can safely check its
    801     # hash before overwriting the original.
    802     with temporary_file(path=str(fname.parent)) as tmp:
--> 803         downloader(url, tmp, pooch)
    804         hash_matches(tmp, known_hash, strict=True, source=str(fname.name))
    805         shutil.move(tmp, str(fname))

File ~/miniconda3/envs/climatematch/lib/python3.10/site-packages/pooch/downloaders.py:226, in HTTPDownloader.__call__(self, url, output_file, pooch, check_only)
    224     progress = self.progressbar
    225     progress.total = total
--> 226 for chunk in content:
    227     if chunk:
    228         output_file.write(chunk)

File ~/miniconda3/envs/climatematch/lib/python3.10/site-packages/requests/models.py:816, in Response.iter_content.<locals>.generate()
    814 if hasattr(self.raw, "stream"):
    815     try:
--> 816         yield from self.raw.stream(chunk_size, decode_content=True)
    817     except ProtocolError as e:
    818         raise ChunkedEncodingError(e)

File ~/miniconda3/envs/climatematch/lib/python3.10/site-packages/urllib3/response.py:628, in HTTPResponse.stream(self, amt, decode_content)
    626 else:
    627     while not is_fp_closed(self._fp):
--> 628         data = self.read(amt=amt, decode_content=decode_content)
    630         if data:
    631             yield data

File ~/miniconda3/envs/climatematch/lib/python3.10/site-packages/urllib3/response.py:567, in HTTPResponse.read(self, amt, decode_content, cache_content)
    564 fp_closed = getattr(self._fp, "closed", False)
    566 with self._error_catcher():
--> 567     data = self._fp_read(amt) if not fp_closed else b""
    568     if amt is None:
    569         flush_decoder = True

File ~/miniconda3/envs/climatematch/lib/python3.10/site-packages/urllib3/response.py:533, in HTTPResponse._fp_read(self, amt)
    530     return buffer.getvalue()
    531 else:
    532     # StringIO doesn't like amt=None
--> 533     return self._fp.read(amt) if amt is not None else self._fp.read()

File ~/miniconda3/envs/climatematch/lib/python3.10/http/client.py:466, in HTTPResponse.read(self, amt)
    463 if self.length is not None and amt > self.length:
    464     # clip the read to the "end of response"
    465     amt = self.length
--> 466 s = self.fp.read(amt)
    467 if not s and amt:
    468     # Ideally, we would raise IncompleteRead if the content-length
    469     # wasn't satisfied, but it might break compatibility.
    470     self._close_conn()

File ~/miniconda3/envs/climatematch/lib/python3.10/socket.py:705, in SocketIO.readinto(self, b)
    703 while True:
    704     try:
--> 705         return self._sock.recv_into(b)
    706     except timeout:
    707         self._timeout_occurred = True

File ~/miniconda3/envs/climatematch/lib/python3.10/ssl.py:1274, in SSLSocket.recv_into(self, buffer, nbytes, flags)
   1270     if flags != 0:
   1271         raise ValueError(
   1272           "non-zero flags not allowed in calls to recv_into() on %s" %
   1273           self.__class__)
-> 1274     return self.read(nbytes, buffer)
   1275 else:
   1276     return super().recv_into(buffer, nbytes, flags)

File ~/miniconda3/envs/climatematch/lib/python3.10/ssl.py:1130, in SSLSocket.read(self, len, buffer)
   1128 try:
   1129     if buffer is not None:
-> 1130         return self._sslobj.read(len, buffer)
   1131     else:
   1132         return self._sslobj.read(len)

KeyboardInterrupt: 

Section 2: Radiative-Convective Equilibrium

From the plot you just made, one of the largest differences between observations and the pure radiation model with all gases lies in the lower atmosphere, where the surface air temperature is 20 degrees too warm and the 200 hPa pressure surface is 40 degrees too cold. What could be the issue?

One thing we have not included in our model yet is dynamics (motion of the air). The model’s temperature profile is what’s known as statically unstable (note this definition of stability is different than that used in the previous tutorials).

Here static means not due to wind, and unstable means the atmosphere wants to adjust to a different state because the surface air is relatively light and wants to rise into upper layers of the atmosphere. As the air rises, it creates convective turbulence (similar to boiling water, where convective circulation is introduced from heating water from below). The rising air, and the resultant turbulence, mixes the atmospheric column. This mixing often occurs across the troposphere, which is roughly the lowest 10km of the atmosphere. Most of the weather we experience lies in the troposphere.

When air rises adiabatically, it expands and cools due to the lower pressure. The rate of cooling depends on whether the air is saturated with water vapor. When rising air is unsaturated, it cools following the dry adiabats. If the air saturates, it cools at a lesser rate denoted by the moist adiabats (we did not have time to discuss these moisture effects in the mini-lecture).

To identify unstable atmospheric layers, let’s take another look at the SkewT plot, but this time we will plot the dry adiabats. We can then compare the rates of cooling of our model to these adiabats.

skew = make_skewT()
add_profile(skew, radmodel)
skew.ax.set_title("Pure Radiative Equilibrium", fontsize=18)

Near the surface, the reanalysis temperature profile is steeper than the dry adiabats. In these layers vertical motion is inhibited, and the surface conditions are considered stable. However, the model profile is shallower than the dry adiabats. In the model, the surface air is statically unstable, which would lead to convective mixing if this physical process was included in the model (the model currently only includes radiative processes). In this tutorial we will see whether including convective mixing in our model can bring the model closer to the reanalysis temperature profile.

To build a radiative-convective model we can take the radiative model we have already made and couple it to a convective model. Here the term couple implies there is communication between the models such that their effects are both incorporated into the final product, which in our case is the temperature profile.

# restate the model here for ease of coding

# make a model on same vertical domain as the water vapor data
mystate = climlab.column_state(lev=Qglobal.lev, water_depth=2.5)

# create the radiation model
rad = climlab.radiation.RRTMG(
    name="Radiation (net)",
    state=mystate,
    specific_humidity=Qglobal.values,
    timestep=climlab.constants.seconds_per_day,
    albedo=0.25,  # surface albedo, tuned to give reasonable ASR for reference cloud-free model
)

# create the convection model
conv = climlab.convection.ConvectiveAdjustment(
    name="Convection",
    state=mystate,
    adj_lapse_rate=6.5,  # the adiabatic lapse rate of the atmopshere
    timestep=rad.timestep,  # same timestep as radiation model
)
# couple the two components
rcm = climlab.couple([rad, conv], name="Radiative-Convective Model")

Now let’s run the radiation part to equilibrium, which should give us the same profile as in the previous section. Once we get this radiative equilibrium profile we will add convective mixing physics to the model and continue running it until it reaches radiative-convective equilibrium.

The new model does not actually resolve the actual vertical motion and mixing that occurs in convection. Instead, the model includes a parameterization for convection which automatically mixes regions of the atmospheric column that are statically unstable.

# run JUST the radiative component to equilibrium
for n in range(1000):
    rcm.subprocess["Radiation (net)"].step_forward()

# compute diagnostics
rcm.compute_diagnostics()

#  plot the resulting profile (our initial condition once we turn on the physics)
fig, lines = initial_figure(rcm)

# this animation can take a while
animation.FuncAnimation(fig, animate, 50, fargs=(rcm, lines))

Adding convective mixing to the initially unstable temperature profile leads to an instant mixing of air throughout the lower atmosphere, moving the profile towards the observations. The balance at play is between radiative process that warm the surface and cool the troposphere (lower atmosphere) as well as convection which moves heat away from the surface, leading to a colder surface and warmer troposphere. Note the differences in surface versus tropospheric temperatures in the new versus old equilibrium profiles.

Questions 2: Climate Connection

1. The static instability was removed in the first time step! In reality, which process do you think changes temperature faster, convection or radiation?

2. What do you think the next step would be to move towards a more realistic climate model?

Coding Exercises

1. Recreate the video above except using an isothermal atmosphere (uniform temperature profile) set to the surface temprature (rcm.state.Ts). Does the equilbrium profile look any different than the one you made before?

# set initial temperature profile to be the surface temperature (isothermal)
rcm.state.Tatm[:] = rcm.state.Ts

# compute diagnostics
_ = ...

#  plot initial data
fig, lines = initial_figure(rcm)

# make animation - this animation can take a while
# animation.FuncAnimation(fig, animate, 50, fargs=(rcm, lines))

Summary

In this tutorial, you explored the concept of radiative-convective equilibrium, learning to implement a Python model using the climlab package. You learned the limitations of a pure radiation model and the necessity to include atmospheric dynamics such as convection. Through coupling a radiation model with a convective model, you have simulated a more realistic atmospheric temperature profile. Even though the model does not directly resolve the vertical motion during convection, it incorporates a parameterization to automatically account for mixing in statically unstable regions. By the end of the tutorial, you could comprehend how this model expands on the previous energy balance models from tutorials 1-4.

Resources

Data from this tutorial can be accessed for specific humidity here and reanalysis temperature here.