Skip to article frontmatterSkip to article content

Modeling non-scattering radiative transfer

University at Albany (SUNY)

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


1. The two-stream Schwarschild equations


Here we are going to introduce the governing equations of radiative transfer to put what we’ve been doing on a more solid theoretical footing.

Our derivations here will also serve as a coherent documentation for how some of the radiation solvers are implemented in climlab.

Optical thickness

The optical thickness of a layer of absorbers is Δτν\Delta \tau_\nu. It is the emissivity and absorptivity of a layer of atmosphere.

Passing to the limit of very thin layers, we define τν\tau_\nu through

dτνdp=1gκν\frac{d \tau_\nu}{dp} = -\frac{1}{g} \kappa_\nu

where κν\kappa_\nu is an absorption cross-section per unit mass at frequency ν\nu. It has units m2^2 kg1^{-1}. κ\kappa is a measure of the area taken out of the incident beam by absorbers in a unit mass of atmosphere.

In general τν\tau_\nu depends on the frequency of the radiation, as indicated here by the subscript ν\nu.

Using optical depth as vertical coordinate

Since pressure decreases with altitude, τν\tau_\nu increases with altitude.

The equations of radiative transfer can be simplified by using τν\tau_\nu as vertical coordinate instead of pressure.

Absorption cross-section

The specific absorption cross section κ\kappa depends on the number of molecules of each greenhouse gas encountered by the beam and the absorption properties characteristic to each kind of greenhouse gas molecule. Letting qkq_k be the mass-specific concentration of greenhouse gas kk, we may write in general

κ(ν,p,T)=k=1nκk(ν,p,T)qk(p)\kappa(\nu, p, T) = \sum_{k=1}^n \kappa_k\big(\nu, p, T \big) q_k(p)

For a well-mixed greenhouse gas, qkq_k is a constant; for a non-well-mixed gas like water vapor we need to account for the vertical distribution of the gas through qk(p)q_k(p).

The dependence of κk\kappa_k on temperature and pressure arises from certain aspects of the physics of molecular absorption.

Two-stream Schwarzschild equations

For climate modeling we almost always seperate the total flux into two beams: upward and downward.

This involves taking integrals of the full angular dependence of the flux. We’ll skip the details here.

Let UνU_\nu be the upward beam, and DνD_\nu be the downward beam. The governing equations for these beams are the Schwarzschild equations:

dUνdτν=Uν+E(ν,T(τν))dDνdτν=DνE(ν,T(τν))\begin{align} \frac{d U_\nu}{d \tau_\nu} &= -U_\nu + E\big( \nu, T(\tau_\nu) \big) \\ \frac{d D_\nu}{d \tau_\nu} &= D_\nu - E\big( \nu, T(\tau_\nu) \big) \end{align}

where EE is the blackbody emission (both up and down), which in general depends on both frequency and temperature. We have written temperature as a function of the vertical coordinate (optical depth).

The emissions are governed by the Planck function:

E=π B(ν,T)B(ν,T)=2hν3c21exp(hνkT)1\begin{align} E &= \pi~ B\big( \nu, T \big) \\ B\big( \nu, T \big) &= \frac{2 h \nu^3}{c^2} \frac{1}{\exp \left( \frac{h \nu}{k T} \right) -1} \end{align}

with these fundamental physical constants:

  • h=6.626×1034 J sh = 6.626 \times 10^{-34} ~\text{J s} is Planck’s constant
  • c=3.00×108 m s1c = 3.00 \times 10^8 ~\text{m s}^{-1} is the speed of light
  • k=1.38×1023 J K1k = 1.38 \times 10^{-23} ~\text{J K}^{-1} is the Boltzmann Thermodynamic Constant

The two-stream equations basically say the beam is attenuated by absorption (first term) and augmented by emission (second term) in each thin layer of gas.

These equations are valid for beam that are not affected by scattering. We may come back to that later.


2. The Grey Gas Model


The absorption properties κi\kappa_i of most atmospheric gases varies enormously with the frequency of the radiation. Hence the optical thickness τν\tau_\nu also has an intricate dependence on wavenumber or frequency.

Our job as climate modelers is to calculate and understand the net atmospheric absorption and transmission of radiation. To do this thoroughly and accurately, the fluxes must be solved for individually on a very dense grid of wavenumbers, and then the results integrated over all wavenumbers.

Actual radiative transfer codes (as used in GCMs) apply a lot of tricks and shortcuts to simplify this brute-force approach, but lead to sets of equations that are difficult to understand.

