Defining the Atmosphere

In archNEMESIS, the atmospheric structure and composition is defined using the Atmosphere class. In this notebook, we aim to provide a tutorial of the several functionalities of this class, including reading/writing the input files for running retrievals and forward models or just performing useful calculations for the post-processing of the retrievals.

import archnemesis as ans
import numpy as np
import matplotlib.pyplot as plt

1. Reading and writing the input files

The methods within the Atmosphere class can be used to read/write the information from/to the NEMESIS and archNEMESIS input files. In particular, in a NEMESIS run the atmospheric information is defined in the .ref and aerosol.ref files, while in an archNEMESIS run the atmospheric information is defined in the common HDF5 file.

NOTE: For historical reasons, the units of the aerosols in NEMESIS are defined in particles per gram of atmosphere. On the other hand, the Atmosphere class works using SI units and the aerosols are defined in particles per cubic metre. However, in order to keep the file convention consistent, the .ref and aerosol.ref files are read with the units in NEMESIS, and a unit conversion is performed internally when using the reading/writing routines. Nevertheless, since there is a density factor between the two, note that the atmospheric density needs to be defined prior to reading the aerosol file.

Reading the information from the NEMESIS files

First of all, we are going to read the information from the NEMESIS .ref and aerosol.ref files and make a summary plot based on the information in these files.

#Initialising the Atmosphere class
Atmosphere = ans.Atmosphere_0(runname='jupiter')

#Reading the .ref file
Atmosphere.read_ref()

#Making summary plot of the gaseous atmosphere
Atmosphere.plot_Atm()

#Reading the aerosol.ref file (units in file in particles per gram of atm, but units in class of particles m-3)
Atmosphere.read_aerosol()

#Making summary plot of the aerosols in the atmosphere
Atmosphere.plot_Dust()

Writing the information into the NEMESIS files

Similarly, we can use the methods in the Atmosphere class to write the information into files with the format NEMESIS requires (i.e., the specific format of the .ref and aerosol.ref files).

NOTE: Following the explanation in the previous section, there is also an internal conversion of the units from particles per m\(^3\) (archNEMESIS units) to particles per gram of atmosphere (NEMESIS units), so that the aerosol.ref file has the units it requires. Similarly, the units of pressure are also converted from Pa (archNEMESIS units) to atm (NEMESIS units) when writing the .ref file.

#Writing the .ref file
Atmosphere.write_ref()

#Wiring the aerosol.ref file (units in file in particles per gram of atm, but units in class of particles m-3)
Atmosphere.write_aerosol()

Reading and writing the information into the archNEMESIS HDF5 file

In archNEMESIS, the format of the input files was updated with respect to NEMESIS so that a unique file is required for running a retrieval or forward model. In particular, the information is included in an HDF5 file where one of the sections of the files (the Atmosphere section) corresponds to the information in the Atmosphere class. We can use the built-in functions within this class to read/write the information into this HDF5 file.

#Initialising Atmosphere class
Atmosphere = ans.Atmosphere_0()

#Reading the input HDF5 file
Atmosphere.read_hdf5('mars')

#Printing some summary information about the class to the screen
Atmosphere.summary_info()

#Writing the file into HDF5 file
Atmosphere.write_hdf5('mars')

INFO :: summary_info :: Atmosphere_0.py-339 :: Planet :: Mars
INFO :: summary_info :: Atmosphere_0.py-340 :: Number of profiles ::  1
INFO :: summary_info :: Atmosphere_0.py-341 :: Latitude of profiles ::  80.833
INFO :: summary_info :: Atmosphere_0.py-342 :: Number of altitude points ::  72
INFO :: summary_info :: Atmosphere_0.py-343 :: Minimum/maximum heights (km) ::  (np.float64(3.424995893212533), np.float64(130.01110323793642))
INFO :: summary_info :: Atmosphere_0.py-344 :: Maximum/minimum pressure (atm) ::  (np.float64(0.003644059654403374), np.float64(6.430470874101519e-10))
INFO :: summary_info :: Atmosphere_0.py-345 :: Maximum/minimum temperature (K) (np.float64(184.94036865234375), np.float64(131.26361083984375))
INFO :: summary_info :: Atmosphere_0.py-347 :: Maximum/minimum gravity (m/s2) ::  (np.float64(3.73), np.float64(3.47))
INFO :: summary_info :: Atmosphere_0.py-349 :: Maximum/minimum molecular weight ::  (np.float64(0.04365860339841797), np.float64(0.04249212487703371))
INFO :: summary_info :: Atmosphere_0.py-350 :: Number of gaseous species ::  15
INFO :: summary_info :: Atmosphere_0.py-357 :: Gaseous species ::  ['CO2', 'O3', 'CO', 'O2', 'N2', 'H2', 'O', 'H', 'Ar', 'H2O (1)', 'H2O (2)', 'H2O (3)', 'H2O (4)', 'H2O (5)', 'HO2']
INFO :: summary_info :: Atmosphere_0.py-359 :: Number of aerosol populations ::  1

2. Modifying the gaseous Atmosphere

In addition to having methods relevant to read/write the input files, there are also some built-in methods to modify the parameters in the atmosphere. These methods make sure that the dimensions of the arrays and other things are all consistent, giving some errors back otherwise, but note that it is also easy to modify the information in the class just by calling it as Atmosphere.VARIABLE.

Add a new gas in the atmosphere

#Add OH in the atmosphere
ID_oh = 13                                  #Radtran ID for OH
ISO_oh = 0                                  #Radtran isotope ID (0 for all isotopes)
vmr_oh = np.ones(Atmosphere.NP) * 1.0e-9    #1 ppbv of OH (constant with altitude)

