50. Integrating the Schwartzchild equation

To calculate the radiative flux passing through an atmospheric layer we need to know the optical depth as a function of wavelength, height and direction and the temperature as a function of height. That means we need to do three separate integrations with respect to \(\lambda\), \(\tau\), and \(\theta\), assuming that the radiance is independent of azimuth \(\phi\).

50.1. Integration over wavelength

Typically this integration is done over a set of intervals (perhaps 30-40) with different values of the absorption coefficient \(k_\lambda\) for each of the absorbing gasses (assumed independent). We also need to integrate the Planck function over wavelength. From Assignment 6 we know that:

\[ \int_0^\infty B_{\lambda} d\lambda = \frac{\sigma}{\pi} T^4 \]

Below we’ll also define \(\overline{\tau}\) as the average over a wavelength interval \(\Delta \lambda\)

\[ \overline{\tau} = \frac{\int_{\Delta \lambda} \tau_\lambda d\lambda}{\Delta \lambda} \]

50.2. Integration over optical depth

For slant paths start by rewriting (29.6) in The Schwartzchild Equation to account for \(\mu = \cos \theta\)

(50.1)\[ dL_\lambda= \left ( -L_\lambda + B_\lambda (T_{layer} ) \right ) \frac{d\tau_\lambda}{\mu} \]

For constant temperature we know how to solve this. Since \(\mu\) isn’t related to temperature we can just repeat the change of variables trick in Change of variables and get the zenith-angle dependent version.

(50.2)\[ L_\lambda = L_{\lambda 0} \exp( -\tau_{\lambda}/\mu ) + B_\lambda (T_{layer})(1- \exp( -\tau_{\lambda} /\mu)) \]

For height dependent temperature in r (29.13) we used an integrating factor of \(\exp(\tau)\). You should convince yourself that if we do the same thing for (50.1) but with an integrating factor of \(\exp(\tau/\mu)\) we get:

(50.3)\[ L_\lambda = B_\lambda(T_{skin}) \exp(-\tau/\mu) + \int_0^{\tau} \exp\left( - (\tau -\tau^\prime)/\,\mu \right ) B_\lambda(T) d \tau^\prime / \mu \]

and defining the slant trasmission:

\[ t_s = \exp\left( - (\tau -\tau^\prime)/\,\mu \right ) \]

(50.3) simplifies to

(50.4)\[ L_\lambda = B_\lambda(T_{skin}) t_s(\tau) + \int_0^{\tau} B_\lambda(T) d ts(\tau,\tau^\prime) / \mu \]

just like (29.21).

50.3. Integration over zenith angle: no atmosphere

We can’t integrate (50.3) over \(\mu\) analytically, but the integrals are easy to do numerically. As we’ll show in a Calculating the diffuse flux: python code, the following “diffusivity” approximation is quite good:

(50.5)\[ t_f = \int_0^1 \mu \exp \left ( \frac{(\tau - \tau^\prime)}{\mu} \right ) d\mu = \exp \left (-1.66 (\tau - \tau^\prime) \right ) \]

So that the transmissivity for flux is the same as the transmissivity for radiance except that the effective thickness of the atmosphere is increased by a factor of 1.666.

Where does this factor of 1.66 come from?

• In the Finding the flux given the radiance notes we turned blackbody isotropic radiance into a flux by taking the normal component and integrating over the hemisphere:

  (50.6)\[\begin{split} \begin{align} F_\lambda&= \int_0^{2 \pi} \int_0^{\pi/2} \cos \theta\, L_\lambda \sin \theta d\theta d\phi \\ &= 2 \pi \int_0^1 \mu \, L_\lambda d\mu \end{align} \end{split}\]

  assuming no dependence of \(L_\lambda\) on \(\phi\) and \(\theta\) we can take \(L_\lambda\) out of the integral and substituting \(\mu= \cos \theta\) write:

  (50.7)\[ F_\lambda= 2 \pi L_\lambda \int_0^1 \mu d\mu = 2 \pi \frac{\mu^2}{2} \Bigg \rvert_0^1 = 2 \pi L_\lambda \frac{ 1}{ 2} = \pi L_\lambda \]

50.3.1. now add an atmosphere

Once we add the atmopsheric transmissivity, we have a much more difficult integral:

Getting the flux using \(L_\lambda\) from (50.3) looks like this:

(50.8)\[\begin{split} F_\uparrow &= \int_0^{2 \pi} \int_0^{\pi/2} \cos \theta\, L_\lambda \sin \theta d\theta d\phi \\ &= 2\pi \int_0^1 \mu L_\lambda d \mu \notag\\ &= 2 \pi B_{\lambda 0}(T_s) \int_0^1 \mu \exp(-\tau/\mu) d\mu \nonumber\\ &+ 2 \pi \int_0^1 \int_0^{\tau} \exp\left( -\frac{(\tau - \tau^\prime)}{\mu} \right ) B_\lambda(T)\,d\tau^\prime d\mu \end{split}\]

Comparing these two terms with (50.3), note that the \(1/\mu\) has been cancelled from the second term in (50.8), but the structure is otherwise the same for both flux and radiance.

