Earth as an Exoplanet. III. Using Empirical Thermal Emission Spectra as an Input for Atmospheric Retrieval of an Earth-twin Exoplanet

From space, Earth's appearance is dominated by oceans, deserts, vegetation, ice, and clouds. Earth's surface is dominated by oceans (≈70% of surface), and the land-to-ocean ratio differs between the hemispheres (Northern Hemisphere ≈2/3, Southern Hemisphere ≈1/4; Pidwirny 2006). The contribution of different surface types and climate zones to a disk-integrated Earth spectrum (and its seasonal variability) depends on their thermal properties, their fractional contributions, and positions on the observed hemisphere.

In general, in the thermal emission spectrum of Earth, land-dominated views show not only higher flux readings but also larger flux variations over one full orbit than ocean-dominated views (e.g., Hearty et al. 2009; Gómez-Leal et al. 2012; Mettler et al. 2023). Specifically, from Table 2 in Mettler et al. (2023), we see that at Earth's peak emission wavelength (≈10.2 μm) the disk-integrated Northern Hemisphere pole-on view (NP) and the Africa-centered equatorial view (EqA) show annual flux variations of 33% and 22%, respectively. In contrast, the ocean-dominated Southern Hemisphere pole-on view (SP) and the Pacific-centered equatorial view (EqP) show smaller annual variations (≈11%) due to the large thermal inertia of oceans.

Another distinctive characteristic of Earth is its patchy cloud cover (see also Appendix A). Earth's patchy cloud coverage is unique among the three terrestrial planets with significant atmospheres in the Solar System (Venus is completely covered in clouds; Mars has negligible cloud coverage). Using nearly a decade of satellite data, King et al. (2013) show that roughly 67% of Earth's surface is covered by clouds at all times. The cloud fraction over land is approximately 55% and shows a distinct seasonal cycle. Over oceans, cloudiness is significantly higher (≈72%) and shows smaller seasonal variations. In addition, the cloud fraction is nearly identical during day and night, with only modest diurnal variation. Clouds are particularly abundant in the mid-latitudes (latitudes of ≈ ± 60°), and infrequent at latitudes from ±15° to ±30° (often characterized by arid desert conditions). Thus, there are three bands with a high cloud fraction in Earth's atmosphere: a narrow band at the equator and two wider mid-latitude bands.

Atmospheric clouds can significantly impact both the reflected light and thermal emission spectrum of a planet and can reduce or eliminate spectral features (particularly in the UV/O/NIR; e.g., Des Marais et al. 2002; Lu 2023). Parameters such as cloud fraction, composition, particle size, and altitude as well as multilayered cloud coverage and cloud seasonality all affect the resulting spectrum significantly (e.g., Des Marais et al. 2002; Tinetti et al. 2006a, 2006b; Hearty et al. 2009; Kitzmann et al. 2011; Rugheimer et al. 2013; Vasquez et al. 2013; Komacek et al. 2020). Konrad et al. (2023) ran retrievals on simulated MIR thermal emission spectra of a Venus-twin exoplanet. They showed that the presence of clouds can be inferred and requires a minimal spectral resolution of 50 and a signal-to-noise ratio (S/N) of 20. Further, clouds inhibit the accurate retrieval of surface conditions, and inadequate cloud treatment in retrievals (i.e., choosing too complex/simple cloud model given the quality of the input spectrum) can bias the estimates for important planetary parameters (e.g., planet radius, equilibrium temperature, and Bond albedo). However, despite recent efforts to understand how patchy clouds could alter the spectra of terrestrial exoplanets (e.g., May et al. 2021; Windsor et al. 2023), it remains unclear how they affect the characterization of terrestrial HZ exoplanets through MIR retrievals.

Habitability refers to the degree to which a global environment can support life, and depends on a myriad of factors (Meadows & Barnes 2018). The characteristics of a planet and its atmosphere, the architecture of the planetary system, the host star, and the galactic environment all affect habitability (for an extensive list, see, e.g., Meadows & Barnes 2018). For exoplanets, which can only be observed via remote sensing, we require observable characteristics to assess their habitability.

Analyzing MIR thermal emission spectra of exoplanets with atmospheric retrievals (see, e.g., Section 3; Madhusudhan 2018) and/or climate models yields constraints on the planet's atmospheric structure and composition. Such constraints yield valuable insights into a planet's habitability and could be used to infer the presence of a biosphere. In the following, we list observable signatures of habitability and biospheres in ascending order of difficulty to observe:

• 1.  
  Planetary energy budget. A planet's effective temperature and Bond albedo can be calculated from its thermal emission spectrum.
• 2.  
  Water and other molecules. Important atmospheric species, such as water (H₂O), carbon dioxide (CO₂), or ozone (O₃), have strong spectral MIR features.
• 3.  
  Atmospheric P–T structure. The P–T structure can be constrained in retrievals and provides vital information about the atmospheric state.
• 4.  
  Surface conditions. If not fully obscured by clouds, thermal emission spectra contain information about a planet's surface temperature and pressure.
• 5.  
  Molecular biosignatures. Important biogenic gases, such as methane (CH₄) or nitrous oxide (N₂O), have MIR features. The presence of a biosphere can be inferred if abiotic sources can be ruled out.
• 6.  
  Atmospheric seasonality. Seasonal periodicities in molecular abundances that are attributable to life are small for Earth (Mettler et al. 2023) and thus challenging to detect in the MIR. However, they could be strong indicators for biological activity (Olson et al. 2018).

For an in-depth review about evaluating planetary habitability and detectable signs of life, we refer to Schwieterman et al. (2018) and references therein.

In a previous study (Mettler et al. 2020), we analyzed 15 yr of thermal emission Earth observation data for five spatially resolved locations. We investigated flux levels and variations as a function of wavelength range and surface type (i.e., climate zone and surface thermal properties) and looked for periodic signals. From the spatially resolved single-surface-type measurements, we found that typically strong absorption bands from CO₂ (15 μm) and O₃ (9.65 μm) are significantly less pronounced and partially absent in polar regions. This implies that estimating correct abundance levels for these molecules might not be representative of the bulk abundances in these viewing geometries. Additionally, the time-resolved thermal emission spectrum provided insights into seasons/planetary obliquity, but its significance depended on the viewing geometry and spectral band.

In a follow-up study (Mettler et al. 2023), we expanded our analyses from spatially resolved locations to disk-integrated Earth views. We presented an exclusive data set consisting of 2690 disk-integrated MIR thermal emission spectra (3.75–15.4 μm, resolution R ≈ 1200). The spectra were derived from remote sensing observations for four different viewing geometries at a high temporal resolution. Using this data set, we investigated how Earth's MIR spectral appearance changes as a function of viewing geometry, seasons, and phase angles and quantified the atmospheric seasonality of different bioindicators. We found that a representative, disk-integrated thermal emission spectrum of Earth does not exist. Instead, both the thermal emission spectrum and the strength of biosignature absorption features show seasonal variability and depend strongly on viewing geometry.

In this paper, we treat Earth as a directly imaged exoplanet to assess the detectability of its characteristics from MIR observations with LIFE. For the first time, we perform a systematic retrieval analysis of real disk- and time-averaged Earth spectra. We investigate how the retrieval characterization depends on the viewing geometry and the season. Uniquely, in this study, we do have access not only to the real Earth spectra but also to ground truth data from remote sensing satellites. Hence, for the first time, we can compare retrieval results from real spectra to ground truth values. This allows us to evaluate the accuracy of the retrieved constraints and thereby validate our retrieval approach. Despite providing a unique opportunity to validate retrieval frameworks and their underlying assumptions, comparable retrieval studies on solar system observations are rare (e.g., Tinetti et al. 2006a; Lustig-Yaeger et al. 2023b; Robinson & Salvador 2023). However, such studies are indispensable to obtain a correct characterization of terrestrial exoplanets in the future.

In Section 2, we introduce the disk-integrated MIR thermal emission data set and the level 3 (L3) satellite products used to derive the ground truths. We introduce our atmospheric retrieval routine and the used atmospheric model in Section 3. In Sections 4 and 5, we present and discuss our retrieval results. We contextualize these results by discussing implications for characterizing terrestrial HZ exoplanets in Section 6. Finally, in Section 7, we summarize our findings and draw conclusions for future observations.

While there are several methods to study Earth from afar, such as Earth-shine measurements or spacecraft flybys (for a recent review, see, e.g., Robinson & Reinhard 2018, and references therein), we chose a remote sensing approach. This approach offers the extensive temporal, spatial, and spectral coverage needed to investigate the effect of observing geometries on disk-integrated thermal emission spectra and time-varying signals. However, for Earth-orbiting spacecraft, it is impossible to view the full disk of Earth, and the spatially resolved satellite data sets have to be combined into a spatially resolved, global map of Earth, which can then be disk-integrated (e.g., Tinetti et al. 2006a; Hearty et al. 2009; Gómez-Leal et al. 2012). Furthermore, due to the swath geometry of satellites, daily remote sensing data contain gores, which are regions with no data points, between orbit passes near the equator. In the case of Aqua/AIRS, these regions are filled within 48 hr as the satellite continues scanning Earth while orbiting it.

For our analysis, we defined four specific Earth observing geometries as shown in Figure 1: North (NP) and South Pole (SP), as well as Africa- (EqA) and Pacific-centered (EqP) equatorial views. For each viewing geometry, we mapped, calibrated, and geolocated radiances onto the globe and calculated the disk-integrated MIR thermal emission spectra. The spectra cover the 3.75–15.4 μm wavelength range (with a gap between 4.6 and 6.2 μm) at a nominal resolution of R ≈ 1200 and comprise 2378 spectral channels. The radiances originate from an AIRS infrared (IR) level 1C product (V6.7) called AIRICRAD ⁵ and are given in physical units of W m⁻² μm⁻¹ sr⁻¹ (Manning et al. 2019). The total data set contains 2690 disk-integrated thermal emission spectra for 4 consecutive yr (2016–2019) at a high temporal resolution for the four full-disk observing geometries (for an overview, see Table 1 in Mettler et al. 2023).

Zoom In Zoom Out Reset image size

The viewing geometries as portrayed in Figure 1 evolve throughout the year for a distant observer due to Earth's nonzero obliquity. Whereas the equatorial view blends seasons and has a diurnal cycle, the polar views show one season but blend day and night. Over the expected integration time of future direct imaging missions, the spectral appearance and characteristics of a planet change as it rotates around its spin axis and as spatial differences from clear and cloudy regions, contributions from different surface types as well as from different hemispheres vary with time. In accordance with the preliminary minimum LIFE requirements motivated in Konrad et al. (2022), we adopt a typical integration time of 30 days, which is significantly longer than Earth's rotation period. Hence, we average over the EqA and EqP views and denote the resulting data set EqC.

