3.1.1 Control simulation

We first run the weathering model using ERA5 reanalysis of climate fields (1981–2019 climatology average) as a control simulation. The resulting weathering rates (averaged over the selected parameterizations) are presented on a map in Fig. 2b. The total CO₂ consumption by weathering is 5.3 Tmol yr⁻¹. This control weathering field is very similar to the one published in Park et al. (2020, Fig. S5 of their contribution), who used the same model, parameter combinations and slope fields, but different climate fields (CRU TS v.4.03 for the temperature, Harris et al., 2014, and UNH/GRDC Composite Runoff Fields V1.0, Fekete et al., 1999).

The contribution of the different regions of Earth to this control weathering flux is shown in Fig. 2c (blue bars). The weathering flux from the Maritime Continent is about 11 % of the total flux, which is consistent with previous estimates (e.g., Molnar and Cronin, 2015). The main mountain belts (broader Himalaya, the Andes, and North and Central American ranges, including the Rocky Mountains) together contribute ∼22 % of the total. Another significant weathering contributor is the East African Rift (∼8 % of the total) owing to relatively high runoff rates, the substantial amount mafic silicate rock exposed and the relatively steep topography that fosters physical erosion – and hence, elevated chemical weathering fluxes. Taiwan and New Zealand, which are known hotspots of physical erosion rates, do not contribute significantly to the total chemical weathering flux because of the limited areal extent of their orogens.

3.1.2 Climatologies of El Niño and La Niña years

We then recompute the weathering rates with the same dataset, but keeping only the years of the climate time series identified as “El Niño” or “La Niña” (average of May to April, see Methods, Sect. 2.2). The anomaly of Earth-integrated weathering flux is $-\mathrm{90}\pm \mathrm{15.4}$ Gmol yr⁻¹ for El Niño years and $+\mathrm{138}\pm \mathrm{30.1}$ Gmol yr⁻¹ for La Niña (with ± denoting the standard deviation of the weathering anomaly over the selected parameterizations). These calculations assume that these climatologies are imposed over a long timescale such that the chemical weathering profiles adjust to result in these fluxes. Assuming that CO₂ degassing is constant, these weathering anomalies would in reality be compensated for by global warming or cooling so that the net silicate weathering flux anomaly is zero. If we extrapolate from the silicate weathering feedback computed with CESM climate simulations, indicating a global weathering flux increase of ∼0.4 Tmol yr^−1 ∘C⁻¹ of global warming by CO₂ (Fig. 8b, solid black curve), the weathering anomalies should then be compensated for by global mean temperature anomalies of +0.23 ^∘C for long-term El Niño conditions and −0.35 ^∘C for long-term La Niña conditions. These relatively small values of weathering anomaly (1.7 % and 2.6 % of the total control flux, respectively) are the consequences of a local increase and decrease in weathering rates that compensate for each other (Fig. 3). Nevertheless, the differences from 0 are robust among all the selected parameterizations (the absolute anomaly of weathering flux is ∼4 times larger than the standard deviation due to the parameter uncertainties) and approximately symmetrical between El Niño and La Niña conditions. These results support our hypothesis of climate cooling caused by the progressive onset of La Niña-like mean climatic conditions, though the magnitude of this cooling may not have been substantial.

Figure 3(a, b) Anomaly of modeled weathering rate with respect to control weathering (i.e., the field shown in Fig. 2b) for climatology in El Niño years (a) and La Niña years (b). The weathering rate is averaged over the parameter combinations (as in Fig. 2b). The color bar scale is a “symmetric logarithm”, and any value lower than 5×10³ mol km⁻² yr⁻¹ (in absolute value) is approximated as 0. Map projection and graticules are as in Fig. 1. (c) Bar plot of the anomaly of weathering flux integrated over the geographic regions shown in Fig. 2a with respect to ERA5 control weathering (blue bars in Fig. 2c) for climatology in El Niño or La Niña years. The average (main colored bars) and standard deviation (error bars) refer to the parameterizations, as in Fig. 2c. There are many regional changes of note as discussed in the text. For example, there is enhanced weathering over the Maritime Continent during La Niña years relative to El Niño years. The change in weathering rate is most dramatic over the Maritime Continent due to increased precipitation over land during La Niña years compared to being over the Pacific Ocean during El Niño years.

