Open In Colab   Open in Kaggle

Bonus Tutorial 5: Vertical Structure of the Atmosphere#

Week 1, Day 5, Introduction to Climate Modeling

Content creators: Jenna Pearson and Brian E. J. Rose

Content reviewers: Mujeeb Abdulfatai, Nkongho Ayuketang Arreyndip, Jeffrey N. A. Aryee, Younkap Nina Duplex, Will Gregory, Paul Heubel, Zahra Khodakaramimaghsoud, Peter Ohue, Abel Shibu, Derick Temfack, Yunlong Xu, Peizhen Yang, Ohad Zivan, Chi Zhang

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

Production editors: Wesley Banfield, Paul Heubel, Jenna Pearson, Konstantine Tsafatinos, Chi Zhang, Ohad Zivan

Our 2024 Sponsors: CMIP, NFDI4Earth

Tutorial Objectives#

Estimated timing of tutorial: 30 minutes

Note that the previous, the present, and the subsequent tutorials contain solely Bonus content, that serves experienced participants and for later investigation of 0- and 1-dimensional models after the course. Please check out these tutorials after finishing Tutorial 7-8 successfully.

In this tutorial, you will run a one-dimensional radiative equilibrium model that predicts the global mean atmospheric temperature as a function of height. Much of the code shown here was taken from The Climate Laboratory by Brian E. J. Rose. You are encouraged to visit this website for more tutorials and background on these models.

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

  • Implement a 1-D model that predicts atmospheric temperature as a function of height using the python package climlab.

  • Understand how this model builds off of the energy balance models developed in the previous tutorials.

Setup#

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

# !pip install -q climlab cftime metpy
# !pip install -q git+https://github.com/climlab/climlab-rrtmg
# !pip install -q git+https://github.com/climlab/climlab-cam3-radiation
# !pip install -q git+https://github.com/climlab/climlab-emanuel-convection
# 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 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
Downloading data from 'http://www.atmos.albany.edu/facstaff/brose/resources/climlab_data/absorptivity/abs_ems_factors_fastvx.c030508.nc' to file '/home/runner/.cache/pooch/61707eb16bd77eb262c3c7db3706b35d-abs_ems_factors_fastvx.c030508.nc'.
/home/runner/micromamba/envs/climatematch/lib/python3.11/site-packages/climlab/radiation/cam3.py:45: UserWarning: Cannot import and initialize compiled Fortran extension, CAM3 module will not be functional.
  warnings.warn('Cannot import and initialize compiled Fortran extension, CAM3 module will not be functional.')

Install and import feedback gadget#

Hide code cell source
# @title Install and import feedback gadget

!pip3 install vibecheck datatops --quiet

from vibecheck import DatatopsContentReviewContainer
def content_review(notebook_section: str):
    return DatatopsContentReviewContainer(
        "",  # No text prompt
        notebook_section,
        {
            "url": "https://pmyvdlilci.execute-api.us-east-1.amazonaws.com/klab",
            "name": "comptools_4clim",
            "user_key": "l5jpxuee",
        },
    ).render()


feedback_prefix = "W1D5_T5"

Figure settings#

Hide code cell source
# @title Figure settings
import ipywidgets as widgets  # interactive display

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

Helper functions#

Hide code cell source
# @title Helper functions


def pooch_load(filelocation=None, filename=None, processor=None):
    # this is different for each day
    shared_location = "/home/jovyan/shared/Data/tutorials/W1D5_IntroductiontoClimateModeling"
    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

Video 1: Vertical Structure of the Atmosphere#

Submit your feedback#

Hide code cell source
# @title Submit your feedback
content_review(f"{feedback_prefix}_Vertical_Structure_of_the_Atmosphere_Video")
If you want to download the slides: https://osf.io/download/dxthu/

Submit your feedback#

Hide code cell source
# @title Submit your feedback
content_review(f"{feedback_prefix}_Vertical_Structure_of_the_Atmosphere_Slides")

Section 1: Setting up the Radiative Equilibrium Model Using Climlab#

The energy balance model we used earlier today was zero-dimensional, yielding only the global mean surface temperature. We might ask, is it possible to construct a similar, one-dimensional, model for an atmospheric column to estimate the global mean temperature profile (i.e., including the height/\(z\) dimension). Additionally, can we explicitly include the effects of different gases in this model, rather than just parametrizing their collective effects through a single parameter \(\tau\)? The answer is yes, we can!