However, there is a lot we can understand about the basics of radiative transfer and the greenhouse effect by ignoring the spectral dependence of the flux.

Specifically, we make the approximation

κ(ν,p,T)=κ(p)\kappa(\nu, p, T) = \kappa(p)

so that the optical depth τ\tau is now independent of frequency.

This is known as the grey gas approximation.

Grey gas versions of the Schwarzschild equations

Integrating the Planck function

If we assume τν\tau_\nu is independent of frequency, we can then integrate the two-stream equations over all frequencies.

The integral of the Planck function gives our familiar Stefan-Boltzmann blackbody radiation law

E(T)=0πB(ν,T)dν=σT4E(T) = \int_0^{\infty} \pi B\big( \nu, T \big) d\nu = \sigma T^4

with

σ=2π5k415c2h3=5.67×108 W m2 K4\sigma = \frac{2 \pi^5 k^4}{15 c^2 h^3} = 5.67 \times 10^{-8} ~\text{W m}^{-2}~\text{K}^{-4}
# climlab has these constants available, and actually calculates sigma from the above formula
import numpy as np
import climlab 

sigma = ((2*np.pi**5 * climlab.constants.kBoltzmann**4) / 
         (15 * climlab.constants.c_light**2 * climlab.constants.hPlanck**3) )
print( sigma)
sigma == climlab.constants.sigma
5.6703726225913323e-08
True

The grey gas equations

This gives the governing equations for the grey gas model:

dUdτ=U+σT(τ)4dDdτ=DσT(τ)4\begin{align} \frac{d U}{d \tau} &= -U + \sigma T(\tau)^4 \\ \frac{d D}{d \tau} &= D - \sigma T(\tau)^4 \end{align}

These equations now say that the beam is diminished by absorption in a thin layer (first term) and augmented by blackbody emission, where both processes are assumed to be independent of frequency.

How do the fluxes change across a finite layer of absorbers?

The above equations are linear, first order ODEs, and they are uncoupled from each other (because we neglected scattering).

Consider a layer of atmosphere from optical level τ0\tau_{0} to τ1\tau_{1}. The optical thickness of the layer is Δτ=τ1τ0\Delta \tau = \tau_{1} - \tau_{0}.

The incident upwelling beam from below is denoted U0U_{0}; the upwelling beam that leaves the top of our layer is denoted U1U_{1}. We want to calculate U1U_{1}, which we can get by integrating the equation for dU/dτdU/d\tau over our layer.

The result is

U1=U0exp(Δτ)+τ0τ1Eexp((τ1τ))dτU_{1} = U_{0} \exp(-\Delta \tau) + \int_{\tau_{0}}^{\tau_{1}} E \exp \big( -(\tau_{1} - \tau) \big) d\tau

where τ\tau is now a dummy variable of integration, and we’ve written the blackbody emissions as E=σT4E = \sigma T^4.

The change in the downwelling beam is similar:

D0=D1exp(Δτ)+τ0τ1Eexp((ττ0))dτD_{0} = D_{1} \exp(-\Delta \tau ) + \int_{\tau_0}^{\tau_1} E \exp\big( -(\tau - \tau_{0}) \big) d\tau

3. Discretizing the Grey-gas equations on a finite grid


One particularly important reason to think about these changes over finite layers is that we will typically solve the radiative transfer equations in a numerical model where the temperature is represented on a discrete grid (often using pressure coordinates).

If we assume that temperature is uniform everywhere in our layer then the blackbody emission EE is also a constant across the layer, which we will denote E0E_{0}.

It can come out the integral and the expressions simplify to

U1=U0 exp(Δτ)+E0 (1exp(Δτ))D0=D1 exp(Δτ)+E0 (1exp(Δτ))\begin{align} U_1 &= U_0 ~ \exp(-\Delta \tau ) + E_{0} ~\Big( 1 - \exp(-\Delta \tau) \Big) \\ D_0 &= D_1 ~ \exp(-\Delta \tau ) + E_{0} ~\Big( 1 - \exp(-\Delta \tau) \Big) \end{align}

Transmissitivity

The first term is the transmission of radition from the bottom to the top of the layer (or vice-versa). We can define the transmissivity of the layer (denoted t0t_{0}) as the fraction of the incident beam that is passed on to the next layer:

t0=U0exp(Δτ)U0=exp(Δτ)\begin{align} t_{0} &= \frac{U_{0} \exp(-\Delta \tau)}{U_{0}} \\ &= \exp(-\Delta \tau)\\ \end{align}

Emissivity

The second term is the net change in the beam due to emissions in the layer.

