AN ANALYTIC RADIATIVE–CONVECTIVE MODEL FOR PLANETARY ATMOSPHERES

The flux profiles of thermal radiation and the atmospheric temperature profile are key parameters in a radiative–convective model, and these quantities are inter-related. The forms of these profiles are different in the regions of the atmosphere that are convection-dominated versus radiation-dominated. We start by writing the general equations that describe the relationships between these key quantities.

In a one-dimensional, plane-parallel atmosphere, the two-stream Schwarzschild equations for the upwelling and downwelling thermal radiative fluxes (F⁺ and F⁻, respectively) are (Andrews 2010, p. 84)

$$\begin{equation} \frac{{dF^ + }}{{d\tau }} = D({F^ + - \pi B} ) \end{equation} \tag{ 1 }$$

$$\begin{equation} \frac{{dF^ - }}{{d\tau }} = - D({F^ - - \pi B} ), \end{equation} \tag{ 2 }$$

where τ is the gray infrared optical depth (0 at the top of the atmosphere and increasing toward larger pressures), which is used as a vertical coordinate. For brevity, we will simply refer to these fluxes as "thermal fluxes" for the remainder of this paper. The variable B is the integrated Planck function, with

$$\begin{equation} \pi B\left(\tau \right) = \sigma T^4 \left(\tau \right), \end{equation} \tag{ 3 }$$

where σ is the Stefan–Boltzmann constant (5.67 × 10⁻⁸ W m⁻² K⁻⁴) and T(τ) is the atmospheric temperature profile. The parameter D is the so-called diffusivity factor, which accounts for the integration of the radiances over a hemisphere. The value of D is often taken as 1.66, which compares well with numerical results (Rodgers & Walshaw 1966; Armstrong 1968). Others commonly set D = 3/2 (Weaver & Ramanathan 1995), or D = 2, which is the hemi-isotropic approximation. For clarity, we note that others (e.g., Andrews 2010 or Pierrehumbert 2010) sometimes absorb the value of D into their definition of τ. We choose to leave it in our expressions so that the diffusivity approximation is not hidden. As we shall see in Section 2.3, Equations (1) and (2) allow us to solve for the upwelling and downwelling thermal flux profiles if T(τ) and a boundary condition are specified.

The net thermal flux, F_net, is given by

$$\begin{equation} F_{{\rm net}} = F^ + - F^ - \end{equation} \tag{ 4 }$$

so that we can combine Equations (1) and (2) to yield (see Schuster 1905)

$$\begin{equation} \frac{{d^2 F_{{\rm net}} }}{{d\tau ^2 }} - D^2 F_{{\rm net}} = - 2\pi D\frac{{dB}}{{d\tau }}, \end{equation} \tag{ 5 }$$

which is a differential equation that can be used to solve for B(τ), and thus the atmospheric temperature profile, if F_net(τ) and a boundary condition is known.

While the vertical coordinate of the equations governing the transfer of thermal radiation is optical depth, we must relate this to atmospheric pressure, which is the natural physical vertical coordinate of planetary atmospheres. A vertical pressure coordinate is particularly useful and unifying between atmospheres when we consider lapse rates in tropospheres (Section 2.3).

We take the relation between gray thermal optical depth and atmospheric pressure to be given by a power law in the form

$$\begin{equation} \tau = \tau _0 \left({\frac{p}{{p_0 }}} \right)^n , \end{equation} \tag{ 6 }$$

where p is atmospheric pressure, τ₀ is the gray infrared optical depth integrated down from the top of the atmosphere to the reference pressure p₀, and n is a parameter that controls the strength of the scaling. In the simplest case, when an absorbing gas is well mixed and the opacity does not depend strongly on pressure, we will have n = 1. This can physically correspond to Doppler broadening. A common scenario is to have a well-mixed gas providing collision-induced opacity (e.g., H₂ in gas giant atmospheres in the solar system) or pressure-broadened opacity so that n = 2. Typically, n takes a value between 1 and 2, as has been parameterized into some radiative models (Heng et al. 2012). In some cases, where the mixing ratio of the absorbing gas depends strongly on pressure, larger values of n have been proposed, such as for water vapor in Earth's lower troposphere (Satoh 2004, p. 373; Frierson et al. 2006).

In the convection-dominated region of a planetary atmosphere, we take the temperature–pressure profile to be similar to a dry adiabat, although somewhat less steep because of, for example, the latent heat released by condensation of volatiles during convective uplift. By specifying T(p), we can then derive the expressions for the upwelling and downwelling thermal flux profiles, which we can join to the profiles from the radiatively dominated region of the atmosphere.

The dry adiabatic temperature variation in the lower, convective part of a troposphere is given by Poisson's adiabatic state equation (e.g., Wallace & Hobbs 2006, p. 78):

$$\begin{equation} T = T_0 \left({\frac{p}{{p_0 }}} \right)^{\left({\gamma - 1} \right)/\gamma }. \end{equation} \tag{ 7 }$$

Here T₀ is a reference temperature at a reference pressure p₀, and γ is the ratio of specific heats at constant pressure (c_p) and volume (c_v):

$$\begin{equation} \gamma = \frac{{c_p }}{{c_v }}. \end{equation} \tag{ 8 }$$

Kinetic theory also allows us to relate the ratio of specific heats to the degrees of freedom, N, for the primary atmospheric constituent(s), where

$$\begin{equation} \gamma = 1 + \frac{2}{N}.\end{equation} \tag{ 9 }$$

