A Greenhouse Model You Can Check by Hand, and Where It Stops Being Honest

Thirty-Three Degrees, and a Model That Cannot Reach Them

The layer model of the greenhouse effect is the best piece of physics teaching in the whole of climate science. It needs no calculus and no computer. A single line of algebra turns sunlight into a surface temperature, and the line is short enough to write on the back of a receipt. Generations of students have met the greenhouse effect this way [10][12][17].

We set it a task. Earth's surface averages about 288 K. Strip the atmosphere away in the infrared and the same sunlight supports only 255 K. How many layers of the model does it take to cover the difference?

Nought point six two two.

Zero layers is 32.79 K too cold. One layer is 15.50 K too hot. There is no third option, since layers come in whole numbers, and the number the planet wants is between them and closer to one than to zero. That gap is not a rounding problem. Fifteen and a half kelvin is nearly half of everything the model was built to account for, and it is roughly five times the warming the world has argued about for thirty years.

This article does three things. It builds the model carefully and checks it against its own exact solution. It repairs the overshoot by letting the layer be an imperfect absorber, which is the standard fix and works. Then it spends its second half finding out precisely which physical facts the repaired model is still getting wrong, and by how much, in kelvin and in watts per square metre. We like this model. We would teach it tomorrow. That is exactly why the second half matters.

Building It From One Sentence

Start with a sphere in sunlight. The Sun delivers \(S = 1361\) W/m² to a disc of the planet's cross-section [12]. The planet reflects a fraction \(\alpha = 0.293\) of that straight back [19], and it spreads what it keeps over a sphere four times the area of the disc. So the absorbed flux, averaged over the whole planet, is

In equilibrium the planet must radiate that back. If it radiates as a blackbody at a single temperature \(T_e\), then \(\sigma T_e^4 = F\) and \(T_e = 255.212\) K. That is the whole zero-layer model.

Now add the atmosphere, in the crudest form that captures anything. Take \(N\) slabs stacked above the surface. Each slab is perfectly transparent to sunlight, so the ground still absorbs all of \(F\). Each slab is a perfect absorber and therefore a perfect emitter in the infrared, which is Kirchhoff's law doing the work. The idea of treating an atmosphere as a stack of grey emitting layers in radiative balance goes back to Schwarzschild's 1906 work on the solar atmosphere [2], which is where the two-stream picture underneath all of this began. Each slab has one temperature and radiates \(\sigma T^4\) upward and the same amount downward. Space, above the top slab, sends nothing back.

Write \(x_j = \sigma T_j^4\) for the surface (\(j=0\)) and each slab (\(j=1\) to \(N\)). Then the energy balance at every level is a linear equation in the \(x_j\):

That is a tridiagonal system, and the club's code assembles it and hands it to numpy.linalg.solve for any \(N\) you like. The reason for solving numerically something with a known closed form is the next section.

The closed form, for the record, is that \(x_j = (N+1-j)F\). Every level differs from the one below it by exactly \(F\) in units of flux. The surface sits \(N\) steps above the top, so

Figure 2 draws the column that comes out of that. Look at what the five stacks have in common. Whatever \(N\) is, the top layer always sits at \(T_e\), because the top layer is the thing space actually sees and space is owed exactly \(F\). Adding layers does not change the top of the atmosphere. It changes how far below the top the ground has been pushed.

That is the physical content of the model, and it is genuinely correct physics: the greenhouse effect works by raising the altitude from which the planet radiates to space, and then letting the lapse rate do the rest [11][17]. The slabs are a caricature of a lapse rate. A bad one, as we will see. But the mechanism they caricature is the real one.

The Arithmetic

Everything above reduces to numbers. Here they are, done slowly.

Absorbed sunlight. \(1361 \times (1 - 0.293) = 962.227\), divided by 4 is \(240.557\) W/m².

Effective temperature. \(240.557 / 5.670374419 \times 10^{-8} = 4.2424 \times 10^{9}\). The fourth root of that is \(255.212\) K. In Celsius, −17.94.

One layer. \(2^{1/4} = 1.189207\). Times 255.212 is 303.500 K, or 30.35 °C. The real planet is 15.00 °C.

