2.1 Model Description

All simulations were performed with the CLIMBER-3α Earth System Model of Intermediate Complexity (EMIC) (Montoya et al., 2005), consisting of a modified version of the ocean general circulation model MOM3 (Hofmann & Morales Maqueda, 2006; Pacanowski & Griffies, 2000), run at a horizontal resolution of 3.75° × 3.75° with 24 vertical levels, a dynamic/thermodynamic sea-ice model (Fichefet & Maqueda, 1997), and a fast statistical-dynamical atmosphere model (Petoukhov et al., 2000) with a coarse spatial resolution of 22.5° in longitude and 7.5° in latitude. The main difference between CLIMBER-3α and a coupled atmosphere-ocean general circulation model (AOGCM) is the simplified atmosphere model which describes the dynamics of large-scale circulation patterns, but statistically parametrizes the effects of synoptic-scale systems rather than solving the fundamental equations. This makes the model computationally efficient and allows running large ensembles of simulations. Although CLIMBER-3α has evolved over the years, differences between simulated and observed climatic variables under present-day boundary conditions are still broadly similar to the original version described in Montoya et al. (2005). Specifically, surface air temperatures are overestimated in mountain regions like the Andes or the Tibetan plateau. Furthermore, simulated temperatures are generally too high over present-day Antarctica, in particular during Southern winter, and too low over the Northern Atlantic during Northern winter. The model has been successfully used in a number of paleoclimate studies (Brugger et al., 2017, 2019; Feulner, 2017; Landwehrs et al., 2020), and compares well with other models in model intercomparison projects (Eby et al., 2013; Zickfeld et al., 2013).

2.2 Boundary Conditions

Climate simulations were conducted for a set of 40 timeslices with a 5 Myr spacing from 255 to 60 Ma (denoted as T_5Myr), using paleotopographies/-bathymetries based on the reconstruction of Scotese and Wright (2018) (also see Figure S1). These were adjusted regarding model requirements, for example, by avoiding narrow ocean channels and isolated ocean cells on the model grid resolution (see supporting information section 1). Vegetation cover was prescribed for each timeslice based on the proxy-based estimates of climatic zonation by Boucot et al. (2013) and Scotese (2016) (also see Figure S1). Their five climate types were represented by certain fractions of tree, grass/shrub, and bare soil cover for tropical/boreotropical (75%, 20%, and 5%), arid (15%, 35%, and 60%), warm temperate (70%, 20%, and 10%), cold temperate (65%, 15%, and 20%), and polar (10%, 25%, and 65%) climates, using Bonan et al. (2002), Poulter et al. (2015), and Pfadenhauer and Klötzli (2015) for orientation. From these spatially heterogeneous vegetation patterns (VegHet), homogeneous vegetation patterns with average compositions over the entire land area were defined for sensitivity experiments: VegHom reflects the average composition for each respective timeslice, while VegFix has a constant average composition (51% trees, 25% grass/shrub, and 24% bare soil) over all timeslices (see Table  S1 and Figure S3). Runoff is routed along topographic gradients, with necessary modifications, for example, where internal basins occur. The solar constant (S₀) increases approximately linearly from 1,332.7 to 1,354.7 W/m² (∼97.9% to 99.5%; from 255 to 60 Ma) according to standard solar evolution (Bahcall et al., 2001) relative to its present-day value of 1,361 W/m² (Kopp & Lean, 2011).

Different pCO₂ pathways are tested here to account for the large uncertainties associated with existing pCO₂ reconstructions for the Mesozoic: The P_{pCO2_proxy} pathway corresponds to the fit to compiled proxy data provided by Foster et al. (2017) (see Figure 1b). This is contrasted with P_{pCO2_COPSE}, which reflects the evolution of pCO₂ obtained with the COPSE carbon cycle model (Mills et al., 2019). This model estimate was included because it correlates reasonably with proxy-based global temperature trends (Mills et al., 2019) and helps to illustrate the spread of different proxy- and model-based pCO₂ reconstruction approaches. For the purpose of this study, the pCO₂ values of both reconstruction pathways were rounded to the nearest hundred at each timeslice. A fixed value of pCO₂ = 1,000 ppm (P_{pCO2_1000ppm}) appears a reasonable median for most of the Mesozoic, except the Middle and Late Jurassic and the latest Cretaceous where both considered reconstructions agree on lower values (see Figure 1b). Additional simulations were performed for a subset T_25Myr (250, 225, 200, 175, 150, 125, 100, 75, and 65 Ma) of the timeslices, to further test the sensitivity to different boundary conditions: P_{pCO2_const} for pCO₂ fixed at (250, 500, 1,500, 2,000, 4,000) ppm, P_S0ini for a fixed S₀ of 1,332.7 W/m² as well as P_VegHom and P_VegFix for homogeneous vegetation patterns (see Table 1).

