Non-isothermal model of atmospheric density variation as altitude increases

The Earth's atmosphere is divided into layers that vary in temperature, density, pressure, and altitude. Today we will focus on one of these variables: air density; it is defined as the ratio between mass and volume $(\frac{Kg}{m^2})$ :

\rho = \frac{m}{V}

The higher the altitude and the further from the surface, the lower the density. But why?

We can attribute this decrease to three factors:

As altitude increases, there is less "air" above us, pressing down on us, decreasing the pressure and thus the density;
Because gravity decreases as altitude increases, the effective weight of the air decreases, also decreasing the density;
At higher altitudes, the air temperature becomes colder, thus becoming less dense.

Therefore, it is clear that understanding the atmospheric division is not simple, given the many factors involved. Thus, starting from a model in which the atmosphere is a single isothermal layer, that is, it has no temperature variation, and in which there is an exponential decrease in density with altitude, we will complicate our study.

The original simplified function would be:

\rho = \rho_0 e^{\frac{-h}{H}}

For $\rho_0 = 1225 \frac{Kg}{m^3}$ as the density at sea level. We're going to create code that modifies this initial idea and shows us the decrease in density considering temperature changes!

First, we establish two constants: gravity and the gas constant. These will be necessary to calculate density using the ideal gas law:

\rho = \frac{P}{RT}

g0 = 9.80665      # gravity (m/s²)
R = 287.05        # gas constant (J/kg·K)

So we established the atmospheric layers up to 80 km, according to the international standard. Each layer will have a base altitude (m), base temperature, a temperature gradient (K/m or K/km), and the final altitude. We will also assume a range for the altitude, with 800 values between 0 and 80 km.

# Layers: (base_alt, base_temp, lapse_rate, top_alt)
layers = [
   (0,     288.15, -0.0065, 11000),   # Troposfera
   (11000, 216.65,  0.0,    20000),   # Tropopausa
   (20000, 216.65,  0.001,  32000),   # Estratosfera 1
   (32000, 228.65,  0.0028, 47000),   # Estratosfera 2
   (47000, 270.65,  0.0,    51000),   # Estratopausa
   (51000, 270.65, -0.0028, 71000),   # Mesosfera 1
   (71000, 214.65, -0.002,  80000),   # Mesosfera 2
]


# Altitudes and inicialization
altitudes = np.linspace(0, 80000, 800)
densities = np.zeros_like(altitudes)

We defined the initial variables for the first layer, which will be updated as we increase the altitude, and created a loop for each layer of the atmosphere. This loop will be responsible for calculating the density in that layer, depending on whether there is a temperature gradient per kilometer or not.

# Inicial conditions
p_base = 101325  # Pa
T_base = 288.15  # K
h_base = 0       # m


alt_pointer = 0
# Loop layers
for layer in layers: #Iteramos
   h_start, T_start, L, h_end = layer
   while alt_pointer < len(altitudes) and altitudes[alt_pointer] <= h_end: # Dentro de cada camada
       h = altitudes[alt_pointer]
       delta_h = h - h_base
# Density calculation depends on whether or not there is a gradient in that layer.
       if L == 0: 
           # Isotérmica
           T = T_base
           p = p_base * np.exp(-g0 * delta_h / (R * T))
       else:
           # Com gradiente
           T = T_base + L * delta_h
           p = p_base * (T / T_base) ** (-g0 / (R * L))


       rho = p / (R * T) # Para calcular a densidade

The T_base and p_base data are updated at the end of each layer to serve as a baseline for the next layer. So, with the model established, let's now plot these relationships on a graph, with colored bands indicating the studied atmospheric layers:

These are simulated data, not based on experimental measurements. This is simply a representation of the direct proportionality between the variables.

Therefore, we can see how density decays in this more complex model. The decrease is more subtle and sensitive to the characteristics of each layer. However, the decrease is still observed and proven.

Code used for the simulations:

import numpy as np
import matplotlib.pyplot as plt
import imageio
import os