Most appreciable atmospheres in the solar system are dominated by linear diatomic gases, such as H₂ (Jupiter, Saturn, Uranus, and Neptune), N₂ (Titan), or an N₂–O₂ mixture (Earth). These molecules have three translational and two rotational degrees of freedom, so that N = 5. Thus, γ = 7/5 = 1.4 is the same for these worlds. In the case of Venus and Mars, whose atmospheres are dominated by CO₂, there are 3 translational, 2 rotational, and ∼2 vibrational degrees of freedom (weakly depending on temperature; Bent 1965), so N = 7 and γ = 1.3. Thus, in Equation (7), the dry adiabatic temperature T varies as p^0.3 and p^0.2 for diatomic- and CO₂-dominated atmospheres, respectively. These T–p scalings will also apply to exoplanet tropospheres.

Planetary tropospheres do not follow a true dry adiabat, though, because of the condensation of volatiles during the convection process, such as the effects of water on Earth or methane on Titan, or because the degrees of freedom for the primary atmospheric constituents vary with altitude. Thus, following Sagan (1962), we modify the temperature structure in the convection-dominated region of a planetary atmosphere as

$$\begin{equation} T = T_0 \left({\frac{p}{{p_0 }}} \right)^{\alpha \left({\gamma - 1} \right)/\gamma }. \end{equation} \tag{ 10 }$$

Here, α is a factor, typically around 0.6–0.9, which accounts for deviations from the dry adiabatic lapse rate, primarily because of latent heat release. Physically, α represents the average ratio of the true lapse rate in the planet's convective region to the dry adiabatic lapse rate. While this parameterized representation of the temperature–pressure profile works well for all of the worlds of the solar system with thick atmospheres, it breaks down for very moist atmospheres (Pierrehumbert 2010, p. 108), and for one-component condensable atmospheres, which will have temperature varying as 1/ln (p/p₀) according to the Clausius–Clapeyron relation (Satoh 2004, p. 254).

By solving Equation (6) for p/p₀ and inserting this into Equation (10), we can rewrite the temperature profile with the gray infrared optical depth as the vertical coordinate

$$\begin{equation} T = T_0 \left({\frac{\tau }{{\tau _0 }}} \right)^{\alpha \left({\gamma - 1} \right)/n\gamma } = T_0 \left({\frac{\tau }{{\tau _0 }}} \right)^{\beta /n} , \end{equation} \tag{ 11 }$$

where we define β = α(γ − 1)/γ. Inserting this temperature profile into Equation (1) and solving (see the Appendix), we obtain an expression for the upwelling thermal flux in the convective portion of the atmosphere (see also Mitchell et al. 2006; Mitchell 2007; Caballero et al. 2008),

$$\begin{equation} F^ + \left(\tau \right) = \sigma T_0^4 e^{ - D\left({\tau _0 - \tau } \right)} + D\sigma T_0^4 \int_\tau ^{\tau _0 } {\left({\frac{{\tau ^{\prime}}}{{\tau _0 }}} \right)^{4\beta /n} } e^{ - D\left({\tau ^{\prime} - \tau } \right)} d\tau ^{\prime}, \end{equation} \tag{ 12 }$$

where we have used the boundary condition F⁺(τ₀) = σT⁴₀, which assumes that the reference location in the atmosphere is either a solid surface or sufficiently deep so that the opacity is large enough to drive the upwelling and downwelling thermal fluxes toward that of a blackbody radiating at the reference temperature. Using the upper incomplete gamma function Γ (defined in the Appendix), we can write Equation (12) as

$$\begin{eqnarray} F^ + (\tau ) &=& \sigma T_0^4 e^{D\tau } \bigg[ e^{ - D\tau _0 } + \frac{1}{({D\tau _0 })^{4\beta /n} }\bigg(\Gamma \bigg({1 + \frac{{4\beta }}{n},D\tau } \bigg)\nonumber\\ && - \Gamma \bigg({1 + \frac{{4\beta }}{n},D\tau _0 } \bigg) \bigg) \bigg] , \end{eqnarray} \tag{ 13 }$$

which provides an analytic expression for the upwelling thermal flux in the convective portion of the atmosphere. The incomplete gamma function ought to be considered no different from a sine or cosine, in the sense that it is evaluated with a single function in modern programming languages, such as MATLAB, IDL, or Python, so that the flux can be computed readily, if all input parameters are specified. Note that Equations (12) and (13) only apply in the region τ_rc ⩽ τ ⩽ τ₀, where τ_rc is the optical depth at the boundary between the convection-dominated region and the radiation-dominated region. Also, in deriving these expressions, we have assumed that n, γ, and α are constant parameters through the convection-dominated region. Physically, the first expanded term, $\sigma T_0^4 e^{ - D(\tau _0 - \tau)}$, in Equation (13) is the upwardly attenuated contribution to the upwelling thermal flux from the reference level at temperature T₀. The second expanded term is the flux contribution from the atmosphere above the reference level.

A similar line of argument (see the Appendix) allows us to write an expression for the downwelling thermal flux in the convective region of the atmosphere

$$\begin{eqnarray} F^ - \left(\tau \right) &=& F^ - \left({\tau _{{\rm rc}} } \right)e^{ - D\left({\tau - \tau _{{\rm rc}} } \right)} \nonumber\\ &&+ D\sigma T_0^4 \int_{\tau _{{\rm rc}} }^\tau {\left({\frac{{\tau ^{\prime}}}{{\tau _0 }}} \right)^{4\beta /n} } e^{ - D\left({\tau - \tau ^{\prime}} \right)} d\tau ^{\prime}, \end{eqnarray} \tag{ 14 }$$

where the boundary condition is that the downwelling flux at the top of the convection-dominated region is equal to the downwelling flux coming from the radiation-dominated region above. While this equation is not required in order to obtain the solution to the radiative–convective thermal structure, it does offer a means to compute the downwelling thermal flux in the convective region.

Above the convection-dominated region of the planet's atmosphere, the temperature profile tends toward radiative equilibrium. In this region, emitted thermal radiation, absorbed stellar radiation, and any source of energy from the planet's interior are all in balance. Thus, by writing an expression for the profile of absorbed stellar flux, we can solve for the thermal flux profiles and the temperature profile using the results from Section 2.1.

The planetary atmospheres of the solar system exhibit absorption of solar radiation deep in the atmosphere (in the troposphere or at the surface) as well as in the stratosphere. The latter can lead to an inversion in the atmospheric temperature profile, while the former can cause the upper part of the troposphere to be radiatively dominated to considerable depth, as in the case of Titan. Consequently, for generality it is necessary to distribute the stellar flux across two short-wave channels, and write the net absorbed stellar flux, F^☉_net, as

$$\begin{equation} F_{{\rm net}}^ \odot \left(\tau \right) = F_1^ \odot e^{ - k_1 \tau } + F_2^ \odot e^{ - k_2 \tau } , \end{equation} \tag{ 15 }$$

where F^☉₁ and F^☉₂are the top-of-atmosphere net absorbed stellar fluxes in the two channels, and k₁ and k₂ parameterize the strength of attenuation in these two channels, and are a ratio of the visible optical depth to the gray infrared optical depth. (Note that k₁ and k₂ effectively incorporate a spatial and time mean zenith angle.) This description of the absorbed stellar flux is similar to that of McKay et al. (1999) except that they only allowed for the attenuation of sunlight in one of their channels. At the top of the atmosphere (i.e., where τ = 0), the net absorbed stellar flux must obey

$$\begin{equation} F_{{\rm net}}^ \odot \left(0 \right) = (1 - A)\frac{{F^ \odot }}{4}, \end{equation} \tag{ 16 }$$

where A is the planet's Bond albedo, F^☉ is the stellar flux incident on the top of the planet's atmosphere, and the factor of four accounts for an averaging over both the day and night sides as well as over the illuminated hemisphere.

A temperature profile results from a balance of fluxes subject to boundary conditions. Balancing the net thermal flux with the absorbed stellar flux and any internal energy source requires that

$$\begin{equation} F_{{\rm net}} \left(\tau \right) = F_{{\rm net}}^ \odot \left(\tau \right) + F_i , \end{equation} \tag{ 17 }$$

where F_i is the energy flux from the planet's interior (e.g., ∼5 W m⁻² for Jupiter; Hanel et al. 1981), which is assumed to be independent of pressure. By combining Equations (5), (15), and (17), and using the boundary condition that there is no downwelling thermal radiation at the top of the atmosphere (i.e., where τ = 0), we obtain the temperature profile in the radiation-dominated region (i.e., where 0 ⩽ τ ⩽ τ_rc)

$$\begin{eqnarray} \sigma T^4 \left(\tau \right) &=& \frac{{F_1^ \odot }}{2}\left[ {1 + \frac{D}{{k_1 }} + \left({\frac{{k_1 }}{D} - \frac{D}{{k_1 }}} \right)e^{ - k_1 \tau } } \right]\nonumber\\ && + \frac{{F_2^ \odot }}{2}\left[ {1 + \frac{D}{{k_2 }} + \left({\frac{{k_2 }}{D} - \frac{D}{{k_2 }}} \right)e^{ - k_2 \tau } } \right]\nonumber\\ && + \frac{{F_i }}{2}\left({1 + D\tau } \right). \end{eqnarray} \tag{ 18 }$$

The upwelling and downwelling thermal flux profiles in the region where 0 ⩽ τ ⩽ τ_rc are then given by

$$\begin{eqnarray} F^ + \left(\tau \right) &=& \frac{{F_1^ \odot }}{2}\left[ {1 + \frac{D}{{k_1 }} + \left({1 - \frac{D}{{k_1 }}} \right)e^{ - k_1 \tau } } \right]\nonumber\\ && + \frac{{F_2^ \odot }}{2}\left[ {1 + \frac{D}{{k_2 }} + \left({1 - \frac{D}{{k_2 }}} \right)e^{ - k_2 \tau } } \right] + \frac{{F_i }}{2}\left({2 + D\tau } \right)\nonumber\\ \end{eqnarray} \tag{ 19 }$$

$$\begin{eqnarray} F^ - \left(\tau \right) &=& \frac{{F_1^ \odot }}{2}\left[ {1 + \frac{D}{{k_1 }} - \left({1 + \frac{D}{{k_1 }}} \right)e^{ - k_1 \tau } } \right] \nonumber\\ &&+ \frac{{F_2^ \odot }}{2}\left[ \!{1 + \frac{D}{{k_2 }} - \left({1 + \frac{D}{{k_2 }}} \right)e^{ - k_2 \tau } }\! \right] + \frac{{F_i }}{2}D\tau.\nonumber\\ \end{eqnarray} \tag{ 20 }$$

The formulae in the preceding sections describe our simple, one-dimensional radiative–convective climate model. In applying the model, 10 parameters are often fixed (p₀, T₀, n, γ, α, F^☉₁, F^☉₂, F_i, k₁, and k₂), which leaves two variables,τ₀ and τ_rc, that the model provides in a solution. However, if the total optical depth τ₀ is specified then one of the aforementioned parameters (typically T₀) can become a variable. The temperature, upwelling flux, and downwelling flux profiles follow Equations (18), (19), and (20), respectively, in the radiative region from the top of the atmosphere down to the optical depth at the radiative–convective boundary, τ_rc. For optical depths (or pressures) below this level in the atmosphere, the temperature profile follows the adiabat described by Equation (11), and the upwelling and downwelling flux profiles follow Equations (13) and (14), respectively.

The requirement that the temperature and upwelling flux profiles be continuous at the radiation–convection boundary will implicitly solve for two variables in our model (the downwelling flux profile is guaranteed to be continuous due to the boundary condition applied in Equation (14)). Thus, at τ_rc, the upward thermal fluxes given by Equation (13) (for the convective formulation) and Equation (19) (for the radiative formulation) must be equal. We also require the convective temperature (Equation (11)) and radiative temperature (Equation (18)) to be equal, so that

$$\begin{eqnarray} \sigma T^4 _0 \left({\frac{{\tau _{{\rm rc}} }}{{\tau _0 }}} \right)^{4\beta /n} &= & {\frac{{F_1^ \odot }}{2}\left[ {1 + \frac{D}{{k_1 }} + \left({\frac{{k_1 }}{D} - \frac{D}{{k_1 }}} \right)e^{ - k_1 \tau _{{\rm rc}} } } \right]} \nonumber\\ && + \frac{{F_2^ \odot }}{2}\left[ {1 + \frac{D}{{k_2 }} + \left({\frac{{k_2 }}{D} - \frac{D}{{k_2 }}} \right)e^{ - k_2 \tau _{{\rm rc}} } } \right] \nonumber\\ &&+ \frac{{F_i }}{2}\left({1 + D\tau _{{\rm rc}} } \right). \end{eqnarray} \tag{ 21 }$$

These two equalities result in analytic expressions. But these expressions contain transcendental functions and so must be solved numerically, for example, using Newton's scheme. Model parameters that remain after applying the requirements above must be specified or fit. Greater simplicity may be justified by the transparency of an atmosphere, so that one of the stellar channels can be removed (eliminating F^☉₂ and k₂). Also, an internal energy source is irrelevant for many worlds (e.g., Earth, Titan, and Venus), setting F_i = 0 in many cases. As mentioned before, the value of γ is set by the degrees of freedom associated with the primary atmospheric constituent, and α is around 0.6–0.9 for worlds in the solar system.

Using the thermal fluxes from Equations (13) and (14), as well as the parameterized net stellar flux from Equation (15), the convective flux, F_conv, can be easily computed according to

$$\begin{equation} F_{{\rm conv}} \left(\tau \right) = F_i + F_{{\rm net}}^ \odot \left(\tau \right) - ({F^ + \left(\tau \right) - F^ - (\tau )}) \end{equation} \tag{ 22 }$$

in the region τ_rc ⩽ τ ⩽ τ₀.

It is worth considering the simplest form of our radiative–convective model wherein there is no attenuation of sunlight (k₁ = k₂ = 0), so that

$$\begin{equation} F_{{\rm net}}^ \odot = F_1^ \odot + F_2^ \odot = (1 - A)\frac{{F_{}^ \odot }}{4} \end{equation} \tag{ 23 }$$

at all locations in the atmosphere, which reduces our radiative equilibrium expressions to those of a single short-wave channel. In the radiative region, the temperature and thermal flux profiles then simplify to

$$\begin{equation} \sigma T^4 \left(\tau \right) = \frac{{F_{{\rm net}}^ \odot + F_i }}{2}\left({1 + D\tau } \right) \end{equation} \tag{ 24 }$$

$$\begin{equation} F^ + \left(\tau \right) = \frac{{F_{{\rm net}}^ \odot + F_i }}{2}\left({2 + D\tau } \right) \end{equation} \tag{ 25 }$$

$$\begin{equation} F^ - \left(\tau \right) = \frac{{F_{{\rm net}}^ \odot + F_i }}{2}D\tau. \end{equation} \tag{ 26 }$$

For dealing with the convective region, the equalities that ensure temperature and upwelling flux continuity now become

$$\begin{equation} \sigma T^4 _0 \left({\frac{{\tau _{{\rm rc}} }}{{\tau _0 }}} \right)^{4\beta /n} = \frac{{F_{{\rm net}}^ \odot + F_i }}{2}\left({1 + D\tau _{{\rm rc}} } \right) \end{equation} \tag{ 27 }$$

and

$$\begin{eqnarray} &&\sigma T_0^4 e^{D\tau _{\rm rc} } \bigg[ e^{ - D\tau _0 } + \frac{1}{({D\tau _0 } )^{4\beta /n} }\bigg(\Gamma \bigg(1 + \frac{{4\beta }}{n},D\tau _{{\rm rc}} \bigg) \nonumber \\ &&\quad - \Gamma \bigg(1 + \frac{{4\beta }}{n},D\tau _0 \bigg) \bigg) \bigg] = \frac{{F_{{\rm net}}^ \odot + F_i }}{2}\bigg({2 + D\tau _{{\rm rc}} } \bigg). \end{eqnarray} \tag{ 28 }$$

The dependence on the Bond albedo, which is incorporated into F^☉_net, can be removed by scaling the temperatures to the equilibrium temperature, T_eq, which is defined by

$$\begin{equation} \sigma T_{{\rm eq}}^4 = (1 - A)\frac{{F^ \odot }}{4} + F_i. \end{equation} \tag{ 29 }$$

Note that the simple model described in this section is similar to an exercise discussed by Pierrehumbert (2010; see Problem 4.28, p. 311). However, we (1) generalized the temperature–pressure relationship in the convection-dominated region (Equation (10)), (2) generalized the relationship between optical depth and pressure (Equation (6)), and (3) use the upper incomplete gamma function to provide an analytic evaluation of the upwelling thermal radiative fluxes in the convection-dominated region.

Our simplest radiative–convective model (Section 2.6) allows us to make straightforward deductions about the location of radiative–convective boundary over a wide range of conditions. We can combine Equations (27) and (28) to yield an expression for τ_rc:

$$\begin{eqnarray} &&\left({\frac{{\tau _0 }}{{\tau _{{\rm rc}} }}} \right)^{4\beta /n} e^{ - D\left({\tau _0 - \tau _{{\rm rc}} } \right)} \bigg[ 1 + \frac{{e^{D\tau _0 } }}{{\left({D\tau _0 } \right)^{4\beta /n} }}\bigg(\Gamma \left({1 + \frac{{4\beta }}{n},D\tau _{{\rm rc}} } \right) \nonumber\\ &&\quad - \Gamma \left({1 + \frac{{4\beta }}{n},D\tau _0 } \right) \bigg) \bigg] = \frac{{2 + D\tau _{{\rm rc}} }}{{1 + D\tau _{{\rm rc}} }}, \end{eqnarray} \tag{ 30 }$$

which is independent of the net solar flux or the internal energy flux, and only depends on τ₀ and the parameter 4β/n. Figure 1 shows contours of solutions to Equation (30) for Dτ_rc over a range of values for τ₀ and 4β/n. Since typical values of 4β/n are between 0.3 and 0.5 for worlds in the solar system (taking n = 2), the shaded region in this figure demonstrates that Dτ_rc will typically be less than unity for realistic values of 4β/n, placing the "radiating level" or "emission level" of Dτ = 1 in the convective region of the atmosphere. As was discussed in Sagan (1969), worlds with Dτ_rc < 1 will have "deep" tropospheres and "shallow" stratospheres, whereas worlds with Dτ_rc > 1 will have "deep" stratospheres and "shallow" tropospheres.

Zoom In Zoom Out Reset image size

Figure 1 also shows that, in the limit that τ₀ ≫ 1 and τ₀ ≫ τ_rc, the value of τ_rc becomes independent of the value of τ₀. In this limit, the upwelling flux at the radiative–convective boundary is no longer sensitive to the radiation coming from the deepest atmospheric layers, and Equation (30) simplifies to

$$\begin{equation} \frac{{\Gamma \big({1 + \frac{{4\beta }}{n},D\tau _{{\rm rc}} } \big)}}{{\left({D\tau _{{\rm rc}} } \right)^{4\beta /n} e^{ - D\tau _{{\rm rc}} } }} = \frac{{2 + D\tau _{{\rm rc}} }}{{1 + D\tau _{{\rm rc}} }}, \end{equation} \tag{ 31 }$$

which only depends on the value of 4β/n. Figure 2 shows the solution to Equation (31), computed over a range of values for 4β/n, which, again, shows that, for common values of 4β/n, worlds will have deep tropospheres.

Zoom In Zoom Out Reset image size

Several previous authors (e.g., Sagan 1969; Weaver & Ramanathan 1995) have shown that the requirement for convective instability in a gray radiative equilibrium model without solar attenuation is

$$\begin{equation} \frac{{d\log T}}{{d\log p}} = \frac{{Dn\tau }}{{4\left({1 + D\tau } \right)}} > \frac{{\gamma - 1}}{\gamma } = \left({\frac{{d\log T}}{{d\log p}}} \right)_{{\rm ad}} , \end{equation} \tag{ 32 }$$

where the subscript "ad" indicates the (dry) adiabatic value. As mentioned earlier, γ typically has a value of 1.3–1.4, so that the right-hand side of this expression is typically between 0.23 and 0.29. Sagan (1969) took the solution to Equation (32) as defining the radiative–convective boundary, which is also shown in Figure 2 (taking α = 1). Sagan's solution is not the same as a true radiative–convective solution as he only ensures continuity in temperature across the radiative–convective boundary, whereas a realistic radiative–convective model must also maintain upwelling flux continuity across this boundary, which places our radiative–convective boundary more correctly higher in the atmosphere. At large values of Dτ_rc Sagan's solution agrees with ours because, in the optically thick limit, the upwelling flux approaches the local blackbody flux, so that temperature continuity and upwelling flux continuity are equivalent. Note, however, that high values of Dτ_rc are only achieved for values of 4β/n that are far larger than values in the solar system.

The simplest model without a stellar channel is unstable to convection, as shown in Figure 3. Here, the radiative equilibrium temperature profile (labeled as "no atten."), for an atmosphere with n = 2, has a portion that is unstable to convection (for γ = 1.4). We can increase the generality of our simplest model by allowing for a single stellar channel with attenuation, obtained from Equation (18) by eliminating terms in F^☉₂ and F_i, and by dropping the subscripts on the remaining stellar channel. Figure 3 demonstrates example temperature profiles (taking n = 2) for such a model for different values of k, which is the ratio of the stellar optical depth in the single channel to the gray thermal optical depth. The logarithmic temperature gradient, or lapse rate, in such a model is

$$\begin{equation} \frac{{d\log T}}{{d\log p}} = \frac{{nk\tau }}{4}\left[ {\frac{{({D^2 - k^2 } )e^{ - k\tau } }}{{kD + D^2 + ({k^2 - D^2 })e^{ - k\tau } }}} \right]. \end{equation} \tag{ 33 }$$

Note that, for k > D, this expression is strictly negative, and, thus, everywhere stable against convection, and the profile has a temperature inversion, as shown in Figure 3. This is similar to an argument in Pierrehumbert (2010, p. 212), albeit in different notation. For k = D, the temperature profile is isothermal and is stable against convection for all values of n. However, for a given n value, it is possible to have a value of k < D but still be stable against convection, although in this case there is no temperature inversion.

Zoom In Zoom Out Reset image size

We investigated this other threshold by examining the maximum value of the lapse rate. Figure 4 shows the maximum value of the lapse rate, which is only a function of k and n. This figure demonstrates that, for a given value of n and γ (which defines the dry adiabatic lapse rate), there is a threshold value of the relative stellar absorption, represented by k, above which the single stellar channel model is everywhere stable against convection. The thick lines in this figure indicate the values of k for which profiles would be unstable to convection as compared to dry adiabats in a CO₂-dominated atmosphere and an atmosphere dominated by a diatomic gas. For example, for the n = 2 case, the threshold value of k in a CO₂-dominated atmosphere is 0.2 and is 0.1 for an atmosphere dominated by a diatomic gas.

Zoom In Zoom Out Reset image size

In general, the gray infrared optical depth of the radiative–convective boundary in a model with only a single attenuated stellar channel and no internal heat flux depends on τ₀, k, and 4β/n, and can be computed from the expression

$$\begin{eqnarray} && \left({\frac{{\tau _0 }}{{\tau _{{\rm rc}} }}} \right)^{4\beta /n} e^{ - D\left({\tau _0 - \tau _{{\rm rc}} } \right)} \bigg[ 1 + \frac{{e^{D\tau _0 } }}{{\left({D\tau _0 } \right)^{4\beta /n} }}\bigg(\Gamma \left({1 + \frac{{4\beta }}{n},D\tau _{{\rm rc}} } \right)\nonumber\\ &&\quad - \Gamma \left({1 + \frac{{4\beta }}{n},D\tau _0 } \right) \bigg) \bigg] = \frac{{1 + D/k + \left({1 - D/k} \right)e^{ - k\tau _{{\rm rc}} } }}{{1 + D/k + \left({k/D - D/k} \right)e^{ - k\tau _{{\rm rc}} } }}, \nonumber\\ \end{eqnarray} \tag{ 34 }$$

which comes from combining the single-channel versions of Equations (18) and (19) with Equations (13) and (11). Note that this expression does not depend on the absorbed stellar flux. Contours of τ_rc as a function of τ₀ and k are shown in Figures 5(a) and (b) for two different values of 4β/n (appropriate for a CO₂-dominated atmosphere and for an atmosphere dominated by a diatomic gas, respectively).

Zoom In Zoom Out Reset image size

Figures 5(a) and (b) demonstrate that, for small values of k, we have that Dτ_rc < 1, which corresponds to a deep troposphere. However, for values of k larger than about 0.1, we have Dτ_rc > 1, so that the troposphere is shallow because much of the lower atmosphere is stabilized against convection by the absorption of stellar energy throughout the deep atmosphere, as opposed to depositing this energy abruptly at p₀, which will happen in the limit of small values of k. As we shall see later, cases corresponding to Figures 5(a) and (b) are applicable, respectively, to Venus and Titan. For certain combinations of τ₀ and k, there exists two values of τ_rc which satisfy Equation (34), but the larger of the two always represents an unphysical solution where the lapse rate in the radiative regime (from Equation (33)) exceeds that for the adiabat in the convective regime at a range of pressures above the radiative–convective boundary.

Gray radiative equilibrium models with a single stellar channel and an internal energy source have been used to understand certain properties of hot Jupiters (e.g., Hansen 2008; Guillot 2010). We can derive a similar model by dropping all terms in F^☉₂ in Equation (18) (and by dropping the subscripts on the remaining stellar channel). The resulting radiative temperature profile only depends on n, k, and the ratio of the absorbed stellar flux to the internal energy flux, F^☉/F_i. Figure 6 shows example temperature profiles for this model for two different values of k where thickened lines show regions unstable to convection. Both models assume strongly irradiated (F^☉/F_i = 10⁴) deep atmospheres with n = 2. In Figure 6, the deep convectively unstable zone at p/p₀ > 100 is caused by the internal flux of the giant planet rather than the stellar flux. In the case of the curve with k/D = 0.1, the stellar flux is absorbed around the base of a detached, convectively unstable zone where the atmosphere is optically thick in the infrared, and this gives rise to the separate, detached convective zone.

Zoom In Zoom Out Reset image size

The lapse rate in the radiative portions of such a model is given by

$$\begin{eqnarray} &&\hbox{\fontsize{8.7}{10.7}\selectfont{$\dsty\frac{{d\log T}}{{d\log p}} = \dsty\frac{{nk\tau }}{4}\left[ {\dsty\frac{{D^2 + ({F^ \odot /F_i } )({D^2 - k^2 } )e^{ - k\tau } }}{{kD\left({1 + D\tau } \right) + \left({F^ \odot /F_i } \right)({kD + D^2 + ({k^2 - D^2 } )e^{ - k\tau } } )}}} \right],$}}\nonumber\\ \end{eqnarray} \tag{ 35 }$$

which, depending on the values of n, k, and F^☉/F_i, can either be (1) everywhere stable against convection, (2) only unstable to convection deep in the atmosphere, or (3) unstable to convection in both deep layers and in a detached convective zone, as has been noticed in more complex, numerical models (Fortney et al. 2007).

To investigate the range of parameter space in which detached convective zones can form, we plotted contours of the gray infrared optical depth where Equation (35) first goes unstable to convection (assuming γ = 1.4, appropriate for an atmosphere dominated by a diatomic gas) for a range of values for k and F^☉/F_i (Figure 7). We see that, for large values of F^☉/F_i and for k/D larger than about 0.1, the profile only becomes stable to convection deep in the atmosphere, where the temperature profile is dominated by the internal energy flux.

Zoom In Zoom Out Reset image size

For larger values of k/D, the stellar flux is absorbed higher in the atmosphere where Dτ < 1 so that the energy can be radiated to space and there is no detached convective zone. However, for smaller values of k/D, the stellar flux is absorbed deeper down where Dτ > 1, so that a steep, super-adiabatic temperature gradient can cause a detached convective region. We find that for k/D smaller than about 0.1, a convective zone separate from that caused by the internal heat flux is possible. This limiting value of k can be seen in Figure 4 where the n = 2 curve crosses the dry adiabatic lapse rate for γ = 1.4.

The shaded region of Figure 7 shows the portion of parameter space where a detached convective region can form, which has a dependency on the internal heat flux. Note that such a region can only form in our model for F^☉/F_i larger than about 3, which is consistent with the results shown in Fortney et al. (2007). For a fixed absorbed stellar flux, F^☉, larger internal fluxes from a giant planet cause the temperature at depth to increase and the convective region extends higher and higher upward until it joins the detached convective zone, making a continuous region of convection.

Our simplified radiative–convective model can be applied to Venus, which, due to a lack of a stratospheric inversion combined with an opaque atmosphere (at thermal wavelengths), means that we can take k₁ and k₂ to be roughly equal to zero. Taking Venus' Bond albedo to be 0.76 (Moroz et al. 1985), the mean solar flux absorbed by Venus is about 160 W m⁻². Venus' atmosphere is primarily CO₂, so we take γ = 1.3, and comparing the average lapse rate in Venus' lower atmosphere to an average dry adiabat gives α = 0.8. We take T₀ and p₀ to be the surface temperature and pressure (730 K and 92 bar, respectively), set F_i = 0, and explore two cases where either n = 1 or n = 2, which bounds the two extremes of n that we would expect.

Figure 8 shows a comparison between Venus' observed temperature–pressure profile and our model-generated profiles for these cases. Our Venus model solves for τ_rc and τ₀, and finds these to be 1 and 400, respectively, for the n = 1 case. These values put the boundary between convection-dominated and radiation-dominated, p_rc, to be 0.2 bar. For the n = 2 case, we find τ_rc = 0.1, τ₀ = 2 × 10⁵, and p_rc = 0.07. Our values for p_rc agree with the observed boundary, which is located at 0.1–0.3 bar, and varies with latitude (Tellmann et al. 2009). Performing a correlation analysis on the data and the n = 1 and n = 2 models yields a square of the correlation coefficient, r², of 0.99 for both models, indicating a very good fit of the model to the observed temperature–pressure profile. Note that extensive overlap of absorption lines at high pressures may cause the pressure dependence of the optical depth to be less steep than the n = 2 case deep in the Venusian atmosphere.

Zoom In Zoom Out Reset image size

In Figure 9, we compare our model-generated thermal fluxes (from our n = 1 case) to those from a line-by-line radiative transfer model that has seen extensive application to Venus (Crisp 1986, 1989; Meadows & Crisp 1996), the Spectral Mapping Atmospheric Radiative Transfer (SMART) model. The agreement is quite good (considering we are using an analytic gray model), with differences tending to be at or below about 10%. Also shown is the curve representing σT(p)⁴ for our model-generated temperature–pressure profile, which demonstrates the expected result that our model thermal fluxes approach the local blackbody flux at high pressure where there are large optical depths.

Zoom In Zoom Out Reset image size

Jupiter provides a case in which we compare a hierarchy of approaches to analytically computing atmospheric temperature–pressure profiles. We apply three different versions of our analytic model to Jupiter: a radiative equilibrium model without solar attenuation, a radiative–convective model without solar attenuation, and a radiative–convective model with solar attenuation. Steps up this chain of models incorporate additional physics and, thus, require larger numbers of parameters.

For our purely radiative equilibrium model without solar attenuation, we use Equation (24) to compute the temperature profile, with a modification to include an internal energy source, taking A = 0.34, F^☉ = 50 W m⁻², and F_i = 5.4 W m⁻² (Hanel et al. 1981). In addition, we take p₀ = 1.1 bar as a reference level, n = 2, and τ₀ = 6. The 1.1 bar reference level was selected to coincide with a value reported in the observed profile, and the value of τ₀ was computed from the broadband Rosseland mean opacity tables from Freedman et al. (2008) by interpolating to Jupiter's metallicity. A Rosseland mean opacity was chosen (rather than a Planck mean) because it should apply in the optically thick, deep atmosphere where we set p₀.

Our radiative–convective model without solar attenuation builds on the previous model. In the convective regime, we specify the ratio of specific heats as γ = 7/5 and the average ratio of the lapse rate to dry adiabatic lapse rate as α = 0.85. We then add solar attenuation to this model, which develops a stratospheric inversion. Using outputs from the model atmospheres of Fortney et al. (2011), we set F₁ = 1.3 W m⁻² (taken as the solar flux absorbed above 0.1 bar) and F₂ = 7 W m⁻² (taken as the solar flux absorbed below 0.1 bar). About 2 W m⁻² are absorbed between 0.1 bar and 1.1 bar, so that k₂ = −ln ((7 − 2)/7)/τ₀ = 0.06, and the value of k₁ is found to be 100 by comparing the model temperature profile to the known temperature at the top of Jupiter's stratosphere (165 K; see Equation (18)).

Figure 10 shows the temperature–pressure profiles from these models as compared to a measured profile synthesized from a number of observations (Moses et al. 2005). The radiative equilibrium model with solar attenuation (model "a") fits through the upper troposphere, but is super-adiabatic at deeper levels. This is corrected in the radiative–convective model without solar attenuation (model "b"), but both aforementioned models cannot match the stratospheric temperature profile since they ignore solar absorption. The radiative–convective model with solar attenuation (model "c") more closely matches the observed temperature–pressure profile, even through the stratosphere. This model finds p_rc = 0.25 bar, while the radiative–convective model without solar attenuation finds p_rc = 0.22 bar, which can be compared to the ∼0.5 bar boundary found in non-gray numerical radiative–convective models (Appleby & Hogan 1984). Performing a correlation analysis on the data and the radiative–convective model with solar attenuation yields a correlation coefficient of r² = 0.92, indicating a good fit to the observed temperature–pressure profile, despite having only a few parameters in our model.

Zoom In Zoom Out Reset image size

In Figure 11 we show the flux profiles from our radiative–convective model with solar attenuation. The difference between the net solar flux and the net thermal flux at the top of the atmosphere is equal to the internal energy flux (5.4 W m⁻²), and, as is the case for our Titan models (see below), the net solar flux profile shows some structure that is related to its parameterization. An important characteristic of our model is that it can easily compute a convective flux, which is also shown in Figure 11.

Zoom In Zoom Out Reset image size

McKay et al. (1999) applied simple, analytic radiative equilibrium models to Titan. They were able to fit Titan's temperature–pressure curve using a model similar to our Equation (18), except that attenuation was ignored in one of their channels (i.e., extinction of sunlight was ignored in Titan's troposphere, setting k₂ = 0). These authors took n = 4/3, which is intermediate between the dependence in the troposphere (n = 2) and high stratosphere (n = 1). Other model parameters were determined by comparing the temperatures in the analytic model to the known temperature–pressure curve (Lellouch et al. 1989) at the surface, tropopause, and the top of the atmosphere, giving τ₀ = 3 and k₁ = 160.

For our model, we take p₀ = 1.5 bar and $T_0 = 94{\kern 1pt} {\rm K}$ (the known surface pressure and temperature), α = 0.77 (the average ratio of the lapse rate in Titan's troposphere to the dry adiabatic lapse rate), and γ = 7/5. We also take F^☉₁ = 1.5 W m⁻² and F^☉₂ = 1.1 W m⁻², which represent the solar flux absorbed in Titan's stratosphere and troposphere, respectively. The division of these fluxes between the stratosphere and troposphere follows McKay et al. (1991), while the sum of these fluxes matches the net solar flux absorbed by Titan, and is in agreement with Titan's Bond albedo of roughly 0.27 (Neff et al. 1985) and mean insolation of about 3.6 W m⁻². Since Titan has no appreciable internal heat flux, we set F_i = 0.

Of the 1.1 W m⁻² that is absorbed in the troposphere, only about 0.35 W m⁻² is absorbed at the surface (McKay et al. 1991). This gives an optical depth of roughly unity to sunlight in this channel, so that k₂ ≈ 1/τ₀, and k₁ is then determined by comparing the model temperature profile to the known temperature at the top of Titan's stratosphere (175 K, see Equation (18)). As was mentioned earlier, the requirement that the temperature and upwelling flux remain continuous across the radiation–convection boundary allows us to fit for τ_rc and τ₀.

Figure 12 shows a comparison between our radiative–convective model and the best-fit model from McKay et al. (1999). We show two radiative–convective models: one with n = 4/3 and the other in which, following the optical depth profile from McKay et al. (1999) and the parameterizations from Frierson et al. (2006) and Heng et al. (2012), the power of the τ–p relationship varies smoothly from ∼2 at the top of the convective region to ∼1 at the top of the atmosphere according to

$$\begin{equation} \tau \left(p \right) = \tau _{{\rm rc}} \left[ {f\frac{p}{{p_{{\rm rc}} }} + \left({1 - f} \right)\left({\frac{p}{{p_{{\rm rc}} }}} \right)^n } \right]. \end{equation} \tag{ 36 }$$

Here f controls the extent of the region of the atmosphere that is dominated by collision-induced absorption (n = 2), the top of which we place at roughly 0.2 bar.

For our n = 4/3 model, we find τ_rc = 4.8, τ₀ = 5.3, k₁ = 120, and p_rc = 1.4 bar. Our value of τ₀ is larger than that of McKay et al. (1999) because we consider attenuation of sunlight in Titan's troposphere, which then requires a larger greenhouse effect (i.e., more infrared-opaque atmosphere) to maintain the same surface temperature. This also causes our value of k₁ to be slightly smaller, since it is inversely proportional to τ₀.

Zoom In Zoom Out Reset image size

Our variable n model finds τ_rc = 4.0, τ₀ = 5.4, k₁ = 120, and p_rc = 1.3 bar. Note that our low values for p_rc are in agreement with observations (McKay et al. 1997; Fulchignoni et al. 2005) and results from models (Flasar 1983, 1998). We discuss the shallowness of Titan's convective region below in Section 5 in the context of general model behaviors presented in Section 3. Both our variable n model and the n = 4/3 model produce good fits to the observed temperature–pressure profile, with the variable n model (dashed line in Figure 12) performing quite well in both the lower and upper atmosphere (yielding a correlation coefficient, r², of 0.99 when compared to the data). Figure 13 shows the thermal fluxes from our variable n model, and compares the net thermal flux from our model to that from a validated Titan radiative transfer model (Tomasko et al. 2008). The agreement is reasonably good, with discrepancies typically less than 30%, and additional structure appearing in our model due to our parameterization of the net solar flux (Equation (15)).

Zoom In Zoom Out Reset image size