Cálculo numérico em R da pressão atmosférica na região da troposfera - método das diferenças finitas

Este artigo destina-se ao estudo numérico do comportamento da pressão na camada troposférica até seu limite superior, a tropopausa. Essa aproximação é válida apenas neste limite, uma vez que tanto na estratosfera quanto na termosfera a temperatura atmosférica não segue o mesmo comportamento aproximadamente linear como na troposfera e na mesosfera. O estudo aqui discutido levará em consideração um decrescimento que é aproximadamente linear na troposfera, o que está de acordo com o que é observado. Eu escolhi este tema para este tópico por ele estar alinhado com uma simulação computacional simples proposta por um professor do Departamento de Física Aplicada, da USP de São Paulo.

A equação termodinâmica de estado para a atmosfera, que é derivada considerando-se apenas ar seco ($\textit{dry air}$) e vapor, uma vez que eles representam 99,96% da concentração de gases da atmosfera. Assim: $P_a = P_d + P_v$, onde $P_d$ é a pressão para o ar seco e $P_v$ é a pressão de vapor. Para gases ideais (uma suposição bastante razoável é a de que os gases na atmosfera comportam-se como ideais), a equação é $Pv = nRT$, que pode ser reescrita em termos da densidade $\rho$ do gás:

$$P = \frac{m_a}{V}n \frac{R}{m_a}T = \rho R' T$$