# Constants
g0 = 9.80665      # gravidade (m/s²)
R = 287.05        # constante dos gases (J/kg·K)


# Layers: (base_alt, base_temp, lapse_rate, top_alt)
layers = [
   (0,     288.15, -0.0065, 11000),   # Troposfera
   (11000, 216.65,  0.0,    20000),   # Tropopausa
   (20000, 216.65,  0.001,  32000),   # Estratosfera 1
   (32000, 228.65,  0.0028, 47000),   # Estratosfera 2
   (47000, 270.65,  0.0,    51000),   # Estratopausa
   (51000, 270.65, -0.0028, 71000),   # Mesosfera 1
   (71000, 214.65, -0.002,  80000),   # Mesosfera 2
]


# Altitudes and inicialization 
altitudes = np.linspace(0, 80000, 800)
densities = np.zeros_like(altitudes)


# Inicial conditions
p_base = 101325  # Pa
T_base = 288.15  # K
h_base = 0       # m
 

alt_pointer = 0


for layer in layers:
   h_start, T_start, L, h_end = layer
   while alt_pointer < len(altitudes) and altitudes[alt_pointer] <= h_end:
       h = altitudes[alt_pointer]
       delta_h = h - h_base


       if L == 0:
           # Isotérmica
           T = T_base
           p = p_base * np.exp(-g0 * delta_h / (R * T))
       else:
           # Com gradiente
           T = T_base + L * delta_h
           p = p_base * (T / T_base) ** (-g0 / (R * L))


       rho = p / (R * T)
       densities[alt_pointer] = rho
       alt_pointer += 1


   # Uptades base for next layer
   h_base = h_end
   delta_h_layer = h_base - h_start
   if L == 0:
       T_base = T_start
       p_base = p_base * np.exp(-g0 * (h_base - h_start) / (R * T_base))
   else:
       T_base = T_start + L * delta_h_layer
       p_base = p_base * (T_base / T_start) ** (-g0 / (R * L))


# Layers to show in the graphic
visual_layers = [
   {"name": "Troposfera",    "start": 0,   "end": 11,  "color": "#d0f0c0"},
   {"name": "Tropopausa",    "start": 11,  "end": 20,  "color": "#96dee9"},
   {"name": "Estratosfera",  "start": 20,  "end": 50,  "color": "#ffcc9a"},
   {"name": "Estratopausa",  "start": 50,  "end": 51,  "color": "#ff9a9a"},
   {"name": "Mesosfera",     "start": 51,  "end": 80,  "color": "#d79aff"},
]


# Create directory for the frames
if not os.path.exists("frames"):
   os.makedirs("frames")


# Create frames
for i in range(1, len(altitudes) + 1):
   plt.figure(figsize=(10, 6))


   for layer in visual_layers:
       plt.axvspan(layer["start"], layer["end"], color=layer["color"], alpha=0.8, zorder=1)
       plt.text(
           (layer["start"] + layer["end"]) / 2,
           0.3,
           layer["name"],
           ha='center',
           va='center',
           fontsize=9,
           rotation=90,
           color='black',
           zorder=2
       )


   plt.plot(altitudes[:i] / 1000, densities[:i], label='Densidade Atmosférica', color='black', zorder=10)
   plt.xlabel('Altitude (km)')
   plt.ylabel('Densidade (kg/m³)')
   plt.title('Densidade Atmosférica vs. Altitude (Modelo Estratificado)')
   plt.yscale('log')
   plt.grid(True, which='both', linestyle='--', linewidth=0.5)
   plt.ylim(1e-5, 1.5)
   plt.xlim(0, 80)
   plt.legend()


   filename = f"frames/frame_{i:03d}.png"
   plt.savefig(filename)
   plt.close()


# Create GIF
frames = [imageio.imread(f"frames/frame_{i:03d}.png") for i in range(1, len(altitudes) + 1)]
imageio.mimsave("atmospheric_density_ni.gif", frames, fps=10)


print("GIF com modelo atmosférico corrigido criado com sucesso!")