To test the systematic effect of changing sea levels on global climates, further experiments were performed for raised or lowered sea levels and subsequently modified land–sea distributions for the 200 Ma (sea-level offset ΔSL = [+40, +200] m), 150 Ma ([−40, +40] m) and 100 Ma ([−200, −40, +40] m) timeslices. The ΔSL = ±200 m offsets were chosen to represent potential effects of the long-term sea-level rise from the Late Triassic to the mid-Cretaceous that has been reconstructed by various methods (e.g., Wright et al., 2020, Figure 16). The ΔSL = ±40 m offsets were tested to represent the variability of sea levels on the million-year scale captured by sedimentary regression-transgression sequences throughout the Mesozoic (Haq et al., 1988; Miller et al., 2005; Ray et al., 2019; Sames et al., 2020) and because 40 m is the vertical resolution of the elevation models from Scotese and Wright (2018). In the latter, coastlines were inferred from plate-tectonic modeling and lithofacies and should represent the average best-guess paleogeography for the covered time interval (Scotese & Wright, 2018). Independently from the determination of the coastlines, the authors assume land-ice where glacial deposits occur (Boucot et al., 2013), but for the Mesozoic only minor ice volumes at 140–120 Ma have been inferred by Scotese (2018). For the purpose of this study, continental ice-sheets were thus assumed to be absent in all simulations.

In most of the simulations, Earth's orbital parameters were fixed at a configuration representing a circular orbit (eccentricity e = 0, precession angle ω = 0°) with an intermediate obliquity (ϵ = 23.5°). To represent effects of lower-/upper-end obliquities and of a strongly eccentric orbit, four additional configurations with ϵ = 22.0° and ϵ = 24.5°, e = 0.06, ω = [90, 180, 270] ° were tested for the T_25Myr subset of timeslices. All 221 simulation runs were initialized with the same idealized modern ocean temperature and salinity profile without sea ice and were integrated for 5,000 model years (see also supporting information section 2). All analyses are based on averages over the last 500 years. For the plots and analyses in Sections 6-9, atmospheric model data were interpolated on the ocean model grid to enable investigating the influence of changing land and sea areas. Table 1 provides an overview over the described ensemble of boundary conditions and an additional table with information on all individual simulation runs is included in Table S1.

3.1 Global Mean Temperature Evolution

The simulated global mean surface air temperatures (GMST) for all tested pCO₂ pathways are summarized in Figure 1a, and for all other sensitivity experiments in Figure 2a. Overall, the GMSTs at pCO₂ = 1,000 ppm (P_{pCO2_1000ppm}, black dots and lines in Figures 1a and 2a) exhibit a slight cooling through the Triassic to a minimum of 17.7°C in the Late Triassic and then a +3.5°C warming to up to 21.2°C at 80 Ma in the Late Cretaceous. The two considered pCO₂ reconstructions (P_{pCO2_proxy} and P_{pCO2_COPSE}, Figure 1b) do, however, suggest deviations of up to ∼5°C from this median pathway (Figure 1a). To decompose these temperature trends, we further examine the P_{pCO2_1000ppm} pathway, which reflects the effects of changing paleogeography, vegetation patterns, and the increasing solar constant S₀. The latter elevates the GMST at 65 Ma by +1.9°C (difference between P_{pCO2_1000ppm} and P_S0ini in Figure 2a, also see Table S1), which indicates the order of warming through the Mesozoic caused by the increasing solar luminosity if everything else would have remained the same.