This model is too complex to construct from scratch, as we did in the previous tutorials. Instead, we will use a model already available within the python package climlab.

The model we will first use is a radiative equilibrium model. Radiative equilibrium models consider different layers of the atmosphere. Each of these layers absorbs and emits radiation depending on its constituent gases, allowing the model to calculate the radiation budget for each layer as radiative energy is transferred between atmospheric layers, the Earth’s surface, and space. Radiative equilibrium is reached when each layer gains energy at the same rate as it loses energy. In this tutorial, you will analyze the temperature profile of this new model once it has reached equilibrium.

To set up this model, we will need information about some of the mean properties of the atmosphere. We are going to download water vapor data from the Community Earth System Model (CESM), a global climate model that we will go into detail on in the next tutorial, to use a variable called specific humidity. Specific humidity is the mass of water vapor per mass of a unit block of air. This is useful because water vapor is an important greenhouse gas.

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
ds
Downloading data from 'https://osf.io/c6q4j/download/' to file '/tmp/cpl_1850_f19-Q-gw-only.cam.h0.nc'.
SHA256 hash of downloaded file: 006bce4cbf3616123ec4bd67403f73e68fa78e8e3c9e2684cb6cae6124b7ff0f
Use this value as the 'known_hash' argument of 'pooch.retrieve' to ensure that the file hasn't changed if it is downloaded again in the future.
<xarray.Dataset> Size: 345MB
Dimensions:  (time: 240, lev: 26, lat: 96, lon: 144)
Coordinates:
  * lev      (lev) float64 208B 3.545 7.389 13.97 23.94 ... 929.6 970.6 992.6
  * time     (time) object 2kB 0001-02-01 00:00:00 ... 0021-01-01 00:00:00
  * lat      (lat) float64 768B -90.0 -88.11 -86.21 -84.32 ... 86.21 88.11 90.0
  * lon      (lon) float64 1kB 0.0 2.5 5.0 7.5 10.0 ... 350.0 352.5 355.0 357.5
Data variables:
    Q        (time, lev, lat, lon) float32 345MB ...
    gw       (lat) float64 768B ...
Attributes: (12/13)
    Conventions:               CF-1.0
    source:                    CAM
    case:                      cpl_1850_f19
    title:                     UNSET
    logname:                   br546577
    host:                      snow-30.rit.alba
    ...                        ...
    revision_Id:               $Id$
    initial_file:              b40.1850.track1.2deg.003.cam.i.0501-01-01-0000...
    topography_file:           /data/rose_scr/cesm_inputdata/atm/cam/topo/USG...
    history:                   Tue Feb 26 17:17:15 2019: ncrcat atm/hist/cpl_...
    NCO:                       4.6.8
    nco_openmp_thread_number:  1
# the specific humidity is stored in a variable called Q
ds.Q
<xarray.DataArray 'Q' (time: 240, lev: 26, lat: 96, lon: 144)> Size: 345MB
[86261760 values with dtype=float32]
Coordinates:
  * lev      (lev) float64 208B 3.545 7.389 13.97 23.94 ... 929.6 970.6 992.6
  * time     (time) object 2kB 0001-02-01 00:00:00 ... 0021-01-01 00:00:00
  * lat      (lat) float64 768B -90.0 -88.11 -86.21 -84.32 ... 86.21 88.11 90.0
  * lon      (lon) float64 1kB 0.0 2.5 5.0 7.5 10.0 ... 350.0 352.5 355.0 357.5
Attributes:
    mdims:         1
    units:         kg/kg
    long_name:     Specific humidity
    cell_methods:  time: mean
ds.time
<xarray.DataArray 'time' (time: 240)> Size: 2kB
array([cftime.DatetimeNoLeap(1, 2, 1, 0, 0, 0, 0, has_year_zero=True),
       cftime.DatetimeNoLeap(1, 3, 1, 0, 0, 0, 0, has_year_zero=True),
       cftime.DatetimeNoLeap(1, 4, 1, 0, 0, 0, 0, has_year_zero=True), ...,
       cftime.DatetimeNoLeap(20, 11, 1, 0, 0, 0, 0, has_year_zero=True),
       cftime.DatetimeNoLeap(20, 12, 1, 0, 0, 0, 0, has_year_zero=True),
       cftime.DatetimeNoLeap(21, 1, 1, 0, 0, 0, 0, has_year_zero=True)],
      dtype=object)
