Terran World Spectral Simulator

The planet specification follows from a set of "tiles", each representing a spatial section of the planet surface. This tileset represents the discretization of the exoplanet into planar facets, allowing a plane-parallel RTM to represent the light scattering from a combination of atmosphere and surface. The software stores the tiles as a one-dimensional list. In principle any discretization of the sphere could be used, with the only requirement that the observation angles and the tile areas must be readily computed.

Each tile must then be described as a fractional mixture of possible surface and atmosphere types. While there is no limit imposed on the number of types, in practice these must be kept small (less than 30) otherwise the number of pre-computed RTM outputs becomes large. The different surface types would be selected based on the surface composition of the exoplanet. The typical atmosphere specification would include one or more cloud types, with a fixed molecular composition and thermodynamic profile (meaning the profiles can vary in the vertical dimension, but are constant in the horizontal dimension). In principle, spatial variations in the atmospheric composition could be included with additional classes. For example, slowly varying composition gradients (e.g., water vapor in Earth's atmosphere) may be approximated with a small number of discrete atmospheric composition classes.

The RTM pre-computes the spectral reflectance associated with the surface and atmosphere types used in the planet specification. The required spectral reflectance quantity is defined as the bidirectional reflectance factor (BRF) of the surface and atmosphere (Schaepman-Strub et al. 2006). This quantity is the ratio between the simulated reflected radiance to the reflected radiance from a purely reflective Lambertian surface at the same orientation. The BRF is often called "reflectance" for simplicity, and uses the symbol ρ. This can be expressed as a relationship between the reflected radiance (L_r) and the incident irradiance (E_i):

$$\begin{eqnarray}&&\rho ({\theta }_{i},{\theta }_{r},{\phi }_{i},{\phi }_{r})=\displaystyle \frac{{L}_{r}}{{L}_{\mathrm{Lambertian}}}=\displaystyle \frac{{L}_{r}}{\cos {\theta }_{r}{E}_{i}/\pi },\end{eqnarray} \tag{ 1 }$$

where the reflectance is a function of the incident zenith and azimuth angles (θ_i and ϕ_i, respectively) and reflected angles (θ_r and ϕ_r). The spectral dependence (either per wavelength or per frequency) is implied, and omitted for brevity. Note that the reflectance will generally not be Lambertian even if the surface is assumed to be, due to the effect of atmospheric scattering and absorption.

The RTM must be run individually for each combination of surface and atmosphere type referenced by the planet specifications. For example, a simple desert world, with a homogeneous surface and one type of cloud, would require two RTM model output databases: one for the clear condition and one for the cloudy condition, over the spectral reflectance appropriate for a desert surface. Each of these model databases requires a set of spectral reflectances that span the range of possible observation angles. For a Lambertian surface model the system is azimuthally symmetric, so only three angles are needed to describe the geometry between the star and observer: the incident zenith angle from the star (θ_i), the reflection zenith angle to the observer (θ_r), and the relative azimuth angle (Δϕ = ϕ_i − ϕ_r). Complicated non-Lambertian surfaces can have azimuthally asymmetric surface reflectances. Examples are geological land forms such as parallel ridges. These are likely not strong contributions to the global disk-integrated reflectance, so our framework currently assumes the reflectance database is azimuthally symmetric and only requires the three angles to describe the radiation geometry. The pre-computed RTM outputs are stored in database files for later use. The primary output data are the four-dimensional reflectance array and the grid coordinates for each dimension.

The combination of the planet specification and the RTM output database is combined in a python object class, exoplanet. The main method exposed by the exoplanet class computes the per-tile spectral reflectance for a given set of per-tile observation angles. In addition, several visualization tools are included, which render maps of the exoplanet surfaces under a fixed illumination to each tile. This is not intended to simulate true stellar illumination (since it illuminates the entire surface), but rather to visualize the full global map of the exoplanet.

The planet tile specifications and the associated RTM output databases represent the pre-computed information needed for the reflectance integration over the observed exoplanet disk. For a particular observation geometry we must: (1) determine the per-tile observation geometry, (2) interpolate the RTM database to the per-tile observation angles, and (3) integrate over all tiles visible to the observer, applying weighting for each tile's projected area. Figure 1 shows a schematic view of the observation angle interpolation at a particular tile. The final disk-integrated reflectance, R can be expressed as

$$\begin{eqnarray}&&R=\displaystyle \frac{{\displaystyle \sum }_{k=1}^{N}{\rho }_{i}({\theta }_{i,k},{\theta }_{r,k},{\rm{\Delta }}{\phi }_{i,k}){T}_{k}\cos ({\theta }_{r,k}){V}_{k}}{{\displaystyle \sum }_{k=1}^{N}{T}_{k}},\end{eqnarray} \tag{ 2 }$$

where k is the index over the set of N exoplanet tiles, T_k represents the tile area, and V_k is a visibility function (with a value of 1 when the tile is illuminated and visible to the observer, and 0 if not). In this definition, R is equivalent to the product of the exoplanet's geometric albedo, A, and the planetary phase function, Φ(α). The planetary phase angle (α) does not appear directly in Equation (2), as all the relevant scattering angles must be defined relative to individual tiles' coordinates. Finally, note that both the pre-computed RTM database and the calculations used in Equation (2) refer to unpolarized reflectances. In principle, these calculations could be done in a similar way with polarized reflectance (via the Stokes parameters), but the current implementation assumes these are unpolarized ρ and R.

Zoom In Zoom Out Reset image size

To solve for the observation geometry, an exosystem object class is implemented which combines the exoplanet object with an orbit object, which contains an orbital geometry solver. The orbit object contains the classical orbital parameters that define the orbit shape and orientation in three-dimensional space (eccentricity, argument of periastron, longitude of ascending node, inclination angle) as well as the exoplanet's rotational axis (obliquity and longitude of the equinox). A simple iterative solution is used to solve for the mean anomaly as a function of orbital time (Prussing & Conway 1993). The orbit object can also compute the phase angle for the given observation. Since the exosystem class has the full information available to compute the observation geometry, it has visualization functions based on the Python MayaVI package to render a three-dimensional image of the viewed exoplanet. Figure 2 shows an example color rendering of a cloud-free Earth model at an arbitrary observation geometry.

Zoom In Zoom Out Reset image size

To illustrate the TSS we collapse the full range of Earth surface conditions into four types: ocean, snow, sand, and vegetation; see Table 1. These four types use Lambertian surface spectral reflectances from the Advanced Spaceborne Thermal Emission and Reflection Radiometer (ASTER) database (Baldridge et al. 2009), and represent the primary surface types that affect the global reflectance signature of the Earth. The surface spectral reflectances for these four types are shown in Figure 3. Several important properties are apparent: the "red edge" in the vegetation (grass) spectrum around 700 nm, the very low reflectance of water at all wavelengths, the relatively smooth reflectance for sand, the very high reflectance for snow at short wavelengths (<1000 nm), and the low reflectances in the snow and vegetation spectra at 1500 nm and 2000 nm due to water absorption. Since we are using the unpolarized library spectra from ASTER, the radiative transfer calculations will not include the highly non-Lambertian scattering from ocean surfaces that produces the "glint" signature.

Zoom In Zoom Out Reset image size

For the atmosphere specification, we use a simple atmosphere model which has a globally constant vertical profile of temperature and molecular composition. This fixed vertical profile is the US Standard atmosphere profile, and we include the concentration profiles of seven important absorbing molecular species: H₂O, CO₂, O₃, N₂O, CO, CH₄, and O₂. This atmosphere specification is then applied to three conditions: a clear atmopshere, a low-altitude water cloud, and a high-altitude ice cloud. We do not use an overlapping cloud type that includes both water and ice clouds. These atmosphere types are then combined with each of four surface types: ocean, snow, sand, and vegetation. The combination of three atmosphere and four surface types yields a total of 12 possible combinations that can exist within each tile.

For clouds, each of the two classes (low-altitude water cloud and high-altitude ice cloud) has a fixed optical optical thickness, altitude, and particle effective diameter (see Table 2 for a summary). The water clouds use single-scattering properties from Mie calculations and a particle size distribution from the gamma distribution form used in Hansen & Travis (1974). The ice clouds use single-scattering properties from the general crystal habit mixture in Baum et al. (2014). The two effective diameters are intended to be reasonable global average values, in the range of remotely sensed particle sizes from contemporary Earth observation satellites (Platnick et al. 2017). The atmosphere composition is assumed to be spatially homogeneous over the entire globe, specified by the US Standard atmosphere profile. Adding further classes beyond the 12 described above would add significant computational complexity that may not make the simulated reflectances more realistic. The distribution of surface types is taken from the International Geosphere–Biosphere Program (IGBP) land cover classification data set (Loveland et al. 2009). The 17 IGBP classification types are mapped into weights among the four ASTER types (given in the table in Appendix A.1). These weights are mostly ad hoc selections motivated by the detailed IGBP-type descriptions. For example, the "open shrublands" is mapped to a mixture of sand and vegetation at a 70% to 30% ratio. The IGBP data are at 0.25° resolution, while the simulations are most commonly done at 3° resolution to match the cloud fraction maps (explained below). The weights assigned at the IGBP resolution are therefore spatially averaged over 12 × 12 sample blocks to reduce the spatial resolution. Because the IGBP classification is a single static data set, the TSS tile fraction map based on the IGBP data does not contain any seasonality in the coverage of snow, ice, or vegetation.

The cloud fraction maps are taken from European Centre for Medium-Range Weather Forecasts Reanalysis-Interim (ERA-I) data products (Dee et al. 2011).⁴ These data contain gridded fields of cloud water and ice concentration, and cloud fraction. There is no straightforward relationship between the reanalysis data values and the cloud class weighting fractions required for the planet specification, so we used an ad hoc relationship that produced qualitatively reasonable renderings of Earth. The ERA-I data are available at high resolution for this application (0.75° blocks in latitude and longitude), so the original data are averaged across 4 × 4 samples to produce a final resolution of 3° × 3°.

A similar level of cloud model complexity was used in Robinson et al. (2011). In that study, the authors concluded that four cloud classes were required: two optical depths (5 and 15) for the water and ice clouds. However, we have found in our own tests that it is difficult to distinguish finer detail in the cloud types from differences in cloud fraction. In the Robinson et al. (2011) study, the cloud fractions were fixed using Earth observation data (from the Moderate Resolution Imaging Spectroradiometer). It is likely that to fully represent the EPOXI spectra with observed cloud fractions, two different cloud optical depths were needed. The ERA-I data used in this study to define the spatial cloud distribution are less precise and add a degree of freedom (in the mapping from ERA-I variables to tile fractional weights) that was not present in the Robinson et al. (2011) study.

Finally, in the last presented case study we use data from the TraCE-21ka project (He 2011)⁵ to specify the Earth surface. This paleoclimate-focused project investigated different mechanisms involved in Earth's climate evolution from the LGM, 21 kyr before the present, to the (pre-industrial) present. The Community Climate System Model (CCSM version 3)⁶ (Collins et al. 2006), a coupled atmosphere–ocean–land model, was used to simulate the Earth through this period. We use the simulation output to represent an Ice Age Earth by taking the model output fields of ice cover (on land and ocean), vegetation, and clouds, and compare it to the present-day Earth using the same simulation output.

The input surface reflectance and cloud properties are the same as described above. The TraCE-21ka CCSM output was used to define the coverage of the ice and vegetation surface fractions, as well as the cloud cover. Since the CCSM in this case is run at a 3.75° resolution, we use the gridded output data at the native resolution. The key feature of interest in these simulations is the surface ice, which changes significantly between the glacial maximum and the present. In addition, since the TraCE-21ka data were stored at monthly time resolution, the simulation does contain seasonality in the ice and vegetation coverage. Therefore, the tile map derived from TraCE-21ka will contain seasonal variations, unlike the static IGBP-derived tile map.

The most accurate translation of CCSM data to the TSS simulation would use instantaneous fields (for example, the model state at each internal time step); however, only the monthly means for TraCE-21ka were stored. This decision was made because of the complexity of the model run, which made the full time resolution over the 21,000 yr simulation period too large to store. So, our translation of the TraCE-21ka could introduce biases due to the nonlinear relationships between the cloud fractions, surface properties and the resulting spectral reflectance. However, in the example presented here, we are focused on the change in the spectrum between the two time periods, so to first order the reflectance biases cancel.

The monthly mean fields are mapped into fractional tile composition for the TSS, using the same four classes as described above in Section 3.1: ocean, vegetation, snow, and sand. The snow category is the sum of the land ice and sea ice categories, and a fractional scaling of the snow depth on land. The vegetation category is the sum of all vegetation types used in the CCSM land model. The residual fraction for each tile is then assigned to either the ocean or sand class, depending on whether the tile exists within the land mask defined in the CCSM output. All tile data are used at the native resolution of the model (3.75° × 3.75°). The atmospheric gas concentration and cloud scattering properties are the same as those described in Section 3.1.

We use the LBLDIS (LBLRTM-DISORT) RTM (Turner 2005) to compute the tile reflectances. This model combines tables of gas absorption optical depths from the Line-By-Line Radiative Transfer Model (LBLRTM), a high-accuracy line-by-line model (Clough et al. 2005), with the Discrete Ordinates Radiative Transfer code (DISORT) (Stamnes et al. 1988), to enable multiple scattering calculations at high spectral resolution. The high spectral resolution allows for the detailed absorption band structures to be accurately modeled for arbitrary instrument spectral resolutions or spectral filter shapes. While the LBLDIS was developed primarily for calculations in the thermal infrared, the most recent release (v3.0) has included a Rayleigh scattering component according to the approximation from Hansen & Travis (1974). The advantage of LBLDIS for this purpose is its great flexibility in terms of wavelength coverage and atmospheric properties. The primary limitations are the lower accuracy at large zenith angles, the computational complexity of the monochromatic calculations over large wavelength ranges, especially at visible wavelengths, and the lack of polarization in the DISORT scattering algorithm. Finally, the current implementation of LBLDIS does not allow for non-Lambertian surfaces.

The first limitation is inherent to plane-parallel models of multiple scattering, which lose accuracy at large zenith angles. For LBLDIS runs, we compute over a grid of angles up to 80 degrees. When the integration over an exoplanet disk involves larger zenith angles, the reflectance at 80 is reused. This is implemented by simply repeating the 80 degree result within the reflectance database with a "pseudo-observation" at 90 degrees.

The second limitation is solely related to computational complexity. In practice, we must reduce the number of monochromatic wavelengths needed to simulate the entire wavelength range. In the visible to near-infrared wavelengths, the Earth's atmospheric absorption is mostly confined to absorption bands (primary molecular oxygen, ozone, and water vapor), separated by highly transmissive window regions. Within the window regions, the cloud properties, Rayleigh scattering, and ozone absorption are important, but the spectral variation for all these components is much slower and requires fewer wavelengths.

To thin the set of required monochromatic wavelengths, we take an ad hoc approach that selects a subset of the finely spaced monochromatic grid produced by LBLRTM. We start with a coarse, evenly spaced grid of monochromatic wavelengths. The atmospheric transmission is computed at zenith on this coarse grid and then linearly interpolated back to the fine monochromatic transmission grid. The difference between these two transmission calculations is the error incurred by attempting to reconstruct the original fine monochromatic transmission from the coarse subsample. Wavelengths with large reconstruction error are identified and added to the coarse grid. Since these newly added wavelengths can be at arbitrary locations, the new coarse grid becomes a sparse and unevenly spaced subset of the original fine grid. This process is iterated a number of times until the reconstruction error is within a tolerance (0.1%) over the entire wavelength range. Finally, in order to accurately capture the effect of clouds on the depths of absorption features, we repeat this process for the transmission of the entire atmosphere, as well as the partial atmosphere from each cloud top (see Table 2). The final list of selected monochromatic frequencies is the union of these three sets. Figure 4 shows an example of a resulting wavelength selection around the oxygen A-band absorption feature. In the center of the R-branch at 759.6–762.0 nm, the sampling is very dense in order to correctly capture the shapes of the various absorption lines. In the wing of the band at 759.0 nm, the spacing increases rapidly, since the high transmission and slow spectral variation do not require dense sampling.

Zoom In Zoom Out Reset image size

The third limitation is the fact that the DISORT solver is unpolarized. Because molecular scattering has a strong polarization, this can lead to errors in the simulation of the unpolarized reflectance (Mishchenko et al. 1994). Polarization effects are strongest for the second order of scattering, so the relative error is highest for moderate optical depth from molecular scattering over absorbing surfaces in cloud-free atmospheres. The resulting reflectance error can be up to 10% for Earth's atmosphere for specific scattering directions (the unpolarized radiance is too low near forward and backward scattering, and too high near 90° scattering angles). Since the error is primarily related to the scattering angle (and thus the planetary phase angle), the resulting impact on the disk-integrated reflectance is a bias as a function of phase angle. The upper bound on this error is about 10% (see Appendix A.2 for detailed information), computed for a total molecular scattering optical depth of 0.5 over a surface albedo of 0. For the standard Earth atmosphere, this optical depth corresponds to a wavelength of approximately 370 nm, where the contribution of molecular scattering to the Earth reflectance is largest. This bias is an upper bound, since additional scattering from the surface (ignoring ocean glint) and clouds will be nearly unpolarized and reduce the relative impact of the polarization-related bias.

Since the LBLDIS cannot use non-Lambertian surfaces, the ocean surface model is the Lambertian spectral reflectance from the ASTER model (see Section 3.1). Thus, the current implementation will not simulate ocean glint, which could be a potentially important signature for detecting liquid water on exoplanets through direct detection. However, the ocean glint signature is only significant at large phase angles ("crescent" phases) (Robinson et al. 2010; Lustig-Yaeger et al. 2018). For quadrature or smaller phase angles, the relative contribution from ocean glint is not significant (Lustig-Yaeger et al. 2018). All the case studies presented here are at quadrature or smaller phase angles, so lack of ocean glint should not impact the results.

In Ford et al. (2001), rotational light curves were computed for Earth models with various configurations using a Monte Carlo-based radiative transfer tool. Although they used different input data sets for surface reflectance, different cloud scattering properties, and utilized non-Lambertian surface scattering models, the first-order features of the light curves are driven by the variation of the visible surface features which are very similar in both frameworks. So, despite the differences in the detailed simulation aspects, we can expect to recreate many of the primary features in the light curves.

Figure 5 shows the results for three Earth models (no atmosphere, clear atmosphere, and cloudy atmosphere) at several wavelengths (450, 550, 750, 950 nm). The observation geometry places the observer over the Equator with the Sun illumination at 90 degrees phase angle (quadrature phase). The cloudy atmosphere simulation includes cloud fields from ERA reanalysis at three times, separated by two days, to minimize the spatial correlation between cloud fields. The difference between the three light curves for the cloudy atmospheres is an estimate of variability in the disk-integrated signal due to changes in the cloud distribution. The cloud fields are fixed within each simulation, so there is no realistic diurnal variability of the cloud field.

Zoom In Zoom Out Reset image size

Our results show qualitatively similar light curves to those presented by Ford et al. (2001), especially at wavelengths with small extinction optical thickness from the atmosphere. Unlike the Ford et al. (2001) study, we can also examine the impact of the molecular absorption and scattering on the rotational light curve, as our model is based on full multiple-scattering radiative transfer simulations including molecular absorption bands. At the shortest wavelength, 450 nm, the reflectance is dominated by atmospheric scattering. The no-atmosphere model has very low reflectance, due to the relatively low surface reflectance of natural materials at the near-ultraviolet wavelength. The clear atmosphere has much higher reflectance from molecular scattering. Clouds increase the reflectance even more, and a relative maximum appears to occur in the range of 0.1–0.2 days. The view of Earth at this time is over the West Pacific, suggesting that this maximum in reflectance is due to the cloudy region in the Tropical West Pacific Warm Pool (Matus & L'Ecuyer 2017). At the blue visible wavelength, 550 nm, the signatures are qualitatively similar, except for a relative increase in contribution from the land surfaces and relative decrease in the molecular scattering. The African continent rotates into the illuminated portion at approximately 0.5–0.6 days, so the light curve shape now shows behavior consistent with contrast between land and ocean surfaces. At the 750 nm wavelength, the atmosphere has a very low extinction optical depth, meaning there is little contribution from molecular scattering. The shapes of the light curves for both the no and clear atmospheres are now dominated by the land–ocean contrast. Some of these features are still visible in the cloudy atmosphere, but the clouds reduce the relative amplitude. Finally, at 950 nm, absorption from water vapor reduces the total reflectance substantially in the clear and cloudy atmospheres. The cloudy atmosphere now shows a stronger relative maximum at 0.15 days, which is consistent with this reflectance signature's relation to high-altitude clouds in the Tropical Western Pacific. The water vapor column increases rapidly near the surface, and the high-altitude cirrus and anvil clouds associated with convection will be above the atmospheric layers with strong water vapor absorption.

NASA's EPOXI mission consisted of several observation programs using the Deep Impact spacecraft after the primary mission was completed. One component of the EPOXI mission was the Extrasolar Planetary Observation and Characterization (EPOCh) investigation. This investigation included observations of Earth, as well as observations of stellar systems with known transiting exoplanets. The Earth observation data set contains time-resolved spectrophotometry of Earth from the visible to infrared. The spacecraft distance from Earth ranged from 0.18 to 0.34 AU during the Earth observations in 2008, allowing for observation of the whole Earth disk in an analogous manner to extrasolar planets. The primary instruments relevant to EPOCh are the high-resolution instrument visible camera (HRIVIS) and the high-resolution instrument infrared spectrometer (HRIIR). The HRIVIS measures photometry in broadband short-wavelength channels (350–950 nm) and the HRIIR measures moderate-resolution spectra in the short-wavelength infrared (1000–4500 nm). These data were collected at three observation geometries in a small range of phase angles (58–77 degrees). Details on the Earth observations are documented in Livengood et al. (2011). These data have been used to validate other three-dimensional exoplanet simulations (Robinson et al. 2011). To provide an initial evaluation of the TSS, we simulated the "EarthObs1" and "EarthObs5" from Livengood et al. (2011), from 2008 March 18 and June 4. (The "EarthObs4" was not used in this comparison, due to the fact that the Moon was also within the field of view.) We used the diurnally averaged photometry for comparison, as we are most interested in evaluating the framework's ability to represent the spectral variation of Earth's reflectance. The HRIVIS radiometric values are taken from Livengood et al. (2010, Table 2), which reports the diurnally averaged at-aperture detected irradiance of Earth in W m^–2 μm^–1 scaled to 1 AU distance. The irradiance is converted to equivalent disk-integrated radiance by dividing by the solid angle of Earth at 1 AU. The HRIIR spectra were processed using the methods described in Livengood et al. (2011), and were converted to disk-integrated radiance by applying a correction to rescale the values relative to the solid angle of the full Earth disk. The TSS output spectra, in reflectance, were converted to radiance by multiplying by the solar irradiance spectrum at Earth. The solar irradiance spectrum is provided by the AER solar source function program⁶ , based on the solar spectrum provided by Kurucz (1992).

Figure 6 shows the results of using the model directly with no rescaling. The simulated values are typically within 10% of the HRIVIS measurements, with the largest difference at the 950 nm absorption band where the simulation is about 15% too low. In the comparison to the HRIIR data, the simulated radiance appears 20%–25% too large in the atmospheric window regions (1250, 1700, 2200 nm) and 10%–40% too small in the absorption features (1400, 1800, 2550 nm). The absorption features are all due to water vapor, implying the simulation's water vapor absorption is too strong. Since a single atmospheric water vapor concentration was used (the US Standard atmosphere), it is likely that the true water vapor in polar regions is lower, and increasing the total disk-integrated reflectance. The reason for the high bias in the simulated radiance in the window regions is not clear, but could be due to the simplistic choice of surface reflectance. The ASTER library's sand reflectance may be too bright at wavelengths larger than 1000 nm, where deserts tend to be much brighter than other Earth surface types.

Zoom In Zoom Out Reset image size

The model differences in Figure 6 compare favorably to the results of Robinson et al. (2011). In their Figure 5, similar data–model spectral comparisons are shown for individual EPOXI observations. The overall differences between the HRIIR data and model have a similar pattern as the TSS comparisons presented here, in that the radiance in atmospheric windows is often too large (up to 15%) and too low (10%–30%) in the absorption bands. The smaller data–model differences are expected, due to the higher complexity of the input data sets used by Robinson et al. (2011) in the construction of the Earth model. Specifically, the TSS used a two-class cloud model, with a globally constant atmospheric composition equal to the US Standard atmosphere, while the Robinson et al. (2011) model used a four-class cloud model, along with spatially varying atmospheric composition data derived from global Earth observations. In particular, the three-dimensional water vapor fields used by Robinson et al. (2011) would improve the modeled radiance in the absorption bands at 950, 1400, and 1800 nm.

The final case study illustrates an application of the TSS to simulate an exoplanet other than a clone of modern Earth, in this case corresponding to a past epoch in Earth's history. This case also highlights how the TSS can be applied to climate simulation output. As described in Section 3.1, we translate the monthly mean output from the TraCE-21ka climate simulation into TSS planet specifications. Through the late Quaternary period of the past one million years, the Earth has experienced approximately regular ice age cycles once every 100,000 years (Hays et al. 1976), with every ∼90,000 years of cold glacial period followed by ∼10,000 years of warm interglacial period. The modern Earth is in the Holocene, the latest interglacial stage of glacial–interglacial cycles. The TraCE-21ka simulation covers the last 21,000 years, covering the most recent transition into an interglacial period. While the simulation represents major climactic change within the late Quaternary, viewed over the entire Earth history, these changes are relatively minor compared to extreme conditions like the "snowball" earth (with global ice coverage) or the Paleocene–Eocene thermal maximum (with likely no surface ice) (Zachos et al. 2001; Kasting & Ono 2006).

Figure 7 shows the resulting disk-integrated spectra from 400–1000 nm, and associated three-dimensional renderings of the planets at the start and end of the simulation, specifically the LGM and present day. All simulations were performed for 90° phase angle observations, and in order to emphasize the surface properties, we show simulations of cloud-free atmospheres. In boreal summer (upper row), the north polar regions are the primary surface feature, and the reflectance spectra are relatively flat. The influence of the ice sheet is apparent, both in the three-dimensional rendering and in the shape of the spectrum. During the LGM, an ice sheet has replaced the northern boreal forests that exist in the present-day Earth. The result is the vegetation "red edge" spectral slope across 650–750 nm is suppressed in the LGM relative to the present-day Earth. For a simple estimate of the magnitude of the red edge, Figure 7 includes the ratio of the reflectances at two wavelengths (680 and 755 nm) that bracket the region where the vegetation reflectance sharply increases. Larger values of this ratio indicate a stronger red edge signature. In boreal winter (lower row), most of the reflectance arises from lower-latitude surfaces which show less change from the LGM to the present day. The southern part of the North American ice sheet is still visible and contributes to the total disk-integrated reflectance. The resulting spectral reflectance shows a relatively constant increase in reflectance over the entire spectral range, since the ice surface model has very high (>0.8) and smoothly varying reflectance over the entire wavelength range. The constant increased reflectance reduces the relative amplitude of the vegetation signature, but there is still clearly a strong "red edge" due to the influence of mid-latitude and equatorial regions.

Zoom In Zoom Out Reset image size

The TSS has similar capabilities to several exoplanet spectral modeling tools, but with important differences in the model assumptions and the intended applications.

As metioned in Section 2.1, the three-dimensional exoplanet model definition, with independent geometry calculations to discretized tiles in the spherical model, is very similar to the approach used in the Virtual Planetary Laboratory (VPL) Earth model (Robinson et al. 2011; Lustig-Yaeger et al. 2018). The TSS is intended to be a simpler and more modular approach, by splitting the radiative transfer and geometric calculations. The VPL model performs the high spectral resolution calculations for each tile, which is the most accurate method, but has a large computational cost. The TSS does interpolation of the radiative transfer output at the tile geometries, which is less accurate but allows for higher simulation throughput. Thus, the TSS would more easily allow for long time period simulations or ensembles of simulations over a wider range in observation geometries and planet specifications.

Another simulation tool is available as post-processing software for the Resolving Orbital and Climate Keys of Earth and Extraterrestrial Environments with Dynamics (ROCKE-3D) general circulation model (Way et al. 2017). The ROCKE-3D method performs similar geometric calculations over the exoplanet, discretized according to the ROCKE-3D model grid. The model's spatially varying atmospheric and surface components will naturally be included in these calculations. This method makes a simplifying assumption that the radiation field is isotropic and directly coupled to the model's internally computed radiative flux at the top of the atmosphere. The calculation would be limited to the available spectral band specifications available within ROCKE-3D, so high spectral resolution calculations, in particular absorption features, may be impractical. The TSS is relatively more complex, since it can use output of radiative transfer codes that more accurately model non-isotropic scattering from clouds.

Finally, the Planetary Spectrum Generator (PSG; Villanueva et al. 2018) is a suite of simulation tools intended to be used for a variety of planetary bodies and observation scenarios. The PSG is extremely flexible, and includes a variety of state-of-the-art radiative transfer model physics. The PSG is applicable to a very wide range of planetary bodies (comets, asteroids, planets, and exoplanets), for a wide range of observation platforms and spectral ranges. The PSG does not directly compute the radiation field across the exoplanet disk with varying tile geometries. Instead, it uses a single geometry that is expanded across the visible exoplanet disk by assuming an isotropic radiation field and weighting the isotropic expansion on a fixed grid over the disk. A more detailed calculation similar to the TSS could be performed by a suite of runs of the PSG, using the programmatic interface, and using finer grained controls over the model geometry. This technique would incur the large computational expense of individual RTM runs for each tile and each geometry. It may be possible to use the PSG to generate an RTM database that could then be used within the TSS. This could allow for more efficient calculation of time series but still include desired aspects of the PSG's available RTM physics.

While our specific implementation for the RTM in the presented case studies used a Lambertian surface, the TSS could immediately utilize more sophisticated RTM model calculations that included non-Lambertian surface models. For example, an RTM could include specular water reflection (ocean glint). Users must ensure that enough observation angle grid points are performed in order to accurately resolve the reflectance variation around the specular reflection geometry. We plan to develop new RTM databases that include a non-Lambertian ocean surface model to enable more accurate simulations at large phase angles.

A polarized RTM could also be used to generate the reflectance, which would eliminate the errors discussed in Section 3.2 and Appendix A.2. A polarized RTM may be needed for terrestrial exoplanets with thicker atmospheres and thus increased molecular scattering optical thicknesses. On the other hand, a fully polarized disk integration will require upgrades to the TSS, since the individual tiles' Stokes vectors would require additional geometric calculations.

The software is hosted on github at https://github.com/aronnem/TerranSpecSim.

Table 3 contains the mappings between the IGBP surface type classifications and the four ASTER surface spectral reflectance types used in the radiative transfer model runs. See Section 3.1 for details.

Molecular scattering is strongly polarizing, and neglecting the polarization effects in the radiative transfer solver will result in errors in the total reflectance. As discussed in Section 3.2 and Mishchenko et al. (1994), the error depends primarily on the scattering angle. The error is largest for near-forward and backward scattering (low bias in radiance computed with an unpolarized solver) and near-90° scattering (high bias with an unpolarized solver). The relative error is also largest for a molecular scattering layer over a dark surface. In order to estimate an upper bound on the disk-integrated reflectance, we used the TSS framework to perform numeric integration of the estimated error for a model planet with only molecular scattering. To derive the error magnitudes for specific observation angles, we used the published reflectance tables in Natraj et al. (2009), which contain very accurate reference values for a pure molecular scattering layer over a Lambertian surface. These reference values were computed from numerical solutions to the fully polarized radiative transfer equations. The Natraj et al. (2009) tables are computed for a set of optical depth from 0.02 to 1.0, and the optical depths are converted to the wavelengths that produce the same total depth for the standard Earth atmosphere column density. We then computed an equivalent set of reflectance values at the same grid points (observation angles and optical depths) using LBLDIS. Two TSS runs were peformed for a full phase curve of a 90° inclination ("edge-on") orbit, using each of the two RTM databases. This yielded a modeled phase curve for a non-absorbing Earth atmosphere over a Lambertian surface for a set of wavelengths covering the near-ultraviolet to near-infrared.

Figures 8 and 9 show the resulting error estimates, plotted as a phase curve from 0 to 153°. In Figure 8, the error is plotted as a percentage difference, which shows the larger impact over the non reflecting (albedo = 0) surface, of up to 10% at zero phase angle and the shortest wavelength. The relative impact is reduced to 3% over the bright (albedo = 0.8) surface, due to the relatively larger contribution of the unpolarized surface scattering. In Figure 9 the error is displayed in apparent albedo, where the reflectances are divided by the reflectance of a perfectly reflecting Lambertian sphere. This figure shows the magnitude of the error is nearly independent of albedo, and the smaller relative errors for the albedo = 0.8 surface in Figure 8 are simply due to the decreased relative contribution of the molecular scattering to the total reflectance. Finally, we again note that these values represent upper bounds, since realistic atmospheres that include cloud and aerosol scattering will tend to reduce these errors.

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size