A Discrete Model Atmosphere

(A few additional comments since this post seems to have attracted certain attention. Although it might look odd, in reality it describes an ordinary diffusion process. The corresponding differential equation is described on the post "A new attempt". Here, the thermal conductivity k is normalized to 1 but if you want to vary it all you have to do is to divide the forcing vector F by k. The derivation of this I leave to the reader.)


This post is a bit technical but if you bear with me I hope it might shed some light on certain important topics. The reason why I have constructed a discrete model is that by doing so I have avoided the tricky problem of specifying the boundary conditions, which you are invariably faced with if you formulate the model in terms of differential equations. The model which I am going to present admittedly bear some SUPERFICIAL ressemblance with the greenhouse effect, but one of the questions I am asking is if it ACTUALLY reproduces the greenhouse effect. Take a look at the diagram below:

The atmosphere has been discretized into a finite number of layers each having its own "weight" w1, w2 and so on. The arrows may represent either a "real" or a "formal/mathematical" energy flow. What do I mean by this? There is actually a point in not specifying this at the moment. It could be interpreted as if the "weight" of the atmospheric layer is simply a function of the amount of greenhouse gases present and that the arrows represent upward or downward radiation. On the other hand, if we allow a certain "flexibility" in our thinking, the "weight" could be interpreted as a function of the actual weight (mass) of the layer and the arrows merely illustrate a "formal" energy conserving representation of the obstruction to cooling caused by these layers. Although that may sound like mumbo-jumbo in your ears, somehow we have to account for the eventual energy loss to outer space, and I haven't found any better way to do it than this. Now for the technicalities:

The amount of energy absorbed per unit time by each layer is the amount of incident energy times its wight. (Hence, the wight must be less than 1). Each layer loses energy per unit time proportional to its temperature AND weight. This loss is equally distributed upwards and downwards. Keep in mind that the energy that is not absorbed by the layer is free to pass on to the next and so on. (This is accounted for by inserting a factor (1-w) for each intermediate layer). Below is a Mathematica worksheet where the energy balance equations for the STATIONARY state of a model with four layers are written down and solved. Staionary means that all of the layers have reached energy balance, hence, the amount of energy lost per unit time equals the amount of energy absorbed per unit time. The left hand sides represents the energy lost and the right hand sides contains the absorbed energy.

This model is linear in temperature, unlike mainstream climatology where one is instead obsessed with the SB-law. A linear dependence on temperature may very well be physically motivated, especially if we are dealing with ideal gases, though admittedly, here it is primarily employed for computational efficiency, things would become very much more cumbersome with a T4-dependence. If the wight of the layers are assumed to decline exponentially with altitude, from the analytical and numerical solutions we see not only the formation of a stratosphere but we can also infer with almost certainty (although I have no strict proof for this) that the temperature approaches the finite value F/2 at the top of the atmosphere, where F is the solar forcing at the surface. This can be verified also with a numerical algorithm solving the equations for a model atmosphere discretized in 60 layers as is shown below (The Python-script for this is appended at the end)

What is then the main conclusion to be drawn from this. The most obvious one is that this model does not predict a cooling of the stratosphere upon an increases in the weight-function. On the contrary it either produces a warming of the entire atmosphere upon an increase in the weight or solar forcing, or conversely a cooling of the entire atmosphere upon a decrease in solar forcing or weight. In other words, the temperature gradient does not pivot, and this holds true even if the weight is interpreted as the amount of greenhouse gases present. (Admittedly it also somewhat contradicts some of my previous statements/guesses, which the reader can verify on his/her own.) 

How the cooling of the stratosphere comes about in the mainstream models is very obscure. From the most simple equations I can see that this follows from foolishly applying a Dirichlet boundary condition at the surface rather than at the top of the atmosphere, but how it emerges in their time-stepping algorithms remains a mystery to me. One of the lessons to be learned, however, is that there is much mathematical subtlety involved and that any "wordy" explanations from people like Roy Spencer must now be replaced by reproducible algorithms/equations.

Moving on, one could now ask if the model presented above could actually be useful as a model atmosphere. My answer to that is: Yes maybe, if the weight is interpreted as a function of the TOTAL MASS. The reason why I believe this is that:

1. It predicts an elevation in surface temperature as a function of atmospheric mass

2. It predicts the formation of a stratosphere whose temperature is a function of solar forcing alone

3. It can easily be extended so as to reproduce a thermosphere higher up if you add solar forcing terms to the uppermost layers.


Below is the Python code for your enjoyment:

***

import numpy as np
import matplotlib.pyplot as plt

interval = 6
meshsize = 0.1
N = int (interval/meshsize)

weight = np.zeros(N)
weight[0] = 1


for idx in range(1,N):

    weight[idx] =  np.exp(-idx*meshsize)*meshsize
    


A = np.zeros([N,N])
    
for idx1 in range(N):
    for idx2 in range(N):
        if idx1 == idx2: 
            if idx1 == 0:
                A[idx1,idx2] = -1
            else:
                A[idx1,idx2] = -2
            
        else:
            A[idx1,idx2] = weight[idx2]

            if idx2>idx1:

                for idx3 in range(idx1+1, idx2):
                    A[idx1,idx2] = A[idx1,idx2]*(1-weight[idx3])

            else:
                for idx3 in range(idx2+1, idx1):
                    A[idx1,idx2] = A[idx1,idx2]*(1-weight[idx3])    
            
F = np.zeros(N)
F[0] = 1



temp = np.linalg.solve(A,-F)
print temp

x = np.arange(0,interval,meshsize)

plt.plot(x, temp)
plt.ylabel('Temperature')
plt.xlabel('Position')
plt.show()

***