12.4: Correction of the Temperature Distribution and Radiative Equilibrium

Having created an accurate representation of the radiation field from the previously obtained physical structure, we must see how well the solution conforms to the condition of radiative equilibrium. Departures of the radiation field from that required to satisfy radiative equilibrium will form the basis for correcting the temperature distribution throughout the atmosphere. We have developed the concept of radiative equilibrium several times in this book and most recently in Chapter 10 [equations \ref{10.4.4} and \ref{10.4.5}] as \[\begin{aligned}
\int_0^{\infty} \kappa_\nu \rho\left\{J_\nu\left(\tau_0\right)-B_\nu\left[T_0\left(\tau_0\right)\right]\right\} d\nu & =0 \\
\mathbf{F} & =\int_0^{\infty} F_\nu d\nu=\frac{\sigma T_e^4}{\pi}
\end{aligned}\label{12.4.1}\]

Even though these two conditions are logically equivalent, their utilization for generating a temperature correction scheme will be quite different. Although a substantial number of temperature correction schemes have been developed during the last 40 years, we describe only two. The first is chosen for its simplicity and historical interest while the second represents the most widely used method in contemporary use.

a. Lambda Iteration Scheme

The first of equations \ref{12.4.1} is obtained by setting the total flux derivative to zero. In general, the radiation field obtained from our approximate structure will not satisfy this expression. If we assume that the reason for this is that the temperature used to evaluate the local Planck function is incorrect, we can replace the temperature with a first-order Taylor series expansion about the current temperature. Thus, \[\begin{aligned}
\int_0^{\infty} \kappa_\nu J_\nu\left(\tau_0\right) d\nu & =\int_0^{\infty} \kappa_\nu\left\{B_\nu\left[T_0\left(\tau_0\right)+\delta T\right] d\nu\right. \\
& =\int_0^{\infty} \kappa_\nu B_\nu\left[T_0\left(\tau_0\right)\right] d\nu+\delta T \int_0^{\infty} \kappa_\nu \frac{\partial B_\nu}{\partial T} d\nu
\end{aligned}\label{12.4.2}\]

or solving for the temperature correction we have \[\delta T\left(\tau_0\right)=\frac{\int_0^{\infty} \kappa_\nu\left[J_\nu\left(\tau_0\right)-B_\nu\left(\tau_0\right)\right] d\nu}{\int_0^{\infty} \kappa_\nu\left[\partial B_\nu\left(\tau_0\right) / \partial T\right] d\nu}\label{12.4.3}\]

This is known as the Λ iteration scheme since \(\mathrm{J}_{\mathrm{n}}\left(\mathrm{t}_0\right)=\Lambda\left[B_{\nu}\left(\tau_0\right)\right]\) [see equation \ref{10.1.16}]. The method yields suitable corrections to the temperature distribution as long as the opacity \(\kappa_\nu\) is decidedly nongray. However, as one moves deeper and deeper in the atmosphere, \(\mathrm{J}_\nu\rightarrow\mathrm{B}_\nu\) and the integrand vanishes for all frequencies. Thus, this method relies on the departure of the source function from the value it would have in statistical equilibrium to provide corrections to the local temperature. So while the method may produce a useful temperature correction near the surface, the correction will become smaller and smaller as one descends into the atmosphere. This fact will be reflected in the rate at which the atmosphere converges to a self-consistent value. Indeed, it may become difficult to even know when meaningful convergence has been achieved. To make matters worse, equation \ref{12.4.3} guarantees − in principle − that a self-consistent atmosphere with zero total flux derivative can be calculated. However, it may not have the desired flux, \(\sigma T_{\mathrm{e}}^4 / \pi\). Thus, we should look for a method for correcting the temperature that employs the second of equations \ref{12.4.1} as well as the first. Such a scheme is due to E. Avrett and M. Krook³ although it is more lucidly described by D. Mihalas¹ (pp. 35 - 39).

b. Avrett-Krook Temperature Correction Scheme

Since the temperature correction scheme is to form the basis for an iteration algorithm, it is not essential that it produce the correct temperature the first time it is applied. However, repeated application should produce a series of temperature distributions which approach the one that is correct for radiative equilibrium. Thus, all temperature corrections must vanish asymptotically as the sequence approaches radiative equilibrium. This is the only essential criterion for an iteration scheme. Therefore, it is not necessary to justify all assumptions made in establishing the iterative equations as long as the final result converges to a temperature distribution that is consistent with radiative equilibrium.

The beauty of the Avrett-Krook scheme is that it simultaneously uses both expressions of radiative equilibrium as given in equations \ref{12.4.1}. There are two ways of correcting the temperature distribution. The first is the obvious one of simply changing the value of the temperature at some given value of the independent variable \(\tau_0\). This is the approach taken by the Λ-iteration scheme. A second way to find an improved temperature distribution is to find the value of the independent variable, in this case the reference optical depth, for which the given temperature is the correct temperature. This approach amounts to inverting the problem and treating the temperature as the independent variable and perturbing \(\tau_0\).

The Avrett-Krook scheme does both, using one statement of radiative equilibrium to calculate a temperature correction and the other condition of radiative equilibrium to find a new value of the optical depth at which the corrected temperature is to be applied. Thus, both temperature and optical depth become independent variables in the perturbation calculation. The perturbation equations for the temperature are very similar to the Λ-iteration equations and therefore provide good corrections near the surface. The perturbation equations for the optical depth yield small corrections near the surface, but become significant at larger optical depths where the Λ-iteration scheme is ineffective. Thus, the combination yields a temperature correction scheme which converges fairly quickly throughout the entire atmosphere. Unfortunately, the resulting temperature distribution will not directly give the corrected temperatures at the reference optical depths. However, the appropriate temperatures at the reference optical depths can be obtained from the new temperature distribution by interpolation.

The basic approach is to express both the correct temperature and the optical depth in terms of the given values and a first order correction to them, namely, \[\tau_0=t+\lambda \tau^{(1)}(t) \quad T\left(\tau_0\right)=T^{(0)}(t)+\lambda T^{(1)}(t)\label{12.4.4}\]

The parameter λ simply measures the order of significance for the particular term and will eventually be set to unity. Substitution of these expressions into the equation of transfer will produce similar corrections in the parameters that describe the radiation field so that \[\begin{aligned}
I_\nu\left(\mu, \tau_0\right) & =I_\nu^{(0)}(t, \mu)+\lambda I_\nu^{(1)}(t, \mu) \\
J_\nu\left(\tau_0\right) & =J_\nu^{(0)}(t)+\lambda J_\nu^{(1)}(t) \\
F_\nu\left(\tau_0\right) & =F_\nu^{(0)}(t)+\lambda F_\nu^{(1)}(t)
\end{aligned}\label{12.4.5}\]

We can expand the normalized opacity [equation \ref{12.2.3}] and the Planck function in a Taylor Series in \(t\) and \(T\), respectively, and get \[\begin{aligned}
k_\nu\left(\tau_0\right) & =k_\nu^{(0)}(t)+\lambda \tau^{(1)}(t) \frac{d k_\nu(t)}{d t} \\
B_\nu\left[T\left(\tau_0\right)\right] & =B_\nu\left[T^{(0)}(t)\right]+\lambda T^{(1)}(t) \frac{d B_\nu\left[T^{(0)}(t)\right]}{d T}
\end{aligned}\label{12.4.6}\]

For simplicity, from now on we denote differentiation with respect to optical depth and temperature by \[\frac{d}{d t}=^{\prime} \quad \frac{d}{d T}=\cdot\label{12.4.7}\]

In addition, for clarity we ignore scattering and treat the problem of pure absorption only. Later we give the perturbation equations appropriate for a source function that includes scattering, justifying the results on physical grounds alone.

Perturbed Equation of Radiative Transfer The general nongray equation of transfer for a plane-parallel atmosphere for the case of pure absorption is \[\mu \frac{d I\left(\mu, \tau_0\right)}{d \tau_0}=k_\nu\left\{I_\nu\left(\mu, \tau_0\right)-B_\nu\left[T\left(\tau_0\right)\right]\right\}\label{12.4.8}\]

If we insert the expansions given by equations \ref{12.4.4} and \ref{12.4.5} into this equation and ignore all second order terms (i.e., terms involving λ²), we get \[\begin{aligned}
\mu \frac{d I_\nu^{(0)}}{d t}-\lambda \mu \tau’^{(1)} \frac{d I_\nu^{(0)}}{d t}+\mu \lambda \frac{d I_\nu^{(1)}}{d t}= & \left(k_\nu^{(0)}+\lambda \tau^{(1)} k’_{v^{(0)}}\right)\left(I_\nu^{(0)}-B_\nu^{(0)}\right) \\
& +\lambda k_\nu^{(0)}\left(I_\nu^{(1)}+T^{(1)} \dot{B}_\nu^{(0)}\right)
\end{aligned}\label{12.4.9}\]

Since this equation must hold for any value of λ, we can separate the zeroth- and first-order terms. The zeroth-order equation is then \[\mu \frac{d I_\nu^{(0)}}{d t}=k_\nu^{(0)}\left(I_\nu^{(0)}-B_\nu^{(0)}\right)\label{12.4.10}\]

We can use this result to eliminate \(\mathrm{dI}_\nu{ }^{(0)} / \mathrm{dt}\) from the first-order equation so that it becomes

\[\mu \frac{d I_\nu^{(1)}}{d t}=k_\nu^{(0)}\left(I_\nu^{(1)}+T^{(1)} \dot{B}_\nu^{(0)}\right)+ \left(\tau^{\prime(1)} k_\nu^{(0)}+\tau^{(1)} k_\nu^{\prime(0)}\right)\left(I_\nu^{(0)}-B_\nu^{(0)}\right)\label{12.4.11}\]

This equation can be solved by using the Eddington approximation to moments of the equation in a manner that should be familiar by now.

Forming the first two moments of equation \ref{12.4.11} (i.e., just integrating over all µ to obtain the first and multiplying by µ and integrating to get the second), we obtain \[\begin{aligned}
& F’^{(1)}=k_\nu^{(0)}\left(J_\nu^{(1)}+T^{(1)} \dot{B}_\nu^{(0)}\right)+\left(\tau’^{(1)} k_\nu^{(0)}+\tau^{(1)} k_\nu^{\prime(0)}\right)\left(J_\nu^{(0)}-B_\nu^{(0)}\right) \\
& K_\nu^{\prime(1)}=\frac{J_\nu^{\prime(1)}}{3}=k_\nu^{(0)} F_\nu^{(1)}+\left(\tau’^{(1)} k_\nu^{(0)}+\tau^{(1)} k_\nu^{\prime(0)}\right) F_\nu^{\prime(0)}
\end{aligned}\label{12.4.12}\]

In the second equation, we have already assumed that the Eddington approximation can be applied to the first-order perturbations as it is to the entire radiation field.

Tau Perturbation Equation Now we integrate the second of equations \ref{12.4.12} over all frequencies and get \[\frac{1}{3} \int_0^{\infty} \frac{J_\nu^{\prime(1)}}{k_\nu^{(0)}} d\nu=\int_0^{\infty} F_\nu^{(1)} d\nu+\tau^{\prime(1)} \int_0^{\infty} F_\nu^{(0)} d\nu+\tau^{(1)} \int_0^{\infty} \frac{k_\nu^{\prime(0)}}{k_\nu^{(0)}} F_\nu^{(0)} d\nu\label{12.4.13}\]

Requiring that \[J_\nu^{\prime(1)}(t)=0\label{12.4.14}\]

guarantees that the left-hand side of equation \ref{12.4.13} will vanish. The assumption stated by equation \ref{12.4.14} is justified by expediency alone. However, it is an assumption concerning the perturbation only and therefore can affect only the rate of convergence. There may be some instances where this approximation should be replaced. However, to do so, we must know something additional about the problem.

The first term on the right-hand side of equation \ref{12.4.13} is just the integrated flux error so that \[\int_0^{\infty} F_\nu^{(1)}(t) d\nu=\mathbf{F}-F^{(0)}\label{12.4.15}\]

where \[\mathbf{F}=\frac{\sigma T_e^4}{\pi}\label{12.4.16}\]

With this, we can rewrite equation \ref{12.4.12} as a first-order linear differential equation for the perturbed optical depth \[\tau^{\prime(1)}+\tau^{(1)} \int_0^{\infty} \frac{k_\nu^{\prime(0)}}{k_\nu^{(0)}} \frac{F_\nu^{(0)}}{F^{(0)}} d\nu=1-\frac{\mathbf{F}}{F^{(0)}}\label{12.4.17}\]

All that remains is to specify a boundary condition for the solution of the equation. An appropriate condition is \[\tau^{(1)}(0)=0\label{12.4.18}\]

While this condition appears to be arbitrary, it anticipates the result for the \(T\) perturbation which will provide the majority of the correction at the surface. The boundary condition given in equation \ref{12.4.18} will ensure that the tau corrections are small near the surface and thus will not compete heavily with the \(T\) corrections.

Temperature Perturbation Equation To obtain the \(T\) perturbation equation, we begin with the first of equations \ref{12.4.12}. Since we required that the derivative of the perturbed mean intensity \(\mathrm{J}_\nu^{(1)}\) be zero at all frequencies and depths [equation \ref{12.4.14}], we may get the last term on the right-hand side of the first of equations \ref{12.4.12} from the second equation, so that \[\tau^{\prime(1)} k_\nu^{(0)}+\tau^{(1)} k_\nu^{\prime(0)}=-\frac{k_\nu^{(0)} F_\nu^{(1)}}{F_\nu^{(0)}}\label{12.4.19}\]

That same assumption on the derivative of the perturbed mean intensity will require that \[J_\nu^{(1)}(t)=\text { const }=J_\nu^{(1)}(0)=a F_\nu^{(1)}(0)\label{12.4.20}\]

The last term implies the Eddington approximation; so that \(a\) is usually taken to be ½. However, some authors use somewhat different values for a based on empirical work. As with any iteration scheme, one that works is a good one. Remembering that we have assumed a boundary condition on the \(\tau^{(1)}\) - equation of \(\tau^{(1)}=0\), we see that equations \ref{12.4.17} and \ref{12.4.19} give \[\tau’^{(1)}(0)=-\frac{F_\nu^{(1)}(0)}{F_\nu^{(0)}(0)}=1-\frac{\mathrm{F}}{F^{(0)}(0)}\label{12.4.21}\]

Thus, we may obtain the perturbed\nualue for J as \[J_\nu^{(1)}(t)=-a F_\nu^{(0)}(0)\left[1-\frac{\mathbf{F}}{F^{(0)}(0)}\right]=\text { const }\label{12.4.22}\]

Inserting this result and equation \ref{12.4.19} into the first of equations \ref{12.4.12} we get \[\begin{aligned}
F_\nu^{\prime(1)}(t) & =k_\nu^{(0)}(t)\left\{a F_\nu^{(0)}(0)\left[1-\frac{\mathbf{F}}{F^{(0)}(0)}\right]+T^{(1)}(t) \dot{B}_\nu^{(0)}\right\} \\
& =-\frac{k_\nu^{(0)} F_\nu^{(1)}}{F_\nu^{(0)}}\left[J_\nu^{(0)}(t)-B_\nu^{(0)}\right]
\end{aligned}\label{12.4.23}\]

From the definition of \(\mathrm{F}_\nu^{\prime(1)}\) we know that \[F_\nu^{\prime(1)}(t)=\frac{d}{d t}\left[\mathrm{~F}-F_\nu^{(0)}(t)\right]=-F_\nu^{\prime(0)}\label{12.4.24}\]

Incorporating this into equation \ref{12.4.23}, integrating over all frequencies, and remembering that the condition of radiative equilibrium applies to the zeroth-order equations, we finally get the perturbation equation for the temperature as \[T^{(1)}(t)=\frac{\mathbf{F} \int_0^{\infty} k_\nu^{(0)}(t)\left\{a\left[1-F_\nu^{(0)}(0) / \mathbf{F}\right]+\left[J_\nu^{(0)}(t)-B_\nu^{(0)}(t)\right] / F_\nu^{(0)}(t)\right\} d\nu}{\int_0^{\infty} \dot{B}_\nu^{(0)}(t) k_\nu^{(0)}(t) d\nu}\label{12.4.25}\]

We now have expressions for the temperature corrections \(T^{(1)}(t)\) and the corrected\nualues of the optical depth \(t+\tau^{(1)}\), to which they are to apply. Interpolation of this temperature distribution back onto the original optical depth scale completes the temperature correction procedure. A comparison of equations \ref{12.4.25} and \ref{12.4.3} shows that the Avrett-Krook temperature correction equation is indeed\nuery close to the Λ-iteration equation. However, an additive constant appears in the Avrett-Krook equation which ensures that the corrections will converge to the correct flux F.

Perturbation Equations Including Scattering The inclusion of scattering significantly complicates the algebra of deriving the perturbation equations, but not the concept. However, the essence of the problem can be seen without suffering through the algebra of the derivation. Consider a very general source function such as that given in equation \ref{10.1.7}. The parameter \(\varepsilon_\nu\) is a measure of the fraction of photon interactions that can be viewed as pure absorptions. Thus, \(1-\varepsilon_\nu\) is the relative fraction of scatterings. Since at the microscopic level scattering is a fully conservative process, we should expect it to have no influence on the physical structure of the atmosphere. Scattering decouples the radiation field from the physical domain of the gas. Thus, any temperature correction procedure will become less well defined for an atmosphere where the opacity becomes more nearly gray.

To carry out the perturbation analysis, we must add a perturbation equation for \(\varepsilon_\nu(\tau_0)\) similar to equations \ref{12.4.6}. It could take the form \[\epsilon_\nu\left(\tau_0\right)=\epsilon_\nu^{(0)}(t)+\lambda \tau^{(1)}(t) \epsilon_\nu^{\prime(0)}(t)\label{12.4.26}\]

As with the opacity, an assumption is made that the derivatives with respect to optical depth are more important than the derivatives with respect to temperature. The appropriate equation of radiative transfer analogous to equation \ref{12.4.8} is then \[\mu \frac{d I\left(\mu, \tau_0\right)}{d \tau_0}=k_\nu\left\{I_\nu\left(\mu, \tau_0\right)-\epsilon_\nu\left(\tau_0\right) B_\nu\left[T\left(\tau_0\right)\right]-\left[1-\epsilon_\nu\left(\tau_0\right)\right] J_\nu\left(\tau_0\right)\right\}\label{12.4.27}\]

where \(\kappa_\nu\) is now defined by \[k_\nu\left(\tau_0\right) \equiv \frac{\kappa_\nu\left(\tau_0\right)+\sigma_\nu\left(\tau_0\right)}{\kappa_0\left(\tau_0\right)+\sigma_0\left(\tau_0\right)}\label{12.4.28}\]

Development of the two moment equations analogous to equations \ref{12.4.12} will show that the second is unchanged by the presence of scattering. This leads to the happy result that the tau perturbation equation is also unchanged, so that equation \ref{12.4.17} and its solution are correct for the more general case including scattering.

The presence of scattering does modify the first moment equation of equations \ref{12.4.12}. This yields a somewhat different temperature perturbation equation from equation \ref{12.4.25}. With scattering, it takes the form \[T^{(1)}(t)=\frac{\mathbf{F} \int_0^{\infty} k_\nu^{(0)}(t)\left\{a \epsilon_\nu^{(0)}(t)\left[1-F_\nu^{(0)}(0) / \mathbf{F}\right]+\left[J_\nu^{(0)}(t)-B_\nu^{(0)}(t)\right]\right.{\left.\times\left[\tau^{(1)}(t) \epsilon_\nu^{\prime(0)}(t)+\epsilon_\nu^{(0)}(t) \mathbf{F} / F_\nu^{(0)}(t)\right] / \mathbf{F}\right\} d\nu}}{\int_0^{\infty} \dot{B}_\nu^{(0)}(t) k_\nu^{(0)}(t) \epsilon_\nu^{(0)}(t) d\nu}\label{12.4.29}\]

In the limit of pure absorption where \(\varepsilon_\nu\rightarrow1\), we recover immediately equation \ref{12.4.25}. As we approach the limit of a pure scattering atmosphere \(\varepsilon_\nu\rightarrow0\). All terms in the numerator of equation \ref{12.4.29} clearly vanish. Unfortunately so does the denominator, leaving the asymptotic behavior of T⁽¹⁾ in doubt. An application of L'Hospital's rule shows that the temperature correction terms indeed formally go to zero for the case of pure scattering. However, many of the terms of equation \ref{12.4.29} are difficult to calculate numerically so that the practical result of increased scattering will be to at first slow the rate of convergence of the iteration procedure. The iteration procedure will become unstable as the amount of scattering becomes very large. This is not surprising since the instability merely reflects the decoupling of the radiation field from the physical structure of the atmosphere.

Equations \ref{12.4.17} and \ref{12.4.29} provide the mechanism by which departures from radiative equilibrium can be translated to an improved temperature distribution. With this temperature distribution, we may return to the beginning of this chapter and re-compute the structure and improved radiation field of the atmosphere. The entire process can be iterated until radiative equilibrium is satisfied at the appropriate level.