Coordinates:
  * time     (time) object 2kB 0001-02-01 00:00:00 ... 0021-01-01 00:00:00
Attributes:
    long_name:  time
    bounds:     time_bnds

However, we want an annual average profile:

# take global, annual average using a weighting (ds.gw) that is calculated based on the model grid - and is similar, but not identical, to a cosine(latitude) weighting

weight_factor = ds.gw / ds.gw.mean(dim="lat")
Qglobal = (ds.Q * weight_factor).mean(dim=("lat", "lon", "time"))
# print specific humidity profile
Qglobal
<xarray.DataArray (lev: 26)> Size: 208B
array([2.16104904e-06, 2.15879387e-06, 2.15121262e-06, 2.13630949e-06,
       2.12163684e-06, 2.11168002e-06, 2.09396914e-06, 2.10589390e-06,
       2.42166155e-06, 3.12595653e-06, 5.01369691e-06, 9.60746488e-06,
       2.08907654e-05, 4.78823747e-05, 1.05492451e-04, 2.11889055e-04,
       3.94176751e-04, 7.10734458e-04, 1.34192099e-03, 2.05153261e-03,
       3.16844784e-03, 4.96883408e-03, 6.62218037e-03, 8.38350326e-03,
       9.38620899e-03, 9.65030544e-03])
Coordinates:
  * lev      (lev) float64 208B 3.545 7.389 13.97 23.94 ... 929.6 970.6 992.6

Now that we have a global mean water vapor profile, we can define a model that has the same vertical levels as this water vapor data.

# 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)
mystate
AttrDict({'Ts': Field([288.]), 'Tatm': Field([200.  , 203.12, 206.24, 209.36, 212.48, 215.6 , 218.72, 221.84,
       224.96, 228.08, 231.2 , 234.32, 237.44, 240.56, 243.68, 246.8 ,
       249.92, 253.04, 256.16, 259.28, 262.4 , 265.52, 268.64, 271.76,
       274.88, 278.  ])})

To model the absorption and emission of different gases within each atmospheric layer, we use the Rapid Radiative Transfer Model, which is contained within the RRTMG module. We must first initialize our model using the water vapor.

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,                       # set the SURFACE shortwave albedo
    # set the timestep to one day (measured in seconds)
    timestep=climlab.constants.seconds_per_day,
)
radmodel
Downloading data from 'http://www.atmos.albany.edu/facstaff/brose/resources/climlab_data/ozone/apeozone_cam3_5_54.nc' to file '/home/runner/.cache/pooch/0e7ac27a9b97c3a843dcfdd7cd89a9b7-apeozone_cam3_5_54.nc'.
---------------------------------------------------------------------------
TimeoutError                              Traceback (most recent call last)
File ~/micromamba/envs/climatematch/lib/python3.11/site-packages/urllib3/connection.py:239, in HTTPConnection._new_conn(self)
    238 try:
--> 239     sock = connection.create_connection(
    240         (self._dns_host, self.port),
    241         self.timeout,
    242         source_address=self.source_address,
    243         socket_options=self.socket_options,
    244     )
    245 except socket.gaierror as e:

File ~/micromamba/envs/climatematch/lib/python3.11/site-packages/urllib3/util/connection.py:85, in create_connection(address, timeout, source_address, socket_options)
     84 try:
---> 85     raise err
     86 finally:
     87     # Break explicitly a reference cycle

File ~/micromamba/envs/climatematch/lib/python3.11/site-packages/urllib3/util/connection.py:73, in create_connection(address, timeout, source_address, socket_options)
     72     sock.bind(source_address)
---> 73 sock.connect(sa)
     74 # Break explicitly a reference cycle

TimeoutError: timed out

The above exception was the direct cause of the following exception:

ConnectTimeoutError                       Traceback (most recent call last)
File ~/micromamba/envs/climatematch/lib/python3.11/site-packages/urllib3/connectionpool.py:793, in HTTPConnectionPool.urlopen(self, method, url, body, headers, retries, redirect, assert_same_host, timeout, pool_timeout, release_conn, chunked, body_pos, preload_content, decode_content, **response_kw)
    792 # Make the request on the HTTPConnection object