Atmosphere.add_gas(ID_oh,ISO_oh,vmr_oh)

Remove gas from the atmosphere

#Remove the fifth isotope of water
ID_h2o = 1
ISO_h2o = 5

Atmosphere.remove_gas(ID_h2o,ISO_h2o)

Update volume mixing ratio of gas in the atmosphere

#We want to update the VMR of Ar
ID_Ar = 76                                  #Radtran ID for Ar
ISO_Ar = 0                                  #Radtran isotope ID (0 for all isotopes)
vmr_Ar = np.ones(Atmosphere.NP) * 0.01

Atmosphere.update_gas(ID_Ar,ISO_Ar,vmr_Ar)

Re-adjust altitudes based on hydrostatic equilibrium

Using the function Atmosphere.adjust_hydrostatH() we will keep the pressure levels in the atmosphere fixed and re-compute the altitudes based on the hydrostatic equation. The altitudes are re-computed with respect to the altitude level that is closest to 0 km.

#Saving array before the re-adjustment
h_before = Atmosphere.H/1.0e3

#Re-calculate altitude levels based on hydrostatic equilibrium equation (fixing pressure levels)
Atmosphere.adjust_hydrostatH()

#Saving results after re-adjustment
h_after = Atmosphere.H/1.0e3

#Printing the results
print(f"P (bar), {('H_before (km),','H_after (km)')}")
#for i in range(Atmosphere.NP):
for i in range(10):
    print(Atmosphere.P[i]/1.0e5,h_before[i]/1.0e3,h_after[i]/1.0e3)

P (bar), ('H_before (km),', 'H_after (km)')
0.003692343444824219 0.0034249958932125327 0.0034249958932125327
0.00247947998046875 0.006422610045284344 0.0064134805288371625
0.001654933319091797 0.009428527090137006 0.009408455428165846
0.0013440388488769532 0.010950178536527432 0.010926622599227689
0.001088509826660156 0.012471829982917857 0.012447284370977868
0.0008802003479003906 0.013994307152503492 0.013962412720921613
0.0007102473449707032 0.015515958598893919 0.015477438567657134
0.0005719057846069336 0.017037610045284347 0.01699368564809638
0.0003315441131591797 0.020813390107650834 0.02076086444406526
0.00019087339401245118 0.024593344446822107 0.024531071268190235

Re-adjust the pressures based on the hydrostatic equilibrium equation

Using the function Atmosphere.adjust_hydrostatP() we can fix the altitude levels in the atmosphere and re-compute the pressure levels. The reference pressure at a given tangent height must be defined, and the rest of the pressure levels are re-computed with respect to this.

#Saving the pressure array before computations
p_before = Atmosphere.P

#Performing the calculations
htan = 12442.    #Altitude at which to fix the pressure (m)
ptan = 1.0e+02   #Pressure at htan (Pa)
Atmosphere.adjust_hydrostatP(htan,ptan)
p_after = Atmosphere.P

#Printing the results
print(f"H (km), {('P_before (Pa),','P_after (Pa)')}")
#for i in range(Atmosphere.NP):
for i in range(10):
    print(Atmosphere.H[i]/1.0e3,p_before[i],p_after[i])

H (km), ('P_before (Pa),', 'P_after (Pa)')
3.424995893212533 369.2343444824219 338.9658145719878
6.413480528837162 247.947998046875 227.62185299033686
9.408455428165846 165.4933319091797 151.92608530847525
10.926622599227688 134.4038848876953 123.3850419043048
12.447284370977867 108.85098266601561 99.92674766208638
13.962412720921613 88.02003479003906 80.8033578162203
15.477438567657135 71.02473449707031 65.2011970315601
16.99368564809638 57.19057846069336 52.501071501440165
20.76086444406526 33.15441131591797 30.435351180890827
24.531071268190235 19.087339401245117 17.521623460907847

3. Perform useful calculations

Some of the methods in the Atmosphere class also include functions that are normally performed as internal calculations within archNEMESIS, but that can be useful from the user end too. Here, we show a few examples of these calculations. In particular, we are going to calculate:

• Molecular weight (kg/mol).

• Gravity (m/s\(^2\)).

• Number density (molec m\(^{-3}\)).

• Mass density (kg m\(^{-3}\)).

#Performing calculations
#####################################################################

#Molecular weight (based on the molecular weight and vmr of each of the gases)
Atmosphere.calc_molwt()

#Gravity
Atmosphere.calc_grav()

#Number density
numdens = Atmosphere.calc_numdens()

#Specific density
rho = Atmosphere.calc_rho()


#Making summary plot
######################################################################

fig,(ax1,ax2,ax3,ax4) = plt.subplots(1,4,figsize=(10,3),sharey=True)

ax1.plot(Atmosphere.MOLWT,Atmosphere.H/1.0e3)
ax1.set_xlabel('Molecular weight (kg mol$^{-1}$)')
ax1.set_ylabel('Altitude (km)')
ax1.set_facecolor('lightgray')

ax2.plot(Atmosphere.GRAV,Atmosphere.H/1.0e3)
ax2.set_xlabel('Gravity (m s$^{-2}$)')
ax2.set_facecolor('lightgray')

ax3.plot(numdens,Atmosphere.H/1.0e3)
ax3.set_xlabel('Number density (m$^{-3}$)')
ax3.set_xscale('log')
ax3.set_facecolor('lightgray')

ax4.plot(rho,Atmosphere.H/1.0e3)
ax4.set_xlabel('Specific density (kg m$^{-3}$)')
ax4.set_xscale('log')
ax4.set_facecolor('lightgray')

plt.tight_layout()