The fraction. We want \(T_e (N+1)^{1/4} = 288\), so \(N+1 = (288/255.212)^4\). The ratio is 1.128473. Squared, 1.273451. Squared again, 1.621676. Subtract one: \(N^\ast = 0.6217\).

The greenhouse effect in flux rather than temperature. \(\sigma \times 288^4 = 390.11\) W/m² leaving the ground, against 240.557 W/m² leaving the top, a difference of 149.55 W/m². Observation puts the same quantity at \(398.2 - 239.0 = 159.2\) W/m² [16]. We are 9.65 W/m² low, or 6.1%. That discrepancy is real and we come back to it.

Venus, for scale. Its fact sheet gives \(S = 2601.3\) W/m² and \(\alpha = 0.770\), so \(T_e = 226.627\) K, and its mean temperature is 737 K [23]. Then \(N^\ast = (737/226.627)^4 - 1 = 110.85\). Venus needs a hundred and eleven layers. Mars, on the same arithmetic, needs 0.082.

Checking It Against Itself

A model with an exact solution is a gift, because it lets you find out whether your code is computing the thing you think it is computing before you ask it anything you do not already know. We ran three checks and printed all three.

Check one, against the textbook. The zero-layer case should give the famous 255 K. Ours gives 255.212 K with the CERES albedo of 0.293 [19], and 254.578 K with the round 0.30 that most textbooks use. Both round to 255. Notice that the spread between the two albedo choices, 0.634 K, is bigger than either one's distance from the quoted value, which is a small warning about how much precision the famous number can carry.

Check two, against the analytic solution. For every \(N\) from 0 to 10 the code prints the numerically solved surface temperature beside \(T_e(N+1)^{1/4}\). Nine of the eleven rows agree bit for bit. The worst row, \(N = 10\), differs by \(5.684 \times 10^{-14}\) K out of 464.78, a relative error of \(1.223 \times 10^{-16}\). Double precision has an epsilon of \(2.220 \times 10^{-16}\), so the largest disagreement anywhere in the table is 0.55 of one representable step. There is no physics in that number. It is the last bit of a floating-point division.

We pushed it further than the article needs, to \(N = 500\), where the condition number of the matrix has climbed to \(4.1 \times 10^{5}\) and the relative error with it, to \(2.4 \times 10^{-15}\). Still eleven orders of magnitude below anything physical. We print it because that is where "exact" turns into "exact enough", and a reader deserves to be shown the boundary rather than told it is far away.

Check three, against a closed-form series. The two-band work later needs the fraction of a blackbody's emission that falls between two wavelengths, which we get by Simpson's rule. The blackbody fraction function also has an exact series solution, so we computed both. Across the five bands the model actually uses, the largest disagreement is \(8.3 \times 10^{-12}\). Then, because a 401-point grid that is generous for a 4 µm interval is not generous for a 96 µm one, we printed the refinement table for a deliberately wide band and watched the error fall by a factor of about 256 per fourfold refinement, which is the fourth-order convergence Simpson's rule promises. Nothing in the study integrates a range wider than 5.5 µm.

Working Notes From the Club Table

Half a Layer Short

Figure 1 is the answer to the question in the title, and it is a slightly awkward picture. The measured planet sits between the first two points on the curve.

Read the table underneath it and the awkwardness gets worse rather than better. The fourth-root law means the steps shrink fast. Going from zero layers to one buys 48.29 K. Going from one to two buys 32.38 K. From nine to ten, 10.94 K. So the model's resolution is coarsest exactly where Earth lives, and the single interval a reader most wants to subdivide is the widest one in the whole table.

The last column is the punchline of the next section, written early. It asks what single-layer emissivity would give the same surface temperature as each \(N\). For \(N = 1\) the answer is 1, as it must be. For every \(N\) beyond that the answer is greater than 1, which is physically impossible for one layer, and the numbers march sedately towards 2 without ever getting there. A stack of layers is not equivalent to a thicker single layer. It is a different object.

Letting the Layer Be Imperfect