To capture the largest variability between observations, our analyses focus on observing Earth at its extremes in January and July. This choice is motivated by the measured relative flux change for these months at Earth's peaking wavelength in the disk-integrated thermal emission signal (Mettler et al. 2023). Although a Pacific-dominated view shows comparable variability to the South Pole view, Earth's rotation causes Africa and the Pacific to rotate in and out of the field of view. This rotation impacts the observed seasonal variability due to the different surface characteristics of EqA and EqP. The study details are summarized in Table 1.

The disk-integrated thermal emission spectra for this study are derived from our previously published data set. Since Earth's MIR spectrum exhibited negligible differences between consecutive years for a fixed viewing geometry (e.g., Mettler et al. 2020, 2023), we randomly chose the year 2017 and used the data of that year in order to calculate the monthly averages for January and July for the three viewing geometries: NP, SP, and EqC. Blending day and night data to simulate the phase of Earth at its orbital position was unnecessary for polar views due to Earth's obliquity, so they naturally include data of both types. However, in the case of the EqC view, we blended day and night data to simulate a rotating Earth at quadrature. This orbital position is preferred for the direct imaging of exoplanets due to the large apparent angular separation between the exoplanet and its host star.

AIRS spectra exhibit a gap between 4.6 and 6.2 μm due to dead instrument channels. This gap lies in a H₂O absorption feature centered at 6.2 μm (e.g., Catling et al. 2018). Due to concerns that the partially missing H₂O feature might deteriorate our retrieval results, we sourced level 1C data ⁶ for the year 2017 from the IASI instrument on board the MetOP satellite. We applied the same data reduction steps as for the AIRS data set described in Section 2.1 and Section 2 of Mettler et al. (2023). Covering the 3.62–15.50 μm wavelength regime with 8461 channels, IASI delivers a continuous spectrum comparable to that of AIRS, which makes it a suitable alternative instrument (see Figure 2). However, test retrievals showed no significant discrepancies between the retrieval results obtained for the gapped AIRS and continuous IASI spectra. The lack of discrepancies can be attributed to LIFE's noise level at these lower MIR wavelengths (e.g., Figure 4). Thus, since no significant differences were observed and the fact that our ground truth data introduced in Section 2.3 is based on Aqua/AIRS L3 monthly standard physical retrievals, we opted to use AIRS spectra for this study for consistency.

Zoom In Zoom Out Reset image size

In Section 4, we compare the retrieval outputs to an L3 satellite product comprising the P–T profile and the trace-gas abundances. Specifically, we have used the Aqua/AIRS L3 Monthly Standard Physical Retrieval (AIRS-only) 1° x 1° V7.0 (AIRS3STM) product (AIRS Project 2020), from which we extracted the surface temperature (land and sea surface) as well as the P–T profile. From the trace-gas parameters, we extracted the total integrated column burdens and vertical profiles (mass mixing ratios) of H₂O, CO, CH₄, and O₃. Both the P–T profile and trace-gas abundances are reported on 24 standard pressure levels ranging from 1000 to 1.0 hPa, which are roughly matched to the instrument's vertical resolution (Tian et al. 2020). The H₂O profile is an exception, as it is only provided at twelve layers ranging from 1000 to 100 hPa, spanning from the surface to the tropopause.

Since the AIRS3STM product did not contain any CO₂ abundances, we sourced the corresponding ground truth from a gridded monthly CO₂ assimilated data set ⁷ based on observations from the Orbiting Carbon Observatory 2 (OCO-2). The OCO-2 mission provides the highest quality space-based XCO2 retrievals to date, where the L3 data are produced by ingesting OCO-2 L2 retrievals every 6 hr with GEOS CoDAS, a modeling and data assimilation system maintained by NASAs Global Modeling and Assimilation Office (GMAO; NASA/GSFC/GMAO Carbon Group 2021). The data assimilation (or state estimation) technique is employed in order to estimate missing values based on the scientific understanding of Earth's carbon cycle and atmospheric transport. The missing values are mainly the result of the instrument's narrow 10 km ground track and limited ability to penetrate through clouds and dense aerosols.

Following the data reduction of the radiances in Section 2.1, the P–T profile and trace-gas abundances were mapped onto the globe for the different viewing geometries and then disk-integrated at each pressure level. For consistency, we also applied the empirical limb/weighting function to the ground truths. The uncertainties of the retrieved parameters from the AIRS L3 standard product were error propagated, and the resulting error bars are displayed for each data point. The results obtained for July and January are shown in Figure 3 and Appendix B, respectively.

Zoom In Zoom Out Reset image size

As input for our atmospheric retrievals, we use reduced-resolution versions of the disk-integrated ARIS spectra from Section 2.1 (NP, SP, and EqC viewing geometries for January and July). All spectra cover the 3.8–15.3 μm wavelength range, with a gap between 4.6 and 6.2 μm.

Based on the preliminary minimal LIFE requirements presented in Konrad et al. (2022, 2023), Alei et al. (2022) (R = 50, S/N = 10), we consider two resolution cases (R = 50, 100) and two S/Ns (S/N = 10, 20) for each of the six disk-integrated spectra. We define R as λ/Δλ, with the width of a wavelength bin Δλ and the wavelength at the bin center λ. Further, the S/N value corresponds to the S/N in the 11.2 μm wavelength bin. We choose the 11.2 μm bin because it does not coincide with any strong spectral features. In Figure 4, we show the six R = 50 input spectra together with the two different noise levels.

Zoom In Zoom Out Reset image size

We model the wavelength-dependent S/N expected for LIFE with LIFEsim (Dannert et al. 2022), which accounts for astrophysical noise sources (photon noise of planet emission, stellar leakage, and local dust emission as well as exozodiacal dust emission). ⁸ To estimate the LIFEsim noise, we put Earth on a 1 au orbit around a G2V star located 10 pc from the observer. The exozodiacal dust emission of the system was assumed to reach 3 times the local zodiacal level. ⁹

In our retrievals, we interpret the noise as uncertainty to the points of the disk-integrated spectra. Thus, the spectral points correspond to the true flux values and are not randomized according to the LIFEsim S/N. While randomized spectra would provide a more accurate simulated observation, a retrieval study based on a single noise realization will yield biased parameter estimates. Ideally, we would run retrievals for multiple (≳10) noise realizations of each spectrum. However, the number of retrievals required makes such a study computationally unfeasible. Yet, Konrad et al. (2022) motivate that the results from retrievals on unrandomized spectra provide reliable estimates for the average expected retrieval performance on randomized spectra.

For this study, we utilized the Bayesian retrieval routine introduced in Konrad et al. (2022). The initial routine was improved and modified in Alei et al. (2022), Konrad et al. (2023). We provide a brief summary of the routine here, and refer to the original publications for an in-depth description.

Our retrieval framework uses the radiative transfer code petitRADTRANS (Mollière et al. 2019, 2020; Alei et al. 2022) to calculate the theoretical emission spectrum of a 1D plane-parallel atmosphere model. petitRADTRANS assumes a blackbody (BB) spectrum at the surface and models the interaction of each atmospheric layer with the radiation to calculate the spectrum at the top of the atmosphere. The model atmosphere is defined via a set of forward model parameters (see Section 3.3 for our forward model). In a retrieval, we search the space spanned by the prior probability distributions (or priors) of the forward model parameters for the parameter combination that best reproduces the input spectrum. To efficiently search the prior volume, we use the pyMultiNest (Buchner et al. 2014) package, which uses the MultiNest (Feroz et al. 2009) implementation of the nested sampling algorithm (Skilling 2006). Here, we ran all retrievals using 700 live points and a sampling efficiency of 0.3. ¹⁰

The retrieval yields the posterior probability distribution (or posterior) for the model parameters. The posterior estimates how likely a certain combination of model parameter values is given the observed spectrum. Further, our routine estimates the Bayesian evidence ${ \mathcal Z }$, which is a measure for how well the used forward model fits the input spectrum and can be used for model comparison (see Appendix C).

As in Konrad et al. (2022, 2023), Alei et al. (2022), we characterize each layer of the model atmosphere by its temperature, pressure, and the opacity sources present. We provide a list of all model parameters in Table 2. A comparison between different forward models to justify our choice is provided in Appendix C.

Table 2. Parameters of the Retrieval Forward Model

Parameter	Description	Prior
a4	P–T parameter (degree 4)	${ \mathcal U }(0,10)$
a3	P–T parameter (degree 3)	${ \mathcal U }(0,100)$
a2	P–T parameter (degree 2)	${ \mathcal U }(0,500)$
a1	P–T parameter (degree 1)	${ \mathcal U }(0,500)$
a0	P–T parameter (degree 0)	${ \mathcal U }(0,1000)$
${\mathrm{log}}_{10}({P}_{0})$	${\mathrm{log}}_{10}(\mathrm{Surfacepressure}\left[\mathrm{bar}\right])$	${ \mathcal U }(-4,2)$
Rpl	Planet radius $\left[{R}_{\oplus }\right]$	${ \mathcal G }(1.0,0.2)$
${\mathrm{log}}_{10}({M}_{\mathrm{pl}})$	${\mathrm{log}}_{10}(\mathrm{Planetmass}\left[{{\rm{M}}}_{\oplus }\right])$	${ \mathcal G }(0.0,0.4)$
${\mathrm{log}}_{10}({{\rm{N}}}_{2})$	${\mathrm{log}}_{10}({{\rm{N}}}_{2})$ mass fraction	${ \mathcal U }(-10,0)$
${\mathrm{log}}_{10}({{\rm{O}}}_{2})$	${\mathrm{log}}_{10}({{\rm{O}}}_{2})$ mass fraction	${ \mathcal U }(-10,0)$
${\mathrm{log}}_{10}({\mathrm{CO}}_{2})$	${\mathrm{log}}_{10}({\mathrm{CO}}_{2})$ mass fraction	${ \mathcal U }(-10,0)$
${\mathrm{log}}_{10}({{\rm{H}}}_{2}{\rm{O}})$	${\mathrm{log}}_{10}({{\rm{H}}}_{2})$ O mass fraction	${ \mathcal U }(-10,0)$
${\mathrm{log}}_{10}({{\rm{O}}}_{3})$	${\mathrm{log}}_{10}({{\rm{O}}}_{3})$ mass fraction	${ \mathcal U }(-10,0)$
${\mathrm{log}}_{10}({\mathrm{CH}}_{4})$	${\mathrm{log}}_{10}({\mathrm{CH}}_{4})$ mass fraction	${ \mathcal U }(-10,0)$

Note. The third column lists the priors assumed in the retrievals. We denote a boxcar prior with lower threshold x and upper threshold y as ${ \mathcal U }(x,y);$ for a Gaussian prior with mean μ and standard deviation σ, we write ${ \mathcal G }(\mu ,\sigma )$.

Download table as:  ASCIITypeset image

In our forward model, we parameterized the atmospheric P–T profile using a fourth-order polynomial:

$$\begin{eqnarray}&&T(P)=\displaystyle \sum _{i=0}^{4}{a}_{i}{P}^{i}.\end{eqnarray} \tag{ 1 }$$