Local weathering anomalies are largely driven by changes in runoff rates. Many regions display dipolar patterns under El Niño conditions (Fig. 3a). For example, weathering fluxes increase on the Ethiopian Traps associated with El Niño conditions but decrease on the southern part of the East African Rift. They increase in the Andes near the Equator but decrease in the northernmost part of the range, and they decrease on the central Himalaya but increase in the western part. Hence, the anomaly of weathering flux integrated over those regions is close to zero (Fig. 3c). Most of the negative anomaly of El Niño weathering comes from the Maritime Continent. The weathering integrated over this region is almost equal to the total anomaly (Fig. 3c, purple bars). This result highlights the key role this area has potentially played in cooling Earth's climate over the last few million years. The signal from the Maritime Continent is expanded by a weathering decrease on the Deccan Traps and the Guiana highlands associated with El Niño conditions and partially offset by an increase on the Brazilian highlands (Fig. 3).

It is interesting to note that although the sign of weathering anomalies under La Niña conditions is the opposite of the ones under El Niño conditions nearly everywhere, the magnitude of those anomalies is not perfectly symmetrical (Fig. 3c).

3.2.1 Weathering and equilibrium temperature with reduced zonal gradient

In this section, we examine the simulations where the Pliocene SST field is applied in the 10^∘ S–10^∘ N band (10^∘ SN Pliocene SST) with the slab ocean version of CESM.

The standard CO₂ experiment for this configuration (299.4 ppmv, see Appendix B, Sect. B2) has a global mean 2 m temperature of 13.67 ^∘C, which is 0.25 ^∘C warmer than the pre-industrial control (at 284.7 ppmv). The precipitation anomaly (Fig. 5a) exhibits a large increase in the eastern Pacific, similar to El Niño events. This increase does not extend southward as much as for observed El Niño conditions (Fig. 1b), likely because we imposed SST perturbation only in the 10^∘ S–10^∘ N band. Instead of drier conditions around the Maritime Continent and wetter conditions in the western Indian Ocean characteristic of El Niño, precipitation increases along the Equator ($\sim \pm \mathrm{5}$^∘) across the Indian ocean and extends up to the eastern coast of Borneo over the Maritime Continent. The islands of Borneo, Sulawesi and New Guinea show a similar pattern of wetter conditions except on their western edge. The tropical Atlantic also shows significant differences with El Niño because the imposed SST in the Atlantic is warmer almost uniformly along the Equator (Fig. 4 b and d), which favors atmospheric convection. The runoff anomaly (Fig. 5a) reveals some similarities with the El Niño case (Fig. 1c): a decrease over India and Southeast Asia, an increase over China and the eastern side of the Himalaya, an increase in the Andes around the Equator, a decrease in the rest of equatorial South America (though the enhanced convection in the Atlantic generates wetter conditions along the Atlantic coast), an increase in the southeast part of South America, and an increase over the Congo basin. The largest discrepancies are found on the Maritime Continent and the East African Rift.