The repair is old and standard [5][14][17]. Keep one layer. Stop insisting it absorbs everything. Give it a longwave emissivity \(\varepsilon\) between 0 and 1, which by Kirchhoff is also its absorptivity, and redo the two balances. The layer's own balance is

and the \(\varepsilon\) cancels. That deserves a pause. The temperature of the layer does not depend on how absorbing it is. A gauzy veil and an opaque slab sit at exactly the same temperature, \(T_a = T_s/2^{1/4}\). They differ only in how much of the ground they can see, and therefore in how much they send back down. Put that into the surface balance and

Set \(T_s = 288\) and invert. We did it twice, once by algebra and once by bisection on the forward function, and the two agreed to \(6.7\times10^{-16}\):

The number is satisfying. It is comfortably below 1, so the fix is not straining. It sits where the curve is steep enough to be informative: at \(\varepsilon^\ast\) the surface moves 58.4 K per unit of emissivity, so one percentage point of emissivity is worth 0.58 K. And the layer temperature it implies, 242.18 K, corresponds to an emission altitude of about 7.0 km at a standard lapse rate, which is the right order for the real effective emission level of around 5 km and about 41% too high.

How well is \(\varepsilon^\ast\) determined? We ran 200,000 Monte Carlo trials, drawing the solar constant from Normal(1361.0, 0.5) [12], the albedo from Normal(0.293, 0.005) [19] and the observed surface temperature from Normal(288.0, 0.5). Figure 5 shows the running mean settling.

The answer is \(\varepsilon^\ast = 0.7667\) with a standard deviation of 0.0122 and a standard error of \(2.7 \times 10^{-5}\), a 95% interval of 0.743 to 0.790. The albedo and the surface temperature contribute almost equally, moving \(\varepsilon^\ast\) by 0.0087 and 0.0085 per standard deviation, while the solar constant is worth only 0.0005 and might as well be exact. The Monte Carlo mean sits \(4.2\times10^{-5}\) below the point estimate, which is 1.53 standard errors, and is nonlinearity rather than error: the mean of a nonlinear function of the inputs is not the function of the mean of the inputs. We print it instead of rounding it away.

Where It Starts Lying: Two Constraints That Do Not Cross

Now the second half. Everything so far assumed the atmosphere is grey, meaning it absorbs equally at every infrared wavelength. It does not. Water vapour and carbon dioxide absorb in bands, and between roughly 8 and 12 µm there is a region where the clear atmosphere is comparatively transparent. This is the atmospheric window, and it is not a minor feature.

Integrating the Planck function, a 288 K surface puts 25.28% of its emission into that 8 to 12 µm interval. The Wien peak at 288 K is at 10.06 µm. The window sits directly on top of the emission maximum, not out in a tail where it could be ignored.

So we built the two-band version. Same single layer, same geometry, but the spectrum is split. The window has optical depth \(\tau_w\) and the rest has \(\tau_b\), with emissivity \(1 - e^{-\tau}\) in each, and the band fractions come from Planck integration at whatever temperature the surface has reached, so the whole thing is solved by fixed-point iteration.

The first result is a disappointment dressed as a success. Collapse the algebra and the two-band model presents the surface with an effective grey absorptivity \(A = \varepsilon_w f + \varepsilon_b (1 - f)\), where \(f\) is the window fraction, and then the surface temperature obeys exactly the same formula as before with \(A\) in place of \(\varepsilon\). Calibrate \(\tau_w = 0.30\) and solve for \(\tau_b = 2.7878\), and the two models agree on \(A = \varepsilon^\ast = 0.766708\) to the last printed digit. For the present-day state, splitting the spectrum bought us nothing at all. That is worth knowing, and it is the reason the grey model survives in textbooks.

The second result is not a disappointment.

Two constraints that do not cross

Ask the calibrated two-band model how much surface emission reaches space without ever being absorbed. In the window it is \(f e^{-\tau_w}\sigma T_s^4 = 73.07\) W/m², plus another 17.94 through the leaky absorbing band, for 91.01 W/m² in total.

The measured value is 22 W/m² [13].

