Types of atmospheres
PSG handles two types of atmospheres: hydrostatic equilibrium (typical for planets)
and expanding coma (typical of comets and small bodies). In both cases, a set of three
parameters can be used to describe the overall structure of the atmosphere: temperature (K), mean molecular weight and
surface pressure (for equilibrated atmospheres) or production rate (for expanding ones).
PSG permits providing detailed vertical information of molecular abundances and temperatures,
and atmospheric templates (vertical profiles of temperature and
abundances) are available for the main atmospheres (Venus, Earth, Mars,
Titan, Neptune, Uranus), while general atmospheric and surface parameters
are available for the other bodies. For expanding atmospheres, PSG assumes isotropic outflow,
constant expanding velocity (following vexp
= 0.8 Rh-0.5
a dust outgassing model.
The approximate brightness for the object (and for the parent star in the case of exoplanets) is displayed in the surface section,
which is computed based on the object size, albedo, distances and stellar temperature.
For cometary coma, three brightness metrics are calculated and displayed next to the production rate field:
the expected visual magnitude (mv
), the infrared Figure-Of-Merit (FOMIR
the radio Figure-Of-Merit (FOMRadio
). These parameters are calculated based on
the heliocentric distance (RH, [AU]), the geocentric distance (Delta, [AU]), and the gas production rate (Q, [mol/s]),
employing these approximations: mv
= (30.675 - log(Q))/0.2453 + 5⋅log(Delta) [Jorda et al. 2008]; FOMIR
= Q⋅1E-28/Delta; and where Q = QAU
|Types of atmospheres:
the user can select between two types of atmospheres: hydrostatic equilibrium (typical for planets, profiles adapted from Robinson & Catling 2014)
and expanding coma (typical of comets and small bodies). For the main planets, vertical profiles are available,
while the user can also load any arbitrary vertical strucutre. For expanding atmospheres, PSG assumes isotropic outgassing,
a constant temperature across the coma, and an outgassing velocity established by the heliocentric distance.
Defining the profiles and abundances
Two important steps are addressed in this section: 1) the structure and composition of the atmosphere/surface, 2) the databases
that are being used for each component. PSG ingests a broad range of spectroscopic information, with the following types:
0:REFL Reflectance: the reflecting properities of the material are described as a scaling factor (0-1).
1:OPTC Optical constants: the spectroscopy is tabulated as n-k, where n is the refractive index and k is the extinction coefficient.
2:ALPH Alpha parameter: the α index indicates the extinction coefficient per slab width.
3:XSEC Cross sections: the molecular absorptions are described as a table of cross-sections [cm2/molecule].
4:SCAL Scattering function (Legendre): the scattering function is described as a summation of Legendre polynomials.
5:SCAT Scattering function (Henyey-Greenstein): the scattering function is defined as a H-G function.
6:LBLN Line-by-line database: each line is described individually, requiring expensive line-by-line RT calculations.
7:CKTB Correlated-k tables: pre-computed corr-k opacity tables are available in PSG for efficient RT calculations.
The user can provide an arbitrary profile of P/T and of molecular abundances by selecting "File Template".
The format of this file is a text file with all the entries as in the configuration file for the sections "ATMOSPHERE" and "SURFACE"
The pressures are provided in 'bars', temperatures in 'K' and abundances in volume mixing ratio -
[molecules / molecules] for gases (e.g., water vapor, methane) and [kg / kg] for hazes (e.g., water ice clouds, methane ice). Particle sizes in [m] are provided by adding the keyword '_size'.
The base (molecules or kg) is always wet air (all gases, including water vapor) without hazes. For Earth, the standard profiles (US-Standard,
Tropical, mid-latitude-winter, midlatitude-summer, subarctic-summer, subarctic-winter) are available by selecting "Template", while
for Mars an array of climatological templates computed using the LMD General Circulation Model (Millour et al. 2008) are available.
For exoplanets, PSG will select from a list of pre-computed templates (considering the planet's mass and density), and
will compute P/T profiles for giants based on a non-gray analytical model (Parmentier et al. 2014, A&A, Vol. 562, A133) according to the planet's
equilibrium temperature (based on distance to star, stellar radius and stellar temperature) and the planetary gravity.
Abundances are computed layer-by-layer considering the equation-of-states (EOS) computed by Kempton et al. (PASP, Vol. 129, Issue 974, pp. 044402, 2017; Mbarek & Kempton, ApJ, Vol. 827, Issue 2, id. 121, pp. 10, 2016) for the
T-P grid of 100 to 3000 K and 1000 to 1E-9 bar.
When no vertical profile is provided, the temperature (T) is assumed to be constant across the atmosphere, and the pressure (P)
decreases with altitude (z) following the scale-height: P = Psurf
where g is the gravity (defined in the object section) and R is the gas constant (8.3144598 J / K / mol).
In order to change the abundance of gases (or their isotopologues), first select the gas in the left column (press "Search" to
search for available species within the catalogues), then select the isotopologue / linelist in the "Type" section.
For instance, the HITRAN database is described by "HIT" (first number is mol ID, second is isotopic numer), the JPL Molecular Spectroscopy database with "JPL", the Cologne Database of Molecular Spectroscopy by "CDMS",
and the GSFC Fluorescence database by "GSFC".
: most atmospheres are a mix of gases and of aerosols (dust/ice/clouds/hazes),
with the latter scattering the light in a complex pattern depending on the viewing angle,
aerosol size/shape and composition. Inclusion and realistic treatment of this phenomenon in the radiative transfer
analysis is extremely computationally expensive, yet several approaches to attack this problem do exist.
An efficient scattering package is available in PSG, which is based on a Martian scattering model.
It incorporates the latest radiative transfer numerical methods, and it is parameterized
for LTE (Local-Thermodynamic-Equilibrium) calculations.
(Smith, M. D. et al., JGR, Vol. 114, E00D03, 2009; Villanueva, G. L. et al., Science, Volume 348, Issue 6231, pp. 218-221, 2015).
Aerosols abundances can be provided in terms of total aerosol column (kg/m2),
a scaling factor (scl) of the provided vertical profile, or a fixed value across alitude (%, ppm, ppb, ppt).
When an aerosol profile is available, PSG will employ this vertical distribution, but will scale the column
to match the defined total column(kg/m2). The natural scaling unit of PSG is kg/m2.
Molecular/atomic databases (1259 species)
The key parameter to be entered in the "Atmosphere" section is the definition of the molecular species and
the corresponding linelist to be used for these molecules. HITRAN is generally very complete (for IR, optical and UV at low temperatures)
for most typical planetary atmosphere and has become the main repository of line information. At radio
wavelengths, the JPL Molecular Spectroscopy and Cologne Database of Molecular Spectroscopy (CDMS) are generally
more complete and have a better description of the rotational spectrum of complex molecules. NASA-Goddard currently
holds the main repository for non-LTE fluorescence linelists, suitable when synthesizing cometary spectra in the UV/optical/IR range.
Wavelength range: 0.3 μm to radio|
Number of lines: 5,399,562
Number of molecules: 50
Number of isotopologues: 126
Number of CIA spectra: 554
Number of cross-section spectra: 987
Number of aerosols: 98
|Gordon, I.E. et al; "The HITRAN2016 molecular spectroscopic database"; Journal of Quantitative Spectroscopy and Radiative Transfer (2017)|
Wavelength range: 0.2 to 100,000 um
Number of temperatures: 20 (40 to 2000 K)
Number of pressures: 17 (1E-6 to 100 bar)
Number of species: 21
|Absorption coefficients computed with PUMAS, and assuming wings of 25 cm-1 and a fine core of 1 cm-1 where maximum resolution calculations are applied.|
Molecules: H2O, CO2, O3, N2O, CO, CH4, O2, SO2, NO2, NH3, HCl, OCS, H2CO, N2, HCN, C2H2, C2H4, PH3, H2S, C2H4, H2
|Line-shape parameters (line width, pressure shift, and temperature dependence of these factors)|
Air: H2O CO2 O3 N2O CO CH4 O2 NO SO2 NO2 NH3 HNO3 OH HF HCl HBr HI ClO OCS H2CO HOCl N2 HCN CH3Cl H2O2 C2H2 C2H6 PH3 COF2 SF6 H2S HCOOH HO2 O ClONO2 NO+ HOBr C2H4 CH3OH CH3Br CH3CN CF4 C4H2 HC3N H2 CS SO3 C2N2 COCl2 He
Self: H2O CO2 O3 N2O CO CH4 O2 NO SO2 NO2 NH3 HNO3 HF HCl HBr HI ClO OCS H2CO N2 HCN CH3Cl C2H2 C2H6 PH3 COF2 H2S HCOOH HO2 O NO+ C2H4 CH3OH CH3Br CH3CN C4H2 HC3N H2 CS SO3 C2N2 COCl2
H2: H2 CO SO2 HF HCl OCS C2H2 (soon CO2 N2O H2CO HCN H2S OH)
He: He CO SO2 HF HCl OCS C2H2 (soon CO2 N2O H2CO HCN H2S OH)
CO2: H2O CO2 CO SO2 HF HCl OCS C2H2 (soon N2O H2CO HCN H2S OH)
|GSFC Fluorescence Database||
Wavelength range: shorter than 10 μm|
Number of lines: 530,281
Based on lines (non-LTE): 2 billions
Number of species: 24
For daughter species, rotational populations are heavily affected by photodissociation, and for simplicity in PSG, we assume
an elevated Trot of 600K for OH,CN,CH,NH and 3200K for C2.
|Villanueva, G. L., Mumma, M. J., DiSanti, M. A., Bonev, B. P., Gibb, E. L., Magee-Sauer, K., Blake, G. A., Salyk, C., "The molecular composition of Comet C/2007 W1 (Boattini): Evidence of a peculiar outgassing and a rich chemistry". Icarus, Volume 216, Issue 1, p. 227-240. (2011)|
Wavelength range: 0.3 μm to radio|
Number of lines: 5,023,277
Number of species: 52
Number of isotopologues: 118
|Jacquinet-Husson, N. et al; "The 2015 edition of the GEISA spectroscopic database"; Journal of Molecular Spectroscopy, Volume 327, Pages 31-72 (2016)|
|JPL Molecular Spectroscopy||
Wavelength range: 2.65 μm to radio|
Number of lines: 888,113
Number of species: 383
|Pickett, H. M.; R. L. Poynter, E. A. Cohen, M. L. Delitsky, J. C. Pearson, and H. S. P. Muller, "Submillimeter, Millimeter, and Microwave Spectral Line Catalog"; J. Quant. Spectrosc. and Rad. Transfer 60, 883-890 (1998)|
|CDMS Cologne Database for Molecular Spectroscopy||
Wavelength range: 1.81 μm to radio|
Number of lines: 1,612,154
Number of species: 792
|Muller, Holger S. P.; Schloder, Frank; Stutzki, Jurgen; Winnewisser, Gisbert; "The Cologne Database for Molecular Spectroscopy, CDMS: a useful tool for astronomers and spectroscopists"; Journal of Molecular Structure, Volume 742, Issue 1-3, p. 215-227 (2005)|
Wavelength range: 0.1 to 170 um|
Number of spectral points: 7454
Number of temperatures: 30 (100 to 3000 K)
Number of pressures: 13 (1E-9 to 1000 bar)
Number of species: 30
|Kempton, E. M-R, Lupu, R., Owusu-Asare, A., Slough, P., Cale, B., "Exo-Transmit: An Open-Source Code for Calculating Transmission Spectra for Exoplanet Atmospheres of Varied Composition", PASP 129 044402 (2017)|
|The MPI-Mainz UV/VIS Spectral Atlas||
Wavelength range: 0.01 to 1 um|
Number of cross-sections: 52 (selection)
Number of species: 22 (selection)
Keller-Rudek, H., Moortgat G.K., Sander R., Sörensen R., "The MPI-Mainz UV/VIS Spectral Atlas of Gaseous Molecules of Atmospheric Interest", Earth Syst. Sci. Data, 5, 365–373 (2013)|
Venot et al., A&A, 609-A34 (2018)
Serdyuchenko et al.,Atmos. Meas. Tech., 7, 625-636, 2014
|The CFA/Harvard Kurucz atomic database
Wavelength range: 0.001 to 1000 um|
Number of species: 80 elements and their ions
Number of lines: 2,309,497
Kurucz, Smith, Heise, Esmon (gfall08oct17), "The CFA/Harvard Kurucz Atomic spectral line database" (2017)'|
Modeling: the model of the lineshapes includes a Voigt analysis (Thermal, Natural, van-der-Waals) within the impact theory region and
an extended wing region beyond for the statistical regime. The transition between these regimes is defined by the detuning frequency (Nefedov+1999, Burrows+2000, Iro+2005).|
Einstein Aul [s-1] = gf / (1.499E-14 ⋅ (2J+1) ⋅ λ2) where gf is the tabulated oscillator strength, J is the upper-state rotational number and λ is the wavelength [nm].
Natural halfwidth [cm-1] = γrad / (4πc) where γrad [s-1] is the radiative damping constant and c the speed of light [cm/s].
Collissional halfwidth [cm-1/atm] (T/296)-0.7 = γ6 ⋅ 8.62489e+18 / (4πc) where γ6 [s-1] is the van der Waals damping constant (neutral hyrogen) and T is the temperature [K].
Detuning frequency [cm-1] = 25 ⋅ [T/(ma⋅250)]0.5 in this approximation, ma is the molar mass [amu] of the atmosphere.
Extended lineshape model = (v-v0)-3/2⋅exp(-h(v-v0)/kT) where v0 is the line center [cm-1].
The presence of aerosol particles (dust/ice/clouds/hazes) in a planetary atmosphere has a significant impact on the intensity and morphology of planetary spectra, and can severely hinder the ability to infer atmospheric properties such as abundances from gaseous absorption signatures and atmospheric temperatures (e.g., Sing et al. 2013; Vasquez et al. 2013a; 2013b). In addition, without explicitly including the radiative effects of aerosols (including their particular microphysical properties), the inconsistencies between the predicted spectra and the subsequently inferred equation of state are likely to be considerable, as has been demonstrated for the well-studied case of Mars (e.g., Madeleine et al. 2011; 2012).
The general approach of the community has involved the coupling of a highly developed molecular transmission code (i.e., LBLRTM, GENLN3) with a multiple scattering radiative transfer algorithm (i.e., DISORT, Stamnes et al. 1988; Lin et al. 2015) such as has been done by Turner and collaborators with LBLDIS (Turner 2005). However, while these codes are publicly available and include documentation by the authors, their use as a tool has generally been limited to radiative transfer specialists and the availability of non-trivial computational resources. In PSG, we provide this capability for specific set of modes that balance spectral resolution and range with computation efficiency.
A treatment of multiple scattering from atmospheric aerosols is enabled within PSG using the discrete ordinates method (e.g., Goody and Yung, 1989; Thomas and Stamnes, 1999; Smith et al., 2013). The radiation field is approximated by a discrete number of streams distributed in angle with respect to the plane-parallel normal. The number of stream pairs (pairs of corresponding upward and downward radiation streams) can be set as high as necessary to accurately model the angular dependence of the aerosol scattering phase function while maintaining computational feasibility. The “two-stream approximation” often used when modeling planetary atmospheres is an example of the discrete ordinates method using one stream pair.
The angular dependence of the scattering phase function for a particular aerosol is described in terms of an expansion in terms of Legendre Polynomials, typically with the number of expansion terms equal to the number of stream pairs. As implemented in PSG, the Legendre expansion coefficients are pre-computed using an assumed particle size distribution for each available aerosol type and using either Mie scattering (e.g., Wiscombe, 1980) or T-matrix (e.g., Mishchenko et al., 1996) codes as specified in the associated information files for each aerosol type. The underlying indices of refraction for aerosols are empirically derived from spacecraft observations in the case of Mars dust and water ice aerosols (Wolff et al., 2009), or from the HRI (HITRAN Refractory Index, Massie and Hervig 2013 and references therein) database. For the case of the HRI constants, we calculate scattering coefficients employing a Mie implementation (Bohren and Huffman 1983) that derives Henyey-Greenstein scattering g-factors. Internally, PSG converts these g-factors into Legendre expansion coefficients in order to be ingested by the radiative-transfer suite.
The discrete ordinates formulation computes the diffuse radiation field for a plane-parallel atmosphere. When spherical geometry is important (e.g., limb geometry observations), the pseudo-spherical approximation (e.g., Spurr, 2002; Thomas and Stamnes, 1999) is used for computational efficiency. In this scheme, the source functions computed using the diffuse field from the discrete ordinates plane-parallel geometry are integrated along an equivalent curved path through the model layers. This curved path is defined by computing the correct emission angle for the path at the boundary of each layer. The pseudo-spherical approximation is accurate over a wide range of conditions and is orders of magnitude faster than an “exact” Monte Carlo code (Smith et al., 2013).
Converting aerosol abundances
: PSG ingests aerosol abundances X in [Kg/Kg] or [g/g], yet many references report aerosols
in several different units. A common unit is N [particles/cm3], to determine X [Kg/Kg] for N=1 [particles/cm3
], P=5e-3 [bar], T=340 [K], matm
=18 [g/mol], particle radius (not diameter) rhaze
=1e-5 [cm], ρhaze
= 1E-1 ⋅ P / (KB [1.38064852e-23] ⋅ T) = 1.065e+17 [molecules/cm3
= π (4/3) ⋅ rhaze3
= 4.2e-15 [cm3
] Volume of each haze particle
= 5.7e-15 [g]
= N ⋅ Mhaze
⋅ AVOG[6.022140857e23] / (ρatm
) = 1.8E-9 or 1.8 [ppb]
Defining the surface
The surface defines one of the "boundary" conditions of the radiative-transfer analysis. Currently, PSG describes the surfaces as Lambertian (with the corresponding albedo) and emitting as a black-body (with the corresponding emissivity).
We have developed a versatile surface module that combines a realistic Hapke scattering model (Protopapa et al. 2017) and the capability to ingest a broad range of optical constants,
permitting PSG to accurately compute surface reflectances and emissitivities. The surface model considers areal mixing and a one-term Henyey-Greenstein phase funtion.
Cometary/exospheres dust particles:
Dust and icy particles in cometary comae are in many cases the main source of continuum/broad radiation
in small bodies and in exospheres. In PSG, the surface properties of the nucleus and the dust particles are assumed to be same,
with the exception of the phase function. Small particles have a very different response to the
solid nuclear body, with a strong forward scattering peak and a less shallow scattering phase function. PSG employs the Halley-Marcus (H-M) phase function
compiled by Schleicher et al. 2011 to characterize the differing scattering efficiencies with phase (download table
The intensity of the dust continuum is related to the ammount of dust particles, and different quantities exist to quantify this:
- Dust/gas ratio: in this approach the dust particles are treated as behaving like the surrounding gas
and a dust/gas mass ratio of 1.0 provides consistent continuum fluxes to the brightness vs. gas-activity relationship identified by Jorda et al. 2008.
The reflected sunlight flux is affected by the H-M phase curve, while the thermal emission is assumed to be isotropic and not affected by phase.
- A(Θ)fρ: is a quantity introduced by A'Hearn et al. in 1984 (AJ 89, 579, 1984) that describes continuum intensity and is generally
independent of the different image-scale and measuring window sizes used in the photometry. Since this quantity intrinsically
includes a phase correction, the A(Θ)fρ reflected fluxes are not corrected by the H-M phase curve, yet the thermal fluxes are
corrected by 1/p(H-M).
Properly modeling spectroscopic features of planetary surfaces over a wide wavelength range requires a comprehensive and inclusive spectroscopic database.
These parameters are used to establish the "boundary" conditions for the radiative transfer calculations.
There is currently no single repository that integrates optical constants and reflectances of solid surfaces and ices over a wide spectral range (only
specialized databases exist). We have identified eleven spectral databases that are applicable to the synthesis of planetary spectra,
and we have developed a program to standardize and homogenize these libraries of spectral constants.
Retrieval of planetary parameters
PSG allows to compare user-provided data to synthetically generated spectra, and to perform physical retrievals of
the parameters included in the model employing the optimal estimation method
The retrieval module also implements an additional regularization parameter that improves the retrieval accuracy
and convergence speed of the inverse scheme. The optimal estimator that PSG uses to get an estimate of the parameters
from observations is obtained by considering a linear Gauss–Newton iterative scheme. The formal retrieval equation
implemented follows this scheme:
x (γSa-1 + KTSy-1K) = y (KTSy-1)
The scheme is based on the a-priori information available for the variability of the input quantities
and the knowledge of the measurement uncertainty; if no error is provided, PSG will assume a 5% uncertainty.
The minimum/maximum limits entered in the form for each parameter as assumed to be 5-sigma, and therefore
the "variance" of each parameter is defined as (maximum-minimum)/10. Both covariance matrices, data and parameters, are assumed
to be diagonal and filled with the squares of the data and parameter variances, respectively. The Jacobians are computed
with forward steps of 0.1% of the maximum value.
The parameter γ
is the extra regularization parameter, and acts as a tradeoff between the a-priori values
given to the parameters to be retrieved (x-space) and the observations (y-space).
Large values of γ will constrain the retrieval scheme more to the a-priori parameter values,
hence a large γ might be recommended in case of reliable a-priori knowledge of retrieval parameters.
As γ approaches 0, the solution scheme tends to a constrained least-square. For γ = 1, the Rodgers’ classical scheme is run.
Together with the optimal estimate for the parameters, the retrieval module implemented in PSG
computes a series of a-posteriori quantities. The uncertainties of the parameters are computed
from the a-posteriori covariance matrix, which is directly derived from the Rodgers’ formalism.
In addition, the averaging kernels matrix
(AK) is computed for each of the retrieved parameters.
AKs provide information about the information content of the data with respect to the retrieved parameters.
The diagonal of the AK matrix contains values between 0 and 1, which indicate the sensitivity
of the retrieval to each parameter. Low values could be indicative of intrinsic poor sensitivity
of the data themselves to the retrieval parameters, or could be due to very stringent constraints
imposed to the variability of the parameters to be retrieved.
In general, the contribution functions matrix
and the AKs are presented to highlight the properties
of the retrieved parameters, and their significance. AKs can also be provided prior to the retrieval,
to allow the user to evaluate the information content of the spectrum, and tune the input parameters.
The retrieval provides also an accurate a-posteriori estimation of the measurement uncertainty.
Beyond fitting the data to the spectroscopic model, PSG can also correct for a broad range of instrumental issues typically affecting data (see the figure below).
By employing pre-computed telluric spectra, PSG can also fit the water and column parameters affecting ground-based data.
For fitting these instrumental effects, PSG employs the Levenberg–Marquardt algorithm, also known as the damped least-squares method,
in solving the non-linear least squares problem. This method interpolates between the Gauss–Newton algorithm and
the method of gradient descent, with the Levenberg–Marquardt method being more robust, meaning that it finds a solution even if
the initial conditions are far from the final solution. The Levenberg-Marquardt algorithm is based on
the MPFIT program (Craig Markwardt), which draws from the robust package called MINPACK-1 (Garbow, Hillstrom and More).
|Optimal Estimation Method (OEM) C/C++ package
Developed by Liuzzi and Villanueva, Sep/2018
Based on methodologies described in Rodgers 2000 and therein.
|Levenberg-Marquardt (LM) minimization C/C++ package
Developed by Craig Markwardt, June/2016
Based on the MINPACK-1 (Garbow, Hillstrom and More).
The input file should be formatted as a text file with two or three columns. The first column should indicate
the frequency/wavelength of each pixel, while the second column should describe the measured flux.
A third column (optional) should indicate the 1-sigma uncertainty in flux, and if no error is provided, PSG
will assume a 5% uncertainty. See below a list of example configuration files (that also include the data):
- Comet IR retrieval: retrieval of production rates and rotational temperatures in comet C/2007 W1 (Boattini) from high-resolution Keck/NIRSPEC data acquired on July/9/2008 (Villanueva et al. 2011, Icarus 216, 227-240).
- Mars retrieval: retrieval of deuterated water from high-resolution Keck/NIRSPEC data acquired in January/2014. This is a relatively complex and expensive retrieval, in which several instrumental issues are removed simultaneously (Villanueva et al. 2015, Science 348-6231).
- Comet radio retrieval: retrieval of cometary methanol production rate and its rotational temperature from high-resolution data of comet C/2002 T7 (LINEAR) acquired with the Heinrich Hertz Submillimeter Telescope (Villanueva et al. 2005, AAS/DPS meeting).
- Earth retrieval: retrieval of methane and carbon dioxide in Earth's atmosphere via laser heterodyne radiometry (LHR). This is part of a ground-based network of stations located across the planet (Wilson et al. 2017).
beyond performing the retrieval of planetary parameters, PSG can also correct for several typical issues affecting spectroscopic data.
The different panels show how these issues could affect the data, and their impact on the residuals. The methods employed to remove these
instrumental effects are described in Villanueva et al. 2013 (Icarus 222, 11–27).