--> 793 response = self._make_request(
    794     conn,
    795     method,
    796     url,
    797     timeout=timeout_obj,
    798     body=body,
    799     headers=headers,
    800     chunked=chunked,
    801     retries=retries,
    802     response_conn=response_conn,
    803     preload_content=preload_content,
    804     decode_content=decode_content,
    805     **response_kw,
    806 )
    808 # Everything went great!

File ~/micromamba/envs/climatematch/lib/python3.11/site-packages/urllib3/connectionpool.py:499, in HTTPConnectionPool._make_request(self, conn, method, url, body, headers, retries, timeout, chunked, response_conn, preload_content, decode_content, enforce_content_length)
    498 try:
--> 499     conn.request(
    500         method,
    501         url,
    502         body=body,
    503         headers=headers,
    504         chunked=chunked,
    505         preload_content=preload_content,
    506         decode_content=decode_content,
    507         enforce_content_length=enforce_content_length,
    508     )
    510 # We are swallowing BrokenPipeError (errno.EPIPE) since the server is
    511 # legitimately able to close the connection after sending a valid response.
    512 # With this behaviour, the received response is still readable.

File ~/micromamba/envs/climatematch/lib/python3.11/site-packages/urllib3/connection.py:567, in HTTPConnection.request(self, method, url, body, headers, chunked, preload_content, decode_content, enforce_content_length)
    566     self.putheader(header, value)
--> 567 self.endheaders()
    569 # If we're given a body we start sending that in chunks.

File ~/micromamba/envs/climatematch/lib/python3.11/http/client.py:1346, in HTTPConnection.endheaders(self, message_body, encode_chunked)
   1345     raise CannotSendHeader()
-> 1346 self._send_output(message_body, encode_chunked=encode_chunked)

File ~/micromamba/envs/climatematch/lib/python3.11/http/client.py:1106, in HTTPConnection._send_output(self, message_body, encode_chunked)
   1105 del self._buffer[:]
-> 1106 self.send(msg)
   1108 if message_body is not None:
   1109 
   1110     # create a consistent interface to message_body

File ~/micromamba/envs/climatematch/lib/python3.11/http/client.py:1044, in HTTPConnection.send(self, data)
   1043 if self.auto_open:
-> 1044     self.connect()
   1045 else:

File ~/micromamba/envs/climatematch/lib/python3.11/site-packages/urllib3/connection.py:398, in HTTPConnection.connect(self)
    397 def connect(self) -> None:
--> 398     self.sock = self._new_conn()
    399     if self._tunnel_host:
    400         # If we're tunneling it means we're connected to our proxy.

File ~/micromamba/envs/climatematch/lib/python3.11/site-packages/urllib3/connection.py:248, in HTTPConnection._new_conn(self)
    247 except SocketTimeout as e:
--> 248     raise ConnectTimeoutError(
    249         self,
    250         f"Connection to {self.host} timed out. (connect timeout={self.timeout})",
    251     ) from e
    253 except OSError as e:

ConnectTimeoutError: (<HTTPConnection(host='www.atmos.albany.edu', port=80) at 0x701f52842710>, 'Connection to www.atmos.albany.edu timed out. (connect timeout=30)')

The above exception was the direct cause of the following exception:

MaxRetryError                             Traceback (most recent call last)
File ~/micromamba/envs/climatematch/lib/python3.11/site-packages/requests/adapters.py:696, in HTTPAdapter.send(self, request, stream, timeout, verify, cert, proxies)
    695 try:
--> 696     resp = conn.urlopen(
    697         method=request.method,
    698         url=url,
    699         body=request.body,  # type: ignore[arg-type]  # urllib3 stubs don't accept Iterable[bytes | str]
    700         headers=request.headers,  # type: ignore[arg-type]  # urllib3#3072
    701         redirect=False,
    702         assert_same_host=False,
    703         preload_content=False,
    704         decode_content=False,
    705         retries=self.max_retries,
    706         timeout=resolved_timeout,
    707         chunked=chunked,
    708     )
    710 except (ProtocolError, OSError) as err:

File ~/micromamba/envs/climatematch/lib/python3.11/site-packages/urllib3/connectionpool.py:847, in HTTPConnectionPool.urlopen(self, method, url, body, headers, retries, redirect, assert_same_host, timeout, pool_timeout, release_conn, chunked, body_pos, preload_content, decode_content, **response_kw)
    845     new_e = ProtocolError("Connection aborted.", new_e)