Figure 5Annual mean climatology of the slab ocean simulation with 10^∘ SN Pliocene SST at 299.4 ppmv anomaly with respect to the pre-industrial control. (a) Total precipitation and (b) continental runoff. Map projection and graticules as in Fig. 1. In each panel, non-significant differences are hatched (p value of Welch's t test > 0.1, see Appendix C).

Download

Figure 6Results for the 10^∘ SN Pliocene SST GEOCLIM simulation at C cycle equilibrium (i.e., with silicate weathering flux balancing control degassing). (a) Atmospheric CO₂ level (logarithmic axis) and (b) anomaly (with respect to pre-industrial simulation) of global mean 2 m temperature. The distribution – shown as a box plot – represents the ensemble of GEOCLIM parameterizations (similarly to Figs. 2c and 3c). The x axis of panel (a) is bounded so that a given CO₂ in panel (a) corresponds to the aligned global temperature in panel (b). This is with the approximation that the CO₂ is continuously logarithmic (from lower to upper bound), whereas GEOCLIM log(CO₂) interpolates climate fields between the CO₂ levels – the ones considered here being 250, 299.4 and 427.1 ppmv. Because of this approximation, the box plots in (a) and (b) are not perfectly aligned, but this misalignment is negligible. The pre-industrial CO₂ level and the 0 ^∘C anomaly are highlighted by the thin grey dashed or dotted vertical lines. Note that those two lines do not align because of the offset in the CO₂–temperature relationship of the 10^∘ SN Pliocene SST simulation with respect to the pre-industrial one.

Download

We then used the 10^∘ SN Pliocene SST simulations at several CO₂ levels to compute the equilibrium CO₂, which is the value at which the silicate weathering flux with 10^∘ SN Pliocene SST balances CO₂ degassing. More precisely, we consider control CO₂ degassing to be equal to the pre-industrial silicate weathering flux. Therefore, for each selected combination of GEOCLIM parameters (see Methods, Sect. 2.1), we imposed the silicate weathering flux with 10^∘ SN Pliocene SST as equal to the silicate weathering flux with the same parameter combination in the pre-industrial control simulation (i.e., computed with climate fields from pre-industrial slab simulation; this “control” weathering is shown in Fig. 2c, orange bars). Figure 6 shows this equilibrium CO₂ (Fig. 6a) and the global mean 2 m temperature corresponding to that CO₂ concentration (Fig. 6b). The Pliocene SST (within 10^∘ S–10^∘ N) generates climatic conditions less favorable for weathering, causing CO₂ to increase and global temperature to rise by ∼0.4 ^∘C, which is about twice that estimated using ERA5 reanalysis of El Niño years. This result is consistent with the fact that, in the 10^∘ SN Pliocene SST simulation, the west-to-east Pacific SST difference is reduced by more (∼2.5^∘) than in the climatology of El Niño years (∼1.5^∘). Note that though the climate sensitivity (global temperature difference for a CO₂ doubling) stays the same, there is an offset in the CO₂–temperature relationship, meaning that CO₂ of 284.7 ppmv does not correspond to a 0 ^∘C temperature anomaly (see also Fig. 8a).

Figure 7(a) Anomaly of the weathering rate of the slab ocean climate simulation with 10^∘ SN Pliocene SST at 299.4 ppmv with respect to pre-industrial control (averaged over the parameter combinations, similar to Fig. 3). (b) Bar plot of the anomaly of weathering flux from panel (a) integrated over the geographic regions (similar to Fig. 3c).

This decrease in weatherability (or weathering efficiency) can be investigated by looking at the anomaly of the silicate weathering rate of 10^∘ SN Pliocene SST at 299.4 ppmv with respect to the pre-industrial control weathering rate (Fig. 7). Though the global temperature is not the same (10^∘ SN Pliocene SST is 0.25^∘ C warmer), it still illustrates the causes of this decrease in weatherability. Similar to what was observed with El Niño and La Niña composites (Fig. 3), several positive and negative regional anomalies of weathering, driven by runoff changes, compensate for each other, leading to a weak warming. Nevertheless, this warming is consistent over the whole ensemble of parameter combinations. Figure 7b shows the integration of the weathering rate anomaly from Fig. 7a over the geographic regions presented in Fig. 2a. Unlike the El Niño case, the weathering anomaly integrated over the Maritime Continent is positive, thereby not contributing to the warming. Indeed, there are more increases than decreases in runoff on the Maritime Continent in the 10^∘ SN Pliocene SST simulation. The East African Rift also has an opposite behavior, with a decrease in weathering instead of a slight increase with El Niño. The last main difference concerns the category “the rest of Asia”, showing a large decrease in weathering. The decrease in runoff in India and associated weathering decrease overwhelm the increase in China, whereas those two balance out under modern El Niño conditions.

The weatherability decrease in the 10^∘ SN Pliocene SST simulation can also be appreciated by looking at the global temperature–total weathering flux relationship (Fig. 8b). The 10^∘ SN Pliocene SST curve (dashed blue) is systemically beneath the pre-industrial one (solid black) for any given global temperature. This means that regardless of the chosen background CO₂ degassing flux, the simulation with 10^∘ SN Pliocene SST will be warmer.

In summary, applying our estimation of Pliocene SST in the 10^∘ S–10^∘ N band – whose main feature is a flattening the tropical Pacific zonal gradient – generates a similar weatherability decrease and global warming to that estimated from El Niño (Sect. 3.1.2), though it is approximately twice as large. However, the mechanisms of that weatherability decrease are different due to differences in the regional patterns of simulated runoff.

3.2.2 Simulation with full Pliocene sea surface temperature

We now examine the slab ocean climate simulations wherein the full Pliocene SST field is applied. The most noticeable feature of the full Pliocene SST simulation is its mean temperature anomaly. As shown in Fig. 8a, for any given CO₂ level, the full Pliocene SST slab ocean simulation is about 2.5 ^∘C warmer than its counterpart with pre-industrial boundary conditions.

Figure 8(a) Global mean 2 m temperature (GMST) plotted against atmospheric CO₂ level for the three series of slab ocean general circulation model (GCM) simulations: pre-industrial boundary conditions, full Pliocene SST and 10^∘ SN Pliocene SST. Each dot represents an actual simulation, and the diamonds correspond to the “standard” CO₂ levels (284.7 ppmv for pre-industrial; see the previous section for Pliocene SST). (b) Total silicate weathering flux anomaly (with respect to pre-industrial) plotted against GMST for the same three series of slab ocean GCM simulations (same color code as panel a). Each circle represents a single simulation at fixed CO₂, with its level indicated by the color scale. The semi-transparent box plots show, for each simulation, the variability of weathering across the parameter combinations. In each panel, null weathering anomaly, pre-industrial GMST and pre-industrial CO₂ level are highlighted by the vertical and horizontal thin grey dotted lines.

Download

These warmer conditions at a given CO₂ level can be explained by analyzing the components of the radiative budget. Figure 9 shows the comparison between the full Pliocene SST simulation and pre-industrial simulations interpolated at 309.4 ppmv (same CO₂ level as Pliocene SST simulation) and at 538.3 ppmv (CO₂ level at which the global mean 2 m temperature equals the one of the Pliocene SST simulation). A significant signal of the full Pliocene SST simulation can be found in the net solar radiative flux at the top of the atmosphere (Fig. 9a and e). There are two peaks of incoming radiation in the tropics (around 5^∘ N and 10^∘ S, Fig. 9e), which do not appear just by raising CO₂ with pre-industrial boundary conditions. Such tropical peaks are caused by reduced cloudiness (Fig. 9d and h) due to less intense convection because of the reduced zonal temperature gradient, lowering Earth's shortwave albedo. Indeed, the cloud contribution to the peaks, estimated by the difference “full sky” minus “clear sky” (dashed lines in Fig. 9e), is almost 100 % of the signal. The SST pattern effect on global cloud radiative feedback has indeed been recognized (e.g., Zhou et al., 2017). Another peak in incoming solar radiation can be seen around 40^∘ N or S Fig. 9e and corresponds to a poleward shift in cloudiness at those latitudes (Fig. 9d), mostly discernable in the Southern Hemisphere. This feature, however, is also found in the pre-industrial simulation at 538.3 ppmv and may be a mere consequence of rising global temperature, though its amplitude is higher in the full Pliocene SST simulation. More solar radiation also enters at high latitude in the Pliocene simulations relative to the pre-industrial control due to the retreat of sea ice, also lowering shortwave albedo. Though it is countered by an increase in cloudiness (Fig. 9d and h), as can be seen in the cloud contribution in Fig. 9e, the net effect is still positive. Similar to the peak at 40^∘ N or S, this phenomenon is also present in the pre-industrial 538.3 ppmv simulation but is strengthened in the Pliocene SST simulation. On global average, the net solar radiation forcing of the full Pliocene simulation (with respect to pre-industrial at same CO₂) is +3.9 W m⁻², the tropical part alone being equivalent to +1.0 W m⁻² (equivalent if spread over the whole Earth). The latitude profile of downwelling longwave radiation at Earth's surface also shows some variations (Fig. 9 b and f), essentially following the anomaly of water vapor (Fig. 9c and g). The reduced meridional temperature gradient in the full Pliocene SST simulation is responsible for less water vapor in the tropics and more water vapor in the extratropics compared to the pre-industrial simulation at the same mean temperature (Fig. 9g), though it is mostly visible in the Southern Hemisphere. The downwelling longwave radiation at the surface exhibits the same behavior. The more significant drops in the tropics (still compared to the pre-industrial simulation at the same mean temperature) are likely due to the reduced Hadley circulation. The averaged downwelling longwave flux anomaly is not so different between the Pliocene SST simulation and the pre-industrial one at 538.3 ppmv (+14.0 and +15.8 W m⁻², respectively); both seem to be a consequence of a warmer atmosphere (and amplifier through a positive feedback loop). In summary, we can conclude that the main warming cause of the full Pliocene SST simulation is the decrease in cloudiness in the tropics due to a weaker convection and Hadley circulation, combined with a decrease in cloudiness around 40^∘ N and S and of sea ice at high latitude, these last two also being part of a positive feedback enhancing the “initial” warming.

Figure 9Comparison of zonal annual mean fluxes and climatology for the three simulations, which are pre-industrial boundary conditions log (CO₂ interpolated at 309.4 ppmv (indigo) or at 538.3 ppmv (teal) and full Pliocene SST (at 309.4 ppmv). (a) Zonally integrated net (incoming minus outgoing) solar (SW) radiative flux at the top of the atmosphere (TOA), (b) zonally integrated downwelling longwave (LW) radiative flux at the Earth's surface (SFC), (c) zonally averaged total precipitable water (vertically integrated) and (d) zonally averaged total cloud fraction (vertically integrated). All variables are annual mean climatology. Panels (e, f, g, h) are similar to (a–d) (respectively) but with 309.4 ppmv pre-industrial subtracted; the thin vertical grey dashed line highlights zero anomaly. The dashed color lines in (e) and (f) show the difference of full sky (SF) minus clear sky (CS), illustrating the contribution from clouds. The solid lines (full sky) are the regular values, as in (a) and (b). The “0 anomaly” in panels (e)–(h) is highlighted by a thin grey dotted vertical line. The color code holds for the entire figure.

Download

In term of precipitation, it is easier to compare the full Pliocene SST simulation with the pre-industrial simulation log(CO₂)-interpolated at the same global temperature to cancel the effect of increased precipitation with global warming. This precipitation (and runoff) difference is shown in Fig. 10. A redistribution of precipitation, with a decrease in the tropics and an increase in the extratropics, can be observed, superimposed on the anomalies already observed in the 10^∘ SN Pliocene SST simulation (Fig. 5): an increase in the eastern tropical Pacific and on the western boundaries of the tropical Indian and Atlantic oceans. The results are different for land precipitation. Despite the reduced meridional temperature gradient, precipitation and runoff mostly increase in tropical Africa and South America, with the exception of the western part of the Amazon basin and the East African Rift, where a drying tendency was already observed in the 10^∘ SN Pliocene SST simulation. This goes against our initial hypothesis based on the “wet gets drier, dry gets wetter” feature from Burls and Fedorov (2017). The Maritime Continent exhibits contrasting behavior, with wetter conditions on the western sides of the most equatorial islands and drier conditions on their eastern side. Such behavior was also present in the 10^∘ SN Pliocene SST simulation, with a larger precipitation increase on the western sides, but is highlighted here because the anomaly is computed with respect to the pre-industrial simulation at 538.3 ppmv instead of 284.7 ppmv: the increase on the eastern sides is less than the roughly uniform increase due to higher CO₂, whereas the increase on the western sides is larger. On global average, precipitation increases by 15 mm yr⁻¹ and continental runoff increases by 16 mm yr⁻¹.

Figure 10Annual mean climatology of the slab ocean simulation with full Pliocene SST at 309.4 ppmv with respect to the pre-industrial control log (CO₂)-interpolated at 358.3 ppmv. This last CO₂ level was chosen so that the two simulations in the difference have the same global mean 2 m temperature. (a) Total precipitation and (b) continental runoff. Map projection and graticules are as in Fig. 1. In each panel, non-significant differences are hatched (p value of Welch's t test > 0.1, see Appendix C).

The offset in global temperature for a given CO₂ level (in full Pliocene SST simulations with respect to pre-industrial ones) poses a technical problem: to go back down to pre-industrial temperature, with full Pliocene SST, CO₂ needs to be lowered to ∼140 ppmv. Such CO₂ levels are unrealistically low and are even lower than those during Pleistocene glacial maxima. Moreover, at this CO₂ level, plant stomates cease to operate normally in the land model of CESM, with huge consequences for the hydrologic cycle. Shut stomate behavior contrasts with the regular tendency for stomates, which is to be more open at lower CO₂, and severely reduces transpiration, eventually leading to increased global runoff, though the mean temperature is lower. Therefore, it is not possible to properly simulate silicate weathering with GEOCLIM. Increased runoff would lead to increased weathering at lower CO₂ (thereby killing the negative feedback). This effect is likely responsible for the “flattening” of the temperature–weathering curve at the lowest tested CO₂ levels (dashed red curve in Fig. 8b). On the other hand, the fact that plants cannot survive would undermine weathering, as plant roots are a major weathering agent. One may expect this second effect to be stronger and hence to see a collapse of total weathering flux with even lower CO₂. For this reason, we cannot perform the inversion to compute the equilibrium CO₂ wherein silicate weathering flux balances pre-industrial degassing. We instead analyze the weathering fluxes at fixed CO₂.

Figure 8b illustrates the silicate weathering fluxes of the pre-industrial (solid black curve) and full Pliocene SST (dashed red curve) simulations. Despite the higher continental runoff compared to the pre-industrial simulation at the same global temperature, the full Pliocene SST temperature–weathering curve and the pre-industrial curve are superimposed, except for very low CO₂ (<200 ppmv) Pliocene simulations. This result means that Pliocene SST does not generate any significant increase or decrease in weatherability – except at low CO₂, but this is likely due to CO₂–plant interaction and stomate closure that generates higher runoff. If the carbon cycle was left to reach balance, with pre-industrial degassing, CO₂ would decrease, but it is not possible to determine what would be the equilibrium temperature. The dashed red curves in Fig. 8b cannot be extrapolated to the lower bound because it would fall outside the range of validity of the weathering model.

The fact that higher mean runoff does not lead to higher weatherability means that the regions where runoff decreases matter more in terms of weathering. Comparing the full Pliocene SST simulation to the pre-industrial interpolated at the same global temperature (538.3 ppmv) reveals weathering decreases in the almost all of Southeast Asia, India, China, the East African Rift, and even the Brazilian highlands and the Himalaya (Fig. 11). The combined decrease in those regions – mostly the first four – is enough to counteract the increase in the rest of the tropical land masses (Africa and South America) and along most of the American Cordillera (from the Andes to the Cascade Range, Fig. 11a) as shown in Fig. 11b.

Figure 11(a) Anomaly of the weathering rate of the slab ocean climate simulation with full Pliocene SST at 309.4 ppmv; anomaly computed with respect to the pre-industrial simulation interpolated at 358.3 ppmv (see Fig. 10); (b) bar plot of the anomaly of weathering flux from panel (a).