We can the define the emissivity ϵ0\epsilon_0 of the layer analogously

ϵ0=1exp(Δτ)=1t0\begin{align} \epsilon_{0} &= 1 - \exp(-\Delta \tau) \\ &= 1 - t_{0} \end{align}

The discrete two-stream equations on a finite pressure grid

Putting this all together gives us two simple equations that govern changes in the upwelling and downwelling beams across a discrete layer of optical depth Δτ\Delta \tau:

U1=t0 U0+ϵ0 E0D0=t0 D1+ϵ0 E0\begin{align} U_{1} &= t_0 ~ U_{0} + \epsilon_{0} ~ E_{0} \\ D_{0} &= t_0 ~ D_{1} + \epsilon_{0} ~ E_{0} \\ \end{align}

Our model will typically be discretized in pressure coordinates. Suppose the thickness of the layer is Δp\Delta p, then the optical depth is

Δτ=κgΔp\Delta \tau = -\frac{\kappa}{g} \Delta p

(the minus sign accounts for the opposite sign conventions of the two coordinates).

Thus the emissivity of the layer is

ϵ0=1exp(κgΔp)\epsilon_{0} = 1 - \exp\big( \frac{\kappa}{g} \Delta p \big)

In the grey gas approximation we ignore all spectral dependence of the flux, so κ\kappa is a constant (we could let it vary in the vertical to represent non-well-mixed greenhouse gases like water vapor), and the blackbody emissions are simply

E0=σT04E_0 = \sigma T_0^4

4. Matrix equations for the grey gas model


climlab implements the discretized two-stream equations as written above.

We will now write out the equations on a discretized pressure grid with NN layers.

Let the upwelling flux be a column vector

U=[U0,U1,...,UN1,UN]T{\bf{U}} = [U_0, U_1, ..., U_{N-1}, U_N]^T

If there are NN levels then U\bf{U} has N+1N+1 elements (i.e. the fluxes are defined at the boundaries between levels). We will number the layers starting from 0 following numpy index conventions.

  • U0U_0 is the upwelling flux from surface to layer 0.
  • U1U_1 is the upwelling flux layer 0 to layer 1, etc.
  • UNU_N is the upwelling flux from layer N-1 (the top level) to space.

Same for the downwelling flux

D=[D0,D1,...,DN]T{\bf{D}} = [D_0, D_1, ..., D_N]^T

So DND_N is the flux down from space and D0D_0 is the back-radiation to the surface.

The temperature and blackbody emissions are defined for each NN pressure level, and are related by

Ei=σTi4E_i = \sigma T_i^4

Emissivity and Transmissivity vectors

Let the vector of absorptivity / emissivity be

ϵ=[ϵ0,ϵ1,...,ϵN1]{\bf{\epsilon}} = [\epsilon_0, \epsilon_1, ..., \epsilon_{N-1}]

where each element is determined by the thickness of the layer and the absorption cross-section for this particular spectral band:

ϵi=1exp(κgΔpi)\epsilon_{i} = 1 - \exp\big( \frac{\kappa}{g} \Delta p_i \big)

And the transmissivity of individual layers is

ti=1ϵit_{i} = 1 - \epsilon_{i}

It is convenient to define a vector of transmissivities t{\bf{t}} with N+1N+1 elements:

t=[1,t0,t1,...,tN1]{\bf{t}} = [1, t_0, t_1, ..., t_{N-1}]

The downwelling beam

For the downwelling beam, we define a column vector of emissions with N+1N+1 elements:

Edown=[E0,E1,...,EN1,EN]T{\bf{E_{down}}} = [E_0, E_1, ..., E_{N-1}, E_N]^{T}

where we define the last element ENE_N as the incident flux at the TOA.

For a longwave model, we would usually set EN=0E_N = 0. For a shortwave model, this would be the incident solar radiation.

We want a matrix Tdown{\bf{T_{down}}} that, when multiplied by ENE_N, gives the downwelling beam at each of the N+1N+1 layer interfaces:

D=Tdown Edown{\bf{D}} = {\bf{T_{down}}} ~ {\bf{E_{down}}}

The upwelling beam

Define a vector of emissions for the upwelling beam thus:

Eup=[upsfc,E0,E1,...,EN1]{\bf{E_{up}}} = [up_{sfc}, E_0, E_1, ..., E_{N-1}]

We need to add the reflected part of the downwelling beam at the surface to any emissions from the surface:

upsfc=Esfc+αsfc D[0]up_{sfc} = E_{sfc} + \alpha_{sfc} ~ D[0]

