Elevation-dependent temperature response in early Eocene using paleoclimate model experiment

To study the elevation-dependent temperature changes during the Eocene hothouse, we investigated CESM1.2 model simulations which is available at the deep-time model intercomparison project framework (DeepMIP, described in [33]). All participating modeling group experiments were forced by proxy-estimated CO₂ levels. The model simulation was performed with CO₂ concentrations from 280 ppm (pre-industrial levels) to 2520 ppm (9 × CO₂). As our ultimate goal was to investigate elevation-dependency on the temperature at various atmospheric CO₂ concentrations, we used CESM1.2 simulations for the following reasons. First, its sensitivity experiments operate with atmospheric CO₂ concentrations ranging from 1 to 9 times the pre-industrial value, under paleogeography with vegetation type conditions, and river routing [42] and preindustrial levels of other gases and orbital forcing. Second, this model shows an overall agreement with proxy records and captured realistic greenhouse warming features during the early Eocene [40, 43]. Even though proxy-estimated data during Eocene are limited, the model simulated zonal mean temperature over land illustrates consistency among corresponding proxy estimates (supplementary figure 1). CESM1.2 configuration used in this study consists of the Community Atmosphere Model 5.3, the Community Land Model 4.0, the Parallel Ocean Program 2, the River Transport Model, the Los Alamos sea ice model 4 and a component coupler [41].

All CESM Eocene simulations were conducted on the horizontal resolution of 1.9° × 2.5° over the atmosphere and land with 30 hybrid sigma-pressure model levels in the atmosphere. The ocean and sea ice at a nominal 1° grid with 60 vertical model levels in the ocean. All CESM model simulations have been integrated for 2000 model years, except for 1 × CO₂, which ran for 2600 model years. To make the model suitable for a paleoclimate simulation under a high CO₂ scenario, the model radiation code has been slightly upgraded to represent missing diffusivity angle specifications [40, 44].

In addition to the different CO₂ experiments (1 × CO₂, 3 × CO₂, 6 × CO₂, and 9 × CO₂), a pre-industrial simulation based on the CMIP5 standard was considered for comparison. The pre-industrial simulation 1 × CO₂ attempted to represent the preindustrial climatic conditions and 4 × CO₂ were conducted to represent the future conditions. However, DeepMIP simulations starting at CO₂ concentration equal to the pre-industrial conditions assumed a planet without a permanent ice sheet or sea ice cover at the poles [45]. Here, we have considered the climatological mean state from the equilibrated last 100 model years of each simulation.

2.2.1. Land amplification

The land amplification factor was calculated [46] using the ratio of land warming per unit of zonally averaged warming of the adjacent oceans (Globally, longitudes from 0° to 360°),

$$\begin{align}A = \frac{{{{\Delta }}{T_{\text{L}}}}}{{{{\Delta }}{T_{\text{o}}}}}\,\,\,\,\,\end{align} \tag{ 1 }$$

where ΔT_L is the change in local air temperature over masked land (latitude from −90° to 90°), and ΔT_o is the zonally averaged change in air temperature over the ocean.

2.2.2. Polar amplification

The polar amplification factor is expressed [47] as ratio of the change in temperature over land to the corresponding change in global mean temperature,

$$\begin{align}A = \frac{{{{\Delta }}{T_{\text{L}}}}}{{{{\Delta }}{T_{{\text{gm}}}}}}\,\,\,\,\,\,\end{align} \tag{ 2 }$$

where ΔT_gm is the change in the global mean temperature.

The model elevation topography was available for each grid resolution used in the model simulation. High mountains were assumed to have an orographic elevation of more than 1 km. Therefore, we screened out orographic regions that exceeded 500 m or more in height (figure 1). For the present analysis, we selected the five major high mountain ranges during the Eocene era, as shown in figure 1(b). The regions were the North American mountain range, the Asian mountain range, the Indian Island mountain, the South African mountain, and the Antarctic mountain range. To derive the elevation dependence of warming, we created vertical profiles using equally sized altitude bins with 500 m intervals. To make a bin, we detected the model elevation at each elevation grid point and clustered them starting from 500 m to the mountain summit and taken average at that mid-class interval (supplementary figure 2). Depending on the size of the respective region and the estimated elevation in the paleogeography, the size and number of bins may vary. The mean temperatures of the grid points in each elevation bin were calculated and changes from the sensitivity experiments relative to the preindustrial simulation are represented by lines in supplementary figure S2. The elevation-dependency of the change in the variables was computed by taking the slope between the variable and elevation over each individual mountain range. Least-squares regression was utilized for calculating the slope.

