= Solution
Use the optical-depth weighting stipulated in the paper. A jump in the <Planck function> at $\tau_1$ gives $dB/d\tau=(B_1-B_0)\delta(\tau-\tau_1)$, so
$$
\boxed{I_0=\pi\left[B_0+(B_1-B_0)e^{-\tau_1}\right]
=\pi B_0(1-e^{-\tau_1})+\pi B_1e^{-\tau_1}.}
$$
The upper isothermal layer emits $\pi B_0(1-e^{-\tau_1})$ and transmits the lower region's blackbody <spectral irradiance> with factor $e^{-\tau_1}$. The lower layer and black surface have the same $B_1$, so their combined radiation entering the upper layer is $\pi B_1$, independent of how its emission is divided between the lower gas and surface.
More generally, integration by parts rewrites the supplied expression as
$$
I_0=\pi B(\tau_s)e^{-\tau_s}+\pi\int_0^{\tau_s}B(\tau)e^{-\tau}\,d\tau.
$$
The terms are attenuated surface emission and the atmosphere's own emission. In the original boundary-plus-gradient form, $\pi B_0$ is the isothermal reference flux and the integral corrects for temperature variation below the top. Optical depths of order unity contribute most strongly. These formulas follow the paper's specified exponential weighting; a full angular <radiative transfer> calculation uses ray-angle-dependent optical paths rather than identifying every ray with a vertical path.
Back to article page
