This week I encountered an interesting simulation stability bug in my climate. All the air in my simulation was steadily traveling from low to high altitude.
I haven’t posted on my climate simulation section yet so here’s a little background on how it works. I hope it’s enough to help this post make sense.
My climate simulation uses three layers: ground, low-altitude air, and high-altitude air. The vertical heat/air exchange (ground to low altitude and low altitude to high altitude) is very simplistic, and the horizontal air exchange within each layer is done through winds generated through Navier-Stokes equations. The simulation is run on a grid of points, and at every point the various properties (temperature, pressure, density, etc.) are calculated. Heat enters the simulation through insolation from the sun to the ground, and heat leaves the simulation through black-body radiation from the ground to space.
On to the bug.
The temperatures in all three layers — ground, low-altitude air, and high-altitude air — were remarkably stable (hurrah!). You can clearly see in the following images that the temperature at a given coordinate is a function of latitude and land vs ocean tile (land is warmer). It took a few hundred iterations for this profile to develop, but once it did, it remained stable. One red flag I should have noticed is that a few hundred iterations is a suspiciously short period of time, since each iteration currently represents only about 20 minutes of simulation time.
The bug became apparent when I started examining the pressure maps, which were not stable. The air at low altitude was reducing in pressure every iteration, and the air at high altitude was increasing in pressure every iteration. This was consistent across the map. I expected the low altitude air to decrease in pressure near the equator and increase in pressure near the poles. This pressure difference is driven by differential heating and, in turn, drives wind patterns.
The pressure map looks good initially, but the low pressure (blue) areas gradually take over the entire layer. I did a second check by summing up the total air mass in each layer at every simulation step. The high-altitude layer was gaining mass steadily and the low-altitude layer was losing mass.
Another red flag I should have noticed was the very slow wind speeds. Most of my wind speeds were below 0.1 m/s. I assumed this was due to magnitude errors somewhere in my simulation; I wasn’t too concerned because the winds were all pointing in reasonable directions. I assumed it would be easy to adjust their magnitudes later.
I hoped that the pressure changes were simply the simulation approaching steady-state. To check, I plotted the pressure changes in one region to see if the graph was asymptotic.
It was immediately clear that there was no asymptotic behavior. That’s about as straight as lines get. The density trends were similar. My low altitude layer was continuously losing a set amount of air to my high altitude layer every iteration.
After investigation, I suspected one possible cause of the instability was that the vertical air transfer between low-altitude and high-altitude air was determined by the temperature difference between the ground and the low-altitude air. Basically, if the ground heats up the air, it causes a local updraft, and if the air is heating the ground, it creates a local downdraft.
Described like this, it seems that the necessary negative feedback loop to keep this section of the simulation stable is missing: the updraft between low and high-altitude air is determined by the relationship between the ground and low-altitude air. Without another connection between the ground and high-altitude air, this is a one-way relationship that can cause unstable behavior.
How could this be fixed?
The obvious brute-force approach would be to scale the updrafts and downdrafts every iteration so that the sum of all vertical air movement is zero. This crudely ensures my air layers are conservative with respect to vertical air movement. It would be effective and wouldn’t be too expensive (one extra iteration through the array to calculate the initial total vertical air movement before scaling), but it’s brute-force and this is a simulation. There have to be more elegant solutions.
I instead tried restructuring the vertical heat/air exchange between the low-altitude and high-altitude layers to exchange only heat. That would guarantee the vertical exchange is mass conservative. I also made the exchange dependent on the temperature difference between the low-altitude and high-altitude air layers (previously this exchange was dependent on the temperature difference between the ground and low-altitude air layers). This introduces a negative feedback loop which should help keep the simulation stable.
After these changes, the density of the high-altitude layer stopped changing. However, the pressure maps continued to be unstable, and the low-altitude map was still losing mass.
Conclusions
It turns out there were two major design flaws causing instabilities, not one. First, more heat was traveling upward than downward. Second, air was being entrained with the heat flowing up. My restructuring fixed the second issue, but the additional negative feedback loop was not enough to fix the first issue. That means the unbalanced heat flow has a deeper root cause than just the vertical exchange between the low-altitude and high-altitude layers. There is also unbalanced heat flow from the ground to the low-altitude air, and from the sun/space to the ground.
It’s likely only one of these steps is driving an imbalance in heat flow in the entire simulation, but identifying that step will require a more holistic review of the flow of heat in the simulation, starting with the insolation and radiation step.