I have been working on neutron Monte Carlo simulations for quite some time, in order to simulate various aspects of nuclear systems. Recently, I realized I can use the exact same approach to model radiative transfer of heat as well.
The overall approach should be surprisingly simple! We will start with a simple, monochromatic radiation spectrum. In reality, the Sun is an approximate black body, the emission spectrum of which can be calculated from the temperature of the solar corona, which is 5778 K.
To do this, we can use Planck’s law, which gives the spectral radiance, defined as the energy per area per frequency interval per steradian. The latter is the SI unit for solid angle, which is the field of view covered.
To simplify things, we will stick to a simpler equation, the Stefan-Boltzmann equation. We can get this from Planck’s law by integrating over the frequencies, the area and the solid angle. The result is the flux in J/m²/s emitted by an object at a particular temperature.
As an aside, we can directly use this to calculate the solar flux on Earth. Taking the surface area of the Sun and multiplying it with the the flux, then dividing by the area of a sphere with the radius of the Earth’s distance from the Sun, we end up with 1369.22 J/m²/s on the Earth’s surface at the equator (neglecting atmospheric effects) — very close to measured values above the atmosphere (1361 J/m²/s), but quite a distance from 1120 J/m²/s measured below the atmosphere.
Clearly, the atmosphere absorbs quite a lot of incoming energy!
Now, we know how the flux comes in at the equator. What about the poles? Imagine taking the dot product between a vector normal to the Earth’s surface with the incoming flux vector, parallel to the equator.
If θ is zero, and the ray is coming in parallel, the area receives the full amount of flux. But what happens if θ is 90°? In that case, the cosine is zero, and no flux is received. Clearly, the polar regions receive less energy.
We can apply this equation to every section of the Earth, but then we are forgetting one essential aspect. We just noticed the serious discrepancy between the flux outside the atmosphere, and the one inside. If we simply take the dot product at the surface, we are ignoring the 100 km of gases between the surface and space (the Kármán line).
The solution is a Monte Carlo simulation. Let’s assume we have a sphere with a frontal surface area of 1 m², with 1000 J/m²/s flying towards it. We divide that 1000 W/m² into 1000 packets, with 1 W/m² each. We then let those packets move straight forward, λ meters at a time.
In all of this, we are assuming the rays are moving straight forward, parallel to the equator. Now, that is not entirely correct, but given the Sun’s radius relative to the distance to the Earth, we can take it as a point source. The same goes for the Earth’s radius. The angles between the two are simply so small that they become irrelevant.
Every step, we check if the packet has intersected with something. If it does, we apply the packet’s energy to that volume. If it hits the sphere’s surface, that heats up a little bit. This approach makes dealing with gases far easier: we can check the distance the packet moves through a particular volume, then calculate the odds of hitting something.
The mathematics for these odds are surprisingly simple. In nuclear physics, we have the macroscopic cross-section Σ. The decrease in intensity after passing a distance dx is given by:
This does pose a problem: if we have a discrete particle, we can’t lose intensity. We can’t end up with a half or quarter or one-hundredth of a neutron. The solution is to make it stochastic.
The probability of the neutron not being absorbed over a distance x is:
As a sanity check: what happens in the limit of x → 0? The exponent becomes zero, and the right-hand side is one: the neutron is not absorbed. The odds of being absorbed is then the complementary probability of this one:
And again, if we reduce x to zero, the probability now also goes to zero: for an infinitely thin sample, it won’t be absorbed, whereas for an infinitely large sample, x → ∞, it will be absorbed.
When moving a distance Δx in a simulation, we just have to calculate this probability p(x) and use an appropriate random number generator to see if the event takes place or not. For a large number of samples, by the law of large numbers, it will converge on the correct values.
After this brief detour into nuclear physics, we can return to light. The exact same mathematics applies here too. For those familiar with chemistry, this is termed the Beert-Lambert law for the attenuation of light.
So — we can now pass packets of light through a gas and calculate the probability of them being absorbed. What happens if they are absorbed? The energy is then released and added to the energy in the gas. This then raises the temperature.
But wait — we just used Stefan-Boltzmann to find that hot objects emit radiation. At that point, not only the Sun will be emitting radiation, but the gases themselves too! And if the light isn’t absorbed by the gas and hits the surface of the Earth, that will emit light too.
This is called Schwarzschild’s equation for radiative transfer. It really is quite intuitive.
To resolve this, we need to split the atmosphere and the Earth’s surface into small volumes, let light pass through them, heat them up, then calculate where the radiation each of them emits goes off to. Since they’re both spheres, the easiest solution for this is a spherical coordinate system, with positions indicated by r, θ and φ.
We let the energy packets themselves move in Cartesian coordinates, then convert those to the spherical coordinates when figuring out which volume they are in.
Over time, more and more of the sphere will heat up. Some radiation will be lost from the outer packets as it travels into space. The regions around the equator will heat up more, for reasons discussed previously. The polar regions will remain cold, but radiation parallel to the surface will heat them up somewhat.
The back of the sphere will not receive any direct radiation, relying entirely on radiative transfer from the front. In reality, the Earth of course rotates. That is the advantage of the spherical coordinate system: we can simply apply a Δθ rotation to the entire thing, shifting all the volume’s coordinates, so that now a new area is facing the Sun.
All of this is rather trivial. What I am curious about is how we could display any of this easily. I suppose, if we zoom in far enough, we can directly project the spherical results on a square: we can construct a manifold of the sphere. A different option would be actually render the results by putting a camera looking straight at the Earth.
Implementing a ray tracing algorithm for this is relatively simple. The intersection between a line with origin o and direction u, and a sphere with center c and radius r is given by the following equation:
The derivation for this is pretty easy, you can represent it with the quadratic equation, then make use of the fact that u is a unit vector to simplify things.
All in all, I’m curious to see how far I can take this — I think I can reuse the majority of my existing Monte Carlo code and directly couple it with an energy balance and one of my existing differential equation solvers.
Once the basic system is working, there are a lot of possibilities: adding data from HITRAN for energy-dependent absorption, simple atmospheric mixing, density variations (I got the 1976 US Standard Atmosphere working in C++ already for a different project!) and so on.