Here, P is the pressure, T is the corresponding temperature, and the a_i is the parameters of the P–T model. As shown in Konrad et al. (2022), a polynomial P–T model allows us to minimize the number of P–T parameters and thereby minimize the retrieval's computational complexity. Learning based P–T models require fewer parameters, but their accuracy for terrestrial planets is currently limited by the availability of sufficient training data (e.g., Gebhard et al. 2023).

We consider various opacities in our forward model. First, we account for the MIR absorption and emission by CO₂, H₂O, O₃, and CH₄ (see Table 3 for line lists, broadening coefficients, and cutoffs). We assume constant vertical abundance profiles for all molecules and discuss potential effects of this simplification in Section 5. Second, we model collision-induced absorption (CIA) and Rayleigh scattering features (CIA-pairs and Rayleigh-species are listed in Table 3).

We neglect scattering and absorption by clouds. The patchy clouds in Earth's atmosphere partially block contributions from high-pressure atmosphere layers and thereby impede the characterization thereof. Konrad et al. (2023) show that neglecting clouds in retrievals can lead to systematic errors in the retrieved surface temperature, surface pressure, and the planet radius. We provide a detailed discussion on potential effects of this simplification in Section 5.

We list the priors assumed for all retrievals in Table 2. The priors on the P–T parameters a_i and the surface pressure P₀ cover a wide range of atmospheric structures (from tenuous Mars-like to thick Venus-like atmospheres). For N₂, O₂, CO₂, H₂O, O₃, and CH₄, we select broad uniform priors that extend significantly below the minimal detectable abundances estimated in Konrad et al. (2022) (≈10⁻⁷ in mass fraction for our R and S/N cases).

As in Konrad et al. (2022, 2023), Alei et al. (2022), we use Gaussian priors for the planet radius R_pl and mass M_pl. The R_pl prior is based on Dannert et al. (2022), who suggest that a planet detection with LIFE yields a constraint on R_pl. ¹¹ The statistical mass–radius relation Forecaster ¹² (Chen & Kipping 2016) is then used to infer the prior on ${\mathrm{log}}_{10}({M}_{\mathrm{pl}})$ from the R_pl prior.

As discussed above, our estimates for the trace-gas abundances are strongly affected by a degeneracy with P₀ and the P–T structure. Further, we expect the trace-gas posteriors to be impacted by a physical degeneracy with the planet's surface gravity g_pl and thus M_pl (see, e.g., Mollière et al. 2015; Feng et al. 2018; Madhusudhan 2018; Alei et al. 2022; Konrad et al. 2022, 2023). ¹³ If the retrieved abundance posteriors of two different trace-gases are affected by these degeneracies in the same way, the biases in our retrieval results can be largely eliminated by considering their point-wise ratio (i.e., divide one posterior by another). Despite not providing information on the absolute trace-gas abundances, such ratios are of interest since they can help identify states of atmospheric chemical disequilibrium, which can indicate biological activity (see Section 6.2 for an extended discussion; Lovelock 1965,1975).

We present the relative abundance posteriors for all trace-gas combinations in Figure 7 (numerical values in Tables E1 to E3). The uncertainties on the relative trace-gas abundances are significantly smaller than on the absolute abundances due to the elimination of the aforementioned M_pl degeneracy. Further, in contrast to the absolute abundance posteriors (Figure 6), all ratios lie within the range of the relative ground truths, indicating that the biases invoked by the P₀ degeneracy are mostly eliminated.

Zoom In Zoom Out Reset image size

Previous studies have evaluated how well LIFE could characterize terrestrial HZ exoplanets. Konrad et al. (2022) find preliminary estimates for LIFE's minimal R and S/N requirements by running retrievals on simulated Earth spectra. Alei et al. (2022) run retrievals on simulated spectra from Rugheimer & Kaltenegger (2018), which represent different stages in Earth's temporal evolution. Konrad et al. (2023) investigate how well LIFE can characterize the atmosphere and clouds of Venus. All studies analyze spectra that were calculated with simplified, temporally constant, 1D atmosphere models. In contrast, we run retrievals on real disk- and time-averaged MIR Earth spectra.

Despite large differences in the complexity of the considered spectra between our study and the previous ones, we confirm the previous findings for the detectability of different trace-gases. Crucially, CH₄, the main driver for the minimum LIFE requirements from Konrad et al. (2022) (R = 50, S/N = 10), remains detectable here. Further, the strengths of the parameter constraints we retrieve here are equivalent to those from the prior studies, which demonstrates their robustness. The retrieved 1σ parameter uncertainties in all studies are < ± 0.5 dex for pressures, < ± 0.1R_⊕ for radii, < ±20 K for temperatures, and < ±1.0 dex for trace-gas abundances.

In Section 4, we state that our R_pl estimates underestimate Earth's true radius. This leads to biased estimates for T_eq and A_B, which are calculated from R_pl (see Appendix D).

We mainly attribute underestimation of R_pl to neglecting Earth's patchy cloud coverage in our forward model (see Section 3.3). Clouds reduce the total MIR emission at the top of the atmosphere by partially absorbing the thermal emission from the warm, high-pressure atmosphere layers below them. By using a first-order approximation (see Appendix G for details), we can demonstrate that the magnitude of the bias on our R_pl estimate can be fully attributed to the missing cloud treatment in our forward model.

Further evidence for links between clouds and biased R_pl estimates is provided by other thermal emission retrieval studies. Konrad et al. (2022), who run retrievals on simulated cloud-free Earth spectra, retrieve bias-free R_pl estimates. In contrast, Alei et al. (2022) run cloud-free retrievals on simulated cloudy Earth spectra and also underestimate R_pl.

As stated in Section 4, the retrieved P₀ and P–T structures are offset to lower pressures relative to the ground truth. These biases are linked to offsets in the retrieved trace-gas abundances, which are degenerate with the pressure-induced line-broadening (see, e.g., Misra et al. 2014; Schwieterman et al. 2015).

We attribute these biases to our assumption of vertically constant trace-gas abundances in our forward model (see Section 3.3). This claim is motivated by comparison with previous LIFE retrieval studies. Konrad et al. (2022) assume constant abundance profiles both to generate their 1D Earth spectra and in their forward model, and retrieve unbiased P₀, P–T, and abundance estimates. In contrast, Alei et al. (2022) assume constant abundance profiles to run retrievals on 1D Earth spectra from Rugheimer & Kaltenegger (2018), which were generated using nonconstant abundance profiles. Their results show offsets in P₀, the P–T structure, and the trace-gas abundances, which are comparable in magnitude to our offsets.

Further, we argue that H₂O is the cause of the observed biases. First, in contrast to CO₂, O₃, and CH₄, H₂O has multiple strong absorption features in Earth's MIR spectrum (see, e.g., Figure 3 in Konrad et al. 2022). Second, the ground truths in Figure 3 show that the variances for H₂O are more than 2 orders of magnitude greater than for the other species. Third, the main H₂O variance occurs in the lowermost atmosphere layers, where H₂O condensation occurs. These layers contribute most strongly to Earth's MIR thermal emission.

In the present study, the ground truth measurements of Earth's atmosphere have allowed us to validate our results. We found important biases in the posteriors, which we attribute to simplifying assumptions made by our forward model. A proposed remedy is to derive quantities that are less affected, such as abundance ratios (see Section 4.1). Also, in a future study, we aim to reduce biases by adding a parameterization for patchy clouds and a vertically nonconstant H₂O profile (motivated by H₂O condensation) to our forward model.

Independent of the success of this future effort, intercomparison efforts (e.g., Barstow et al. 2020) have shown that retrieval results also strongly depend on framework specificities (e.g., parameter estimation algorithms, radiative transfer implementations, and line-lists). To ensure the correct characterization of exoplanets, robust and bias-free retrieval frameworks are required. Thus, community efforts, such as the CUISINES Working Group, ¹⁴ that benchmark, compare, and validate different frameworks on real and simulated spectra with known ground truths are indispensable.

As described in Section 1.1, Earth exhibits an uneven distribution of land and ocean regions. Further, different surface types have different spectral and thermal characteristics (e.g., Hearty et al. 2009; Gómez-Leal et al. 2012; Madden & Kaltenegger 2020). Also, the distribution of life on Earth is nonuniform with a measurable gradient in the abundance and diversity of life, both spatially (e.g., from deserts to rain forests) and temporally (e.g., from seasonal to geological timescales; Méndez et al. 2021). In Mettler et al. (2023), we find that a representative, disk-integrated thermal emission spectrum of Earth does not exist. Instead, the MIR spectrum and the strength of the absorption features show seasonal variations and depend on the viewing geometry. For future observations of HZ terrestrial exoplanets, the viewing geometry will be unknown. Thus, we must understand how the viewing geometry impacts exoplanet characterization, observable habitability markers, and signatures of life.

As we see from Figure 6, most parameter posteriors show no significant dependence on either the viewing geometry or the season (exceptions: T₀, T_eq, and A_B). For the R and S/N levels studied here, both the retrieved R_pl and the trace-gas abundance estimates show no measurable variations with the exoplanet's orientation relative to the observer. Thus, we conclude that their characterization depends on neither the viewing geometry nor the season for an Earth-like exoplanet.

For T₀ and T_eq, the variations in the posteriors are largest for the NP view, where the differences between January and July are robustly detected in all R and S/N scenarios. For the SP and EqC views, the variations in T₀ and T_eq between January and July are much smaller and not confidently detected. This is in agreement with those from Mettler et al. (2023), who find the seasonal disk-integrated thermal emission flux differences for the landmass-dominated NP view to be 33%, as opposed to only 11% for the ocean-dominated SP and EqP (EqP) views. The increased T₀ variance observed for the NP view can be attributed to the large landmass fraction. Also, our results for T₀ and T_eq indicate that the NP Jul, SP, and EqC views cannot be differentiated from one another despite vastly different characteristics like climate zones and landmass fractions. This highlights the strong spectral degeneracy with respect to seasons and viewing geometries, and agrees with other studies (e.g., Gómez-Leal et al. 2012; Mettler et al. 2023).

For A_B, our retrieval results show small differences between the viewing angles and seasons. As for the temperatures, the variations are largest for the NP view (NP, 46%; SP, 17%; EqC, 16%). For the NP and SP views, the retrieved A_B tends to be higher during winter, which agrees with the lower retrieved T₀ and T_eq values. However, due to the uncertainties (±0.05 to ±0.10) and biases on the posteriors, a confident detection of these A_B differences is not possible. Thus, our A_B characterization is independent of viewing geometry and season.

However, as we discuss in Section 4, the accuracy and strength of our constraints for T₀, T_eq, and A_B are limited by our R_pl estimates. As we demonstrate in Appendix F, an accurate and strong R_pl constraint would yield detectable differences in T₀, T_eq, and A_B. In this case, Earth-like seasonal T₀, T_eq, and A_B changes are easily detectable for the NP view with a LIFE-like observatory for all R and S/N cases. Also, for the SP view, detections of seasonal variations are possible (except for the R = 50, S/N = 10 case). For the EqC view, which blends the two hemispheres, the seasonal variations remain undetected.

Earth's MIR spectrum contains features from numerous bioindicator gases. Examples are O₃ (photochemical product of bioindicator O₂), CH₄, and N₂O (see, e.g., for an extensive list, Schwieterman et al. 2018). While N₂O is not detectable at the R and S/N considered, O₃ and CH₄ are (biases ≤ + 1.0 dex, uncertainties ≤ ± 1.0 dex; see Appendix C). However, the sole detection of a bioindicator gases is not sufficient to infer the presence of life, since they can be produced abiotically (see, e.g., Catling et al. 2018; Harman & Domagal-Goldman 2018; Schwieterman et al. 2018). The simultaneous detection of multiple bioindicator gases provides a more robust marker for biological activity.

One detectable multiple bioindicator is the "triple fingerprint," which is the simultaneous detection of atmospheric CO₂, H₂O, and O₃ (see, e.g., Selsis et al. 2002). Depending on the viewing angle, the retrieved ${\mathrm{log}}_{10}({{\rm{H}}}_{2}{\rm{O}}/{\mathrm{CO}}_{2})$ ranges from −0.2 to −0.5 dex, and the ${\mathrm{log}}_{10}({{\rm{O}}}_{3}/{{\rm{H}}}_{2}{\rm{O}})$ ranges from −2.6 to −3.2 dex (uncertainties ≤ ± 0.6 dex). These values lie between the average 1 and 0.1 bar ground truths (${\mathrm{log}}_{10}({{\rm{H}}}_{2}{\rm{O}}/{\mathrm{CO}}_{2})$: 1.2 dex at 1 bar, −2.3 dex at 0.1 bar; ${\mathrm{log}}_{10}({{\rm{O}}}_{3}/{{\rm{H}}}_{2}{\rm{O}})$: −0.3 dex at 1 bar, −5.2 dex at 0.1 bar).

Another promising multiple bioindicator is the simultaneous detection of reducing and oxidizing species in an atmosphere (i.e., a strong chemical disequilibrium). Since the two species will react rapidly with each other, simultaneous presence over large timescales is only possible if both are continually replenished at a high rate by life (Lederberg 1965). One example hereof that we confidently detect is the simultaneous presence of O₂ (or its photochemical product O₃) ¹⁵ and CH₄ (Lovelock 1965; Lippincott et al. 1967). For all but the R = 50, S/N = 10 retrievals, we accurately constrain the ${\mathrm{log}}_{10}({\mathrm{CH}}_{4}/{{\rm{O}}}_{3})$ abundance ratio to 1.1 dex (uncertainty ≤ ± 0.5 dex). In particular, in the context of an Earth-like planet orbiting a Sun-like star, the detection of such an O₂/O₃–CH₄ disequilibrium would represent a strong potential biosignature.

Research on the detectability of exoplanet biosignatures has predominantly focused on a static evidence for life (e.g., the coexistence of O₂ and CH₄). However, the anticipated range of terrestrial planet atmospheres and the potential for both "false positives" and "false negatives" in conventional biosignatures (e.g., Selsis 2002; Meadows 2006; Reinhard et al. 2017; Catling et al. 2018; Krissansen-Totton et al. 2022) underscore the necessity to explore additional life detection strategies. Time-varying signals, such as seasonal variations in atmospheric composition, have been proposed to be strong biosignatures (e.g., Olson et al. 2018), since they are biologically modulated phenomena that arise naturally on Earth and likely also occur on other nonzero obliquity and eccentricity planets. Olson et al. (2018) suggest that atmospheric seasonality as a biosignature avoids many assumptions about specificities of metabolisms. Further, it offers a direct means to quantify biological fluxes, which would allow us to characterize, rather than simply identify, exoplanet biospheres.

To assess the detectability of such time-dependent atmospheric modulations in exoplanets, we consider the retrieved abundance ratios in Figure 7. Abundance ratios are less affected by parameter degeneracies and thus exhibit smaller uncertainties and biases (≤ ± 0.4 dex for the R = 50, S/N = 20, and R = 100 cases). Independent of the viewing geometry, we see no significant differences between the trace-gas ratios retrieved for January and July. Since these months represent Earth's two extreme states, we do not expect differences in the trace-gas ratios to be observable for any two other months. Consequentially, detecting the atmospheric seasonality of trace-gas abundances as a biosignature is not feasible for the studied R and S/N cases.

This agrees with our findings in Mettler et al. (2023), where we studied disk-integrated Earth spectra and quantified the amplitudes of the seasonal variations in absorption strength by measuring the equivalent widths of the biosignature-related absorption features. We detected small seasonal variations for O₃, CO₂, CH₄, and N₂O. For CO₂ and CH₄, the seasonal abundance variations of 1% to 3% are significantly smaller than the uncertainties on our retrieved abundance estimates (≈ ± 0.5 dex), which makes a detection unfeasible.

Significantly higher R or S/N MIR spectra are required to be sensitive to the spectral variations evoked by Earth-like seasonal fluctuations in bioindicator gas abundances. Such observations require either a more sensitive instrument or an integration time greater than the assumed 30 days (see Section 2.1). However, while the magnitude of such spectral variations is unchanged for shorter integration times ¹⁶ (e.g., 10 days), it will decrease for extended integration times (e.g., 90 days). For Earth, significant seasonal changes occur during such extended observations. However, the measured spectrum represents the average state of the observed atmosphere. Thus, the magnitude of the spectral variations evoked by seasonal fluctuations is diminished, which counteracts the sensitivity gain attained via an increase in observation time.