Zoom In Zoom Out Reset image size

The surface energy budget analysis used to investigate the causes of elevation-dependent temperature sensitivity can be written in the following form:

$$\begin{align}Q = {S_ \downarrow } + {S_ \uparrow } + {F_ \downarrow } + \,{F_ \uparrow } + \left( {S + L} \right)\end{align} \tag{ 3 }$$

where Q denotes land surface heat storage and ground heat fluxes, ${S_ \downarrow }$ is the downwelling shortwave radiation, ${F_ \downarrow }$ is the downwelling longwave radiation, ${F_ \uparrow }$ is the longwave radiation from the surface, and (S + L) is the turbulent sensible and latent heat fluxes, respectively. Here all fluxes were described as positive (negative) when upward (downward). The equation (3) can be rewritten as:

$$\begin{align}Q = \left( {1 - \alpha } \right){S_ \downarrow } + {F_ \downarrow } - \varepsilon \sigma {T_{\text{s}}}^4 + \left( {S + L} \right)\end{align} \tag{ 4 }$$

where $\alpha$ is the surface albedo considered as a ratio of upwelling to downwelling shortwave radiation for clear-sky conditions. ${F_ \uparrow }$ can be expressed as σT_s ⁴ as a function of the surface radiative/skin temperature T_s, the approximated surface emissivity ( =1), and the Stefan–Boltzmann constant (σ).

Considering changes in the surface energy budget with respect to changing CO₂ concentrations, following the previous methodology [48, 49],

$$\begin{align}4\varepsilon \sigma {T_{\text{s}}}^3\Delta {T_{\text{s}}} = \Delta [\left( {1 - \alpha } \right){S_ \downarrow }] + \Delta {F_ \downarrow } - \Delta Q- \Delta \left( {S + L} \right).\end{align} \tag{ 5 }$$

In equation (5), $\Delta$ denotes the change in the high CO₂ scenario (3×, 6×, and 9×) with respect to 1 × CO₂.

Lu and Cai [48] defined the cloud radiative effect (CRE) as the difference between the all-sky and clear-sky radiations at the surface,

$$\begin{align}\Delta {\text{CRE}} = \left( {1 - \bar \alpha } \right)\Delta {S_{ \downarrow {} {\text{clr}}}} + \Delta {F_{ \downarrow {} {\text{clr}}}}\end{align} \tag{ 6 }$$

where overbar denotes 1 × CO₂ mean state. $\Delta {S_{{\text{dn,clr}}}}$ and $\Delta {F_{{\text{dn,clr}}}}$ calculated by taking difference between all-sky condition and clear-sky condition radiative fluxes.

Using equation (6), we can rewrite equation (5) and divide it by $4\varepsilon \sigma {T_{\text{s}}}^3\,$ will give us the surface skin temperature changes into partial temperature contributions, resulting in the linearly decomposed form as follows:

$$\begin{align}\Delta {T_{\text{s}}} = \frac{ - \left( {\Delta \alpha } \right)(\overline {{S_ \downarrow }} + \Delta {S_ \downarrow }) + \Delta {\text{CRE}} + (1 - \bar \alpha )\Delta {S_{ \downarrow {} {\text{clr}}}} + \Delta {F_{ \downarrow {} {\text{clr}}}} - \Delta Q - {\mkern 1mu} \Delta \left( {S + L} \right) }{{4\varepsilon \sigma {T_{\text{s}}}^3}}\end{align} \tag{ 7 }$$

Each term on the right-hand side of equation (7) contains the partial temperature change due to the corresponding CO₂ scenario, which are surface albedo feedback (SAF), the CRE, changes in shortwave radiation for clear-sky not related to SAF (SWCS), longwave radiation for clear-sky (LWCS), land heat storage (HSL), and surface turbulent fluxes (FLUX).

We performed the above calculation at each model grid. Furthermore, we illustrated it in the form of elevation dependency as our study mainly emphasizes elevation-dependent temperature response.