onde $n~m_a$ é a massa total, $R^{'} = R/m_a$ é a constante do gás.

Com isso,

$$P_a = P_d + P_v = \rho_d R^{'}_d T + \rho_v R^{'}_v T = T (\rho_v R^{'}_v + \rho_d R^{'}_d)$$

Considere a constante $\epsilon = R_d/R_v = m_v/m_d = 0.622$, então:

$$P_a = R_d T\left( \rho_d + \frac{\rho_v}{\epsilon} \right),$$

que ainda pode ser reescrita em termos da densidade do ar ($\rho_a = \rho_d + \rho_v$) como:

$$P_a = \rho_a R_d T \left( \frac{\rho_d + \frac{\rho_v}{\epsilon}}{\rho_a} \right) = \rho_a R_d T \left( \frac{(\rho_a - \rho_v)}{\rho_a} + \frac{\rho_v}{\rho_a\epsilon}\right).$$

Neste ponto, define-se a umidade específica, $q_v = \frac{\rho_v}{\rho_a}$, de tal forma que:

$$P_a = \rho_a R_d T \left(1 - q_v + \frac{q_v}{\epsilon}\right) =  \rho_a R_d T \left( 1 + \frac{1-\epsilon}{\epsilon}q_v \right),$$

com $(1-\epsilon)/\epsilon = 0.608$. Considera-se, neste momento, uma nova temperatura, a temperatura virtual $T_v$,

$$T_v = T \left( 1 + \frac{1-\epsilon}{\epsilon}q_v \right),$$

que, conceitualmente, é a temeratura que o ar seco teria caso fosse aquecido pelo calor latente armazenado em $q_v$. Algo que é necessário deixar claro é que as grandezas $\rho$, P e $T$ variam verticalmente na atmosfera, independente da equação de estado. Como o ar não está acelerado na direção vertical, isso significa que sua força peso está em equilíbrio com o gradiente de pressão:

$$PA - (P+dP)A = \rho A g dz \\ A dP = \rho A g dz, \\ \frac{dP}{dz} = \rho g,$$

que é a equação hidrostática. Através do método de diferenças finitas, a equação diferencial anterior pode ser resolvida:

$$\frac{P_k - P_{k+1}}{z_k - z_k+1} \approx \frac{\rho_k+\rho_{k+1}}{2} \frac{g_k +g_{k+1}}{2}, \\ \Rightarrow P_k \approx P_{k+1} - \rho_k g (\Delta z_k)$$

em que foi considerado que a densidade varia pouco com a altitude, bem como a constante gravitacional g. Esta última equação pode ser facilmente resolvida em R:

 p2 = c(); p2[1] = 101325; #vetor para a pressao;
rho = 1.3; g = 9.8; 
for (z in 1: 1000){
  
  p2[z+1] = p2[1]-rho*g*z
}

atm2 = p2*9.87e-6

par(mar = c(5, 6, 3, 5))
plot(seq(1:length(atm2)),atm2,type='l', xlab='Altitude above sea level (m)',ylab = 'Air pressure (atm)')
par(new=TRUE)

hbar = p2/100;

plot(hbar,type='l',ylab="",xlab = "",xaxt = "n", yaxt = "n");
axis(side=4);
mtext("Air pressure (hPa)", side = 4, line = 3)

Outra consideração para o cálculo da pressão atmosférica, desta vez com maior acurácia, é através da junção entre a equação de estado e a equação hidrostática. Considerando uma atmosféra isoterma, então:

$$dP = -\rho g dz,\\ P = \rho R_d T_v, \\ \Rightarrow \frac{dP}{P} = - \frac{g}{R_d T_v}dz = -\frac{dz}{H},$$

onde o fator de escala $H = R_d T_v/g$. Resolvendo a equação anterior:

$$P(z) = P_0 e^{-(z-z_0)/H}.$$

 p = c(); p[1] = 101325; #vetor para a pressao;
rho = 1.3; g = 9.8; R = 287.04;
Temp = c(); Temp[1] = 300;
for (z in 1: 15000){
  
  p[z+1] = p[1]*exp(-g*z/(R*Temp[1]));
}

atm = p*9.87e-6
dist2 = seq(1:length(atm))/1000;
par(mar = c(5, 6, 3, 5))
plot(dist2,atm,type='l', xlab='Altitude above sea level (Km)',ylab = 'Air pressure (atm)')
par(new=TRUE)

hbar2 = p/100;

plot(hbar2,type='l',ylab="",xlab = "",xaxt = "n", yaxt = "n");
axis(side=4);
mtext("Air pressure (hPa)", side = 4, line = 3)

Apenas a título de curiosidade, mostro abaixo uma comparação entre a solução da equação hidrostática e a solução da nova equação diferencial, considerando atmosféra isotérmica e equação de estado para até 1000 metros acima do nível do mar:

 plot(hbar2[1:1001],type='l', xlab = "Altitude above sea level (m)", ylab = "Air pressure (hPa)" )   
  points(hbar[1:1001],type='l',col='blue')   
  legend(400,1000,c("State and hidrostatic eq.","Hidrostatic equation"), text.col = c('black','blue'))   

Nota-se que a partir de 200 m de altitute, as soluções de ambas as equações começam a divergir entre si. Agora, considero que há variação da temperatura conforme a altitude aumenta. Especificamente, na troposfera a temperatura diminui com o aumento da temperatura de forma aproximadamente linear, ou seja:

$$T(z) = T_0 - \Gamma z,$$

em que $\Gamma \approx 6.5 K/Km$. Ou seja, a temperatura diminui 6.5 K a cada quilômetro atingido acima do nível do mar. Com essa consideração, a equação diferencial para a pressão fica escrita como:

$$\frac{dP}{P} = -\frac{g}{R_m T}dz = -\frac{g dz}{R_m(T_0 - \Gamma z)}$$

Ao tomar:

$$R_m = R (1 - (1-\epsilon)q_v/\epsilon)$$

como uma grandeza constante, então, integrando os termos correspondentes (0 < z' < z  e Po < P' < P ),:

$$\ln{\left( \frac{P}{P_0} \right)} = \frac{\Gamma g}{R_m T_0 \Gamma}\ln{\left(\left[ \frac{T_0}{\Gamma }\right] - z\right)} \Biggr|_{0}^{z},\\ \Rightarrow \frac{g}{\Gamma R_m} \left[ \ln{\left( \frac{T_0}{\Gamma} - z \right) - \left(\frac{T_0}{\Gamma}  \right)} \right] = \frac{g}{\Gamma R_m} \left[ \ln{\left(\frac{\frac{T_0}{\Gamma}-z}{\frac{T_0}{\Gamma}}\right)} \right], \\ \Rightarrow \ln{\left( \frac{P}{P_0} \right)} = \ln{\left( \frac{T_0 - \Gamma z}{T_0} \right)}\frac{g}{\Gamma R_m}$$

Neste ponto, uso a seguinte propriedade do logaritmo:

$$\ln{A} = B \ln{C} \Rightarrow A = C^B,$$

ou seja,

$$P = P_0 \left( 1 - \frac{\Gamma z}{T_0} \right)^{\frac{g}{\Gamma R_m}}.$$

Esta equação pode ser facilmente resolvida numericamente assim como o caso anterior, com a diferença de que incluiremos a variação da temperatura com a altitude:

p3 = c(); p3[1] = 101325; #vetor para a pressao;
g = 9.8; R = 287.04; gamma = 6.5/1000;
Temp = c(); Temp[1] = 300;

for (z in 1: 15000){
  
  p3[z+1] = p3[1]*(1-gamma*z/Temp[1])^(g/(R*gamma));
}

atm3 = p3*9.87e-6
dist = seq(1:length(atm3))/1000;

plot(dist,atm3, type='l', xlab='Altitude above sea level (Km)',ylab = 'Air pressure (atm)')

Uma comparação desta solução com a do caso anterior, para temperatura constante, é:

 plot(atm, type='l',xlab = "Altitude above sea level (Km)", ylab = 'Air pressure (atm)')   
  lines(atm3,col='green')   
  legend(9500,1,c("Constant Temperature","Variable Temperature"), text.col = c('black','green'))   

mostrando que a solução com temperatura variável tende a ser diferente daquela para temperatura constante a partir de, aproximadamente, 5 Km.