However, terrestrial exoplanets could display seasonality patterns that are very different from that of Earth or other Solar System planets. Given the extensive diversity among exoplanets (e.g., in terms of mass, size, host star type, and orbit), it is likely that some exhibit detectable seasonal variations. Seasonal signals could be amplified by several factors (see, e.g., Section 4.3 in Mettler et al. 2023) such as the following: shorter photochemical lifetimes and/or nonsaturated spectral bands, increased orbital obliquity (leads to greater seasonal contrast due to varying ice and vegetation cover), biological activity promoted by moderately high obliquity (e.g., photosynthetic activity) consequently leading to heightened variations in biosignature gases, and the absence of competing effects from admixed hemispheres (particularly relevant for eccentric planets). The detectability of seasonality depends on both the magnitude of the biogenic signal and the degree to which the observation conditions mute that signal, and is likely maximized for an intermediate obliquity.

In Section F, we reduced the retrieved posterior distributions to fixed values of P₀ = 1 bar, and R_pl = 1R_⊕ by assuming a linear correlation between P₀ or R_pl and the remaining posteriors. Here, we outline the method used to reduce the posterior distributions. A schematic illustration showing both the true and the reduced posteriors can be found in Figure F1.

Zoom In Zoom Out Reset image size

In the following, let us consider one point in the retrieved posterior distribution. We denote the value of P₀ or R_pl (i.e., the parameter we want to reduce over) as θ_red,true and values of the other model parameters as θ_param,true. If we assume a linear correlation between θ_red,true and θ_param,true, we can make a prediction θ_param,pred for θ_param,true using θ_red,true as follows:

$$\begin{eqnarray}&&{\theta }_{\mathrm{param},\mathrm{pred}}=m\cdot {\theta }_{\mathrm{red},\mathrm{true}}+q.\end{eqnarray} \tag{ F1 }$$

Here, m is the slope, and q is the offset with respect to the origin of the linear model. θ_param,pred is the parameter value predicted by the linear model. We search for the best-fit linear model by minimizing the square difference Δ between θ_param,pred and θ_param,true:

$$\begin{eqnarray}&&\begin{array}{l}{\rm{\Delta }}=\displaystyle \sum _{\mathrm{Posterior}}{\left({\theta }_{\mathrm{param},\mathrm{pred}}-{\theta }_{\mathrm{param},\mathrm{true}}\right)}^{2}\\ \,=\,\displaystyle \sum _{\mathrm{Posterior}}{\left(m\cdot {\theta }_{\mathrm{red},\mathrm{true}}+q-{\theta }_{\mathrm{param},\mathrm{true}}\right)}^{2}.\end{array}\end{eqnarray} \tag{ F2 }$$

The best-fit linear models are indicated in Figure F1 as black dashed–dotted lines. From the figure, we see that the correlations between the parameters considered here are well described by our linear model. In the next step, we fix the value of θ_red,true to θ_red,fix and calculate the corresponding reduced posterior values θ_param,red of the other parameters as follows:

$$\begin{eqnarray}&&{\theta }_{\mathrm{param},\mathrm{red}}={\theta }_{\mathrm{param},\mathrm{true}}+m\cdot \left({\theta }_{\mathrm{red},\mathrm{fix}}-{\theta }_{\mathrm{red},\mathrm{true}}\right).\end{eqnarray} \tag{ F3 }$$

This yields the reduced posterior distribution of a parameter, which we plot in Figure F1. This reduction method allows us to remove the effect of one parameter on the posterior distribution, and identify the origin of biases in the retrieval results.

To demonstrate the effects of underestimating P₀, we reduce the abundance posteriors to Earth's true P₀ of 1 bar. The resulting reduced posteriors are shown in the left panel of Figure F2 (numerical values in Tables E1 to E3). The posterior reduction to the true P₀, leads to significantly better estimates for CO₂ and CH₄. For CO₂, the reduced posteriors are perfectly centered on the true value, while the CH₄ abundances are significantly less overestimated. This demonstrates that the shifts in the CO₂ and CH₄ posteriors in Figure 6 are directly linked to the inaccurately retrieved P₀. For O₃ and H₂O, the reduced posteriors are shifted to lower abundances and show a smaller variance between the individual retrievals. These findings suggest that the P₀ reduction also yields improved estimates for atmospheric O₃ and H₂O abundances.

Zoom In Zoom Out Reset image size

In order to investigate how our underestimation of R_pl affects the other posteriors, we reduce the R_pl posterior to 1 R_⊕. We plot the reduced posteriors of P₀, T₀, T_eq, and A_B in the right panel of Figure F2 (numerical values in Tables E1 to E3). First, we observe no significant differences between the reduced and the true P₀ posteriors from Figure 6. Thus, no direct correlation between the R_pl and the P₀ posterior exists. Second, for T₀, which is accurately estimated in all retrievals, the reduced posteriors underestimate the truth by ≥5 K. This finding indicates that Earth's disk-integrated flux is smaller than what is expected for a cloud-free 1 R_⊕ planet with surface temperature T₀. This suggests that patchy clouds, which partially block the emission from the high-pressure atmospheric layers and thereby reduce the total planet flux, are the likely cause of the R_pl biases (see also Appendix G for further evidence). Finally, the reduced T_eq and A_B posteriors are unbiased and provide accurate truth estimates (uncertainties, ≤ ± 2 K for T_eq; ≤ ± 0.1 for A_B), which demonstrates the correlations with R_pl.