A significant difference in the land-sea distribution was evident in the paleogeography [42, 50] compared with a modern topographic map (figure 1). According to the plate tectonics theory [51], continent orientation and ocean bathymetry differ. In the Eocene, mountains were less elevated, and their orographic features were quite different from those in the Anthropocene. Large mountain chains were distributed globally, from the tropics to the poles, and were more concentrated over the mid-latitudes. Figure 1 shows the global mountain ranges in Eocene and pre-industrial simulation. For instance, the highest mountains were located in North America (Rocky Mountains) and the Antarctic during the Eocene, with vegetated land surfaces devoid of a permanent ice sheet [33, 42]. As the Indian subcontinent was not yet connected to Asia, the Tibetan Plateau and Himalayan orogenesis had not occurred. Instead, mountains were distributed on the eastern side of Asia. South Africa had a distinct orography.

The specialty of paleogeography is not just limited to topographic features but also land–sea contracts. Land–sea contrast is a fundamental aspect of the climate system that triggered regional impacts [52, 53]. Figures 2(a) and (b) shows overall land, and polar amplification, respectively, in response to 9 × CO₂, the scenario with the highest CO₂ concentration in the sensitivity experiment (and more information can be found in supplementary table 1). This indicates that there was more prominent warming over land than compared to the ocean, implying that moistening rate with warming over the ocean was greater than that over land (figure 2(a)). The atmosphere gets moister with the ocean warming as the amplification factor is positive, which is nothing but land warming per unit zonal mean warming of the surrounding oceans. The polar amplification response to 9 × CO₂ was also evident with land amplification, with more coverage toward the poles (figure 2(b)). The zonal mean land amplification spread can be seen in the sensitivity experiments. However, the high-elevation warming response (figure 2(e)) was comparatively lower in the context of land amplification, except at the South Pole. This also authenticates that the mean global mountain warming rate is the same in all Paleocene CO₂ responses and nearly consistent with future CO₂ responses, revealing the importance of elevation topography (supplementary table 1). Therefore, we need to carefully analyze the climate response of mountains on a regional scale. Our primary objective is to explore the elevation-dependent temperature response in the early Eocene and assess its sensitivity to CO₂.

Zoom In Zoom Out Reset image size

The elevation-dependent temperature change was evident in both hemispheres (figure 3). However, its sensitivity differed across regions and, with CO₂ concentrations. Here, the slope estimation method indicated the elevation dependency as a linear function, which may be different from that in the real world. The EDW was explicitly amplified at 3 × CO₂ over Asia, Africa, and the Antarctic. The EDW exists over Asia, Africa, and the Antarctic under 3 × CO₂, furthermore heightened under 9 × CO₂ response. In addition, the chaotic shift from EDW to EDC was observed in the Indian plate (figure 3(c)). In contrast, North America shows both flavors of elevation-dependency with less sensitivity under CO₂ concentrations; it shows EDC in 3 × CO₂, but an EDW signal emerged with increased CO₂. Such mixed behavior of elevation-dependent temperature changes among the mountains underscores the complexities across the globe. These changes range from −0.5 °C to 4 °C km⁻¹, slightly negative to strongly positive across the six mountain regions. Interestingly, this temperature gradient seemed more tilted towards the Southern Hemisphere's response to Antarctic amplification. These elevation-dependent changes in surface temperature showed a nonlinear response with more sensitivity to higher CO₂ concentrations.

Zoom In Zoom Out Reset image size