Costa and Shine did this calculation properly, with line-by-line radiative transfer and real atmospheric profiles, and found that about 22 W/m² of the outgoing longwave radiation is surface emission that made it to space untouched, roughly 9% of the total. Our model says 91. It is over by a factor of 4.1.

So close the window. We solved for the window optical depth that reproduces 22 W/m² exactly, with the absorbing band made fully opaque, and got \(\tau_w = 1.704\). The escape is then 22.000 W/m² by construction. The surface temperature is 299.97 K.

Twelve kelvin too hot.

There is a sharper way to put the same thing, and it took us a meeting to see it. Rerun the calibration at several window opacities, re-solving \( au_b\) each time so the surface still lands on 288 K. The window escape falls from 73.07 W/m² at \( au_w = 0.30\) to 22.01 at \( au_w = 1.50\), a factor of 3.3. The total escape reads 91.008 W/m² in every single row. It cannot do anything else: once the model is pinned to 288 K its effective absorptivity is fixed at \(arepsilon^st\), so the untouched fraction leaving the top is \((1-arepsilon^st)\sigma T_s^4 = 91.008\) W/m² however you divide the opacity between the bands. Close the window and the absorbing band is forced open by exactly the compensating amount. You can put the right number in the window column and the model is still leaking four times too much radiation in total.

A single-layer model can match the surface temperature or it can match the window escape. It cannot do both, and the gap between the two is not small in either currency. The reason is structural and no amount of parameter fitting will touch it. One layer has one temperature. In the real atmosphere the window radiation that does get absorbed is absorbed low down, by the water vapour continuum where the air is warm, and re-emitted at nearly the surface temperature, so it costs the planet almost nothing. The band-centre radiation is absorbed and re-emitted high up where the air is cold, which is where the greenhouse effect actually comes from [11][21]. A single slab has to pick one altitude for both. Whichever it picks, the other quantity is wrong.

Stacking more layers fixes this, because a stack has a lapse rate and can put different wavelengths at different heights. Making the layers grey does not fix it, no matter how many you use. The failure here is spectral, not vertical, and it is the one the classroom version never mentions.

What Grey Actually Costs

The present-day state is not where the grey assumption fails, because we forced both models through the same point. It fails in the derivative. Multiply the optical depth by a factor and ask what the surface does.

Double the absorber. The grey model warms by 11.51 K. The two-band model warms by 5.44 K. The grey model overstates the response by a factor of 2.11 on the central settings, and by \(2.02 \pm 0.01\) across 3,677 usable Monte Carlo trials that sampled the window's low edge from 7.5 to 8.5 µm, its high edge from 11.5 to 13.0 µm, and \(\tau_w\) from 0.15 to 0.60. The 2.5th percentile of that distribution is 1.27, so even the most forgiving window we sampled leaves the grey model overstating by a quarter. Another 323 trials were dropped because a wide, clear window made 288 K unreachable, and those were precisely the cases where grey would have looked worst, so 2.02 is conservative.

The reason is visible in the two calibrations. The grey atmosphere at \(\tau = 1.4555\) is only 76.7% opaque, so it can always absorb more. The two-band model's absorbing band is already 93.84% opaque, so doubling it accomplishes very little and the warming has to come out of the window instead.

Push that harder, in the direction carbon dioxide actually pushes. CO₂ works in the band and does almost nothing to the 8 to 12 µm window, so scale \(\tau_b\) alone:

There is a ceiling, and one doubling gets you 94% of the way to it. A hundredfold increase in band opacity is worth 2.64 K and then the band is finished, because a saturated band radiates from its own top and cannot be made to radiate from any higher. The grey model has no ceiling at all. It cannot have one, because a single band with a single optical depth can always be made more opaque. This is the single most important thing the grey assumption hides, and it is part of why real CO₂ forcing goes as the logarithm of concentration rather than linearly [1][9][21].

We should say plainly what this two-band number is not. It is not a climate sensitivity. There are no feedbacks here, no water vapour response [7], no clouds, no lapse-rate adjustment, and the real answer to a CO₂ doubling is around 3 K for reasons this model contains almost none of [22]. Our 5.44 K and our 2.64 K are statements about two toy models compared against each other under identical rules. The comparison is the result. Neither number is a prediction about the planet.

The Strongest Objection We Could Make

Here is the best case against this entire article, and we think it is a good one.

You have not tested the grey assumption. You have tested a single layer. Every failure in §9 and §10 follows from having exactly one temperature available, and you said so yourself. A grey model with fifty layers has a lapse rate, puts its emission at a range of altitudes, and would not show the window pathology at all. So the headline, that grey overstates the doubling response by a factor of two, is really the claim that one-layer grey overstates it relative to one-layer two-band, which is a comparison between two caricatures rather than a measurement of anything.

That objection is largely correct, and we want to grant its correct part fully before arguing with the rest.

It is right that the double bind in §9 is a single-layer failure. A multi-layer grey model would put the window escape and the band emission at different effective altitudes and could match both. We could have built that model. We did not, and a reader is entitled to treat §9 as a statement about single-layer models only, which is how we have labelled it.

It is wrong about the saturation result, and that is the part we would defend. The ceiling in §10 does not come from having one layer. It comes from having a finite band. Add a hundred grey layers and increasing \(\tau\) still buys warming without limit, because in a grey atmosphere every wavelength gets more opaque together and the emission level rises everywhere. Add a hundred two-band layers and the absorbing band's emission level eventually rises out of the troposphere into the stratosphere, where the temperature stops falling with height, and the forcing per doubling flattens for exactly the reason our toy version flattens. Manabe and Wetherald found this with a real radiative-convective column in 1967 [4], and the modern analytic treatment of CO₂ forcing turns on the same mechanism [21]. Our ceiling is a cartoon of a real ceiling. Grey has no ceiling to cartoon.

There is also a second objection we would raise against ourselves, which nobody has yet raised against us: the whole study is in pure radiative equilibrium, and pure radiative equilibrium is not a description of the troposphere. Manabe and Strickler showed in 1964 that a radiative-only column with realistic absorbers produces a lapse rate steeper than the atmosphere can mechanically sustain, and that you must add convective adjustment by hand before anything Earth-like appears [3][5]. Roughly 100 W/m² leaves the real surface as latent and sensible heat rather than radiation [6][8][15]. None of that is in our file. The layer model happens to land on a defensible surface temperature partly because two large omissions point in opposite directions, which is luck dressed as agreement.

Why We Would Still Teach It

A model that overshoots by 15 K with one layer, cannot represent a fractional layer, leaks four times too much radiation through its window, has no convection, no clouds, no water vapour feedback and no spectrum, and whose central parameter depends more on your choice of albedo than on anything anybody measured, sounds like a model to throw away.

We would keep it, and not out of sentiment. It gets the mechanism right. Adding absorber raises the altitude from which the planet radiates to space; the ground, sitting further below that altitude, must be warmer. That sentence is true of the real atmosphere and it is the sentence the layer model teaches, in arithmetic a student can do in their head. Every correct account of the greenhouse effect is a refinement of it [10][11][20]. No amount of spectroscopy replaces it.

What we would change is the framing. The model is usually presented as an explanation that happens to be approximate. It is better understood as an argument that happens to produce numbers, and the numbers are the least reliable thing about it. A student who leaves with "one layer gives 33 degrees" has learned something false. A student who leaves with "raising the emission level warms the ground, and here is the smallest calculation that shows why, and here are four things it is lying about" has learned the physics and learned what a model is, which is two lessons for the price of one afternoon.

We would also teach the checking. The reason this study can say anything at all is that the \(N\)-layer system has an exact solution to be checked against, and we checked, and the agreement was 0.55 of one floating-point epsilon. That agreement proves nothing whatsoever about the atmosphere. It proves the code computes the model. Only then is it worth asking whether the model computes the planet, and the honest answer to that second question, for this model, is partly yes and partly no, in amounts we have now measured.

Where a different choice would have changed the answer

Several decisions in this study could have gone otherwise. Saying which ones move the numbers seems the least we owe a reader.

The albedo. This is the largest single lever and it is not a measurement uncertainty, it is a choice of source. With CERES's 0.293 we get \(\varepsilon^\ast = 0.7667\). With the round textbook 0.30 we get 0.7789. With Stephens and colleagues' 0.29 [15] we get 0.7615. Across the set of defensible choices in the code, including two different baseline periods for the surface temperature, \(\varepsilon^\ast\) runs from 0.7494 to 0.7918, a range of 0.0423. The Monte Carlo standard deviation is 0.0122. So the modelling choices are worth about three and a half times the measurement uncertainty, and anyone quoting \(\varepsilon^\ast\) to three decimals without naming their albedo is quoting noise.

What "no atmosphere" means. The famous 33 K assumes you can delete the greenhouse effect and leave the albedo alone. You cannot: most of Earth's albedo is cloud, and clouds go with the air. Redo it with a surface-only albedo of 0.15 and the airless Earth sits at 267.24 K, making the greenhouse effect 20.76 K rather than 32.79 K. That is a 37% change in the headline number, produced by a modelling decision and not by any new measurement.

Treating the planet as one temperature. Our model radiates \(\sigma \bar{T}^4\). The planet radiates \(\overline{\sigma T^4}\), and Jensen's inequality says the second is larger. At 288 K mean with a spatial and seasonal spread of \(\sigma_T\), the excess is \(\sigma(6\bar{T}^2\sigma_T^2 + 3\sigma_T^4)\). The measured 398.2 W/m² of surface emission [16] is recovered at \(\sigma_T = 16.92\) K, which is an entirely ordinary value. So the 6.1% shortfall we found back in §3 is mostly this, not a broken model. The planet radiates like a 289.48 K surface while averaging 288.00 K.

The Moon makes the same point with the sign reversed and much more violence. Its effective temperature, computed from its own fact-sheet albedo of 0.11, is 270.33 K, and NASA's sheet prints 270.4, so we are using their inputs the way they do [23]. But a body with no atmosphere and no thermal inertia would average 152.92 K over its surface, 117.4 K below the effective temperature that supposedly describes it. No greenhouse gas is involved anywhere in that gap. It is entirely Jensen's inequality, acting on a surface whose equatorial temperature runs from 95 K to 390 K across a lunation [18][23]. The real Moon has thermal inertia and sits between the two figures. Which means the 33 K we all quote for Earth is the gap to a fictitious uniform-temperature airless Earth, not to any airless planet that could exist.

Doubling both bands together. In §10 we scaled \(\tau_w\) and \(\tau_b\) by the same factor, which is the neutral choice and not the realistic one. Scaling the band alone, which is closer to what CO₂ does, gives 2.47 K instead of 5.44 K and makes grey look worse by another factor of two. Scaling the window alone would make grey look better. We reported the neutral choice as the headline and both extremes in the output, because picking the one that flatters your conclusion is the easiest dishonesty available in work like this.

The window's optical depth. The softest number in the study. At \(\tau_w = 0.10\) the grey model overstates the doubling response by 7.75; at \(\tau_w = 0.90\) it overstates it by 1.08. The direction of the conclusion holds everywhere in that range. Its magnitude does not, and a reader should carry the factor of two as an estimate with the window's uncertainty attached rather than as a measured constant.

Reproducing this

Everything in this article comes from one file and one command. You need Python 3.12 and numpy, and nothing else is imported.

python layered-greenhouse.py > layered-greenhouse-output.txt

Expect about nine seconds on a modern laptop. Ours took 8.82 s, on numpy 2.4.2 and Python 3.12.3. Memory never exceeds a few tens of megabytes; the largest array in the file is the 200,000-trial Monte Carlo and the 2,000,001-point grid used for the lunar integral in §11. The master seed is 20260321, hard-coded at the top, and the two Monte Carlo streams are spawned from it through numpy's SeedSequence.spawn, so the output is deterministic to the last digit. The three validation blocks print PASSES or an explicit failure banner, so a broken run announces itself rather than producing plausible nonsense. Numbers differing from ours by more than the printed standard errors mean something is wrong and we would like to hear about it. The interactive model runs both models in your browser, with the same constants, and its default settings reproduce the headline numbers above.

References