1 Introduction

Earth's near-future climate has been extensively explored via the IPCC and associated CMIP studies (e.g., Collins et al., 2013). Earth's ancient climate has also been studied at various levels of detail, including the Cretaceous greenhouse (e.g., Huber et al., 2018), the Neoproterozoic Snowball (Pierrehumbert et al., 2011), and on the supercontinent Pangea (e.g., Dunne et al., 2021; Parrish, 1993). Some authors have explored Earth's deep time future climate by looking at increases in , solar insolation through time (e.g., Sagan & Mullen, 1972) or looking at the future carbon cycle (e.g., Franck et al., 1999). Yet few have investigated climate effects induced by additional changes in topography and land/sea masks (e.g., Davies et al., 2018).

The geological formations on the ever-changing surface of the Earth have a strong influence on our climate. The transition to a cold climate in the Cenozoic, including the glaciation of Antarctica, was induced by the opening of ocean gateways and reduced atmospheric concentrations (Barker, 2001; DeConto & Pollard, 2003; Smith & Pickering, 2003). The development of the Caribbean arc and closing of the Panama Isthmus allowed the Gulf Stream to form, with major consequences for global climate (Montes et al., 2015), whereas the closure of the Strait of Gibraltar led to the Messinian Salinity Crisis (Krijgsman et al., 1999). Furthermore, the Himalayas, a consequence of the India-Eurasia collision, allows for the monsoon (Tada et al., 2016). Recently, Farnsworth et al. (2019) showed that the climate sensitivity for the period 150–35 million years ago is dependent on the continental configuration, particularly ocean area. Schmittner et al. (2011) investigated the effects of mountains on ocean circulation patterns of present day Earth and concluded that the current configuration of mountains and ice sheets determines the relative deep-water formation rates between the Atlantic and the Pacific Oceans.

The continents on Earth aggregate into supercontinents and then disperse on a cycle of 400–600 million years–the supercontinent cycle (Davies et al., 2018; Pastor-Galán et al., 2019; Yoshida, 2016; Yoshida & Santosh, 2018). The latest supercontinent Pangea formed around 310 million years ago and started breaking up around 180 million years ago. The next supercontinent will most likely form in 200–250 million years, meaning Earth is currently about halfway through the scattered phase of the current supercontinent cycle (Davies et al., 2018).

There are obvious and strong links between large-scale tectonics and climate. It would be interesting to know what Earth's climate could be like in the distant future when continental movements will have taken Earth away from the current continental configuration (Davies et al., 2018). Here, we investigate what a climate may look like on Earth in a future supercontinent state. In particular, we focus on changes to land/sea mask, topography, rotation rate and insolation. We do not delve into details of the future carbon cycle or speculations about changes to the Earth's biosphere or atmospheric constituents into the deep future, we keep the latter near modern values. A secondary application of climate modeling of the deep-time future is to create a climate model of an Earth-like exoplanet using the parameters known to sustain habitability and a stable biosphere (Earth). Using the Deep-time future Earth as a basis for exoplanetary climate studies allows us to establish sensitivity ranges for the habitability and climate stability of the future Earth and its distant cousins in our galaxy.

2 Methods

2.3 General Circulation Model Set up

The ROCKE-3D General Circulation Model (GCM) version Planet_1.0 (R3D1) as described in Way et al. (2017) is used for this study. A fully coupled dynamic ocean is utilized. Using data generated via Claire et al. (2012) we use an insolation value of 1361 1.0223 = 1391.3 W for the Amasia simulations (04–06) 200 Myr into the future. We use a value of 1361 1.0260 = 1396.4 W for the Aurica simulations (01–03) 250 My into the future. We do not change the solar spectrum as the changes for such a small leap into the future will be minimal in terms of its effect on the planet's atmosphere.

We use a 50/50 clay/sand mix for the soil given that we have no constraints on what the surface will be like in the deep future and is a value commonly used in the exoplanet community (e.g., Yang et al., 2014; Way et al., 2018). In a 3D-GCM the soil is important for its albedo and water holding capacity, see Section 2 of (Del Genio et al., 2019) for details on the latter. 40 cm of water is initially distributed into each soil grid cell. We use a ground albedo of 0.2 at model start, but the albedo will change via snow deposition (brighter), or from rainfall (darker) as the GCM moves forward in time.

The original topography resolution of 1/ from the tectonic maps discussed in Section 2 is down-sampled to a resolution of in latitude by longitude, which is the default R3D1 resolution. The standard deviation from the down-sampling is used to set the roughness length of the surface in each grid cell. River flow direction is based on the resulting topography and exits to the ocean when possible. Large inland seas (typically less than 15 contiguous grid cells) are defined as lakes rather than ocean grid cells. The GCM allows lakes to expand and contract as dictated by the competition between evaporation and precipitation. The same holds for the possible creation and disappearance of lakes. This allows the model to handle inland surface water in a more sophisticated manner than making all surface water defined as ocean grid cells. This is highly desirable because ocean grid cells cannot be created or destroyed during a model run.

Any ocean grid cell with a depth less than 150 m (from the down-sampled data) was set to have a value of 204 m (the mean depth of ocean model level 6). This is especially important at high latitudes where shallow ocean cells may freeze to the bottom causing the model to crash due to its inability to dynamically change surface types from ocean to land ice.

The down-sampling has a side effect in that the land-sea mask will differ slightly between the three topographic types (CTRL, PD, MTNS). For example, in a case with a collection of ocean or lake grid cells adjacent to a number of high elevation land topography grid cells the down-sampling may change the combined ocean + land grid cells into a land grid cell, or vice-versa if the mean depth of the ocean grid cells is larger than the height of the land grid cells. This is why the land/sea masks differ between CTRL, PD and MTNS in Figure 1, even though their 1/ parents had exactly the same land-sea mask.

One side-effect of having quite distinct land elevations and a lack of oceans in polar regions in the Amasia runs (sims 04–06) is that snow accumulation can result in the growth of ice sheets akin to that of Earth's last glacial maximum (LGM) when the Earth was cooler than present day (Argus et al., 2014; Peltier et al., 2015). The increase in ice sheet height can influence the climate as there may be substantially more snow accumulation at higher elevations, whereas rain would normally fall at lower elevations, due to differences in the lapse rate. Because R3D1 does not have a dynamic ice sheet model we adopt the following approach to deal with these snow accumulations. To accommodate the possibility of such ice sheets we ran models with the original Amasia topography (sims 04–05) and allowed snow to accumulate unhindered. Once these runs reached equilibrium we then used these snow accumulations as the bases for modified production runs. Fifty years climatological averages of snow accumulation (see Figure 2 middle panels) over N. Hemisphere summer months (June, July & August) were used to increase the land elevations where necessary. We choose summer months since those minimum northern hemisphere accumulations work well to allow accumulation in the Fall/Winter months and evaporation in the Spring/Summer months. The same procedure is used in the southern hemisphere with 50 years climatological averages over the months of December, January & February. We then perform small areal averages over the highest latitudes to simulate the effect of ice sheet movement. These summer minima with snow accumulations are then labeled as permanent ice sheets (with appropriate albedo) in the model topography boundary condition files. An offline ice sheet model would be preferred as is typical in LGM studies (Argus et al., 2014; Peltier et al., 2015) but is beyond the scope of the present exploratory work. Figure 3 includes original topography plus snow accumulations (denoted as “with ice sheets” in red dotted lines) versus the original topography (blue solid lines). For comparison purposes Figure 3e over plots the LGM data from Argus et al. (2014); Peltier et al. (2015). Recall that the LGM was at a time of lower solar insolation and differing orbital parameters from our future Earth scenarios. We believe that Figure 3e with the LGM over plotted demonstrates that our approach to dealing with the ice sheets is not unreasonable. The south polar cap is reproduced with high fidelity, while the north polar cap (on average) also mimics the LGM pretty well.

The atmosphere is set to roughly Earth constituents in the year 1850: Nitrogen dominated with 21% Oxygen, 285 ppmv , 0.3 ppmv O, and 0.79 ppmv . No aerosols or Ozone () are included. For the minor species ( and ) this is perhaps the simplest choice given the variability in the past (e.g., Ramstein, 2011), and long-term uncertainties associated with human generated climate change and the subsequent uncertainties associated with the long-term evolution of the carbon cycle (e.g., Franck et al., 1999). For the second most abundant species in Earth's atmosphere () the choice is consistent with recent estimates by Ozaki and Reinhard (2021) who set a 1 limit of the longevity of Earth's 21% oxygenated atmosphere of ≈ years. For comparison purposes with related work (Way et al., 2018) we include a modern Earth-like land/sea mask in Earth #1—#3 (sims 07–09) (Table 1) with these same atmospheric constituents and a bathtub ocean. The Earth-like land/sea mask used in these simulations is described in Way et al. (2018) and shown in Figure 8 of that paper. These changes to the land/sea mask do not greatly affect the mean surface temperature and the bathtub ocean makes the model more resistant to crash conditions often associated with shallow ocean cells freezing to the bottom as would be likely in some of the cases herein. To better understand the possible effects of rotation rate and insolation (given such parameters used in the Aurica & Amasia sims 01–06) we take the same Earth #1 model (sim 07) and slow the rotation rate in Earth #2 sim 08 to be the same as in the Aurica and Amasia sims 01–06, and then increase the insolation in Earth #3 sim 09 to be the same as that of the Aurica sims 01–03 as shown in Table 1 (the higher of the two insolations used at 200 and 250 Myr into the future).

3 Results

Let's first attempt to disentangle any effects of the slower rotation rate. We do this by looking at the modern Earth #1—#2 (sims 07–08). Table 1 shows a minimal difference between the mean surface temperature between our Earth-like world with modern rotation rate (Earth #1 sim 07) and the 24.5 h rotation for Earth #2 sim 08 that is used by our Aurica and Amasia simulations (01–06). Planetary Albedo and snow + ice fraction are also nearly the same. In Figure 4a, visible high latitude regional temperature differences (∼C) are seen between Earth #1 & #2 (sims 07 & 08) even if mean difference is only 0.C.

Looking at Figure 5 (left panels) we see that Earth #1 & #2 (sims 07 & 08) also have very similar atmospheric, ocean and total meridional transport. If one compares the min and max stream functions in the tropics in Figures 6a and 6b (Earth #1 & #2, sims 07 & 08) the differences are small: −9.1 × /–9.2 × 1%, 1.2 × /1.19 × 1% (values are also noted below each figure).

We find very little evidence that the additional 30 min in the length of day has any appreciable effect on the climate dynamics. Work by Showman et al. (2013, Figure 5) has shown that pole to equator temperature differences should decrease as rotation rate slows. There is a marginal difference at high northern latitudes that in fact goes in the opposite direction (Figure 7a). With the slower rotating Sim 08 having a very small increase in equator-to-pole temperature difference. Note that the Showman et al. (2013) result is for much larger changes in rotation rate. Finally in Figure 7b we plot the eddy energy transport fluxes for Earth #1 & #2 (sims 07 & 08). One can see that the mid-latitude eddy energy flux in Earth #1 (sim 07) is slightly larger than that of Earth #2 (sim 08), which would be consistent with that of Showman et al. (2013), but again the differences are marginal.

Next the rotation rate is fixed at 24.5 h, but the insolation is increased from Earth #2 sim 08 (1361 = W ) to Earth #3 sim 09 (1361 1.0260 = 1396.4 W ). The differences are much clearer here with a ∼C difference in the mean surface temperature. The planetary albedo has decreased 0.5% which tracks the decrease in Snow + Ice fraction of 3%.

It should be noted that previous work has shown that some ancient Earth supercontinent phases, which are comparable to our Aurica simulations 01–03, have had more arid interiors where weathering effects and draw down may have been less efficient (e.g., Jellinek et al., 2019). This would increase surface temperatures as the balance of would tend to be larger than present day because volcanic outgassing (sources) would likely remain constant while drawdown (sinks) would decrease. However, there are other climatic effects to consider. For example, the Amasia reconstruction is essentially an arctic supercontinent with an independent and isolated antarctic continent, meaning both poles are covered by land, and much of that is covered by ice. Amasia is thus in essence a shift to consolidate the present day domination of northern latitude land masses even further north.

This increase in land masses at northern latitudes means that there is less ocean heat transport to melt the ice in the northern hemisphere summers as happens on modern Earth. Some of the heating differences can be seen in the middle right panel of Figure 5 where the oceanic meridional transport for the modern Earth #1—#3 simulations (07–09) is lower at lower latitudes than the Amasia simulations (04–06). This is because there are no southern low latitude continents (e.g., S. America or S. Africa) and the northern hemisphere continents are now pushed to higher northern latitudes in the Amasia runs. At the same time, in Figure 8 we see that there are active ocean currents in the modern Earth #3 sim 09 (bottom panels) near the northern polar regions (and in the Aurica sims at high latitudes - top panels), but none are possible in the Amasia PD sim 05 run (middle panels).

The lack of a northern polar ocean means that more ice resides on land and in lakes all year round near the north pole, as we see in present day Antarctica, for the three Amasia simulations (sims 04–06). This is the well known ice-albedo climate feedback and explains why the Amasia simulations tend to be cooler than the Aurica ones. Amasia PD (sim 05) is the coolest of the Amasia simulations. This is because its mean topographic height is higher (especially near the north polar regions) than in Amasia CTRL & MTNS (sims 04 and 06). See Figures 3e vs 3d and 3f. The higher relief means the Amasia PD (sim 05) lapse rate is lower on average and, as discussed in the Methods section above, it is cooler and hence instead of rainfall we tend to get snowfall at high latitudes. This fact is also born out in Figure 2 where grid snow + ice fractional amounts are quite high in the northern hemisphere winter months (top center) and southern hemisphere winter months (bottom center) in comparison with the modern Earth #3 simulation 09 with the same rotation rate and insolation. Note that Earth #3 (sim 09) has ice coverage on Greenland in the northern hemisphere summer. This is because we have not adjusted the height of Greenland assuming it no longer has an ice sheet, so it will accumulate snow and maintain it because of its higher altitude. In reality, it would likely not be snow covered at this higher insolation as its topographic height would surely be far lower, although one would also have to consider the effects of any land rebound height from the removal of the ice sheets.

It is informative to contrast Aurica PD (sim 02) with Amasia PD (sim 05). Aurica PD (sim 02) has land at lower latitudes and uses the same “present day” (PD) topographic height values for inputs as Amasia PD (sim 05) where the landmasses reside at high latitudes. In Table 1, we give their mean surface temperatures, planetary albedo, fractional snow & ice coverage and “Habitable Fraction.” The snow & ice coverage as illustrated in Figure 2 is clearly related to the planetary albedo and mean surface temperatures in Table 1. In Table 1, it is clear that the snow & ice fractions are much higher for the Amasia runs (04–06) compared to the Aurica runs (01–03), and highest for Amasia PD (sim 05) in particular. Amasia PD (sim 05) has the highest snow fraction amount corresponding directly to the lowest mean surface temperature of the Aurica and Amasia simulations (01–06). This coldest of the future climates Amasia PD (sim 05) is nearly C cooler than its corresponding modern Earth #3 simulation (09). We see a lower fractional snow + ice coverage for Earth #3 (sim 09) in Figure 2 versus that of Amasia PD (sim 05). This in turn is related to the fact that Earth #3 (sim 09) maintains open ocean at the northern pole which prevents the year round land ice seen in Amasia PD (sim 05) (see Figure 8). Hence, Amasia PD (sim 05) has 10.2% for the snow + ice versus a mere 6.4% for Earth #3 (sim 09) at the same rotation and insolation.

The general effect of the different land/sea masks between the Aurica (sims 01–03) and Amasia (sims 04–06) simulations and how they compare with the modern Earth #1—#3 simulations (07–09) are seen in Figures 5, 6 and 8. In Figure 5, the largest differences are seen in the oceanic meridional transport between the Aurica & Earth #1—#3 simulations. The weaker values seen for Aurica simulations (01–03) are likely explained by the large low latitude landmass restricting meridional heat transport over a large longitudinal range (left middle panel). In the right middle panel of Figure 5, we see how having a larger low-latitude open-ocean increases the oceanic meridional transport for the Amasia simulations (04–06) versus the modern Earth #1—#3 simulations (07–09). Total (atmosphere + ocean) meridional heat transport is very similar between simulations where the only discernible differences manifest themselves in the larger northern hemisphere transport for Earth #1—#3 versus the Aurica simulations, which are certainly related to the differences in oceanic transport as discussed above.

These general trends are repeated in Figure 6 where we plot the stream function which indicates the strength of the Hadley circulation. The Aurica PD (sim 02) stream function (Figure 6d) is the weaker of the three as we saw in Figure 5 (lower panels). Looking at Amasia PD (sim 05) versus Earth #3 (sim 09) the northern hemisphere values (Figures 6e vs 6c) are very similar, but the southern values differ likely because of the low–mid latitude south American, south African, and Australian continents in Earth #3 (sim 09) that do not exist in Amasia PD (sim 05).