Despite fixed pCO₂ and S₀, the P_S0ini temperatures exhibit variations with time (see Figure 2a) that reflect the changing distribution and extent of land masses and vegetation patterns: The GMST decreases from 18.5°C at 250 Ma to 17.2°C at 200 Ma and subsequently rises to 19.5°C at 75 Ma. These trends correlate with changes in the fraction of Earth's surface occupied by land (Figure 2b, black line), which has a higher average albedo than ocean areas and thus influences Earth's energy balance. The global land fraction in the implemented paleogeographies is highest during the Late Triassic-Early Jurassic (up to 38%) and lowest in the Late Cretaceous (down to 26%) which is related to the change of relative sea levels (Figure 1c, blue line) that are implicitly included in the employed paleogeographic reconstructions (Scotese & Wright, 2018). The effect of raised or lowered sea levels and the resulting flooding/exposure of continental area is also demonstrated for individual timeslices (Figure 2, blue circles). For example, the global land fraction increases by +8% for a 200 m sea-level fall at 100 Ma and the GMST is reduced by −1.2°C. In contrast, a +1.2°C global warming is obtained for a +200 m sea-level rise at 200 Ma. These findings suggest that the lower sea level and larger continental area in the Triassic-Jurassic would have contributed to relatively low temperatures, and that the subsequent sea-level rise provided a major part of the ≈+2°C global warming toward the Late Cretaceous, exhibited by the P_S0ini pathway. To test to which extent these trends are dependent on the uncertain reconstruction of paleotopography, lapse-rate corrected global mean temperatures at sea level were calculated (see Figure S12). On the global scale, they follow the same trend although differences can be noted, especially around 100 Ma where the global mean topography height reaches a maximum. On a regional scale, differences resulting from different paleogeographic reconstructions can potentially be much more pronounced.

The described temperature rise between 200 and 75 Ma is ∼0.5°C larger than for a fixed homogeneous vegetation composition (P_{VegFix_S0ini}, Figure 2a, green line), which yields slightly elevated temperatures in the Triassic–Jurassic. This is related to the large portion of bare land area in the Pangaean low-to mid-latitudes (for VegHet), which contributes to a higher surface albedo in these regions with high insolation and thus slightly lower GMSTs in the Triassic and Jurassic than with a homogeneous vegetation (see difference between P_VegHom and P_{pCO2_1000ppm}, Figure 2a). In our modeling framework, the trend toward denser vegetation, especially in low-latitude land areas, thus contributes to a Late Triassic-Late Cretaceous warming, although on a slightly smaller scale than changing sea levels and solar luminosity. On the P_{VegFix_S0ini} pathway, the distribution of land area is the only varying boundary condition if lapse-rate corrected air temperatures at sea level are considered. We find that global mean sea-level temperatures are indeed strongly negatively correlated with the global land fraction, but also correlate negatively with the average absolute latitude of land area which peaks around 90 Ma (see Figure S13). This indicates that paleogeographic changes contribute to the warming toward the mid-Cretacous both by the decreasing land area and by its distribution to higher latitudes. The simulated GMSTs also vary on the order of ∼0.4°C–0.6°C with the orbital configuration, with low-obliquity configurations (ϵ = 22.0°, e = 0) being generally cooler than those with higher obliquity and eccentricity (ϵ = 24.5°, e = 0.06) (see Figure 2a).

As already noted, the amplitudes of GMST variations resulting from the two considered pCO₂ reconstruction pathways exceed those of the other discussed boundary conditions. The pCO₂ proxy data and the provided fit (Foster et al., 2017) exhibit a large spread in the Triassic, yielding relatively low temperatures in the Early–Middle Triassic but very high temperatures in the Late Triassic to earliest Jurassic (up to 21.1°C at 225 Ma with 2,000 ppm). The COPSE model-derived pCO₂ falls to ∼500 ppm in the Early Jurassic, yielding relatively low temperatures from the latest Triassic to the Middle Jurassic (down to 13.9°C at 200 Ma with 500 ppm). Although discrepancies between both pathways are evident, common relative trends can be identified: Elevated temperatures in the Triassic which decrease into the Jurassic, but subsequently rise again toward the mid-Cretaceous. For both reconstructions, mid-Cretaceous pCO₂ levels do not exceed the relatively high values around the Late Triassic, but the additional warming from the previously discussed drivers, including rising sea levels and solar luminosity, contributes to similar or higher temperatures. The P_{pCO2_proxy} pathway does not capture the prominent Cenomanian–Turonian peak (∼100–90 Ma) suggested by δ¹⁸O proxy data (Friedrich et al., 2012), whereas P_{pCO2_COPSE} yields sustained high temperatures above 20°C from 110 to 75 Ma and generally yields long-term trends similar to the temperature proxies (see Figure 1c), as already noted by Mills et al. (2019). Overall, the simulated GMSTs for both reconstructions are elevated above pre-industrial values of ∼14°C for most of the Mesozoic, and for the mid-Cretaceous (around 110 Ma) both agree on temperatures >20°C. For the latest Cretaceous, they consistently yield a pronounced cooling.