Figure 4 depicts elevation-dependent temperature changes using surface air temperature (T) and skin temperature (T_s). Both temperature changes have notable seasonal variations with an almost identical pattern in the elevation-dependent temperature response, pointing to a close connection between them. The monthly variation increased with higher CO₂ concentrations over each mountain range, and seasonally the reversed response can be noticed at some places. For instance, North America shows EDC during winter, with small values replicated with the CO₂ scenario. However, EDW in summer in higher CO₂ scenarios counteracts the EDC resulting in an altogether EDW in North America (figure 4(a)). The mountain ranges in Asia and Africa experienced EDW in all seasons with different intensities and sensitivities to increasing CO₂ levels, showing very high values and sensitivity in summer, with opposite trends in winter (figures 4(b) and (d)). Similarly, Antarctica showed the most substantial EDW seasonality. In this region, the EDW was the strongest having its highest peaks in the austral summer (December, January, and February), with especially high sensitivity in January. In the austral winter, values are near zero or showed slight EDC (figure 4(e)). In contrast, the Indian plate experienced a decrease in EDW with increased CO₂ concentration (3 × CO₂ and 6 × CO₂) and EDC (9 × CO₂). These effects were reversed in the scenario with the highest CO₂ concentration and weak seasonality (figure 4(c)). This seasonal cycle showed the dependency of surface air temperature on the skin temperature CO₂ forcing. Surface skin temperatures would not respond similarly to surface air temperatures because of the impact of large-scale circulation [46]. The pattern was almost the same for both temperatures; it was more intense in the case of skin temperature. Surface air temperature may be controlled by large-scale circulation; it usually follows the skin temperature, depending on the turbulent fluxes from the surface. These fluxes play an essential role in air–land surface interactions. The air–land surface temperature change and its interaction can be expressed as an imbalance in the surface energy. The skin temperature over land is strongly coupled with sensible heat flux, which is part of the non-radiative feedback.

Zoom In Zoom Out Reset image size

Using the correlation method, previous studies analyzed the driving mechanisms on the interannual-to-decadal scales [8, 14, 19, 54] However, it cannot be applied to the seasonal scale and is also limited in sense of partial contribution to temperature. A previous study [22] specifically investigated the heat storage term using the surface energy budget for the Tibetan Plateau. This study found increased net radiation at the surface due to a decreased total cloud, and surface snow is favorable for EDW response to 4 × CO₂ over the Tibetan Plateau. We attempted to attribute elevation-dependent changes in terms of partial contributions to temperature change and its slope using the model framework.

Figure 5 illustrates the decomposed partial temperature contribution of different terms affecting the energy and heat budget in mountain regions, which can help elucidate the main contributors to the changes in elevation-dependent temperature changes. Such approaches have been extensively used to attribute polar amplification [48, 49]. The terms in equation (7) indicate the partial contributions to CO₂-induced warming and consist of radiative and non-radiative feedback processes. Among these terms, SAF, CRE, SWCS and LWCS represent the radiative processes and HSL and FLUX indicate non-radiative processes. Asia, Africa, and Antarctica have EDW signals and a decrease in EDW or even EDC on the Indian plate, while North America owns mixed signals with increasing CO₂. However, major contributors differed between mountain ranges due to their geographic locations. The LWCS mainly consists of water vapor and cloud feedback [11, 55, 56]. Meanwhile, the SAF is an important contributor to warming in mid- and high latitudes, although it is negligible in low-latitude mountain ranges. Other contributors counteracted this warming trend, including the CRE, SWCS, and LWCS. However, LWCS can seasonally contribute to enhanced EDW. It implies that the topographic process might be enhanced with CO₂ forcing and is consistent with earlier studies [22, 57], which also attributed a similar term [11, 58–60] as LWCS. The CRE was dominant in low-latitude regions, including the Indian plate and Africa (figures 5(c) and (d)). The mountains in South Africa experienced warming mainly because of the LWCS which experienced an increase with increasing CO₂ concentration. The EDW followed this trend, while other radiative terms did not show high sensitivity (figure 5(d)). In contrast, in North America and Asia, the SAF contributed to EDW, and changes in LWCS with increasing CO₂ concentration enhanced this trend. The SWCS imposed EDC. Furthermore, LWCS and CRE did not present clear signs thus contributing to the seasonality of EDW (figures 5(a) and (b)). Over the Antarctic, which has reduced seasonal ice and snow cover in the sensitivity scenarios (3 × CO₂ and 6 × CO₂), the strongest EDW was observed with high seasonality. The SAF followed the annual amplitude and season of EDW which is related to snow cover changes. The strong EDW response interpreted in the 9 × CO₂ scenario might be due to the polar region's ice/snow cover reduction. A recent study [8, 61] also indicated a current decrease in snow depth over Tibetan Plateau corresponding to EDW. In addition, the LWCS contributed to the warming trend with increasing values with increasing CO₂ concentrations (figure 5(e)).

Zoom In Zoom Out Reset image size