Work by Spiegel et al. (2008) uses a metric of “climatic habitability” that defines the amount of surface area of a planet that can host liquid water (e.g., surface temperatures in the range 0 T C) at modern Earth atmospheric pressures. In the right–most column of Table 1 the left values are given using this metric, while the right values utilize a larger temperature range since life on Earth has been found to thrive in temperatures as high as C and as low as –C (e.g., NRC, 2007, Table 3.1). These metrics are calculated from 10 years averages (post-equilibrium) of the ground and sea temperatures. From Table 1 it is clear that the Aurica simulations (01–03) have the largest surface habitable fraction amongst all of the simulations. Since none of our simulations approach the boiling point of water in any region Aurica's high habitability is clearly due to the lack of high-latitude continents found in the Amasia and Earth simulations (04–09) that manifest below freezing temperatures not widely present in the Aurica ones (sims 01–03). Earth #1 & #2 (sims 07 & 08) have large areas with temperatures below freezing–not unexpected given their lower insolations and high latitude land masses compared to the Aurica simulations. The habitable fraction values for Amasia PD (sim 05) are lower than the Earth #1 & #2 simulations (07 & 08) at lower insolation. As noted above, this is attributable to the large ice sheets in the high latitude northern and southern hemispheres. Even though Amasia PD (sim 05) has a higher mean surface temperature than Earth #1 & #2 (sims 07 & 08) the higher global snow fraction appears to influence this metric more than may be expected. However, caution is warranted when using this habitability metric as other work (e.g., Sparrman, 2021) has shown that applying the Spiegel et al. (2008) temperature definition in a 3–D sense reveals little difference in “climatic habitability” between worlds that otherwise appear quite climatically distinct. On Earth life has been found to withstand pressures beyond those of deep sea trenches on Earth (e.g., Sharma et al., 2002; Vanlint et al., 2011), at the bottom of thick ice sheets (e.g., Griffiths et al., 2021) and in extremely deep mines (e.g., Drake et al., 2021; Lollar et al., 2019). Given enough time life has found a way to fill nearly every ecological niche on the modern Earth. While a habitability metric like that used herein may be imperfect it can still provide us a simple way to compare the surface climates of different worlds.

4 Conclusions

The supercontinents of the future can provide us some guidance on how surface temperatures will increase or decrease depending on how the continents are distributed, with implications for exoplanet climate and habitability. But there are other factors to consider related to weathering rates and volcanic outgassing (e.g., Jellinek et al., 2019), not to mention the related role of atmospheric pressure (Gaillard & Scaillet, 2014). We have also used a fixed atmospheric concentration in this paper to avoid introducing a further parameter that can add climate variability and, interesting as it would be, exploring the climate with a dynamic carbon cycle is left for future work.

The 30 min increase in the length of day between simulations 07 and 08 appears to play little to no role in the climate dynamics as there is little discernible difference in the strength or distribution of the Hadley or eddy transport diagnostics. This implies the same for simulations 01–06 with their 30 min longer day lengths than present day Earth.

While we discuss the future climate of Earth we do not touch on the future of life. There are many uncertainties, but recent work provides some guidelines (Mello & Friaça, 2019). The reduced tides during the supercontinent stage (Davies et al., 2020) will lead to reduced vertical mixing rates, that is, a reduced vertical diffusivity in the abyssal ocean (Munk, 1966; Wunsch & Ferrari, 2004). This may have implications for ocean ecosystems, and biodiversity. At the same time, it appears that the formation of Pangea had little effect on the global biodiversity of marine animals (Zaffos et al., 2017) and Pangea was in a very weak tidal state (Green et al., 2017).

It would be interesting to compare the GCM derived climates for the supercontinent at low latitude in the Aurica runs with previous work on Pangea (e.g., Chandler et al., 1992; Chandler, 1994; Fluteau et al., 2001; Gibbs et al., 2002; Roscher et al., 2011). Unfortunately, it is difficult to make a proper comparison for a number of reasons. First, all of these previous works use either atmosphere only GCMs (i.e., no ocean) or shallow mixed layer oceans with either prescribed horizontal heat transport or none at all. Second, unlike Aurica, Pangea spanned not only lower latitudes (like Aurica), but also high southern latitudes where ice/snow forms easily (e.g., Chandler et al., 1992, Figure 5). Finally, there are different reconstructions for different time periods and not all are directly comparable to those we simulate herein. This makes a direct comparison with Pangea complicated and we leave such an analysis for the future.

These new reconstructions may prove useful for exoplanetary researchers who will have a larger library of topographies and land/sea masks to choose from when estimating the probability of surface habitability on neighboring worlds.