3.2 Seasonal Temperature Contrasts

It has been previously suggested, that the climate on the Pangaean supercontinent was characterized by large seasonal contrasts and intense monsoons (e.g., Kutzbach & Gallimore, 1989; Parrish, 1993), but that these conditions decayed with continental breakup through the Mesozoic (e.g., Parrish, 1993). Here, we assess the systematic impact of paleogeography and other boundary conditions on the global-scale amplitude of seasonal temperature contrasts. For this, maximum surface air temperature (SAT) differences between the DJF, MAM, JJA, and SON seasonal averages were calculated. To illustrate the evolution of regional patterns through the Mesozoic, maps of the resulting temperature seasonality are shown for one Late Triassic (225 Ma) and one Late Cretaceous (75 Ma) timeslice of the P_{pCO2_1000ppm} pathway in Figures 3a and 3b. It can be seen that the vast mid- and high-latitude continents of Pangaea experience strong seasonal temperature variations in the early Mesozoic (Figure 3a), while these contrasts are reduced in the fragmented continental configuration of the Late Cretaceous (Figure 3b).

The evolution of the zonal mean continental SAT seasonality obtained from all 40 timeslices of the P_{pCO2_1000ppm} pathway is summarized in Figure 3c (contours). The gray shading indicates the changing latitudinal distribution of land area in the paleogeographies implemented for these simulations. This illustrates the correlation of the amplitude of seasonal temperature changes with the extent of land area, especially in the mid to high latitudes. These amplitudes are highest during the Early Triassic–Early Jurassic on the Southern Hemisphere and during the Late Triassic–Early Jurassic on the Northern Hemisphere and reach up to 40°C. Toward the Late Cretaceous, the mid to high latitude seasonal temperature contrasts are considerably reduced on both hemispheres and only reach ∼20°C–30°C. These trends are also apparent in the tropical latitudes but are strongly dampened and remain below ∼10°C.

On a global scale, the average continental SAT seasonality (Figure 3d) peaks at ∼22°C in the Late Triassic–Early Jurassic and falls down to ∼15°C in the Late Cretaceous. As pCO₂, S₀ and vegetation patterns do not significantly affect these trends, the latter are mainly driven by paleogeographic changes and closely mirror the change of the global land area fraction (see Figure 2b). The major contribution of sea-level rise to reduced seasonality in the Late Cretaceous relative to the Late Triassic-Early Jurassic is demonstrated by the experiments with elevated or lowered sea level at 200 and 100 Ma: with a sea level elevated by 200 m, at 200 Ma, the global SAT seasonality is reduced to levels comparable to those at 100 Ma, and vice versa. Although these results suggest that the long-term seasonality trends were indeed determined by sea-level changes and continental fragmentation, we also find that different orbital configurations modify SAT seasonality to a similar degree (Figure 3d). The global SAT seasonality generally varies by about 4°C–5°C between orbital configurations which yield a lower seasonality (lower obliquity ϵ = 22.0°) and higher seasonality (higher obliquity ϵ = 24.5°). This suggests that orbital cycles would have considerably overprinted Mesozoic long-term seasonality trends.

3.3 Meridional Temperature Contrasts

Simulated annual mean SAT for four selected timeslices at 1,000 ppm (P_{pCO2_1000ppm} pathway) are shown in Figure 4. It can be seen that for these simulations equatorial SATs mostly exceed 25°C, while temperatures below 0°C occur especially where extensive land masses exist in polar latitudes. To systematically assess the evolution of meridional temperature contrasts and their dependence on paleogeographic changes and other boundary conditions, the equator-to-pole thermal gradients where calculated for all conducted simulations. Here, these were defined as the difference between low-latitude (<30° lat.) and high latitude (>60° lat.) annual and zonal mean SATs.