--> 847 retries = retries.increment(
    848     method, url, error=new_e, _pool=self, _stacktrace=sys.exc_info()[2]
    849 )
    850 retries.sleep()

File ~/micromamba/envs/climatematch/lib/python3.11/site-packages/urllib3/util/retry.py:555, in Retry.increment(self, method, url, response, error, _pool, _stacktrace)
    554     reason = error or ResponseError(cause)
--> 555     raise MaxRetryError(_pool, url, reason) from reason  # type: ignore[arg-type]
    557 log.debug("Incremented Retry for (url='%s'): %r", url, new_retry)

MaxRetryError: HTTPConnectionPool(host='www.atmos.albany.edu', port=80): Max retries exceeded with url: /facstaff/brose/resources/climlab_data/ozone/apeozone_cam3_5_54.nc (Caused by ConnectTimeoutError(<HTTPConnection(host='www.atmos.albany.edu', port=80) at 0x701f52842710>, 'Connection to www.atmos.albany.edu timed out. (connect timeout=30)'))

During handling of the above exception, another exception occurred:

ConnectTimeout                            Traceback (most recent call last)
Cell In[15], line 1
----> 1 radmodel = climlab.radiation.RRTMG(
      2     name="Radiation (all gases)",      # give our model a name
      3     state=mystate,                     # give our model an initial condition
      4     specific_humidity=Qglobal.values,  # tell the model how much water vapor there is

File ~/micromamba/envs/climatematch/lib/python3.11/site-packages/climlab/radiation/rrtm/rrtmg.py:125, in RRTMG.__init__(self, icld, irng, idrv, permuteseed_sw, permuteseed_lw, dyofyr, inflgsw, inflglw, iceflgsw, iceflglw, liqflgsw, liqflglw, tauc_sw, tauc_lw, ssac_sw, asmc_sw, fsfc_sw, iaer, tauaer_sw, ssaaer_sw, asmaer_sw, ecaer_sw, tauaer_lw, isolvar, indsolvar, bndsolvar, solcycfrac, **kwargs)
     13 def __init__(self,
     14         # GENERAL, used in both SW and LW
     15         icld = 1,    # Cloud overlap method, 0: Clear only, 1: Random, 2,  Maximum/random] 3: Maximum
   (...)    123 
    124         **kwargs):
--> 125     super(RRTMG, self).__init__(**kwargs)
    127     # Remove specific inputs from kwargs dictionary.
    128     #  We want any changes implemented in the parent __init__ method to be preserved here
    129     remove_list = ['absorber_vmr','cldfrac','clwp','ciwp','r_liq','r_ice',
    130                    'emissivity','aldif','aldir','asdif','asdir','S0','coszen',
    131                    'irradiance_factor','insolation',]

File ~/micromamba/envs/climatematch/lib/python3.11/site-packages/climlab/radiation/radiation.py:216, in _Radiation_SW.__init__(self, albedo, aldif, aldir, asdif, asdir, S0, insolation, coszen, irradiance_factor, **kwargs)
    205 def __init__(self,
    206              albedo = None,
    207              aldif = 0.3,
   (...)    214              irradiance_factor = 1.,  # instantaneous irradiance = S0 * irradiance_factor
    215              **kwargs):
--> 216     super(_Radiation_SW, self).__init__(**kwargs)
    217     self.add_input('S0', S0)
    218     self.add_input('coszen', coszen)

File ~/micromamba/envs/climatematch/lib/python3.11/site-packages/climlab/radiation/radiation.py:273, in _Radiation_LW.__init__(self, emissivity, **kwargs)
    270 def __init__(self,
    271              emissivity = 1.,  # surface emissivity
    272              **kwargs):
--> 273     super(_Radiation_LW, self).__init__(**kwargs)
    274     self.add_input('emissivity', emissivity)
    275     # initialize diagnostics

File ~/micromamba/envs/climatematch/lib/python3.11/site-packages/climlab/radiation/radiation.py:193, in _Radiation.__init__(self, specific_humidity, absorber_vmr, cldfrac, clwp, ciwp, r_liq, r_ice, ozone_file, **kwargs)
    191 self.add_input('specific_humidity', specific_humidity)
    192 if absorber_vmr is None:
--> 193     absorber_vmr = default_absorbers(self.Tatm, ozone_file, self.verbose)
    194 self.add_input('absorber_vmr', absorber_vmr)
    195 self.add_input('cldfrac', cldfrac)

File ~/micromamba/envs/climatematch/lib/python3.11/site-packages/climlab/radiation/radiation.py:137, in default_absorbers(Tatm, ozone_file, verbose)
    135 if ozone_file is not None:
    136     ozonepath_http = _datapath_http + 'ozone/' + ozone_file
--> 137     ozonefilehandle = pooch.retrieve(url=ozonepath_http,
    138         known_hash="bc659bfa129fafa4ed9368bb19278ae15724a5a66599affd317c143ba511ff84")
    139     ozonedata = xr.open_dataset(ozonefilehandle)
    140     ##  zonal and time average

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

File ~/micromamba/envs/climatematch/lib/python3.11/site-packages/pooch/core.py:823, in stream_download(url, fname, known_hash, downloader, pooch, retry_if_failed)
    819 try:
    820     # Stream the file to a temporary so that we can safely check its
    821     # hash before overwriting the original.
    822     with temporary_file(path=str(fname.parent)) as tmp:
--> 823         downloader(url, tmp, pooch)
    824         hash_matches(tmp, known_hash, strict=True, source=str(fname.name))
    825         shutil.move(tmp, str(fname))

File ~/micromamba/envs/climatematch/lib/python3.11/site-packages/pooch/downloaders.py:230, in HTTPDownloader.__call__(self, url, output_file, pooch, check_only)
    228     # pylint: enable=consider-using-with
    229 try:
--> 230     response = requests.get(url, timeout=timeout, **kwargs)
    231     response.raise_for_status()
    232     content = response.iter_content(chunk_size=self.chunk_size)

File ~/micromamba/envs/climatematch/lib/python3.11/site-packages/requests/api.py:87, in get(url, params, **kwargs)
     74 def get(
     75     url: _t.UriType, params: _t.ParamsType = None, **kwargs: Unpack[_t.GetKwargs]
     76 ) -> Response:
     77     r"""Sends a GET request.
     78 
     79     :param url: URL for the new :class:`Request` object.
   (...)     84     :rtype: requests.Response
     85     """
---> 87     return request("get", url, params=params, **kwargs)

File ~/micromamba/envs/climatematch/lib/python3.11/site-packages/requests/api.py:71, in request(method, url, **kwargs)
     67 # By using the 'with' statement we are sure the session is closed, thus we
     68 # avoid leaving sockets open which can trigger a ResourceWarning in some
     69 # cases, and look like a memory leak in others.
     70 with sessions.Session() as session:
---> 71     return session.request(method=method, url=url, **kwargs)

File ~/micromamba/envs/climatematch/lib/python3.11/site-packages/requests/sessions.py:651, in Session.request(self, method, url, params, data, headers, cookies, files, auth, timeout, allow_redirects, proxies, hooks, stream, verify, cert, json)
    646 send_kwargs = {
    647     "timeout": timeout,
    648     "allow_redirects": allow_redirects,
    649 }
    650 send_kwargs.update(settings)
--> 651 resp = self.send(prep, **send_kwargs)
    653 return resp

File ~/micromamba/envs/climatematch/lib/python3.11/site-packages/requests/sessions.py:784, in Session.send(self, request, **kwargs)
    781 start = preferred_clock()
    783 # Send the request
--> 784 r = adapter.send(request, **kwargs)
    786 # Total elapsed time of the request (approximately)
    787 elapsed = preferred_clock() - start

File ~/micromamba/envs/climatematch/lib/python3.11/site-packages/requests/adapters.py:717, in HTTPAdapter.send(self, request, stream, timeout, verify, cert, proxies)
    714 if isinstance(e.reason, ConnectTimeoutError):
    715     # TODO: Remove this in 3.0.0: see #2811
    716     if not isinstance(e.reason, NewConnectionError):
--> 717         raise ConnectTimeout(e, request=request)
    719 if isinstance(e.reason, ResponseError):
    720     raise RetryError(e, request=request)

ConnectTimeout: HTTPConnectionPool(host='www.atmos.albany.edu', port=80): Max retries exceeded with url: /facstaff/brose/resources/climlab_data/ozone/apeozone_cam3_5_54.nc (Caused by ConnectTimeoutError(<HTTPConnection(host='www.atmos.albany.edu', port=80) at 0x701f52842710>, 'Connection to www.atmos.albany.edu timed out. (connect timeout=30)'))

Let’s explore this initial state. Here Ts is the initial global mean surface temperature, and Tatm is the initial global mean air temperature profile.

radmodel.state

One of the perks of using this model is its ability to incorporate the radiative effects of individual greenhouse gases in different parts of the radiation spectrum, rather than using a bulk reduction in transmission of outgoing longwave radiation (as in our previous models).

Let’s display absorber_vmr, which contains the volume mixing ratio’s of each gas used in the radiative transfer model (these are pre-defined; and do not include the water vapor we used as a model input above). The volume mixing ratio describes the fraction of molecules in the air that are a given gas. For example, \(21\%\) of air is oxygen so its volume mixing ratio is \(0.21\).

radmodel.absorber_vmr

To look at carbon dioxide ("CO2") in a more familiar unit, parts per million (by volume), we can convert and print the new value.

radmodel.absorber_vmr["CO2"] * 1e6

We can also look at all the available diagnostics of our model:

diag_ds = climlab.to_xarray(radmodel.diagnostics)
diag_ds

For example to look at \(OLR\):

radmodel.OLR

Note. the \(OLR\) is currently 0 as we have not run the model forward in time, so it has not calculated any radiation components.

Questions 1: Climate Connection#

  1. Why do you think all gases, except ozone and water vapor, are represented by single values in the model?

Click for solution

Submit your feedback#

Hide code cell source
# @title Submit your feedback
content_review(f"{feedback_prefix}_Questions_1")

Coding Exercises 1#

  1. On the same graph, plot the annual mean specific humidity profile and ozone profiles.

fig, ax = plt.subplots()
# multiply Qglobal by 1000 to put in units of grams water vapor per kg of air
_ = ...
# multiply by 1E6 to get units of ppmv = parts per million by volume
_ = ...

# pressure decreases logarithmically with height in the atmosphere
# invert the axis so the largest value of pressure is lowest
ax.invert_yaxis()
# set y axis to a log scale
_ = ...

ax.set_ylabel("Pressure (hPa)")
ax.set_xlabel("Specific humidity (g/kg)")

# turn on the grid lines
_ = ...

# turn on legend
_ = ...
../../../_images/a88a699acb474bc057f1d01399cbe4fdbd2f1f3314cf88b93c7b9dab0be910e0.png

Click for solution

Example output:

Solution hint

Submit your feedback#

Hide code cell source
# @title Submit your feedback
content_review(f"{feedback_prefix}_Coding_Exercises_1")

Section 2: Getting Data to Compare to the Model#

Before we run our model forward, we will download a reanalysis product from NCEP to get a sense of what the real global mean atmospheric temperature profile looks like. We will compare this profile to our model runs later.

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

# this is the long term monthly means (note only 12 time steps)
ncep_air.air
Downloading data from 'https://osf.io/w6cd5/download/' to file '/tmp/air.mon.1981-2010.ltm.nc'.
SHA256 hash of downloaded file: e653932e5153d1b06f04eff66b788014b5d9f4aa1f8fe0f954077e5c79f2c183
Use this value as the 'known_hash' argument of 'pooch.retrieve' to ensure that the file hasn't changed if it is downloaded again in the future.
/home/runner/micromamba/envs/climatematch/lib/python3.11/site-packages/xarray/coding/times.py:995: SerializationWarning: Unable to decode time axis into full numpy.datetime64 objects, continuing using cftime.datetime objects instead, reason: dates out of range
  dtype = _decode_cf_datetime_dtype(data, units, calendar, self.use_cftime)
/home/runner/micromamba/envs/climatematch/lib/python3.11/site-packages/xarray/coding/times.py:170: SerializationWarning: Ambiguous reference date string: 1-1-1 00:00:0.0. The first value is assumed to be the year hence will be padded with zeros to remove the ambiguity (the padded reference date string is: 0001-1-1 00:00:0.0). To remove this message, remove the ambiguity by padding your reference date strings with zeros.
  warnings.warn(warning_msg, SerializationWarning)
/home/runner/micromamba/envs/climatematch/lib/python3.11/site-packages/xarray/core/indexing.py:557: SerializationWarning: Unable to decode time axis into full numpy.datetime64 objects, continuing using cftime.datetime objects instead, reason: dates out of range
  array = array.get_duck_array()
<xarray.DataArray 'air' (time: 12, level: 17, lat: 73, lon: 144)> Size: 9MB
[2144448 values with dtype=float32]
Coordinates:
  * level    (level) float32 68B 1e+03 925.0 850.0 700.0 ... 50.0 30.0 20.0 10.0
  * lat      (lat) float32 292B 90.0 87.5 85.0 82.5 ... -82.5 -85.0 -87.5 -90.0
  * lon      (lon) float32 576B 0.0 2.5 5.0 7.5 10.0 ... 350.0 352.5 355.0 357.5
  * time     (time) object 96B 0001-01-01 00:00:00 ... 0001-12-01 00:00:00
Attributes:
    long_name:     Monthly Long Term Mean of Air temperature
    units:         degC
    precision:     2
    var_desc:      Air Temperature
    level_desc:    Pressure Levels
    statistic:     Long Term Mean
    parent_stat:   Mean
    valid_range:   [-200.  300.]
    actual_range:  [-89.722336  41.616005]
    dataset:       NCEP Reanalysis Derived Products
# 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"))
Tglobal
<xarray.DataArray (level: 17)> Size: 68B
array([ 15.179084  ,  11.207003  ,   7.8383274 ,   0.21994135,
        -6.4483433 , -14.888848  , -25.570469  , -39.36969   ,
       -46.797905  , -53.652245  , -60.56356   , -67.006065  ,
       -65.53293   , -61.48664   , -55.853584  , -51.593952  ,
       -43.21999   ], dtype=float32)
Coordinates:
  * level    (level) float32 68B 1e+03 925.0 850.0 700.0 ... 50.0 30.0 20.0 10.0

Below we will define two helper funcitons to visualize the profiles output from our model with a SkewT plot. This is common way to plot atmospheric temperature in climate science, and the metpy package has a built in function to make this easier.

# 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()
skew = make_skewT()
../../../_images/d97704ed01ce70bfc3b8e5b25399a2a62851bba2395394cd3bb24431e8d2319f.png

SkewT (also known as SkewT-logP) plots are generally used for much more complex reasons than we will use here. However, one of the benefits of this plot that we will utilize is the fact that pressure decreases approximately logarithmically with height. Thus, with a logP axis, we are showing information that is roughly linear in height, making the plots more intuitive.

Section 3: Running the Radiative Equilibrium Model Forward in Time#

We can run this model over many time steps, just like the simple greenhouse model, but now we can examine the behavior of the temperature profile rather than just the surface temperature.

There is no need to write out a function to step our model forward - climlab already has this feature. We will use this function to run our model to equilibrium (i.e. until \(OLR\) is balanced by \(ASR\)).

# 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()
#  check the energy budget to make sure we are really at equilibrium
radmodel.ASR - radmodel.OLR

Now let’s can compare this to observations.

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

Questions 3: Climate Connection#

  1. The profile from our model does not match observations well. Can you think of one component we might be missing?

  2. What effect do you think the individual gases play in determining this profile and why?

Click for solution

Submit your feedback#

Hide code cell source
# @title Submit your feedback
content_review(f"{feedback_prefix}_Questions_3")

Coding Exercises 3#

  1. Create a second model called ‘Radiation (no H20)’ that lacks water vapor. Then re-create the plot above, but add on this extra profile without water vapor.

# make an exact clone of our existing model
radmodel_noH2O = climlab.process_like(radmodel)
# change the name of our new model
radmodel_noH2O.name = ...

# set the water vapor profile to all zeros
radmodel_noH2O.specific_humidity *= 0.0

# run the model to equilibrium
radmodel_noH2O.step_forward()
while np.abs(radmodel_noH2O.ASR - radmodel_noH2O.OLR) > 0.01:
    radmodel_noH2O.step_forward()

# create skewT plot
skew = make_skewT()

# add profiles for both models to plot
for model in [...]:
    ...

Click for solution

Example output:

Solution hint

Submit your feedback#

Hide code cell source
# @title Submit your feedback
content_review(f"{feedback_prefix}_Coding_Exercises_3")

Summary#

In this tutorial, you’ve learned how to use the python package climlab to construct a one-dimensional radiative equilibrium model, and run it forward in time to predict the global mean atmospheric temperature profile. You’ve also visualized these results through SkewT plots.

Resources#

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

Useful links: