[13] We developed a multilayer nongray radiative‐convective transfer computer model to calculate the effect of the addition of artificial greenhouse gases on the surface temperature of Mars. For a given temperature profile, the model computes the upward and downward radiative fluxes throughout the atmosphere, taking into account both CO 2 and greenhouse gas concentrations. The temperature profile is then adjusted until radiative‐convective equilibrium is achieved and the surface satisfies the appropriate convective conditions. In this section we describe the model and the temperature adjustment method used in this computation. Because only temperatures less than ∼260 K are considered, H 2 O opacity was not included in the calculations.

4.1. Radiative‐Convective Calculation

[14] Our computational method uses the hemispheric mean two‐stream approximation in computing the upward and downward fluxes in the atmosphere [ Meador and Weaver, 1980 Toon et al., 1989 F−), is given by I− is the intensity in the downward direction, and μ = cos θ is the projection of the beam onto the z axis. Since we approximate the radiation as uniform in the downward direction, I− is not a function of μ and is the hemispherically averaged incidence angle ( = 0.5 since the radiation is assumed uniform in the downward direction), and B is the Planck function. The first term in McKay et al., 1989 B 0 and B 1 can be solved for by setting the top of the layer (τ = 0) and bottom of the layer (τ = τ i *) temperatures as the boundary conditions. τ i * is the total optical depth of layer i. To find the expression for the downward flux at the bottom of a layer (i.e., at a level) we integrate i * and get F top − is the downward flux at the top of the layer, and the top and bottom subscripts denote the top and bottom edges of the layer, respectively, where down is toward the planetary surface. Our computational method uses the hemispheric mean two‐stream approximation in computing the upward and downward fluxes in the atmosphere []. The atmosphere is divided into discrete layers. In this approach, the thermal radiation is assumed uniform in the upward and in the downward directions. After integrating over the azimuthal angle, the downward flux (), is given bywhereis the intensity in the downward direction, and μ = cos θ is the projection of the beam onto theaxis. Since we approximate the radiation as uniform in the downward direction,is not a function of μ and equation (3) is evaluated toWe use the integral form of the radiative transfer equation,where τ is the optical depth, τ′ is the variable of integration,is the hemispherically averaged incidence angle (= 0.5 since the radiation is assumed uniform in the downward direction), andis the Planck function. The first term in equation (5) represents the energy that is transmitted through the layer, while the second represents the energy that has its source in the layer (that is emitted by gas in the layer) and then is attenuated on the way down through the rest of the layer. Additionally, we approximated the Planck function as a linear function of τ [], so thatandcan be solved for by setting the top of the layer (τ = 0) and bottom of the layer (τ = τ*) temperatures as the boundary conditions. τ* is the total optical depth of layer. To find the expression for the downward flux at the bottom of a layer (i.e., at a level) we integrate equation (5) from 0 to τ* and getwhereis the downward flux at the top of the layer, and theandsubscripts denote the top and bottom edges of the layer, respectively, where down is toward the planetary surface.

[15] For stability reasons, as described later in more detail, the flux in the middle of the layer (τ = τ i */2) must also be calculated and used in converging the model. Integrating i */2, we get the downward flux in the middle of the layer: For stability reasons, as described later in more detail, the flux in the middle of the layer (τ = τ*/2) must also be calculated and used in converging the model. Integrating equation (5) from 0 to τ*/2, we get the downward flux in the middle of the layer:

[16] The same procedure was used in calculating the upward flux: I+ is the upward intensity, given by The same procedure was used in calculating the upward flux:whereis the upward intensity, given byThe Planck function is again approximated as a linear function of τ ( equation (6) ).

[17] The upward flux at the top of the layer (i.e., integrating from 0 to τ i *) is given by i */2) is given by integrating i */2 to τ i *. After simplifying this gives The upward flux at the top of the layer (i.e., integrating from 0 to τ*) is given byThe upward flux in the middle of the layer (τ = τ*/2) is given by integrating equation (10) from τ*/2 to τ*. After simplifying this gives

[18] Our radiative convective calculation is similar in approach to that of Pollack et al. [1987], McKay et al. [1989], and Kasting [1991]. The radiative fluxes were calculated over a grid extending from the surface to an altitude equivalent to 10−9 Pa; the number of layers was set to 30 for all results presented here, but could be varied. It was found that using more than 30 layers resulted in no more than a 0.5% change in the calculated surface temperature change. The layers were about equally spaced in pressure for the bottom 70% of the atmosphere. The albedo in the visible was set to 0.209 in order to give the equilibrium temperature of 212.9 K for the present Mars, as computed by Pollack et al. [1987]. The net incoming solar radiation corresponded to the yearly averaged solar flux and was assumed to be deposited at the ground (the atmosphere was assumed transparent in the visible). The outgoing thermal radiation was resolved into 500 spectral intervals from 2 to 2500 cm−1, each with a width of 5 cm−1 (except for the first interval); band boundaries were rounded to the nearest interval boundary. All greenhouse gases were assumed to be uniformly mixed in the atmosphere and thus the opacity in each wavelength interval and layer was set proportional to the absorption coefficient for that interval and the pressure difference across the layer.

[19] In order to compute the total upward and downward flux for each level, we first computed these fluxes individually for each spectral interval, and for each term in the exponential sum within the spectral interval ( ). The opacity for each exponential sum term was calculated using τ i = kN where τ i is the layer opacity, k is the absorption coefficient for the exponential sum term under consideration [m2 molecule−1], and N is the column density of the absorbing gas [molecules m−2]. The total flux for the spectral interval is then given by . Summing the flux over all spectral intervals gave the total upward and downward flux for that level.

[20] The computation converged on a steady state solution by requiring that the net upward infrared flux at each layer (the total upward flux minus the total downward flux) balance the net incoming solar flux. The temperatures at the levels were changed proportionately to the residual flux (net up minus net incoming), after the method described by Pollack and Ohring [1973], so as to reduce the residual flux to zero. For stability reasons, the flux balance was determined at the mid points of each layer while the temperature was specified at the levels (the layer boundaries) [e.g., McKay et al., 1989]. The computation was iterated upon until there was less than ∼0.5% error in the residual flux, with respect to the incoming flux, above the convective zone. This was possible for all partial greenhouse gas pressures of 0.1 Pa or below. For partial pressures of 1 Pa and 10 Pa, the error was as high as 5%.

[21] A convective adjustment was performed if the temperature lapse rate in the layers near the surface exceeded the adiabatic lapse rate. Because the temperatures considered were always less than ∼260 K, and thus the saturation vapor pressure of water was negligible, the dry adiabatic lapse rate was used in the calculations. If the adiabatic lapse rate was exceeded in the lowermost layer, then the temperature lapse rate for that level was reset to the adiabatic lapse rate and the radiative balance for the rest of the atmosphere was recalculated using the new boundary condition. The procedure was repeated, each time increasing the thickness of the convective zone by one layer, until a solution was found where the adiabatic lapse rate was not exceeded in any layers and all layers above the convective zone were in radiative equilibrium.

[22] Our numerical scheme was tested by running cases with gray opacities and comparing the level fluxes and resulting temperature profile with the known exact solution. Figure 5 shows comparisons of the multilayer model and the theoretical two‐stream approximation results for τ = 0.2, 0.5, 1.0, and 1.5; τ = 0.2 is equivalent to a little less than 0.1 Pa of C 3 F 8 , and τ = 1.5 is equivalent to a little more than 1 Pa of C 3 F 8 . Note that there is no convective layer in this calculation. The results match theory very well, especially for higher optical depths. For lower optical depths, the deviation is due to the ground boundary condition: the theoretical result gives the surface discontinuity, while the model smooths the bottom layers as it tries to balance the incoming and outgoing flux at each layer.

Figure 5 Open in figure viewer PowerPoint Comparing the multilayer radiative‐convective model results (points) with the theoretical results (lines) for a gray atmosphere. The model gives very good results for higher optical depths. For lower optical depths the deviation is the result of the surface discontinuity since the model attempts to balance the fluxes at all levels and create a “smooth” transition.

[23] The results presented here improve on our previous studies of greenhouse warming on Mars using a simplified one‐layer model [Marinova et al., 2000]. The results reported in Table 2 of Marinova et al. [2000], for P C 2 F 6 = 0.01 to 1 Pa, are on average within 27% of the multilayer results listed here. However, this agreement is partly fortuitous, due to an error of a factor of π in the 1‐layer flux computer calculations. Note also the typographical error in equation (3) of Marinova et al. [2000]: the first plus sign should be a minus sign.