The radiative feedback was strong under CO₂ forcing and could explain the EDW pattern and its sensitivity over most of the investigated regions. The results so far describe the radiative feedback effects. Notably, the non-radiative feedback facilitated the radiative contributions and enhanced them. Their negative values (fluxes from land to atmosphere) indicated surface cooling and vice-versa, which ultimately raises the surface air temperature. For instance, the HSL exhibited a similar trend as EDW in North America. The heat stored in the land mass caused the surface temperature to rise, promoting energy release in the form of FLUX (figure 5(a)). The opposite trend can be seen over Africa, where HSL was negative, promoting energy shortage on the surface and triggering positive FLUX values (figure 5(d)). The increased non-radiative feedback (HSL and FLUX), followed by radiative LWCS over the Indian plate under 3 × CO₂ and 6 × CO₂ scenarios, resulted in EWD from non-radiative processes. However, this could not counterbalance the negative effect of CRE with positive non-radiative feedback, resulting EWC in the 9 × CO2 scenario (figure 5(c)).

The partial contribution of the different terms followed a distinct seasonality similar to that of EDW, which changes with rising CO₂ concentration (figures 4 and 5). The total elevation-dependent response and most partial contribution terms to temperature were the strongest in the respective hemisphere's seasonal cycle. Elevation-dependent temperature responses are diverse in each mountains regions with complex local feedback interactions. The SWCS alone did not show this distinct seasonality in any regions and had overall low values. In the summer months, the contribution terms and EDW values were stronger than those in the winter. However, the regions in the Northern Hemisphere did not show seasonality as distinct as those in the Southern Hemisphere. During the summer months in Asia and over the Antarctic CRE and SAF dominated over the radiative terms, respectively, and supported by LWCS, enhanced the seasonal EDW change. The non-radiative terms seemed to counteract these effects (figures 5(b) and (e)). Africa was an exception; here, most terms showed positive values and contributed to warming with elevation. Only slightly negative values were observed for the HSL. The dominant driver for seasonal EDW was LWCS, which was elongated and enhanced by positive CRE (figure 5(d)).

Studies on EDW mainly based on in-situ observation or satellite remote sensing. This is explained by the resolution dependence of micro-physical processes, where resolution dependence of EDW in CESM as thoroughly examined in a previous study using another climate model [19]. Modeling studies supported the assumption, that SAF was relevant in topographic regions where snowline changes [10, 62, 63]. In our paleo-climate setup strongest SAF was indicated for the Antarctic where a melt-down of the entire permanent ice-sheet in the 9 × CO₂ scenario would cause an EDW of more than 10 K km⁻¹ in the local winter months. CRE were also handled as major contributors, where temperatures are sufficient high [19]. CRE was most important for the EDC on Indian plate adding about −6 K km⁻¹ cooling trend in the 6 × CO₂ and 9 × CO₂ scenario, which was not balanced entirely by other processes.

Using present-day-based simulation, all mountain regions (figure 6) reveal EDC in future response to CO₂, which is a future paradox. Here chosen regions are slightly different from the Eocene as geography, and continental configuration is distinct in the present day. As the Indian plate is now part of Asia, we consider the Andes mountain range from South America, shown in figure 6(c). These surprising results suggest that future 4 × CO₂ responses are in contrast to our Eocene framework and the present agreement, which strongly suggests raising CO₂ is responsible for EDW in the future [22, 58, 64]. Even though it is a paradox between the past and the future responses, future EDC is dominantly controlled by LWCS overall mountain regions, underscoring the clear evidence for elevation-dependent water vapor and cloud feedback under future CO₂ warming. Numerous studies (e.g. [58, 64]) suggested that LWCS is a potential contributor to EDW; we accomplished this relationship in the past climate during Eocene. Figure 6 shows typical seasonal amplification in radiative feedback and non-radiative feedback. SWCS is the second most dominant feedback after LWCS over North America, Asia, Andes, and Africa. Increasing CO₂ in the future can lead to more significant warming at low elevations as we have more elevated mountains compared to the Eocene, resulting in EDC. Many studies also underlined that warming is happening to a certain elevation, supporting our argument. For example, a recent study [20] also infers that the Tibetan Plateau below 5 km exhibit EDW; it did not show any EDW above 5 km or in future climate projections. A study [65] demonstrated Antarctic exhibits low elevation warming in response to 2 ×CO₂ forcing. Additionally, this past-future linkage concludes that elevation-dependent temperature changes strongly depend on the geographical orientation of mountains rather than CO₂ concentration.

Zoom In Zoom Out Reset image size