Plots of zonal mean SAT and SST profiles for all timeslices are included in the accompanying data repository (also see Figure S16). The evolution of zonal mean SATs for the P_{pCO2_1000ppm} pathway is summarized in Figure 5a (contours). On the one hand, the contours reflect the previously discussed general warming trend from the Late Triassic-Early Jurassic to the Late Cretaceous. On the other hand, it can be inferred that low polar temperatures and high meridional temperature gradients (Figures 5b and 5c) correlate to first order with the amount of land in polar latitudes. From the end-Permian to the Late Triassic, Pangaea shifted northwards, which leads to a high-latitude warming and a ∼5°C lower thermal gradient (Figure 5c) in the southern hemisphere, whereas the opposite is observed in the northern hemisphere (Figure 5b). The northern hemisphere thermal gradient shows a decrease from the Late Jurassic into the Late Cretaceous in line with the general decrease of land area related to sea-level rise. In contrast, the southern hemisphere meridional temperature gradient increases during this time due to the migration of Antarctica toward the South Pole. Despite the primary control of paleogeography on these trends, our experiments suggest considerable effects of other boundary conditions. For example, the northern hemisphere thermal gradients generally decrease (on the order of ∼5°C) in response to an increase of pCO₂ from 500 to 2,000 ppm and the associated warming. This effect can be related to the snow and sea-ice albedo feedback and is even more pronounced in the southern hemisphere especially before the Middle Jurassic. Relatively low pCO₂ levels suggested by the tested reconstructions for the Jurassic and latest Cretaceous could thus have contributed to higher equator-to-pole thermal gradients relative to the mid-Cretaceous. Through the ice albedo feedback, the warming associated with the increased solar luminosity also contributes to lower thermal gradients in the Late Cretaceous on a smaller, but notable scale.

Periods with lower northern and southern temperature gradients are also found to be times with enhanced meridional heat transport by the ocean in the respective hemisphere (see supporting information section 8 and Figure S18). This also appears to be related to enhanced meridional overturning (see also Figure S19) and ocean mixing (Figures  S20 and S21) and deep water formation. These are enhanced at high northern latitudes in the early Triassic as well as again in the mid- and late Cretaceous, and at high southern latitudes from the Early Triassic to the mid-Cretaceous. However, it will require further work to disentangle these phenomena and evaluate their significance. Ocean model output data are included in the accompanying data repository.

3.4 Zonal Temperature Contrasts

Our simulations suggest the existence of considerable zonal climatic contrasts prior to the breakup of Pangaea, as can be inferred from the simulated SAT patterns and the distribution of arid regions (pink hatches, see Section 9 for explanation) in Figures 4a and 4b: High temperatures and humid conditions are inferred for the tropical latitudes of eastern Pangaea, in contrast to its western low- to mid-latitudes, where temperatures are generally lower and arid regions are more extended. These differences are related to circulation patterns in the Panthalassa ocean where an anticyclonic subtropical gyre in each hemisphere transports warm waters along the equator into the Tethys, but colder water from higher latitudes to the western Pangaean tropics and subtropics. These provide less moisture and precipitation and lower temperatures in contrast to the very warm Tethys Sea. However, the opposite effect occurs in the mid- to high latitudes where cyclonic subpolar gyres warm the western coasts compared to the eastern ones.

To systematically characterize the evolution of zonal temperature contrasts on the continents through the Mesozoic, the deviation of local annual mean SATs from their zonal mean values is calculated. For two selected timeslices of the P_{pCO2_1000ppm} pathway, the resulting patterns are shown in Figures 6a and 6b. As described above, the eastern tropics and the western mid- to high-latitudes of Pangaea are warmer by up to ∼8°C compared to the respective western or eastern regions (Figure 6a). For the fragmented late Cretaceous continental configuration, these east–west contrasts are reduced (Figure 6b). Figure 6c summarizes the zonal mean of the absolute value of these deviations for all considered geologic timeslices of the P_{pCO2_1000ppm} pathway. This indicates that zonal temperature contrasts are amplified in the mid–high latitudes as well as the tropical-subtropical latitudes throughout the Triassic-Jurassic but are clearly reduced after the Middle Jurassic. More uniform zonal temperature patterns were thus established in the mid-Mesozoic while they were considerably perturbed by the supercontinent configuration before. This trend is also indicated by the global average of the absolute SAT deviations from their zonal mean (Figure 6d), which is interpreted as an indicator of global temperature “azonality”. It peaks at ∼2.2°C in the Late Triassic and Early Jurassic and is halved in the Late Cretaceous. Changes in boundary conditions other than paleogeography only lead to relatively small deviations from this trend. Our simulations thus illustrate a general trend toward decreased east–west temperature contrasts through the Mesozoic, mainly driven by continental fragmentation.

3.5 Aridity

Previous modeling work (Chaboureau et al., 2014; Donnadieu et al., 2009) suggested a reduction of the extent of arid climate conditions from the Late Triassic to the Late Cretaceous. To test this, we determine the simulated extent of relatively dry areas using the “radiative index of dryness” of Budyko (1974) which expresses the ratio of potential evaporation (RB/L) to precipitation (PRC): AI = RB/(L ⋅PRC) (Stadler, 2005). Here, RB is the annual mean net surface radiation balance and L is the latent heat of water evaporation and higher values of AI indicate drier conditions. The transition from humid forest climates to dry steppe climates occurs at values around 1.5 today (Budyko, 1974; Mabbutt, 1978). The changing extent of dry land area is indicated in Figure 4 (pink hatches) for four timeslices and in Figure 7a for all timeslices of the P_{pCO2_1000ppm} pathway. It can be seen in Figure 7a that the dry regions have their largest extent in the subtropical to mid-latitudes during the Late Triassic-Early Jurassic. In the low latitudes, dry regions are extended during the Triassic–Middle Jurassic, especially in western Pangaea (see Figure 4), but successively disappear in the further course, which indicates the establishment of more continuous tropical humid conditions. The relative portion of global arid land area on the P_{pCO2_1000ppm} pathway (Figure 7b) peaks in the latest Triassic and continuously decreases into the Late Cretaceous. In agreement with Chaboureau et al. (2014), the fraction of arid land is found to decrease with warming from elevated pCO₂. The relatively low pCO₂ suggested by the COPSE model reconstruction for the latest Triassic to Middle Jurassic would thus have significantly contributed to increased aridity during this time which would have amplified the humidification trend toward the mid-Cretaceous. Also the continuous warming from increasing solar luminosity contributes to more humid conditions through the Mesozoic (Figure 7b, difference between P_{pCO2_1000ppm} and P_S0ini). The cooling in the latest Cretaceous resulting from decreasing pCO₂ in both considered reconstructions is found to increase aridity and to reverse the Early Jurassic to mid-Cretaceous trend toward increasingly humid climates.

Coal and evaporite occurrences are frequently used as paleoclimatic indicators of humid or arid conditions (Markwick, 2007; Ziegler et al., 2003), and a certain agreement can be observed between their distributions and the humid and dry zones simulated here for our baseline pathway P_{pCO2_1000ppm} (Figure 4; proxy data from Boucot et al. (2013) and Cao et al. (2018), locations reconstructed with pyGPlates, see supporting information section 1, for further explanation). For example, many of the coal deposits occur in humid mid- to high-latitude or tropical regions, while evaporites and calcretes are frequent in the relatively dry regions of the western Pangaean continents. This simple scheme can, however, not completely explain these spatial distributions (e.g., evaporites in tropical humid regions). This can be due to a variety of factors, including the range of possible local depositional environments and their variation on shorter (e.g., orbital) timescales. Considering these difficulties and the limitations of the employed climate model, a comprehensive assessment that takes into account the various boundary conditions is beyond the scope of this study. It can still be noted, that Cao et al. (2018) found a statistical Mesozoic–Cenozoic trend from an unimodal to a bimodal zonal distribution in the evaporite occurrence data of Boucot et al. (2013). Our simulations suggest that this could express the establishment of a more consistent tropical humid belt and two arid belts during the fragmentation of Pangaea (see Figure 4d), which could have displaced evaporite formation from equatorial regions especially in the eastern Panthalassa margin. The simulated overall trend toward more humid and less arid continental climates through the Mesozoic also conforms qualitatively with opposite trends of increasing coal versus decreasing evaporite deposition (Figure 7c).