50.3.2. Flux transmissivity

• To make progress, first swap the limits of integration:

(50.9)\[\begin{split} \begin{gather} F_{\lambda \uparrow} = \pi B_{\lambda 0}(T_s)\, 2 \int_0^1 \mu \exp(-\tau/\mu) d\mu + \nonumber\\ \int_0^{\tau} \pi B_\lambda(T)\, 2 \int_0^1 \exp\left( -\frac{(\tau - \tau^\prime)}{\mu} \right ) \, d\mu\, d\tau^\prime \end{gather} \end{split}\]

(Remember that \(T\) does depend on \(z\), and therefore \(B_\lambda(T)\) has to stay inside the \(\tau^\prime\) integral.)

• Look what happens to this equation if we define \(t_f\), the flux transmissivity as:

(50.10)\[ t_f= 2 \int_0^1 \mu \exp(-(\tau - \tau^\prime)/\mu) d\mu \]

and differentiating wrt \(\tau^\prime\):

(50.11)\[ dt_f= 2 \int_0^1 \exp(-(\tau - \tau^\prime)/\mu) d\mu d\tau^\prime \]

Plug these into (50.9) and get:

(50.12)\[ F_\uparrow = \pi B_{\lambda 0}(T_s) \, t_f(0,\tau) + \int_0^{\tau} \pi \, B_\lambda(T)\,d t_f(\tau^\prime,\tau) \]

50.3.3. Exponential integrals

50.3.3.1. calculating \(t_f\)

• But how do we get values for \(t_f\) and \(dt_f\) if we can’t do the integration?

  – In Calculating the diffuse flux: python code we use python to evaluate:

  (50.13)\[ t_f= 2 \int_0^1 \mu \exp(-(\tau - \tau^\prime)/\mu) d\mu = 2 E_3(\tau) \]

• As we show there, it turns out to a very good approximation:

  (50.14)\[ t_f(\tau) = 2 E_3(\tau) \approx \exp \left (- 1.666 \tau \right ) \]

50.3.3.2. calculating \(dt_f\)

It turns out that \(\frac{dE_3}{d\tau^\prime} = -E_2(\tau^\prime)\) (see wikipedia).

In assignment 7 you’ll show that

(50.15)\[ \frac{dt_f(\tau,\tau^\prime)}{d\tau^\prime} = 2E_2(\tau,\tau^\prime) \approx 1.666\exp(-1.666(\tau - \tau^\prime)) \]

In words, that means that the flux sees a layer that is effectively 5/3 times thicker, compared with the radiance case at \(\mu=1\) (straight up/down), for both the flux transmission and the weighting function \(dt_f\). Be sure you understand why this make physical sense.

50.3.4. Matching with Wallace and Hobbs eq. 4.46

Convince yourself that Wallace and Hobbs eq. 4.46

\[ T_\nu^f \simeq e^{-\tau_\nu / \bar{\mu}} \]

where

\[ \frac{1}{\bar{\mu}} \equiv \sec 53^{\circ}=1.66 \]

is identical to (50.14)

and that their equation 4.56:

\[ \left(\frac{d T}{d t}\right)_\nu=-\frac{\pi}{c_p} k_\nu r B_\nu(z) \frac{e^{-\tau_\nu / \bar{\mu}}}{\bar{\mu}} \]

is actually:

\[ \left(\frac{d T}{d t}\right)_\nu=-\frac{\pi}{c_p} k_\nu r B_\nu(z) dt_f \]

50.4. Summary

Bottom line: the longwave radiation code in a climate model:

• Calculates \(\Delta \tau_\lambda\), the optical thickness for every model layer \(n\) and for 30-40 different wavelengths

• Uses the temperature and the Planck function to get \(B_\lambda (T_n)\) in each layer and wavelength

• Starts from the lowest model level for the upward flux, and the highest model level for the downward flux, and finds \(F_{\lambda \uparrow}\) and \(F_{\lambda \downarrow}\) for each layer and wavelength using the angle-integrated Schwartzchild equation. Integrating (50.12) across layer \(n\) assuming constant \(T_{layer}\) we get the upward flux through the top of layer \(n\) as:

  (50.16)\[ F_{\lambda \uparrow n} = F_{\lambda \uparrow n-1} + B_\lambda(T_n)(1 - t_f(\Delta \tau_n)) \]

  and the downward flux (downward negative) though the bottom of layer \(n\) is:

  (50.17)\[ F_{\lambda \downarrow n} = F_{\lambda \downarrow n+1} - B_\lambda(T_n)(1 - t_f(\Delta \tau_n)) \]

• Gets the net flux \(F_{net, n} = F_{\uparrow n} + F_{\downarrow n }\) through every level by summing over all wavelengths

• Gets the heating rate in K/second for each layer

  \[ \frac{dT_n}{dt} = -\frac{1}{\rho c_p} \frac{dF_{net\,n}}{dz} \]

• Update each layer temperature for the timestep (say \(\Delta t = 20\) minutes)

  \[ \frac{dT_{t+1,n}}{dt} = T_{t,n} \Delta t \]