1. Introduction
The rapidly expanding list of confirmed exoplanet detections and accumulating evidence about the histories of planets in our solar system has created an increasing demand for tools that can complement available observations to provide insights about which planets may be habitable or inhabited, now or in their past. To date, studies of the climates and habitability of planets other than modern Earth have been carried out primarily with one-dimensional (1D) radiative-convective models (e.g., Kasting1988; Kasting et al.1993; Pavlov et al.2001; Segura et al.2003; Domagal-Goldman et al.2008,2011; Kitzmann et al.2010; Zsom et al.2012; Kopparapu et al.2013; Rugheimer et al.2013,2015; Grenfell et al.2014; Ramirez et al.2014a,2014b; Meadows et al.2016). These models have the virtue of computational efficiency, permitting exploration of a wide range of parameter space and coupling to complex atmospheric chemistry models. Their limitations include an inability to properly account for the effects of clouds, atmospheric and oceanic heat transports, sea ice and land surface variations, obliquity effects, day-night contrasts, and regional aspects of habitability.
Modeling of terrestrial climate and climate change was initially performed with 1D models as well (e.g., Manabe & Strickler1964; Hansen et al.1981), but soon gave way to three-dimensional (3D) general circulation models (GCMs; sometimes referred to as global climate models), which are lower-resolution versions of the models used for numerical weather prediction. GCMs have evolved from atmosphere-only to coupled atmosphere-ocean-sea ice models, and more recently have added atmospheric and ocean chemistry, land and ocean ecosystem dynamics, and dynamic land ice to create today’s Earth system models (Jakob2014) that are the basis of projections of 21st century anthropogenic climate change.
The first application of a GCM to another planet was the Mars model of Leovy & Mintz (1969), and Joshi et al. (1997) performed the first hypothetical exoplanet GCM simulation. Since these pioneering studies, GCMs have been used to simulate the dynamics and climates of a broad range of rocky planets, past and present, as well as planets with thick H2-He envelopes (e.g., Forget et al.1998; Williams & Pollard2002; Hollingsworth & Kahre2010; Merlis & Schneider2010; Abe et al.2011; Edson et al.2011; Heng et al.2011; Lebonnois et al.2012; Charnay et al.2013; Leconte et al.2013; Shields et al.2013; Yang et al.2013; Hu & Yang2014; Godolt et al.2015; Kaspi & Showman2015; Wolf & Toon2015; Wordsworth et al.2015; Kopparapu et al.2016; Popp et al.2016; Turbet et al.2016; Zalucha2016; Boutle et al.2017). GCMs self-consistently represent all the processes that 1D models cannot, although they have their own limitations: uncertainties in parameterizations of small-scale processes, computational cost that requires radiative transfer and chemistry to be represented in less detail than in 1D models, and a level of detail that cannot be constrained as well by observations for other planets as it can for Earth. Increasingly, GCMs are playing a key role in a “system science” approach that considers planetary climate and habitability in the larger context of the evolution of the solid planet, its parent star, and other planets and planetesimals that affect its evolution.
In this paper, we describe a new planetary and exoplanet GCM, the ROCKE-3D (Resolving Orbital and Climate Keys of Earth and Extraterrestrial Environments with Dynamics) model. ROCKE-3D is developed from its parent Earth climate GCM, the NASA Goddard Institute for Space Studies (GISS) ModelE2 (Schmidt et al.2014). ModelE2 was the GISS GCM version used for the Coupled Model Intercomparison Project Phase 5 (CMIP5), the most recent phase of a protocol by which successive generations of Earth climate model results are made publicly available for systematic analysis by the international community. ROCKE-3D is configured to simulate the present and past atmospheres of rocky solar system planets, as well as rocky exoplanets. Like several other planetary GCMs, ROCKE-3D is adapted from a previously existing Earth GCM (i.e., PlanetWRF, Richardson et al.2007). Unlike most other planetary GCMs, ROCKE-3D is based on the most recent published version of its parent Earth model, developed and used in part by the same people who developed the Earth model, and will evolve in parallel with future generations of the Earth model, thus benefiting from emerging insights from Earth science into physical processes that are also relevant to other planets.
The baseline ROCKE-3D version described in this paper is referred to as Planet 1.0. In the following sections, we discuss the challenges involved in adapting an Earth GCM to simulate other rocky planets, the choices made to make Planet 1.0 as generally applicable as possible, and the remaining limitations that will not be addressed until the next generation of the model has been developed. ROCKE-3D Planet 1.0 has already been used to simulate hypothetical ancient Venus scenarios (Way et al.2016), while simulations of several deep Earth paleoclimate eons, modern Mars, and hypothetical exoplanets are in progress.
In principle, it should be possible to modify an Earth GCM to simulate other planets simply by changing relevant external parameters. In reality, though, terrestrial GCMs are designed with present-day Earth in mind, and are programmed by a large group of people of varying backgrounds and experience whose composition evolves over several decades. At any moment in its history, therefore, a GCM is a mix of modern and obsolete programming approaches, visionary and myopic coding philosophies, and best and worst practices that necessitate new approaches to make the model sufficiently general for planetary applications. Many of those approaches will be discussed herein.
In Section2 below, we discuss the present Planet 1.0 model resolution and possible ocean configurations. In Section3, extensions to the model calendar system are reviewed. These allow for slower or faster rotating worlds (than present day Earth), synchronously rotating worlds, and even retrograde rotation like that of present-day Venus. Section4 discusses the major physics parameterizations in the model, whereas Section5 covers its geophysical properties. Section6 describes several examples of GCM modifications for Planet 1.0 that have fed back to the parent Earth GCM. Section7 covers appropriate uses for ROCKE-3D, and Section8 contains our conclusions. Two appendices provide a description of input and post-processing tools available external to the model.
The Planet 1.0 version of ROCKE-3D is available for community use. It can be downloaded from the NASA GISS ModelE repository6 and documentation for ModelE2/ROCKE-3D is readily available as well.7 Future versions will be made available at the same site. The next generation is expected to be created after GISS releases the first CMIP6 (Eyring et al.2016) version of ModelE2 in late 2017. Its availability will be announced on the GISS Astrobiology website8 upon release. Although we intend ROCKE-3D to parallel the development of its Earth parent, we encourage user implementation of new physics elements that do not yet exist or that cannot be generalized in ModelE2 to expand the range of potential applications.
2. Model Configurations
2.1. Resolution, Dynamics, and Throughput
In describing the physics of ROCKE-3D, we refer to physics from the present operational version of the parent Earth model as derived from “GISS” or “ModelE2,” and new capabilities as those of “ROCKE-3D.” ModelE2 is a Cartesian gridpoint model routinely run at
latitude-longitude atmospheric resolution with 40 vertical layers, and at 1° × 1
25 latitude-longitude ocean resolution with 32 vertical layers. This resolution has been retained for certain deep-Earth paleoclimate simulations, where the higher resolution permits better comparison to geological data, as well as better portrayal of the atmospheric and oceanic dynamics.
GCM atmospheric (as opposed to oceanic) resolution should, at a minimum, be fine enough to crudely resolve the dominant scales of atmospheric motion. Typically, this is assessed using the Rossby radius of deformation (the typical spatial scale of midlatitude low- and high-pressure centers)
, where:N, the BruntVäisälä frequency, is proportional to the static stability;H, the scale height, depends on temperature, gravity, and atmospheric composition; andf, the Coriolis frequency, is proportional to planetary rotation rate. For Earth,
(
Earth’s radius) and
grid boxes are about 200 to 250 km in size, allowing such features to be adequately resolved. For simulations of other planets, most initial studies with Planet 1.0 have been of smaller planets, for which grid boxes at the same resolution are smaller, or more slowly rotating planets, for whichLd is larger than on Earth. For these simulations, it has been possible to run Planet 1.0 at
atmospheric and oceanic horizontal resolution with no loss in accuracy, but at almost an order of magnitude faster speed. This lower-resolution version of Planet 1.0 has 20 atmospheric layers (but with an option for 40 layers) with a model top at
(about
altitude), and in coupled mode, 13 ocean layers with maximum depth up to
.
The ModelE2 dynamical core uses finite differences with atmospheric velocity points on the Arakawa B-grid (Hansen et al.1983). For tracers, nine higher-order moments are carried, as well as the mean tracer amount in each grid box, yielding an effective resolution that is higher than the nominal model resolution (Prather1986). The vertical discretization uses a sigma coordinate from the surface to 150 hPa, and switches to constant pressure layers above.
Planet 1.0 can be run on a capable laptop for modest integrations at this coarser resolution, but the bulk of our research is conducted on the NASA Goddard Space Flight Center Discover cluster of Linux scalable units.9 ROCKE-3D is typically run with an atmosphere and ocean resolution of
, with 20 or 40 atmospheric layers and 13 ocean layers, when run with a fully coupled ocean. The default ROCKE-3D model has 40 atmospheric layers, with SOCRATES in the GA7.0 configuration (see Section4.1) as the radiation scheme. Details on run-times for different configurations are provided in Section4.7.
The parameterized physics in Planet 1.0 is largely the same as that in ModelE2, but several changes that were made after Schmidt et al. (2014) to correct ocean and radiation physics errors have been adopted for Planet 1.0.
2.2. Ocean Models
Oceans are crucial to the accurate 4D portrayal of a planet’s climate system. Energy, moisture, and momentum are exchanged between the atmosphere and oceans, and the transitions between different phases of water drive some of the most significant feedback mechanisms operating in the climate system. The oceans provide the major source of moisture that drives the hydrological cycle, whereas the freezing and melting of surface waters have a major impact on planetary albedo. Together, the transport of heat and these atmosphere-ocean interactions affect the geographic, seasonal, inter-annual, and even geologic-scale variations of a planet’s climate. In Planet 1.0, the oceans differ from other bodies of water (lakes, rivers), in that salinity and temperature combine to alter the 3D density structure, while surface wind stress is allowed to impact movement of water in the upper ocean. Salinity, temperature, and wind stress drive global ocean currents that transport energy on timescales that may exceed the orbital period of the planet by orders of magnitude. Ocean albedo is a function of both water and sea foam reflectance. Water albedo is calculated as a function of the solar zenith angle and wind speed; the sea foam reflectance is derived from Frouin et al. (1996).
Planet 1.0 allows for three different modes of ocean interaction. From simplest to most complex, these are (1) specified sea surface temperature (SST), (2) thermodynamic upper ocean mixed-layer, and (3) coupled dynamic ocean GCM.
2.2.1. Specified Ocean Surface Conditions
Specifying sea surface temperature (SST), including sea ice cover, is a common Earth climate modeling technique where SST observations are used as a surface boundary condition over a range of years or months in order to force an atmospheric GCM (AGCM). The GISS model uses twelve monthly arrays that define the ocean surface temperature and sea-ice distributions. The model interpolates the input into daily values, providing smoother transitions through an annual cycle. Specifying SSTs is the most commonly accepted technique for evaluating the efficacy of AGCM physics parameterizations when surface conditions are well-known (e.g., in performing hindcasts of 20th century climate). It is also used in Earth paleoclimate studies where proxy data can be used to reconstruct past ocean temperature distributions (e.g., MARGO Project Members2009). In this case, the purpose is generally to evaluate the consistency of land- and ocean-based observations, or simply to examine potential states of the atmosphere for various time periods in Earth’s history. Specified SST simulations are also used to collect the atmosphere-ocean flux information to generate the Q-fluxes to run the model in mixed-layer ocean mode.
2.2.2. Mixed-layer (Q-flux) Ocean Model
For other planets with surface oceans, prescribed SSTs are not an option and SSTs must instead be calculated interactively to be consistent with a given planet’s atmosphere and external forcing. The simplest way to do this is to couple the AGCM to a simple, thermodynamically active layer that represents the upper well-mixed layer of the ocean (typically tens to hundreds of meters deep). The temperature of the mixed layer responds to radiative and turbulent (sensible and latent) fluxes of heat across the ocean-atmosphere and ocean-sea ice interfaces, but lateral ocean heat transport is neglected. This approach has been the default choice for most exoplanet GCM studies to date (e.g., Shields et al.2014; Yang et al.2014; Kopparapu et al.2016; Turbet et al.2016). In the literature, this approach is typically referred to as a thermodynamic, mixed layer, slab, or immobile ocean model. The greatest limitations of this method are the assumption of zero horizontal heat transport by ocean currents and the inability to account for deep water formation related to vertical density gradients.
For Earth, where SST observations exist, a variant of the mixed layer approach known as the “Q-flux” method has been commonly applied to simulations of future climates (Miller et al.1983; Russell et al.1985). In the Q-flux approach, a control AGCM run with prescribed SSTs is first conducted to define the radiative and turbulent heat exchanges at the atmosphere-ocean interface that are consistent with the AGCM’s physics parameterizations. The implied horizontal ocean heat transport convergences that would be required to produce the observed SSTs and sea-ice cover in each mixed layer gridbox are then calculated and applied in a second simulation that couples a mixed layer ocean model to an AGCM, as a proxy for the effect of actual ocean heat transports. Sometimes diffusive heat loss through the lower boundary of the mixed layer is also included to mimic exchanges of heat with deeper ocean layers that are otherwise unrepresented in such models. The implied ocean heat transport convergences are themselves fixed, but their presence allows for a more realistic projection of sea-ice changes, and thus ice-albedo feedback, in a changing climate than is possible in a model that completely ignores ocean heat transport. Such models have traditionally been used to define the equilibrium sensitivity of Earth’s climate to a doubling of CO2 concentration, a common benchmark for assessing climate model uncertainty. The Q-flux approach is generally not applied for exoplanet GCM studies, because SST observations are not available to constrain them; hence the use of purely thermodynamic (
) oceans in most such studies without a dynamic ocean. Exceptions to this are the works of Yang et al. (2013), Godolt et al. (2015), and Popp & Eggl (2017), who perform several sensitivity tests with specified ocean heat transports. ROCKE-3D includes a
ocean option, but the error induced by ignoring ocean heat transport must be kept in mind when assessing such studies. An alternative that has been used for sensitivity studies is to prescribe a latitudinal profile of ocean heat transport in a mixed-layer model, with the latitude and magnitude of the peak transport as free parameters that can be varied (e.g., Rose2015). Furthermore, if an existing simulation with a dynamic ocean (see Section2.2.3) is available, the ocean heat transports from this model can, in principle, be used as a specified input to an otherwise thermodynamic ocean model (e.g., Fiorella & Sheldon2017).
2.2.3. Dynamic Coupled Ocean
Given the limitations of Q-flux models, recent generations of Earth climate models have instead coupled more realistic and computationally expensive dynamic ocean GCMs (OGCM) to AGCMs to simulate climate change. Most exoplanet GCM studies have eschewed the use of OGCMs because of the large thermal inertia of the ocean, and consequent long integration times required to reach equilibrium. However, several studies have revealed the importance of interactive ocean heat transport to climates of planets in parameter settings very different from Earth’s (Vallis & Farneti2009; Cullum et al.2014). The most dramatic example of ocean heat transport effects in the exoplanet context is the difference between the concentric “eyeball Earth” open ocean region simulated beneath the substellar point of a synchronously rotating aquaplanet with a thermodynamic ocean (Pierrehumbert2011) and the asymmetric “lobster” ocean pattern produced when a dynamic ocean is used (Hu & Yang2014). The extent of ocean heat transport effects depends on the specific planet configuration used. In particular, modulation of temperature extremes is less efficient on a planet with land masses that produce confined ocean basins than on an aquaplanet (Yang et al.2013).
The exploration of parameter space for saltwater ocean composition differs from that for atmospheric composition, in that the former has a more direct effect on density structure, circulation, and heat transport. The Earth’s thermohaline circulation was recently placed into perspective by Cullum et al. (2016), who demonstrated that an increase in mean salinity can cause the haline component to dominate.
Most ROCKE-3D simulations couple a dynamic ocean to the atmospheric model. The standard configuration uses a
resolution with 13 ocean layers that decreases model throughput by
or less, compared to a thermodynamic ocean. However, it also increases the equilibration time of the climate from decades to centuries of simulated time, with the exact time depending on the assumed ocean depth. Some of our deep Earth paleoclimate studies instead use the same
resolution, 32 layer ocean that is used by ModelE2. Transport by unresolved mesoscale eddies is represented by a unified Redi/GM scheme (Redi1982; Gent & McWilliams1990; Gent et al.1995; Visbeck et al.1997), as in ModelE2. The version used by Schmidt et al. (2014) contained a miscalculation in the isopycnal slopes that led to spurious heat fluxes across the neutral surfaces, resulting in an ocean interior that was generally too warm and southern high latitudes that were too cold. A correction to resolve this problem was implemented for ModelE2, Earth paleoclimate studies (Chandler et al.2013), and ROCKE-3D. The new code uses a mesoscale diffusivity of 600 m2 s−1, although some ROCKE-3D exoplanet simulations have used a value of 1200 m2 s−1 instead. The applicability of mesoscale eddy parameterizations designed for Earth models has not yet been investigated for planets that have different rotation rates, and thus different dominant spatial scales of eddies (Cullum et al.2014).
3. Calendar Changes for Modeling Other Planets
ModelE2 uses a clock and calendar to coordinate model operations that are not active during every time step, and to manage binning/averaging for seasonal and higher-frequency diagnostics. Prior to the development of Planet 1.0, this system made assumptions that were incorrect or inconvenient outside the context of modern Earth. For instance, the number of solar days per month was hardwired for a quasi-Julian 365 day calendar. Note that, unless otherwise specified, “day” refers to thesolar day for the planet being simulated. The original system did permit varying the rotational and orbital periods, as well as other orbital parameters (obliquity, eccentricity, and solar longitude), but provided only limited means to relate these to seasons. Further, a number of model components possessed implicit (hardwired) constants appropriate to the lengths of modern Earth day and year.
To enable the study of exoplanets, the calendar (indeed, the entire time-management system) in ModelE2 has been been redesigned to be extensible and highly encapsulated. The latter was crucial to reduce the likelihood of subsequent developers accidentally reintroducing assumptions about modern Earth into the model. The design of this new time-management system reflects the needs and priorities of climate scientists in several respects. The first priority was to ensure that the default behavior replicates the original behavior for simulations based on modern-day Earth. The other priority was for the new calendar to preserve, as much as possible, correspondence between planetary seasons, months, and days with those of Earth, in terms of basic orbital characteristics. Note that other communities have designed planetary calendars (primarily for Mars) with quite different priorities, such as preserving the number of days per month and seconds per hour (Gangale1986,1997; Allison1997; Allison & McEwen2000; Gangale & Dudley-Rowley2005). Here, the priority is to simplify interpretation of seasonal and diurnal diagnostics—similar to the approach in Richardson et al. (2007). In particular, the new calendar system preserves the intuitive notion of the diurnal cycle being divided into 24 equal “hours,” as well as the seasonal cycle being divided into 12 “months.” Note that a model “hour” will therefore not generally be 3600 s in duration, and months can be significantly longer or shorter than 720 hours. Additional machinery minimally tweaks the orbital period and model timestep to ensure that the simulation has an integral number of time steps per “day” and an integral number of days per year. All times and time intervals are expressed using exact integer arithmetic, to eliminate issues related to numerical round-off. We thus guarantee an exact number of simulation time-steps per day and days per year.
The specific duration of each calendar month is derived as follows. First, the solar longitude
, where
, is computed for the beginning of each month in a reference Earth orbit and calendar. The beginning of each month in the planetary calendar is then determined to have the same solar longitude angle as for the reference month, subject to rounding to ensure an integral number of days in each month. To relate the longitudes to times/durations, the corresponding mean anomalyMi is computed for the planetary orbit for the start of each month. We can deriveMi from the solar longitude via standard Keplerian orbit formulae. The starting day-of-yeardi for each month is then computed by scaling the delta mean anomaly (
by the number of calendar days per radian and rounding to the nearest day:
By default, the system uses the model’s standard Earth-based orbit and pseudo-Julian calendar as the reference. Thus, the planetary “February” will tend to be shorter than average, simply due to the short duration of February in the conventional Earth calendar.
Our basic design is to have the system derive an appropriate calendar directly from the orbital parameters of a given planet. By introducing software abstractions for both the orbit and the calendar, the system provides a natural mechanism for further extension. For example, leap-years and other requested extensions are easy to accommodate.
There is a crucial aspect of Earth’s orbit that is not particularly generic—a large separation of scale between days and years, such that the number of days per year is much, much larger than 1. In terms of the conventional Julian/Gregorian calendars, this permits months to have an integral number of days while simultaneously having roughly uniform duration. It also allows climate models to safely ignore fractional remainders of days that lead to leap-years. However, for extreme orbits, a lack of this separation of scale can can have spectacular consequences. The number of days per year can be less than the number of months, and each day can last a significant fraction of a year (e.g., modern Venus). In such cases, the default for our calendar is to break the correspondence between the calendar day and the solar day and constrain calendar to have at least 120 calendar days, i.e., at least 10 calendar days per month on average.10 The system has runtime switches that can eliminate this constraint, as well as the constraint that the rotational period is commensurate with the orbital period. The latter is crucial to differentiate an orbit such as that of modern Venus from a tidally locked orbit—both of which are of interest to ROCKE-3D modelers. Of course, one must exercise extreme caution when interpreting model diagnostics in such cases. Some months (and even some years!) may have zero solar days. A quantity averaged over one season or even one year may be highly biased, as parts of the planet remain entirely day or entirely night.
For tidally locked planets, it is convenient to have a mechanism to vary the longitude of the subsolar point. For example, Turbet et al. (2016) point out that, for a synchronously rotating world, the continents may be concentrated at either the substellar or anti-stellar point. This variation is supported in our framework by the “hourAngleOffset” parameter, which controls placement of the continents for a synchronously rotating world at any angle with respect to the substellar point. This approach is much simpler than the equivalent shift of all boundary condition data (topography, etc.).
Figure1 demonstrates that the model responds correctly to the calendar modifications for slowly rotating worlds, as the Hadley cells are clearly broadened for the slowly rotating planet versus the rapidly (Earth day length) rotating one.
Figure 1. Left: pressure vs. latitude stream function of the mean meridional circulation of a planet with Earth’s rotation period. Right: as in the left panel, but for a planet with a rotation period 128 days longer than an Earth sidereal day. As expected for a slowly rotating world, the Hadley cells are now much larger in latitudinal extent, due to the decrease in the strength of the Coriolis force at these slow rotation rates. Note that the color bar limits on the 128 day stream function are a little more than a factor of two larger than the one-day stream function, as one might expect.
Download figure:
Standard imageHigh-resolution image
Figure 2. Top row: visible wavelength surface albedo map, spectrally integrated surface albedo map, and snow and ice cover map for an aquaplanet simulation of Proxima Centauri b, using Planet 1.0 with the SOCRATES radiation scheme. Bottom row: corresponding maps for Neoproterozoic Snowball Earth. Snow cover over sea ice or land produces the strongest surface albedo response in both simulations most prominently in the VIS (290 to 690 nm) band.
Download figure:
Standard imageHigh-resolution image
The calendar has also been expanded to handle variable orbital eccentricities in time (Way & Georgakarakos2017). This would be useful in cases where a Jupiter-like planet perturbs the orbital elements of a nearby smaller terrestrial planet (e.g., Georgakarakos et al.2016).
4. Physics Parameterizations
Physical processes that operate on scales smaller than those resolved by a GCM must be parameterized in terms of grid-resolved variables. This section discusses those necessary for ROCKE-3D and how they are accomplished.
4.1. Radiation
4.1.1. The GISS Radiation Scheme
The radiation scheme in ModelE2 was first implemented in Hansen et al. (1983), with more detailed descriptions of the long-wave radiation scheme in Lacis & Oinas (1991) and Oinas et al. (2001), and the short-wave scheme in Lacis & Hansen (1974). Minor updates have been made to improve its accuracy (Schmidt et al.2006,2014; Pincus et al.2015), but its overall structure and parameterizations remain unchanged.
Recently updates have been made to the longwave radiation scheme to improve its accuracy for atmospheres that deviate slightly from that of present-day Earth. The tabulated Planck function has been extended to
, and the gas optical depth table has been updated to enable the major greenhouse gases in Earth’s atmosphere (H2O, CO2 and O3) to be replaced with other gases. The latter enables more accurate calculation of fluxes and heating for cases such as the Archean Earth, which had no O3 or O2, but may have had significantly larger amounts of both CO2 and CH4 than present-day Earth. These updates were recently used in a study of the early climate of Venus (Way et al.2016).
The shortwave radiation scheme uses the doubling and adding method to include the effects of multiple scattering (Peebles & Plesset1951; van de Hulst1963) with two quadrature points (Lacis & Hansen1974). Gaseous absorption is parameterized through analytical expressions for the frequency-integrated absorption as a function of pressure, temperature, and absorber amounts for each gas. The spectrum is divided into 16 gaseous absorption bands, each with one absorbing gas and a corresponding analytical function for the optical depth as a function of pressure, temperature, and absorber amount.
Stellar radiation input to the GCM drives both the planetary energy balance and photochemistry. Stellar spectral irradiance (0.115 to 100μm) to the top-of-the-atmosphere is provided to the model via an input file that can be changed for different stars, to drive the energy balance and provide UV fluxes for ozone calculations and photolysis rates. A software utility provided by GISS can be used to format high-resolution stellar spectra for input to the GCM (see AppendixA). The various modules of the GCM that utilize this spectral irradiance perform different spectral partitioning to suit their functions, such as for surface albedo or photosynthesis by plants and phytoplankton. The dynamic coupled ocean (Section2.2.3) vertically distributes this absorption with a profile that is currently insensitive to changes to the stellar spectrum, but which roughly corresponds to the current-Earth partitioning of visible and NIR:
whereF(z) is the ratio of net solar flux at depth z to the net solar flux at the surface, andf1 = 0.62 is a parameter that fits solar radiation attenuation by typical sea water and roughly represents the red/near-infrared fraction (Paulson & Simpson1977),z1 = 1.5 m,z2 = 20 m. For another stellar spectrum, such as that for M stars with relatively much higher NIR/VIS ratios,f1 should be altered, in lieu of introducing a physical spectral radiative transfer scheme for water absorbance. Some solar spectral radiation assumptions are still hard-coded into the model, so users should consult GISS personnel when interpreting results with alternative stellar spectra. More details are provided below in Section4.4 on the Cryosphere.
We note that, for this radiation scheme to be reasonably accurate, gas concentrations of radiatively active gases, and H2O, CO2, and O3 in particular, should be within a factor of 10 of present-day Earth values throughout the atmosphere. In addition, the short-wave radiation scheme should only be used with stellar spectra that are of the same spectral type as our Sun.
4.1.2. SOCRATES
In ROCKE-3D, we require a radiation scheme that can be applied to a wide variety of planetary atmospheres. Consequently, it is imperative to able to easily change spectral bands, extend the pressure and temperature range of opacity tables, and add and remove absorbers. Unfortunately, the parametrizations used by the radiation scheme in ModelE2 prohibit such a generalization at present. To ease adaptation of ROCKE-3D to different atmospheres, we have coupled it to the Suite of Community Radiative Transfer codes based on Edwards and Slingo11 (SOCRATES, Edwards1996; Edwards & Slingo1996). This radiation scheme is in operational use in the UK Met Office Unified Model, has previously been adapted to hot Jupiters (Amundsen et al.2014,2016), and is available under a BSD three-clause licence.12 Importantly, SOCRATES allows for changing radiation bands, altering pressures and temperatures in the opacity tables, and the inclusion of various combinations of gaseous absorbers with relative ease.
SOCRATES solves the two-stream approximated radiative transfer equation with multiple scattering for both the short- and longwave components. Several different two-stream approximations are available, but we default to using the practical improved flux method version from Zdunkowski & Korb (1985) with a diffusivityD = 1.66 for the longwave component, and the original version of Zdunkowski et al. (1980) with a diffusivityD = 2, for the shortwave component. We note that, unlike the two-stream equations presented in Toon et al. (1989), the two-stream equations of Zdunkowski & Korb (1985) can be applied with a variable diffusivity, enabling improved accuracy compared to the Toon et al. (1989) formulation. In order to improve the representation of the scattering phase functions with strong forward-scattering peaks, delta-rescaling (Thomas & Stamnes2002) is applied for both components.
Gaseous absorption is parameterized using the correlated-k method (Lacis & Oinas1991; Goody et al.1989), withk-coefficients derived via exponential sum fitting of transmissions (Wiscombe & Evans1977), tabulated as a function of pressure and temperature. We use the HITRAN 2012 line list (Rothman et al.2013) to calculate cross-sections line-by-line, using Voigt profiles and the choice of either the MT_CKD 3.0 (Mlawer et al.2012) or the CAVIAR (Ptashnik et al.2011) water vapor continuum. For planets with Earth-like atmospheres orbiting Sun-like stars, we use nine longwave and six shortwave bands, those used by the UK Met Office for global atmosphere configuration 7.0 (GA7.0, D. Walters et al. 2017, in preparation), given in Tables1 and2. We tabulatek-coefficients on 51 pressures equally spaced in
between 10−5 and 1 bar, and 13 temperatures spaced linearly in temperature between 100 and 400 K. The number ofk-coefficients for each gas in each band varies according to the required accuracy (see Amundsen et al.2014; Amundsen2015); for GA7.0, there are up to ∼20k-coefficients per band per gas.
Table 1. The Longwave Bands Adopted for Planets with Earth-like Atmospheres
| Long-wave band | Wavenumber ( | Wavelength ( |
|---|---|---|
| 1 | 1 to 400 | 25 to 10,000 |
| 2 | 400 to 550 | 18.18 to 25 |
| 3 | 550 to 590 and 750 to 800 | 12.5 to 13.33 and 16.95 to 18.18 |
| 4 | 590 to 750 | 13.33 to 16.95 |
| 5 | 800 to 990 and 1120 to 1200 | 8.33 to 8.93 and 10.10 to 12.5 |
| 6 | 990 to 1120 | 8.93 to 10.10 |
| 7 | 1200 to 1330 | 7.52 to 8.33 |
| 8 | 1330 to 1500 | 6.67 to 7.52 |
| 9 | 1500 to 2995 | 3.34 to 6.67 |
Note. These are the bands used by the UK Met Office for Global Atmosphere Configuration 7.0 (GA7.0, D. Walters et al. 2017, in preparation). Note that bands 3 and 5 contain excluded regions.
Download table as: ASCIITypeset image
Overlapping gaseous absorption is treated using equivalent extinction (Edwards1996), although random overlap (Lacis & Oinas1991) is also supported. Both of these methods combinek-coefficients calculated for each gas separately on-the-fly. Equivalent extinction relies on having a major absorber in each band; we use the adaptive equivalent extinction approach described in Amundsen et al. (2017) to determine the major absorber in each band independently for each column, which may also change in time. Pre-mixing of opacities (Goody et al.1989), wherek-coefficients are derived directly for the gas mixture for a given composition, would result in a faster radiation scheme; however, it requires newk-tables to be derived when gas amounts are changed (Amundsen et al.2017).
The GA7.0 configuration of SOCRATES is highly optimized for the present-day Earth atmosphere and irradiation spectrum. Consequently, the bands in Tables1 and2 will need to be changed to improve accuracy for atmospheres with significantly different compositions or stellar irradiation spectra. We are in the process of creating new SOCRATES configuration input files, called “spectral files,” that provide suitable accuracy for other atmosphere-star combinations; for instance, Fujii et al. (2017), where 29 shortwave bands were used to accurately handle absorption of stellar radiation by water vapor at near-IR wavelengths for large specific humidities. We have recently added additional physics to SOCRATES in order to simulate atmospheres with a large amount of CO2. Specifically, self-broadening, non-Voigtian line shapes, and collision-induced absorption (CIA) are now supported. One may reference the ROCKE-3D SOCRATES documentation from the GISS Astrobiology website13 for more information on other possible radiation tables generated for other atmosphere-star combinations of interest.
Table 2. The Shortwave Bands Adopted for Planets with Earth-like Atmospheres Orbiting Sun-like Stars
| Shortwave band | Wavenumber ( | Wavelength ( |
|---|---|---|
| 1 | 31,250 to 50,000 | 200 to 320 |
| 2 | 19,802 to 31,250 | 320 to 505 |
| 3 | 14,493 to 19,802 | 505 to 690 |
| 4 | 8403 to 14,493 | 690 to 1190 |
| 5 | 4202 to 8403 | 1190 to 2380 |
| 6 | 1000 to 4202 | 2480 to 10,000 |
Note. These are the bands used by the UK Met Office for Global Atmosphere Configuration 7.0 (GA7.0, D. Walters et al. 2017, in preparation).
Download table as: ASCIITypeset image
Rayleigh scattering by air is included, and we have also implemented a new Rayleigh scattering formulation that calculates the Rayleigh scattering coefficient consistently with the atmospheric composition in each layer. Water cloud optical properties are derived using Mie scattering, whereas the parameterization of ice crystals, as described in Edwards et al. (2007), is based on the representation of ice aggregates introduced by Baran et al. (2001). Vertical cloud overlap is treated using the mixed maximum-random overlap assumption (clouds in adjacent layers overlap maximally, whereas clouds separated by one or more clear layers overlap randomly).
In order to improve the accuracy of the calculated long- and shortwave fluxes, wavelengths are weighted internally in each band by the thermal Planck function for the longwave component, and by the stellar spectrum for the shortwave component when derivingk-coefficients and aerosol and cloud optical properties. This is important, as our bands are quite broad, particularly for the shortwave component, and the source function varies significantly within the bands. Consequently, changing the stellar spectrum involves computing newk-coefficients and cloud and aerosol optical properties for use in the calculation of short-wave fluxes.
In summary, SOCRATES gives us greatly improved flexibility to model atmospheres with different composition and irradiation than present-day Earth. However, due to the desire to keep the computation time as small as possible while calculating accurate fluxes and heating rates, some adaptation is needed for each planet and star.
4.2. Convection and Clouds
The cumulus parameterization in Planet 1.0 uses a mass flux approach that requires both instability and a trigger based on the buoyancy of moist air lifted a finite distance to initiate convection (version “AR5” in Del Genio et al.2015). In this sense, it is more resistant to convection than parameterizations in other planetary GCMs that require only instability to initiate convection (e.g., Song et al.2013). The mass flux scheme utilizes a cloud model that simulates the thermodynamic, dynamic, and microphysical properties of air rising in convective updrafts. The depth of convection is determined by the distance the updraft penetrates above its level of neutral buoyancy before the diagnosed convective updraft speed goes to zero. The initial mass flux is calculated as that required to produce neutral buoyancy at cloud base, with entrainment increasing the mass flux at higher levels, and detrainment decreasing the mass flux above the level of neutral buoyancy. Simultaneous subsidence of the grid-scale environment that adiabatically warms and dries the gridbox to maintain subsaturated conditions, and convective downdrafts formed from negatively buoyant mixtures of updraft and environmental air, compensate the parameterized updraft mass flux. This approach differs from adjustment schemes that seek to maintain a specified (sometimes saturated) humidity profile and a moist adiabatic lapse rate. The differences in the mass flux and adjustment approaches may affect estimates of the inner edge of the habitable zone (Wolf & Toon2015).
Entrainment of subsaturated environmental air limits convection depth, but the next generation of ROCKE-3D will include stronger entrainment that produces more realistic subseasonal variability, and cools and dries the tropopause region relative to that in Planet 1.0 (Del Genio2016). This may have consequences for estimates of water loss for warm planets. Likewise, the Planet 1.0 version transports condensed water upward too vigorously, relative to the next generation of ModelE2, although this makes little difference to reflected sunlight because these clouds are optically thick (Elsaesser et al.2017). For ROCKE-3D, we have relaxed a limit in the parent Earth ModelE2 that restricts convection top pressures to
.
Stratiform clouds in Planet 1.0 have subgrid cloud fractions that are diagnosed from local relative humidity and stability (Del Genio et al.1996 and updates described in Schmidt et al.2006,2014). This is typical of Earth GCMs that have been adapted to study exoplanets, but it differs from some GCMs developed solely for applications to other planets that require a gridbox to saturate before a cloud forms that fills the gridbox (e.g., the Mars GCM described by Haberle et al.2012 and the Explicit Planetary Isentropic Coordinate Model used for giant planet studies in Palotai & Dowling2008). This is potentially important in estimates of the width of the habitable zone, because the most important cloud feedback (and the one that differs most widely among models) in terrestrial climate change simulations is due to changes in cloud fraction (Zelinka et al.2016). Cloud water mixing ratios evolve prognostically, based on simplified versions of microphysical processes that are not easily generalized to treat cloud processes on planets with different gravity or atmospheric pressure. The next generation model will include a more explicit two-moment microphysics representation (Gettelman & Morrison2015) that can scale more easily to other planets.
ROCKE-3D is adjusted to planetary radiation balance using free cloud parameters that regulate the onset of fractional cloudiness in the free troposphere and boundary layer, and the rate at which small cloud liquid and ice particles are converted to rain and snow that precipitate from the clouds. The specific values of these tuning parameters for Earth are chosen to produce reasonably accurate surface temperatures, but no such observational constraint yet exists for exoplanets. Thus, multiple choices that can bring the model to balance at different temperatures are possible (Way et al.2015, see Figure1). Likewise, more exotic planet configurations (e.g., synchronously rotating, zero obliquity, etc.) may not come into balance at Earth free parameter settings, and are therefore adjusted as needed within reasonable ranges.
Only H2O convection and clouds are represented in the baseline ROCKE-3D model. For the next-generation version, we will allow for the possibility of condensates such as CO2 and CH4 that are important on Mars and Titan, as well as on exoplanets near the outer edge of the habitable zone (see Section5.4).
4.3. Planetary Boundary Layer
The planetary boundary layer (PBL) in Planet 1.0 is described in Schmidt et al. (2006). It is based on nonlocal transport of dry-conserved variables (virtual potential temperature and specific humidity). It includes a diagnosis of the turbulent kinetic energy profile based on large-eddy simulation studies, and uses the resulting profile to define the PBL depth. Cloud top sources of turbulence are not included, although the effect of enhanced mixing at the top of cloudy boundary layers is estimated by the cloud parameterization. The next-generation ROCKE-3D will incorporate a full moist turbulence scheme that transports liquid water potential temperature and total water mixing ratio, and includes cloud-top radiative cooling as a source of turbulence. Boundary layer clouds have largely been absent from discussions of exoplanet habitability to date, but considering that they are the largest source of uncertainty in Earth’s climate sensitivity (Zelinka et al.2016), they warrant more attention in exoplanet studies.
4.4. Cryosphere
The cryosphere in ROCKE-3D—encompassing areas of snow, land ice, and sea ice—has not been significantly modified from the modern Earth version of the GCM (Schmidt et al.2006,2014), so hexagonal ice (ice Ih) is the only natural phase of water ice represented. This means that Planet 1.0 is not yet capable of simulating the physical or spectral characteristics of water ices on worlds with very low surface temperatures (below 75 K; e.g., ice XI), and/or ice under pressures of
or more (e.g., ices II, III, and IX) (Bartels-Rausch et al.2012).
Snow may accumulate on any solid surface, including land ice and sea ice, as long as surface conditions are sufficiently cold. Total snow column depth is divided into two to three layers once the snow depth exceeds 0.15 m; as new snow is added to the uppermost layer, older snow is redistributed to the underlayers. Both heat and water are permitted to pass through the snow column and into the ground (soil) beneath. Areas with snow depths
are considered to have only patchy snow cover. Snow cover over land depends on local topography (Roesch et al.2001) and is never allowed to exceed 95% of the cell. If snow accumulates on top of either land or sea ice to a depth greater than one meter, the bottom of the snow layer is compacted to ice. Note that, under cold global conditions, snow mass may accumulate on land to such an extent that the mass of the ocean is noticeably reduced; however, this ocean mass reduction will not be detectable as a change in the global land/sea fraction.
Land ice has the simplest treatment of the three cryosphere components in ROCKE-3D. Where land ice is distributed as an initial input to the GCM, it consists of two layers that may change thickness in response to mass balance changes (accumulation minus sublimation and melting) induced by snow or rain. However, Planet 1.0 does not have dynamic land ice capabilities that would allow the footprint area or the defined elevation of an ice sheet to grow or shrink in response to forcings, or to affect ocean depth. A glacial melt parameterization permits the return of land ice mass lost as meltwater to the oceans, as well as calving of “icebergs,” into geographic-specific coastal ocean cells.
Sea ice extent and thickness may be prescribed for simulations with a specified SST ocean model (Section2.2.1, or as an initial condition for simulations with mixed-layer oceans (Section2.2.2) or dynamic oceans (Section2.2.3). For the latter two types of simulations, it is also allowed to develop as a consequence of other climate forcings on a world that initially has no sea ice. When sea ice is allowed to respond to other forcings, Planet 1.0 uses the thermodynamic-dynamic sea ice formulation described in Schmidt et al. (2006,2014) to control its formation and transport across the ocean surface. The formulation consists of four layers of variable thickness but fixed fractional height, each of which has prognostic mass, enthalpy, and salt content. Four sea ice surface types are included: bare ice, dry snow on ice, wet snow on ice, and melt ponds (Ebert & Curry1993). Melt pond mass accumulates as a fraction of surface runoff, decays on a timescale that depends on the presence or absence of current melting, and re-freezes when the temperature is below −10 °C. The sea ice thermodynamics in Planet 1.0 have not yet been tested with the internal brine pocket formation capability present in the modern Earth version of the GCM (Bitz & Lipscomb1999; Schmidt et al.2014), which provides a more realistic representation of energy conservation processes, and results in thinner equilibrium ice thickness as compared to non-energy-conserving ice models (Bitz & Lipscomb1999). This capability will be introduced in a future release.
Sea ice dynamics, especially important on synchronously rotating aquaplanets, is treated with a viscous-plastic formulation for the ice rheology that takes into account such factors as the Coriolis force, wind stress, ocean-ice stress, slope of the ocean surface, and internal ice pressure when calculating strain rates and viscosity for ice advection. Frazil ice (spicules or plates of new ice suspended in water) is allowed to form either under existing ice or in open ocean, as long as surface fluxes cool the water to the freezing point, given the local salinity; once formed, the frazil ice advects along with the previously existing sea ice. For a total sea ice thickness of five meters or less, leads (narrow linear fractures) are allowed to form as a result of shearing or divergent stresses. These leads can act as conduits for heat, moisture, and gas fluxes from the ocean below.
As a planet’s climate (and ultimately, its detectability via remote observations) can be highly sensitive to the snow and ice albedo parameterizations used in a GCM, we highlight here the key aspects of Planet 1.0's treatment of albedo in the cryosphere (see also Schmidt et al.2006).
The albedo of any snow or ice surface type is dependent on the zenith angle across all latitudes. The albedo of snow, whether it exists on bare soil, land ice, or sea ice, is also a function of age; on sea ice, whether snow is wet (in the presence of precipitation or melting) or dry, has an additional effect. Wet and aging snow are both less reflective than dry or new snow, respectively (see e.g., Table3). The snow masking depth (i.e., the depth of snow needed to completely counter the albedo properties of the underlying surface type) depends on the underlying surface, though sea ice and land ice are generally considered masked if covered by 0.1 m of snow.
Table 3. Surface Albedos of Various Sea Ice Surface Types in Different Spectral Intervals
| Surface Type | VIS ( | NIR1 ( | NIR2 ( | NIR3 ( | NIR4 ( | NIR5 ( |
|---|---|---|---|---|---|---|
| 330 to 770 | 770 to 860 | 860 to 1250 | 250 to 1500 | 1500 to 2200 | 2200 to 4000 | |
| Bare ice (min) | 0.05 | 0.05 | 0.05 | 0.050 | 0.05 | 0.03 |
| Bare ice (max) | 0.62 | 0.42 | 0.30 | 0.120 | 0.05 | 0.03 |
| Snow (wet) | 0.85 | 0.75 | 0.50 | 0.175 | 0.03 | 0.01 |
| Snow (dry) | 0.90 | 0.85 | 0.65 | 0.450 | 0.10 | 0.10 |
| Melt pond (min) | 0.10 | 0.05 | 0.05 | 0.050 | 0.05 | 0.03 |
Download table as: ASCIITypeset image
The albedo values for sea ice are area-weighted averages for the different surface types, resolved in six spectral intervals (Table3). Melt ponds, which significantly reduce sea ice albedo, are parameterized using a pond fraction and depth that varies as a function of melt pond mass. Bare ice albedo increases between assumed maximum and minimum values as the square root of ice thickness.
Although the surface albedo is resolved spectrally into the six bands, which offers some sensitivity to different spectral irradiance, extinction of radiation with depth through snow, ice, and liquid water currently only distinguishes a VIS (290 to 690 nm) band and one NIR (690 to 1190 nm) band (see Figure2). The extinction for each medium assumes solar-type surface irradiance fractions in these bands and that radiation
is not transmitted. Therefore, this solar assumption will later be revised to capture sensitivity to alternative stellar spectral irradiances. Land ice is assumed to have a spectrally uniform albedo of 0.8. Land ice set as a boundary condition for non-modern Earth simulations may use a different albedo value as a default.
Planet 1.0 does not have the ability to interactively change surface types (e.g., from ocean to land or land to ocean). It therefore cannot treat situations in which sea ice freezes to the ocean bottom along coastlines, or ocean mass decreases/sea level falls as snow accumulates on land, situations that can occur at high latitudes on low obliquity planets, or for low instellations. To avoid this, for ROCKE-3D experiments, ocean bathymetry is deepened along coastlines when needed (e.g., Way & Georgakarakos2017).
4.5. Chemistry
In simulating other planets, including early Earth, certain assumptions that are built into ModelE2 can become invalid, so updated or new parameterizations need to be developed for ROCKE-3D. This is especially the case for reduced atmospheres like those of Archean Earth, Titan, and probably Pluto. Currently, ModelE2 and ROCKE-3D are able to run with interactive gas-phase chemistry and a number of different aerosol configurations, from simple bulk parameterizations to full aerosol microphysics calculations. Simulating ice and gas giant atmospheric chemistry is beyond the capabilities and scope of ROCKE-3D. The implementation of an automated solver of chemistry, which will essentially allow the simulation of any atmospheric composition regardless of its redox state, is underway. This involves the use of the kinetic pre-processor (KPP; Sandu & Sander2006), which will replace the scheme described below, and is expected to gradually become available in ModelE2 and ROCKE-3D in coming years. Its adoption will enable the use of alternate chemical schemes for Earth and planetary science applications, facilitating the easy update and upgrade of the chemical mechanisms currently in the model.
The current chemical scheme in ModelE2 only allows for calculations of O2-bearing atmospheres, and is strongly linked with the expected composition of the present-day atmosphere of Earth. The model uses the CBM-4 chemical mechanism (Gery et al.1989), which explicitly resolves the inorganic chemistry that involves NOx and O3, as well as the chemistry of methane and its oxidation products (Shindell et al.2001). In addition, it resolves the chemistry of higher hydrocarbons via a highly parameterized scheme, based on CBM-4, with only minor changes (Shindell et al.2003), and that of halogens in the stratosphere to account for the ozone hole (Shindell et al.2006). The solution of the chemical system is facilitated by the use of chemical families, which assumes that dynamic equilibrium will be established among species that are very closely related and interchange extremely fast, like the Ox family (O(3P), O(1D), O3), the NOx family (NO and NO2), and the HOx family (OH and HO2), which allows an accurate solution of the system that includes species with lifetimes from sub-seconds to months or even years.
The parameterizations of chemistry involve the solution of the chemical kinetics equations, which is (in principle) independent of conditions. However, some assumptions are made to make the solution of the partial differential equations less stiff, which should not be violated in a different atmospheric composition configuration. One of the most important assumptions is that molecular oxygen is always in excess, so reactions that involve it happen instantaneously. This is the case for the hydrogen radical and all alkyl radicals in the model:
where R is any alkyl radical with more than one carbon atom. These reactions are extremely fast (Burkholder et al.2015), and the assumption that they dominate other competitive loss processes of H, CH3, and R is valid for extremely low O2 levels. For H, the competitive processes would be reactions with O3 or HO2, both of which will go down with reduced levels or O2, whereas for the alkyl radicals, the competitive processes would be reactions with atomic O; O3 would also decrease in a low-O2 atmosphere. Even without the assumption that the levels of the competitive oxidants will go down, the reaction of O2 is still the dominant loss of these radicals for O2 levels as low as 10−6 or present atmospheric levels (PAL), based on their reaction rates alone (Burkholder et al.2015).
We performed ROCKE-3D simulations of present-day Earth under pre-industrial conditions (to minimize the impact of human activities on the atmospheric state), in which we varied the levels of atmospheric O2 from 1 to 10−6 of PAL, to study how chemistry will be impacted, with a focus on O3. We did not allow radiation to be impacted directly by the O2 changes, in order to study the chemical response alone, but the effects of the results in O3 were included. The summary of the simulations is presented in Figure3, which agrees very well with the results of Kasting & Donahue (1979). Interestingly, the calculated vertical profile of O3 for different O2 levels (Figure4) agrees with that of Kasting & Donahue (1979) only for O2 levels down to 10−3 PAL. The model calculates a collapse of the stratosphere for O2 levels below that threshold, because of the colder stratosphere that results from the decrease in O3, whereas the 1D model of Kasting & Donahue (1979) does not.
Figure 3. Global mean O3 column density as a function of O2 concentration.
Download figure:
Standard imageHigh-resolution image
Figure 4. Global mean O3 column profile as a function of O2 concentration.
Download figure:
Standard imageHigh-resolution image
4.6. Aerosols
Simulating aerosols prognostically in other planetary configurations is also possible with the GCM, with few modifications of the original scheme. For Earth applications, the model contains carbonaceous (primary and secondary organic, and black carbon) and non-carbonaceous (sulfate, ammonium, nitrate, sea salt, and dust) aerosols.
The carbonaceous aerosols can be formed either by direct emission in the atmosphere or by the oxidation of precursor volatile organic compounds. In an O2-rich atmosphere, organic hazes like those of the Archean, Titan, and Pluto cannot be simulated in Planet 1.0. This is a limitation of the gas-phase chemistry of the model, which cannot calculate the photochemical formation of condensables in a reduced atmosphere, rather than a limitation of the aerosol scheme. Methane is not able to form aerosols in oxidizing environments, so unless there is life to form higher hydrocarbons (organic aerosol precursors) or any combustible carbonaceous material that can inject organic and black carbon particles in the atmosphere via burning, the carbonaceous aerosols are not needed in non-Earth planetary configurations where O2 is present.
Non-carbonaceous particles are present on both Venus (sulfate aerosols, a SO2 oxidation product) and Mars (dust). Any planet with active volcanism is expected to have some level of sulfate aerosols in their atmosphere, whereas any planet with erodible bare rock is expected to have dust. In addition, any planet with surface saline water and wave breaking is expected to have sea salt aerosols injected into the atmosphere.
During most of Earth’s history, there were salty oceans covering parts of the planet, and the presence of an atmosphere ensures that waves would form, so the presence of sea salt aerosols should be considered ubiquitous from very early on. The GCM can interactively calculate sea salt aerosol sources in the atmosphere using a variety of parameterizations (Tsigaridis et al.2013); the default ModelE2 scheme is that of Gong (2003), which is a function of ocean salinity and surface wind speed over the oceanic grid cells. A parameterization that takes into account sea surface temperature is also available (Jaegle et al.2011), but is tuned toward present-day Earth conditions because there are no physical constraints on the process. Sea salt aerosols are able to run as a standalone component in the model, without requiring the presence of other aerosols, which could save significant computational resources in simulating ocean worlds.
As with sea salt, dust can be calculated interactively in ModelE2 (Miller et al.2006). The source function depends on surface type and orography, as well as wind speed. Topographic depressions tend to be good sources of dust, contrary to mountain tops and steep slopes. For experimental unpublished ROCKE-3D simulations of dust on Mars, we used MOLA14 topography data to construct a map of preferred dust sources based on the topography of the planet, similar to what Ginoux et al. (2001) did for Earth. The fraction of dust available for wind erosion implies that valleys and depressions have accumulated dust, compared to flat basins where dust is more homogeneously distributed. This is represented by calculating the probability to have accumulated sediments,Si, in a
resolution, as a function of the minimum (
) and maximum (
) elevations of the surrounding
topography of the grid boxi, which has an elevation ofzi:
The calculated accumulated sediment probability is shown in Figure5.
Figure 5. Calculated accumulated sediment probability.
Download figure:
Standard imageHigh-resolution image
Active volcanism is another way to form aerosols in the atmosphere of the planet. Volcanic eruptions inject large amounts of ash and SO2 into the atmosphere, among other constituents. Ash, which we are not yet able to simulate in the model, is absorbing and the particles are usually large (and thus have too short a lifetime to be globally important), but SO2 can form sulfate particles in an oxidizing atmosphere, like those of Earth, Venus, or Mars. Although the model is able to simulate both the formation of sulfate from SO2 and the lifetime of sulfate particles in the atmosphere, the sources of SO2 from active volcanism on other planets are virtually unknown, making its interactive simulation difficult. The sulfur cycle on Mars will be studied in the future.
4.7. Model Performance
As stated in Section2, the default version of ROCKE-3D runs at
horizontal resolution for atmosphere and ocean. The vertical resolution of the atmosphere is either 20 or 40 layers, but the number of ocean layers is fixed at 13 layers. It is possible to use fewer ocean layers, if desired. The computation cost per simulated year decreases by a few percent with fewer ocean layers; the major impact of a shallower ocean is the smaller thermal inertia, which reduces the time for the climate to equilibrate.
Using the GISS radiation scheme with a 20-layer atmosphere and 13-layer ocean, one can simulate 100 Earth years in approximately 21 hr of wall-clock time, using 23 cores. These simulations use one node/motherboard each, with two Intel Xeon E5-2697 v3 Haswell
each with 14 cores. The increase from 20 to 40 atmospheric layers increases the computational cost by approximately a factor of two, however, this can be partially offset by increasing the number of cores. Using the 40-layer atmosphere one can simulate 100 Earth years in approximately 24 hr of wall-clock time using 44 cores.15 Using SOCRATES in the GA7.0 configuration, ROCKE-3D is slower. With 40 atmospheric layers, the model is able to simulate approximately 100 Earth years in 37 hr of wall-clock time, using 44 cores.16 Compared to running the model at the same vertical resolution using the GISS radiation scheme, this is an increase in run-time of approximately 50%. Note that, as ROCKE-3D is only parallelized in latitude, the maximum number of cores that can be used with a horizontal resolution of
is 44.
Accurately treating gaseous absorption in atmospheres significantly different from the present-day Earth often requires an increase in the number of spectral bands in the radiation scheme. The speed of SOCRATES decreases with increasing number of bands and absorbers,17 and the run-time for any given planet will therefore depend on the radiation setup adopted.
5. Geophysical Properties
5.1. Land/Ocean Distribution and Topographic Relief
The land/ocean mask used by the GCM can be defined using a fractional or non-fractional method, with the former preferred at coarse resolutions when shorelines are irregular and where ocean gateways may have significant climate impacts. A fractional land mask scheme signifies that an individual grid cell can be defined not only as 100% land or 100% ocean, but also as some percentage of each. This means that coastal grid cells are treated by other routines in the model as some portion ocean and some portion land. Topographic relief then is a weighted average based on the elevation of the land fraction, and zero for the ocean fraction of the cell. Available input data sets for paleoclimate and other planets are described in AppendixA, and users may create their own.
5.2. Continental Drainage (Runoff)
Riverflow and continental drainage redistributes freshwater via the water cycle, and thus impacts soil moisture. This affects land temperatures and precipitation, as well as the ocean salinity distribution, which impacts ocean circulation. For simulations that use modern Earth land/ocean distributions, we use the same riverflow and drainage patterns as defined for modern Earth climate experiments. Drainage directions in the GCM for non-modern-Earth continental configurations must be assigned via a custom input file. We generate the new drainage patterns by examining the topographic elevation boundary condition array and, working inward from continental edges, we calculate the slope of each continental grid cell in eight directions (four sides, four corners). Runoff is then removed from each cell in the direction of maximum slope, tracing a route back to the coast. For coastal grid cells that have more than one border adjacent to an ocean grid cell, runoff crosses the coastal grid cell on the same trajectory as in the adjacent inland grid cell.
Lakes may be treated as pre-defined static features, or permitted to grow and shrink dynamically in response to rainfall. In dynamic mode, lakes may also develop in topographic lows, with drainage developing once the water level rises above the lowest edge of the lake basin. Where a drainage route does not already exist, excess lake water runs off by means of the local drainage patterns defined above.
In ModelE2, glacial ice melts directly from the Greenland and East Antarctic ice sheets, and enters the surface ocean in prescribed cells wherever the edges of the ice sheets coincide with continental edges. However, this can be adjusted for a given topography or other needs by specifying the geographic locations where freshwater (from ice melt) is distributed back into the oceans.
5.3. Ground Hydrology, Albedo, and Surface Vegetation
Ground hydrology employs a six-layer soil heat and water balance scheme (Abramopoulos et al.1988; Rosenzweig & Abramopoulos1997), and the approach for calculating underground runoff is described in Aleinov & Schmidt (2006). The surface energy and water balance algorithm calculates heat and water content on an explicit numerical scheme in the soil layers and vegetation canopies (if present) to predict temperatures and saturation. In the current implementation, the total soil depth is 3.5 m, with the boundary at the bottom impermeable to both heat and water. Surface spectral albedo is partitioned currently into the same six broad bands shown in Table3: (300 to 770 nm, 770 to 860 nm, 860 to 1250 nm, 1250 to 1500 nm, 1500 to 2200 nm, and 2200 to 4000 nm), in an area-weighted average for cover types including vegetation, bare soil (with regard to albedo, see below), snow, and permanent ice. For questions regarding extrasolar planets, the ground hydrology boundary conditions that must be modified include soil texture and albedo maps for bare soil on a lifeless planet, or albedo influenced by vegetation.
Surface fluxes of energy and water for the ground fraction of a grid (i.e., lakes and permanent ice are done separately) are computed as a sum over the surface types, for which subgrid bare soil and/or vegetated fractions can be prescribed. Within these cover fraction columns, fluxes from soil, snow, and vegetation are calculated as components of the column flux. Bare and vegetated surface fractions maintain separate sets of prognostic variables for water and heat content in the soil column and snow. There can be more than one vegetated surface fraction, and although different vegetation types can have significantly different energy and evapotranspiration fluxes, we do not resolve the water and energy balance between different vegetated fractions. Instead, the subgrid canopy conductance is an area-weighted average of that for all vegetated types, which then drives the total vegetated fraction water and energy balance. Therefore, all vegetated fractions see the same temperature and soil moisture. Future development of the GCM is expected to provide partitioned water and energy balances for different vegetated fractions, which will make some difference for paleoclimate and modern Earth studies. In general, the horizontal grid resolution of GCMs is at a scale much larger than that of actual soil heterogeneity, which can be at the sub-meter scale. “Soil moisture” as a quantity in GCMs is highly model-dependent, and not a directly measurable quantity (Koster et al.2009). GCM land surface hydrology models conserve water and energy, and are formulated to capture surface fluxes at the grid scale.
The land albedo is spectrally resolved in the same VIS and NIR bands as described in Table3. The albedo can be calculated in two different ways, which requires different sets of input files:
- 1.The Lambertian albedo scheme from Matthews & Goddard Institute for Space Studies (1984), for which an input file gives grid fractional areas of land cover types, including vegetation types and bare soil. Bare dry soil albedo is specified as a combination of “bright” and “dark” soil of albedo 0.5 and 0.0, respectively, so that their area-weighted averages gives the soil albedo. Soil albedo is spectrally flat and assumed to depend linearly on soil saturation, becoming halved for a completely saturated soil. This is the scheme used from 1984 through to Schmidt et al. (2014). These input files are suitable for Earth vegetation and the solar spectrum. If simulating Earth vegetation under other stars, users should revise the broadband albedos to account for different band irradiances from different stellar types.
- 2.A zenith angle-dependent surface albedo scheme described in Ni-Meister et al. (2010). If vegetation cover is prescribed, then maps of vegetation cover, vegetation height, maximum leaf area index, and soil albedo are separate input files. When ecological dynamics are turned on, vegetation cover does not need to be initialized, because the vegetation will grow and die according to climate interactions. The soil albedo map allows resolution of soil into the six spectral bands of Table3. End member broadband optical properties should be calculated to take into account different stellar spectral types, as described in AppendixA.2.
Surface life, particularly photosynthetic life, can strongly influence a planet’s surface properties like spectral albedo, as we know from the vegetation red-edge on Earth (Tucker & Maxwell1976; Kiang et al.2007b). The Ent Terrestrial Biosphere (Ent TBM) is the Earth dynamic global vegetation model (DGVM) currently coupled to ModelE2 (Schmidt et al.2014; Kim et al.2015), and currently providing land surface albedo, water vapor and CO2 conductance, and seasonal variation in leaf area, while competitive and disturbance dynamics are under development. Although it can be an interesting exercise to subject Earth vegetation18 with full ecological dynamics to conditions on another planet to see what survives, seasonality and physiological differences between plant types are based on close adaptations to the star-planet orbital configuration and climatic regimes, so it would be inappropriate to utilize the current Sun-Earth based plant functional types (PFTs) for extrasolar planets.
As part of proposed work for ROCKE-3D, an Exoplanet Plant Functional Type (ExoPFT) is being introduced to provide a “generic plant” for simulations of alien vegetation influences on exoplanets. This ExoPFT will simply “find the water,” i.e., provide surface life wherever the planet is habitable. This generic plant will be similar to C3 annual grasses currently in the Ent TBM, but will have easily modifiable physiological and optical properties to allow experimentation with the potential distribution of life over a planet’s land surfaces, its impact on the surface energy balance and surface conductance, and its possible detectability. The ExoPFT will be a simple, non-woody, vascular plant with roots to access soil water that simulates the very basic influences of vegetation on climate: surface albedo and water vapor conductance. To “find the water,” the ExoPFT will be parameterized to emerge and senesce according to the presence of water, with broad climatological tolerance, and user-specified leaf spectral albedo that enables one to investigate effects on the climate of photosynthetic pigments adapted to alternative parent star spectral irradiance (e.g., pigment spectral absorbance adaptations similar to those proposed by Tinetti et al.2006 and Kiang et al.2007a). This ExoPFT will be built within the platform of the Ent TBM.
The Ent TBM currently provides the vegetation biophysics and land carbon dynamics to ModelE2 (Schmidt et al.2006,2014). The ExoPFT will utilize the EntTBM scheme for vegetation conductance of water vapor and CO2, and leverage a new canopy radiative transfer model being added to the Ent TBM. In addition, the ExoPFTs phenology (timing of leaf-out and senescence) and growth scheme will introduce its water-seeking parameterization within the Ent TBM framework.
Plant photosynthesis is sensitive to the atmospheric CO2 surface concentration. Leaf conductance of water vapor, which is coupled with photosynthesis, is inversely proportional to surface CO2 concentrations. These sensitivities are represented in the Ent TBM biophysics via the well-accepted Farquhar-von Caemmerer photosynthesis model (Farquhar & von Caemmerer1982), coupled with the leaf stomatal conductance model of Ball et al. (1987) detailed in Kim et al. (2015). This inverse relation to CO2 is infeasible for an atmosphere with zero CO2, which would not be realistic for a planet with photosynthesis. Numerically, in the GCM, the lowest CO2 level recommended is 10 ppm. This is the CO2 compensation point where photosynthesis and respiration just balance each other. This is typical for C4 photosynthesis, a type of photosynthetic carbon fixation pathway that enables uptake of CO2 at lower atmospheric concentrations than the other common pathway, known as C3 photosynthesis. Coupling to the atmosphere currently relies on roughness lengths determined by the ground hydrology scheme for the GCM grid cell scale.
Scaling leaf conductance, as well as optical properties, to the canopy scale for the ExoPFT will be done with the new prognostic vegetation canopy radiative transfer scheme, the Analytical Clumped Two-Stream (ACTS) model (Ni-Meister et al.2010; Yang et al.2010). This model has recently been incorporated in the Ent TBM. The ACTS model depends on zenith angle, direct/diffuse partitioning of radiation, canopy structure,19 and end member spectral optical properties of foliage, soil, and snow. The prior canopy radiative transfer scheme described in Schmidt et al. (2006,2014) relies on prescribed seasonal canopy albedos by vegetation type with fixed seasonal Leaf Area Index (LAI) (Matthews1984), and is not a function of dynamic LAI. Therefore, it is not suitable for use with dynamically changing vegetation.
End member optical properties are summarized into the same six broad bands used for the ground hydrology (see Table3). Alteration of these band albedos must take into account the stellar spectral irradiance, particularly if otherwise investigating the same vegetation optical properties, but with different parent stars. For example, the ACTS Earth vegetation end member broadband spectra (leaf reflectance and transmittance) are derived from convolving hyperspectral leaf data with a solar surface irradiance spectrum at 60 degrees zenith angle (approximating an average over the illuminated face of the planet) with a U.S. standard atmosphere.
Ent TBM vegetation dynamics of phenology (seasonality) and growth have been tested at the site level for several Earth plant functional types (Kim et al.2015). The ExoPFT phenology will be parameterized simply to leaf out and senesce with the availability of water, without mortality or establishment drivers other than water (i.e., insensitive to the plant’s carbon reserves, because this balance is already poorly known for Earth plants). Ecological dynamics involving competition, fire disturbance, and establishment will not be necessary to drive vegetation cover change, because only one ExoPFT will represent vegetation, driven by water availability.
5.4. Variable Atmospheric Mass
Typically, the atmosphere contains one or more components that may condense/evaporate at the surface of the planet or within the atmosphere. One can distinguish three cases: (1) A dilute (small fraction of total atmospheric mass) condensable gas. This is the case for water vapor on modern Earth. (2) A single-component atmosphere that consists of a gas that condenses at typical temperatures and pressures. (3) A non-dilute (significant fraction of total atmospheric mass) condensable gas. In the first case, changes in atmospheric mass and heat capacity due to condensation/evaporation can be neglected—except in the cumulus parameterization, where the effects of water vapor and precipitation loading on parcel buoyancy are non-negligible. The processes at the surface in this case will typically be governed by turbulent diffusion fluxes.
Modern Mars, where CO2 accounts for most (but not all) of the atmospheric mass, is actually an example of case 3. For Planet 1.0, we have taken the first steps toward creating a Mars GCM by ignoring the minor constituents and treating Mars as a pure CO2 atmosphere, corresponding to case 2. In this case, the change in the atmospheric mass over the seasonal cycle may be significant, and cannot be neglected for calculating temperatures and pressure gradients. Also, the amount of condensable substance at the surface is abundant, so the process of condensation/sublimation is governed by the energy balance, rather than by the diffusion fluxes. In the remainder of this section, we present the algorithm we use to model the condensation of a condensable single-component atmosphere at the planet’s surface. The description of similar processes for a dilute condensable component in ModelE2 can be found in Schmidt et al. (2014).
We assume that the condensate is stored in the upper soil layer (recall that ModelE2 has six soil layers). We also assume that the formation of the condensate is controlled purely by energy balance, and the matter is condensed or sublimated as needed to compensate for energy loss or gain by the upper soil layer. Once formed, the condensate is assumed to stay at the condensation temperature
, which depends on the atmospheric pressureps at the planet’s surface. This temperature is described by the Clausius–Clapeyron relation. For the case of CO2 condensation on Mars, it can be expressed approximately via Haberle et al. (1982):
where
is in Kelvin andps is in millibars. We define the latent heat of condensationLc as the amount of heat needed to melt a unit mass of condensate and bring it to surface air temperatureTs
wherecpg andcpc are the specific heat capacities of the condensable substance in gaseous and condensed form, respectively. Here,Lc0 is the latent heat of condensation at some fixed temperatureT0. For CO2 condensation on Mars, we use:
The prognostic variable that defines the ground temperature and the amount of condensate stored in the first layer of soil is the amount of energy per unit area in this soil layerH1. The quantityH1 is defined with respect to some reference temperature
, in the sense that the substance at the temperature
has energy zero. In our model, we set
, because it helps in dealing with freezing/thawing of water in Earth simulations using ModelE2, but one can choose any reference temperature above the condensation point. The ground temperatureTg can be obtained as
where
is the volumetric specific heat capacity of soil and
is the thickness of the upper soil layer. If
then a non-zero amount of condensate is present, and its mass per unit area can be computed as
Because we are dealing with a non-dilute case, the atmospheric pressure is affected by the formation of the condensate, which can be expressed as
whereg is the gravitational acceleration. The heat contentH1 is controlled by the energy balance at the surface
whereRn is net absorbed radiation at the surface,S is the sensible heat flux to the atmosphere,G is the ground heat flux to the lower soil layers, and the last term on the right-hand side is the energy flux due to the gain/loss of the substance from/to the atmosphere (which is assumed to be at temperatureTs).
The algorithm described above is implemented as follows. At each time step,H1 is first updated according to Equation (9), with the assumption that there is no change in the amount of condensate, and the condition in Equation (6) is checked. If true, the system of Equations (3)–(9) is solved iteratively to obtain the new values for
,Tg,ps. The change in the condensate
over the time step is then computed as
where
is the amount of condensate at the end of the previous time step. Next,
is subtracted from the mass of the condensable component of the lower atmosphere layer. The removed/added gas is assumed to be atTs, so the temperature of the lower layer of the atmosphere is updated accordingly. If condition (6) is false, then no condensate is present. If condensate was present at the previous time step, then the corresponding amount of gas should be added to the lower atmospheric layer, and its temperature should be adjusted accordingly.
Figure6 shows the annual cycle of surface pressure on Mars at the location of the Viking 2 lander, and that simulated with the surface condensation routine activated in a version of Planet 1.0 that includes some (but not all) of the physics that affects Mars’ climate (i.e., it uses the GISS radiation scheme, which has limitations in treating atmospheres with composition very different from Earth, and does not yet allow for CO2 clouds or incorporate dust). Despite these limitations, the timing and amplitude of the seasonal variation (Sharman & Ryan1980) are in reasonable qualitative agreement with observations. The site of the Viking 2 lander was chosen for comparison, because it is largely a flat area and can be represented well by a coarse-resolution GCM cell, such as that used in this simulation. Most of the other landing sites have a more complicated terrain and would require higher horizontal resolution for such simulations, which is beyond the scope of our current experiments.
Figure 6. Annual cycle of Mars surface pressure, as measured by the Viking 2 lander (gray crosses) (Hess et al.1977; Tillman1989), and surface pressure simulated by ROCKE-3D (black solid line).
Download figure:
Standard imageHigh-resolution image
In the description of the algorithm above (for simplicity’s sake), we assume that only one (non-dilute) condensable component is present, but our model can also handle the presence of another (dilute) component. This is done by including the latent heat due to the dilute component in Equations (5)–(7), and including the dilute condensate heat capacity into
. Otherwise, the dilute component itself is treated as in ModelE2. The presence of both such components is necessary for a more representative Mars simulation where both dilute (water) and non-dilute (CO2) condensable components are present. We note that Mars’ atmosphere also contains several non-condensing minor constituents, e.g., N2. These are not important for the dynamics, but do affect the ability of CO2 to supersaturate, and thus the occurrence of CO2 cirrus clouds (Colaprete et al.2008). This capability does not yet exist, but will be added in future generations of ROCKE-3D.
Currently, ROCKE-3D does not have the capability to treat case 3, i.e., condensable constituents that represent a significant and variable fraction of the mass of a multi-component atmosphere. This can become important as an Earth-like planet approaches the inner edge of the habitable zone where H2O becomes a non-negligible part of the total atmospheric mass. This will in turn affect pressure gradients and the thermodynamic properties of air, and introduce non-ideal gas behavior. This feature will be added in a future generation of ROCKE-3D.
6. Enhancement to Earth System Modeling (ModelE2) as a Result of ROCKE-3D
Generalizations and extensions to ModelE2 to accommodate the requirements of non-Earth planets can also benefit the Earth model through accelerated implementation of previously planned user-facing improvements to flexibility and accuracy. With a view to future development work and its multi-planet scope, this process also provides an opportunity for restructurings that enhance programmer-facing “code quality.” A visible example of all these trends is the reorganization of the time-management system, discussed in Section3.
Other examples can be found in the modularization of the manner in which the features of a planet are specified by a user. ModelE2 had previously required the presence of input files associated with all surface types (which is inconvenient for desert worlds and aquaplanets), and the time-space distribution of radiatively active constituents that are important for Earth but not other planets (e.g., O3). Ongoing effort to increase the flexibility in the specification of inputs and boundary conditions for Earth runs was extended to cover additional use cases.
Improvements to ModelE2 accuracy can sometimes result from running its modules under conditions sufficiently different from those for Earth to expose inappropriate approximations and/or coding errors. In the first category, the performance of the GISS Long Wave radiation scheme under conditions of extremely low column water vapor (e.g., the Arctic and Antarctic) was improved via better look-up tables generated in response to reports of problems in a cold and dry non-Earth simulation. In the second category, an aquaplanet simulation revealed some oversights in the ocean horizontal diffusion of momentum.
Looking forward, an example of development planned for the Earth model that is also highly convenient for non-Earth simulation includes the option for dynamic surface-type masks due to factors including: sea-level change, sea ice that has thickened to the ocean bottom, and expansion/retreat of glacial ice. Although “transient” simulations of exoplanets in response to imposed time-varying forcings are not a likely near-term objective, and the trajectory followed by a model as it approaches equilibrium for a given set of imposed forcings is typically not of interest either, it is convenient to have a model find its equilibrium in a fully automated manner. A brute-force procedure requiring the user to try a sequence of prescribed land/sea distributions and associated inputs greatly slows the rate at which equilibria can be determined.
Another advance that will benefit the Earth model is the use of the kinetic pre-processor (KPP; Sandu & Sander2006) for interactive chemistry calculations in ROCKE-3D. Its adoption will enable the use of alternate chemical schemes for both Earth and planetary applications, facilitating the easy update and upgrade of the chemical mechanisms currently included in the model.
7. Appropriate Use of ROCKE-3D
Timescale: as this is a GCM that simulates dynamics at time steps of
s, and parameterized physics at time steps of 30 minutes (and less in some submodules), it is best used for scientific questions investigating time-slice equilibrium climate behavior at the scales of decades to centuries. The equilibrium time needed for ocean dynamics can take much longer (some simulations require thousands of years), but the climate characteristics are generally summarized over the last few decades of the run. In some rare instances, simulations tracking secular changes over thousands of years can be accommodated with this GCM (see Way & Georgakarakos2017 for examples). Geological timescale phenomena over millions of years, such as the changes in the carbonate-silicate cycle, cannot be simulated by a GCM, but time-slice atmospheric composition conditions or flux rates could be prescribed.
Atmospheric escape: the ROCKE-3D model top is at
(
for Earth), with 17 layers in the 40 layer model above the tropopause cold trap for Earth-like planets. This is sufficient to resolve the stratospheric general circulation, which becomes important for planets orbiting M stars in which significant shortwave absorption by water vapor occurs at high altitude (Fujii et al.2017). However, this altitude is tens of kilometers below the homopause, where photodissociation of species such as H2O and O2 becomes important. Thus, ROCKE-3D cannot directly simulate atmospheric escape processes; this would require coupling to upper atmospheric models specifically intended to simulate ionization and escape processes (e.g., Gronoff et al.2011). Furthermore, because ROCKE-3D (like all GCMs) can only simulate time slices of hundreds to thousands of years, it cannot be used directly to address problems of atmospheric evolution, such as water loss in moist greenhouse states near the inner edge of the habitable zone. Instead, GCM stratospheric water vapor mixing ratios are traditionally compared to the threshold first estimated for 1D models by Kasting et al. (1993) to characterize planets that may be at risk of significant water loss (e.g., Kopparapu et al.2016). However, more sophisticated approaches (e.g., Wordsworth & Pierrehumbert2013) may be feasible.
8. Discussion
The use of GCMs to study the climate and weather of other planets has increased dramatically in the past few years, in response to increased interest in the past climates of terrestrial solar system planets, the rapidly growing list of rocky and potentially habitable exoplanets, and the promise of more discoveries, as well as atmospheric characterization of exoplanets by upcoming and planned future spacecraft missions. Every GCM has specific strengths and weaknesses in its ability to simulate other planets, and limitations in the range of problems to which it can be applied. The Earth climate modeling community has found that, as a result, a diverse population of GCMs offers advantages over any single model, by revealing robust behaviors that are common to all models and appear to be determined by fundamental well-understood physics, as well as features that differ among models due to differing assumptions in the parameterized physics that highlight more poorly understood processes.
The advantages of ROCKE-3D, relative to other planetary GCMs, are that its physics is identical to the most recent published version of its parent Earth GCM, it will remain current with future generations of the Earth model, and its developers include a subset of the people who develop the Earth model. Thus, it includes much in the way of recent thinking about climate processes that operate to determine Earth’s changing climate, and its coding structure has been generalized to easily allow simulations in parameter settings appropriate to other planets without sacrificing process understanding. It will also be the first exoplanet GCM to represent basic functions of plants that should be generally applicable to any habitable planet (for mock observations based on GCM output, see AppendixB).
That having been said, ROCKE-3D's Earth heritage produces limitations on its use as well. Some of these are structural and cannot easily be modified. The most obvious is that ROCKE-3D is based on a model that is designed to simulate only shallow atmospheres and oceans (i.e., much thinner than the planet radius) with equations of state appropriate to such fluids. Thus, ROCKE-3D can be applied to planet sizes up to the super-Earth range, though not to “waterworld” planets on which water is a significant fraction of the planet mass and a transition occurs from water to ice at high pressure. Likewise, it cannot be used to simulate or predict the transition from super-Earths to sub-Neptunes with thick H2 envelopes, nor can it simulate giant exoplanets.
Other limitations are specific to the Planet 1.0 version of ROCKE-3D, and will disappear as future generations of the model are developed. Planet 1.0 has been applied thus far only to planets with atmospheres composed of constituents found on Earth at pressures equal to or less than that of Earth’s atmosphere and temperatures not too far from those present during Earth’s history, such as snowball Earth periods (Sohl et al.2015) and a hypothetical habitable ancient Venus (Way et al.2016). With the SOCRATES radiation scheme, it is now sufficiently general to handle non-oxygenated atmospheres with prescribed elevated greenhouse gas concentrations, such as Archean Earth and Earth-like planets orbiting M-stars (Del Genio et al.2017; Fujii et al.2017). It has also been run under variable eccentricity (Way & Georgakarakos2017) and rotation periods as slow as 256 days, as well as synchronous rotation. A baseline modern Mars model has also been created. Rotation periods less than Earth’s are also possible, but require the higher horizontal resolution version of the model to accurately capture the dynamics. In its current form, the model cannot simulate atmospheres near the inner edge of the habitable zone, because the radiation does not include information from high-temperature line lists, and also the model does not treat the effects on atmospheric mass, thermodynamics, and dynamics of water vapor concentrations that are a non-negligible fraction of the total atmospheric mass. Atmospheres with compositions fundamentally different from those mentioned above (e.g., H2-dominated) are not yet available, though this is only a matter of developing appropriate radiation tables for such planets. However, even in its current form, ROCKE-3D is well-suited to address a wide range of science questions about habitable and inhabited planets and should be a valuable tool for interpreting near-future spacecraft observations of planets, both within and outside the solar system, and for supporting the planning of a possible, future, direct-imaging exoplanet mission.
We thank James Manners for providing access to SOCRATES and advice on coupling it to ROCKE-3D.
This research was supported by the NASA Astrobiology Program through our participation in the Nexus for Exoplanet System Science, and by the NASA Planetary Atmospheres Program, Exobiology Program, and Habitable Worlds Program. We also acknowledge internal Goddard Space Flight Center Science Task Group funding that triggered the initial development of ROCKE-3D along with substantial assistance from Shawn Domagal-Goldman and Gavin Schmidt. Resources supporting this work were provided by the NASA High-End Computing (HEC) Program through the NASA Center for Climate Simulation (NCCS) at Goddard Space Flight Center. The Viking 2 lander VL1/VL2-M-MET-4-BINNED-P-T-V-V1.0 data set was obtained from the NASA Planetary Data System.
The authors would also like to thank our referee Eric Wolfr for comments that greatly improved the manuscript.
Softwares: SOCRATES is released under a three-clause BSD license. ROCKE-3D is public domain software and under the process of obtaining a NASA Open Source Agreement.
Appendix A: GCM Inputs
Users should consult the GISS GCM Software page for user guides on ModelE2, necessary input data sets, and how to run the model.20 For the Planet 1.0 branch, we offer the following data sets, tools, and guidance for simulating other planets.
A.1. Surface Pressure and Gas Amounts
The surface pressure and amounts of atmospheric constituent gases are provided as input to the model in the form of the total surface pressure
and number/volume gas mixing ratiosri, for each speciesi, respectively. The mean molecular weight of dry air can then be calculated as
where
is the number of atmospheric constituent species, andmi is the molar weight of each gas. As the gases are assumed to be ideal, the number mixing ratios are directly related to the partial pressures
of each gas at the surface:
To calculate the contribution of each gas to the total surface pressure, however, the molar weight of each gas needs to be taken into account. From hydrostatic equilibrium, the total surface pressure is given by
whereM is the total mass of the atmosphere,A is the surface area of the planet,g is the acceleration of gravity, andMi is the total mass of each constituent gas. For well-mixed gases
, where
is the mass mixing ratio of speciesi, and related to the number/volume mixing ratio by
. Consequently, we can write the total surface pressure as
and the contribution to the surface pressure of gasi is therefore given by
By comparing Equation (12) for the surface partial pressure and Equation (15) for the contribution to the surface pressure, it is clear that
only in the special case where
, i.e., where all atmospheric constituents have the same molar weight.
As an example, for an atmosphere composed of
and
, with number/volume mixing ratios
and
, and a total surface pressure
, the partial pressure of each gas at the surface is given by Equation (12):
The mean molecular weight is given by Equation (11):
whereas the contribution of each gas to the total surface pressure is, using Equation (15),
In this example,
and
. In other words, the contribution of N2 to the total surface pressure is smaller than its partial pressure at the surface, whereas the contribution of
to the total surface pressure is larger than its surface partial pressure. This may seem counterintuitive, but is explained by the fact that some of the partial pressure of each gas at the surface results from the weight of all gases.
The above discussion shows that giving gas amounts in units of pressure is ambiguous, and should always be accompanied by a statement specifying whether it is the partial pressures of each gas at the surface or their contributions to the total surface pressure to avoid confusion. For this reason, we prefer to specify gas amounts in terms of total surface pressure and number/volume mixing ratios, as this is both unambiguous and the input required by ROCKE-3D.
A.2. Stellar Spectra
When using the default ModelE2 radiation scheme, instead of SOCRATES, with alternative stellar spectra, a software tool is available to bin and format any high-resolution stellar spectrum for input to the GCM. A Python script is also provided to plot the outputs, with comparisons to the present day Sun. The source spectrum should cover the range 115 to 100,000 nm, and the integral must provide the total stellar flux at any arbitrary known distance from the stellar surface. The software tool partitions the source spectrum into a 190-bin spectrum covering this range. Fluxes outside this range are integrated over wavelength, divided by the bin size of the closest end bin, and added to the flux of this bin. Consequently, the end bins will have a slightly higher flux to account for the tail fluxes outside the bin range. The binned output spectrum is utilized for specification of the stellar constant at the top of the atmosphere of the planet, ultraviolet fluxes to calculate absorption by ozone, and off-line calculation of another input file used for photolysis rate estimation. Additionally, a 16-bin array is output with the fractions of total stellar irradiance per bin for shortwave energy balance calculations. Detailed instructions and the formatting tool can be downloaded from the NASA GISS ROCKE-3D software webpage.21
When using the default ModelE2 radiation scheme, users should be aware that this scheme will have biases with redder stars, and may produce significant errors in fluxes and heating rates for planets with water vapor or other gases that absorb in the near-infrared. In such cases, one would be well-advised to consider using the SOCRATES radiation scheme.
A.3. Topographic Reconstructions
We have made reconstructions of land/ocean distributions and topographic relief for Earth, Venus (Way et al.2016), and Mars, based on the best available resolution digital topographies, with optional ocean coverage dependent on the choice of water depth (see Figure7).
Figure 7. Examples of continental configurations utilized in ROCKE-3D simulations. Top row: modern Earth (left), Cretaceous Earth (100 Mya). Bottom row: Neoproterozoic Earth (715 Mya) (left), paleo-Venus with oceans.
Download figure:
Standard imageHigh-resolution image
Reconstructions of paleo-Earth for the past one billion years use geologic and geophysical information, consisting primarily of paleomagnetic data to determine continental positions, based on reconstructed paleo plate positions. Following plate tectonic reconstruction, sedimentological and paleontological evidence supplies details about shoreline location, which can differ substantially from the continental boundaries, depending on the depth of ocean water. It is worth noting that even minor changes in ocean depth and land/ocean distribution in the GCM can have significant impacts on the planet’s climate, because the opening and closing of ocean gateways, or the orientation of mountains, can have large impacts on the resulting circulation of the oceans and atmosphere, and therefore on the transport of heat and moisture. Topographic relief for paleo-Earth continents is also based on depositional environment reconstructions using sediment and fossil distribution, as well as considering tectonic settings that arise from plate interactions along subduction zones, continental rifts, and continent-continent collision zones. Bathymetry is based on similar evidence, but due to a lack of larger-scale areas of preserved ocean crust in the geologic record, reconstruction of bathymetry (e.g., Xu et al.2006) for time periods older than 180 to 200 million years is not possible.
Appendix B: Post-processing for Exoplanet Mock Observations
B.1. Disk-integrated Spectra and Light Curves
A useful application of the model outputs is the prediction for what the planet would look like if it were an exoplanet. A set of external Python codes are provided for generating the disk-integrated light curves (both reflected and thermal), as if it were observed from an astronomical distance, given the externally specified direction of the observer. The program reads top-of-atmosphere radiation diagnostics (short-wave and long-wave) from the model outputs, and integrates the outgoing top-of-atmosphere fluxes in each radiation band from each pixel over the planetary disk, taking account of the relative configuration of the planet and the observer. Given that GCMs are regularly run with a small number of spectral bands for computational efficiency, users who want to create a higher-resolution spectrum need to run the model for a short amount of time with a larger number of bands after the model reaches an equilibrium state.
An isotropic radiation field is assumed in the current scheme. Thus, any deviation from it (e.g., due to the strongly anisotropic nature of scattering by clouds) is not included. The radiation diagnostics used to calculate these synthetic observations are output as monthly means. Consequently, the effects of temporal variation of cloud cover on timescales shorter than the output frequency are lost.
Figures8 and10 show the annually averaged, disk-integrated, albedo and thermal emission spectra of the Earth, whereas Figures9 and11 present their variations over one orbit. The albedo plotted here is the “apparent albedo,” defined as the observed reflected intensity divided by the reflected intensity of a loss-less Lambert sphere at the same phase. The light curves assume an orbital inclination of 90° (i.e., the planet is on an edge-on orbit), 234 obliquity, and that the vernal equinox for the northern hemisphere coincides with inferior conjunction (e.g., the planet is between the star and the observer). In order to show the diurnal variations together with the yearly variations in one panel, the spin period is artificially set to 100 hr and the orbital period is set at one Earth year in Figures9 and11; diurnal variations with 24 hr periodicity would be smeared out when the horizontal axis spans 365 days (see Figure 3 in Fujii & Kawahara2012 for an example of a 24 hr period). All of the parameters discussed above can be modified to fit a particular planet.
Figure 8. The disk-integrated apparent albedo of the Earth as a function of wavelength. Based on ROCKE-3D outputs simulating present-day Earth.
Download figure:
Standard imageHigh-resolution image
Figure 9. Orbital variation of the disk-integrated scattered light of the Earth based on ROCKE-3D outputs simulating present-day Earth. In order to see the diurnal variations clearly, the spin period is set at 100 hr while the orbital period is one Earth year. Vernal equinox for the northern hemisphere is located at inferior conjunction.
Download figure:
Standard imageHigh-resolution image
Figure 10. The disk-integrated thermal emission spectrum of the Earth as a function of wavelength. Based on ROCKE-3D outputs simulating present-day Earth.
Download figure:
Standard imageHigh-resolution image
Figure 11. Orbital fractional variation of the disk-integrated thermal emission of the Earth, based on ROCKE-3D outputs. In order to see the diurnal variations clearly, the spin period is set at 100 hr and the orbital period is one Earth year.
Download figure:
Standard imageHigh-resolution image
B.2. Transmission Spectra
A second set of external Python codes compute the transmission spectrum of a planet based on output from ROCKE-3D. The code reads the atmospheric columns located near the terminator, interpolates them onto the terminator, and computes the transmission of the starlight based on these profiles. Examples of required input parameters are: location of the GCM outputs, the composition of the background atmosphere, the cross-section tables to use, wavelength resolution, the radius of the star and the planet, and the viewing geometry of the star, planet, and the observer.
We consider polar coordinates
on the spherical plane centered at the planetary center, as illustrated in Figure12. Denote the spectral extinction optical depth along the ray that exits the planetary atmosphere at
by
, the fraction of the intensity absorbed or scattered along the optical path that exits the atmosphere at the altitude betweenb andb +db, and at the angle betweenθ and
,
, is then
Using
, the transit depth,
, is given by
where
is the impact parameter, below which the ray is completely attenuated, and
is the impact parameter, above which the planetary atmosphere may be regarded as transparent. Transmission spectra are also regularly represented by an “effective height”
or “effective radius”
, which is given by
Figure 12. Geometry of transmission spectroscopy.
Download figure:
Standard imageHigh-resolution image
The trajectory of the transmitted ray seen by the observer is traced from the direction of the observer toward the stellar disk, accounting for the refraction due to the planetary atmosphere as described in Misra et al. (2014) and van der Werf (2008). The index of refraction of the atmosphere is computed based on the number density of the atmosphere’s constituents. If the ray that transverses at a certain altitude of the atmosphere does not intersect the stellar disk, it is observed as opaque (Bétrémieux & Kaltenegger2014; Misra et al.2014).
The opacities along the ray trajectories,
in Equation (21), are computed and take into account both gaseous Rayleigh scattering and absorption. The cross-sections of atmospheric molecules are calculated based on HITRAN 2012 (Rothman et al.2013), with both Doppler broadening and pressure broadening by air accounted for in the Voigt profile of each line. The cross-section data are tabulated at temperatures between 100 and 400 K in steps of
and pressures (equally spaced logarithmically in intervals of one order of magnitude in millibar units), that will be interpolated to obtain the absorption coefficient at each location along the trajectory. The cross-section tables for O2, O3, H2O, CO2, CH4, and N2O are provided, whereas those for other molecules can be created from the corresponding HITRAN data with the provided codes. The opacity due to cloud particles may be included in a simplified manner. GCM outputs of the cloud properties are typically averaged over one month, which is significantly longer than the timescale of cloud formation, but it is possible to obtain averages over intervals as small as 30 minutes using the ModelE2 Sub-daily (SUBDD) facility.22
We note that the transmittance based on the time-averaged optical properties (which we calculate) may be somewhat different from the time-averaged transmittance based on the instantaneous optical properties (which is observed).
Figure13 displays examples of the cloud-free transmission spectra based on an Earth GCM simulation, with and without the effect of refraction in the planetary atmosphere. In demonstrating the effect of refraction, we assumed that the planetary center coincided with the center of the stellar disk. These assumed geometrical parameters can easily be modified to fit any particular planet. The results shown in Figure13 are consistent with previous literature results (e.g., Kaltenegger & Traub2009; Bétrémieux & Kaltenegger2014; Misra et al.2014).
Figure 13. Cloud-free transmission spectra based on ROCKE-3D experiments simulating present-day Earth. Results are shown both with the effect of refraction (thick black line) and without (thin gray line). A stellar radius of
and planetary orbital distance of
were used.
Download figure:
Standard imageHigh-resolution image