Now we want a matrix Tup{\bf{T_{up}}} such that the upwelling beam is

U=Tup Eup{\bf{U}} = {\bf{T_{up}}} ~ {\bf{E_{up}}}

The transmissivity matrices

Tup=[100...00t010...00t1t0t11...00t2t1t0t2t1t2...00..................0N1ti1N1ti2N1ti...tN11]{\bf{T_{up}}} = \left[ \begin{array}{cccccc} 1 & 0 & 0 & ... & 0 & 0 \\ t_0 & 1 & 0 & ... & 0 & 0 \\ t_1 t_0 & t_1 & 1 & ... & 0 & 0 \\ t_2 t_1 t_0 & t_2 t_1 & t_2 & ... & 0 & 0 \\ ... & ... & ... & ... & ... & ... \\ \prod_0^{N-1} t_i & \prod_1^{N-1} t_i & \prod_2^{N-1} t_i & ... & t_{N-1} & 1 \end{array} \right]

and

Tdown=TupT{\bf{T_{down}}} = {\bf{T_{up}}}^T

These formulas have been implemented in climlab.radiation.transmissivity.Transmissivity() using vectorized numpy array operations.

# example with N=2 layers and constant absorptivity
#  we construct an array of absorptivities
eps = np.array([0.58, 0.58])
#  and pass these as argument to the Transmissivity class
trans = climlab.radiation.transmissivity.Transmissivity(eps)
print( 'Matrix Tup is ')
print( trans.Tup)
print( 'Matrix Tdown is ')
print( trans.Tdown)
Matrix Tup is 
[[1.     0.42   0.1764]
 [0.     1.     0.42  ]
 [0.     0.     1.    ]]
Matrix Tdown is 
[[1.     0.     0.    ]
 [0.42   1.     0.    ]
 [0.1764 0.42   1.    ]]

5. Band-averaged radiation models in climlab


The Grey Gas model is a useful first step in understanding how radiation shapes the global energy balance. But we quickly run up against its limitations when trying to understand what really determines climate sensitivity.

What’s the next step in the model hierarchy?

Suppose we break up the spectrum into a discrete number MM of spectral bands. The idea is that we find parts of the spectrum in which the absorption characteristics of the important gases are relatively uniform. We then write the band-averaged absorption cross section for gas kk and band jj as

κkj(p,T)=νjκj(ν,p,T)dννjdν\kappa_{kj} \big(p, T \big) = \frac{\int_{\nu_j} \kappa_j \big(\nu, p, T \big) d \nu }{ \int_{\nu_j} d \nu }

where we integrate over whatever part of the spectrum we have chosen to define band jj.

In our band models we will typically ignore any dependence of κ\kappa on temperature. The total absorption cross section for our band is thus (summing over all absorbing gases):

κj(p)=k=1nκkj(p)qk(p)\kappa_j(p) = \sum_{k=1}^n \kappa_{kj}(p) q_k(p)

Notice that once we make this defintion, all of the formulas we wrote down above for the grey gas model can be written nearly identically for the fluxes in each band.

The optical depth in band jj is

Δτj=κjgΔp\Delta \tau_j = -\frac{\kappa_j}{g} \Delta p

from which we can define emissivity and transmissivity for band jj just as above.

The only difference from the Grey Gas formulas is that the blackbody emission in band jj (denoted EjE_j) is now only a fraction of σT4\sigma T^4.

We will denote this fraction as bjb_j.

The fraction bjb_j is temperature-dependent, and can be solved by integrating the Planck function:

Ej(T)=νjπB(ν,T)dνE_j(T) = \int_{\nu_j} \pi B\big( \nu, T \big) d\nu

To simplify our band models, we might choose to fix bjb_j in advance and just assume

Ej(T)=bj σ T4E_j(T) = b_j ~ \sigma ~ T^4

which is sensible if the temperatures don’t vary too much.

Regardless of how we calculate bjb_j, they must add up to one over all the bands in our model:

0M1bj=1\sum_0^{M-1} b_j = 1

Once we’ve figured out this division of the total flux into multiple bands, and we know the absorption cross-sections of each band, we can calculate the upwelling and downwelling fluxes independently for each band, using the same formulas (same code!) as we use in the grey gas model.

To get the total flux, we just need to sum the beams over all bands:

U=0M1UjD=0M1Dj\begin{align} U &= \sum_0^{M-1} U_j \\ D &= \sum_0^{M-1} D_j \end{align}

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.


References
  1. Pierrehumbert, R. T. (2010). Principles of Planetary Climate (p. 680). Cambridge University Press.