Did high Neo-Tethys subduction rates contribute to early Cenozoic warming?

. The 58–51 Ma interval was characterized by a long-term increase of global temperatures ( + 4 to + 6 ◦ C) up to the Early Eocene Climate Optimum (EECO, 52.9– 50.7 Ma), the warmest interval of the Cenozoic. It was recently


Introduction
Based on paleotemperature proxies, a trend of decreasing global temperatures throughout the Late Mesozoic and Cenozoic has long been identified (e.g, Shackelton and Kennett, 1975;Zachos et al., 2001Zachos et al., , 2008;;Cramer et al., 2009;Friedrich et al., 2012).Climatic modeling suggests that this cooling mainly results from decreasing seafloor spreading and subduction rates, as well as increasing CO 2 removal through silicate weathering (Park and Royer, 2011;Godderis et al., 2014;van der Meer et al., 2014).During the Cenozoic, CO 2 consumption was mainly governed by the erosion of the Tethyan orogenic belt, and by continental drift, responsible for the arrival of highly weatherable basaltic provinces in the equatorial belt (Raymo and Ruddiman, 1992;Kent and Muttoni, 2013;Lefebvre et al., 2013).However, global cooling was interrupted by a long-term increase of global temperatures (+4 to +6 • C) and pCO 2 (∼ 450 to ∼ 1000 ppm) from 58 to 50.7 Ma, crowned by the Early Eocene Climate Optimum (EECO, 52.9-50.7 Ma), the warmest interval of the Cenozoic (Zachos et al., 2001;Beerling and Royer, 2011).Because conventional carbon cycle models compute important weathering rates at that time, they fail to reproduce this rise in temperature and atmospheric CO 2 without the addition of excess CO 2 compared to background CO 2 volcanic degassing rates (4-10 × 10 18 molCO 2 Ma −1 at present; Berner, 2004) (Lefebvre et al., 2013;Van der Meer et al., 2014).Carbonates also indicate that from ∼ 58.0 to 52.5 Ma this warming was characterized by a 2 per mil negative shift in marine and terrestrial δ 13 C, referred to as the Late Paleocene-Early Eocene (LPEE) by Komar et al. (2013).This drop in δ 13 C suggests an additional source of depleted CO 2 (i.e enriched in 12 C) or/and decreased net organic carbon burial (Hilting et al, 2008;Komar et al. 2013).In contrast, despite warm temperature, the EECO was associated with a rise in δ 13 C (Cramer et al., 2009), indicative of the addition of heavy CO 2 or/and alternatively by increased net organic carbon burial (e.g., Komar et al., 2013).Various origins of excess CO 2 have been proposed for both periods of the early Cenozoic.Most invoke the activity of large igneous provinces such as the North Atlantic Igneous Province (NAIP), since a mantellic source of CO 2 (δ 13 CO 2 ranging from −3 to −10 ‰) may be compatible with carbon isotope proxies for most of the period of warming (see Reagan et al. 2013 and references therein).Alternatively, Beck et al. (1995), Hilting et al. (2008) and Komar et al. (2013) proposed that large amounts of low-δ 13 C organic carbon were being stored in carbon capacitors separate from the oceanatmosphere-biosphere (e.g., peat, gas hydrates, permafrost) during the Paleocene.They were then massively released during the LPEE warming and progressively vanished during the EECO (Komar et al., 2013).Finally, among several other hypotheses, it was suggested that Neo-Tethys closure may have strongly controlled Cretaceous and early Cenozoic climates, up to the EECO, through the subduction of tropical pelagic carbonates (δ 13 C ∼ 0 ‰) under the Asian plate and their recycling as CO 2 at arc volcanoes (Edmond and Huh, 2003;Kent and Muttoni, 2008;Johnston et al., 2011).These authors argued that the tropical latitudes of the northern Neo-Tethys could have favored deposition of carbonaterich pelagic sediments on the Tethyan seafloor.In detail, Kent and Muttoni (2008) suggested that the Indian plate dominated this "carbonate subduction factory", with a major decrease in CO 2 production as India and Asia collided some 50 Ma ago.However, the same authors recently concluded for low CO 2 outgassing at the Tethyan arc, mainly as a result of low decarbonation during subduction (Kent and Muttoni, 2013).For Kent and Muttoni (2013), high CO 2 could be explained by less efficient weathering close to the EECO, rather than by additional CO 2 production.
In this contribution, we aim to test whether Neo-Tethyan closure, which was obviously associated with widespread arc volcanism, may or may not have had an impact on global warming during the LPEE and the EECO, keeping in mind that this hypothesis hardly conforms to available carbon isotope records during the LPEE.To this end, we first use a simple model that calculates the volume of sediments subducted along with Neo-Tethyan oceanic and Greater Indian margin lithospheres, and computes a range of CO 2 fluxes emitted at active arc volcanoes along the northern Neo-Tethys margin.A coupled climate-carbon cycle model (GEOCLIM) is then used to quantify the impact of CO 2 fluxes obtained from our model and that of Kent and Muttoni (2013), on Paleocene/Eocene pCO 2 and atmospheric temperature.Finally, in light of our results, we discuss the relevance of alternate hypotheses commonly cited to explain the LPEE and the EECO.
Evidence of latest Cretaceous and early Cenozoic subduction-related magmatic activity is widespread along, and restricted to the Eurasian margin.For example, in the Zagros mountains and Turkey (Pontides), widespread arc magmatism occurred during the Mesozoic and the Cenozoic (Sengör et al., 1988;Okay and Şahintürk, 1997;Barrier and Vrielynck, 2008;Agard et al., 2011;Eyuboglu et al., 2011).In southern Tibet, a long-lasting volcanic 'Gangdese' arc was active from Early Cretaceous to Eocene time (Ji et al., 2009), with a short-lived ignimbrite flare-up stage around 50 Ma coinciding with Tibetan Himalaya-Lhasa continental collision (Ji et al., 2009), followed by return of the arc to a background state until the Late Eocene (Sanchez et al., 2013).In Sundaland, Paleocene-Eocene magmatism was likely active since at least ∼ 63 Ma (e.g., McCourt et al., 1996;Bellon et al., 2004).

Volcanic CO 2 release during the LPEE and the EECO by the Carbonate Subduction Factory Model (CSFM)
CSFM is designed to calculate the amount of CO 2 produced during Neo-Tethys closure.It quantifies the Neo-Tethys volcanic arc gas output as a function of subduction flux of oceanic crust, pelagic sediments, and also of Indian margin sediments at the onset of Indian continental lithosphere subduction.Required input parameters are subduction rate, trench length, the thickness, density, carbonate and organic carbon content of sediments and oceanic crust, the decarbonation efficiency of subducted material, and the time-lag to gas emission at the surface.

Subduction rates and trench length estimates
Subduction rates of African, Arabian and Indian plates below Eurasia were calculated from plate motion reconstructions made with GPlates (http://www.gplates.org/)(Boyden et al, 2011) using time steps of 0.5 Ma, between 65 and 35 Ma.Given the controversy regarding the presence or not of continuous subduction in easternmost Neo-Tethys from the Late Cenozoic to the Eocene (Sundaland) (e.g., Whittaker et al., 2007;Hall, 2012), we did not consider Australia-Eurasia convergence and assess the potential role of Neo-Tethys subduction based on the central and western Neo-

G. Hoareau et al.: Early Cenozoic warming
Tethys alone.We used the reconstructed position of three points located on the western, central and eastern syntaxis of each plate, similar to Rosenbaum and Lister (2002), Alvarez (2010)  from India-Asia convergence rates.Given the uncertainties concerning the rate of subduction below widespread ophiolites, and the locations of these subduction zones, we chose to simplify our scenario by assigning all subduction to the zones indicated in Fig. 1.In addition, we assume that there was no active spreading within the Neo-Tethys since 65 Ma.
For each plate, subduction rates at each time step were corrected for convergence obliquity related to the orientation of the subduction trench (Annex 1).Trench lengths were set to 2500 km for Africa, 2900 km for Arabia (from the Levant fault to the Makran) and 2600 km for India (from the Makran to the Indo-Burma range), making a total length of 8000 km (similar to Johnston et al., 2011) (Fig. 1).Three ages, 55, 52.5 and 50 Ma, were tested for the onset of Greater Indian thinned continental lithosphere subduction beneath Eurasia, corresponding to a shift from pelagic to margin sediments on the Indian plate.Scenarios without India-Asia continental subduction were also run, to assess the maximum potential effects of younger collision.

Oceanic crust and pelagic sediments
In the model, all oceanic lithosphere has the same crust and sediment thickness.For the oceanic crust a constant thickness of 7 km and a density of 3 t m −3 were retained.Because subducted Neo-Tethyan crust was older than 40 Ma during the Paleocene/Eocene, it was ascribed a carbonate content of 0.2 wt% (Alt and Teagle, 1999).For pelagic sediments, we adopted a carbonate content of typical deep-sea carbonate oozes (90 wt%) (Kroenke et al., 1991), an organic carbon content of 1 wt%, a thickness of 200 m, and a density of 1.9 t m −3 , similar to uncompacted deep-sea deposits (Sykes, 1996).These values are close to those used by Edmond and Huh (2003) and Johnston et al. (2011) (CaCO 3 = 100 wt%; thickness = 200 m).They can also be compared to those calculated from the model of Kent and Muttoni (2013) (CaCO 3 = 100 wt%; thickness of ∼ 240-270 m from 65 to 50 Ma), who explicitly computed sediment thickness as a function of pelagic carbonate productivity and the timing of residence of the oceanic crust in the Neo-Tethyan equatorial belt zone.Note that sediment thickness may have been locally higher due to the presence of submarine fans or margin deposits such as carbonate platforms.However, paleogeographic reconstructions indicate that north of Greater India, the Paleocene-Eocene Neo-Tethys ocean was deep (Heine et al., 2004), favoring the predominance of pelagic deposition.
For Arabia and Africa, the same conclusion can be drawn from palinspastic reconstructions of Barrier and Vrielynck (2008), who showed that only a small proportion of margin sediments were subducted during the Paleocene and the Eocene, compared to deep sediments.

Continental crust and Indian margin sediments
For subduction of northern Greater Indian passive margin sediments, a simplified passive continental margin geometry consisting of a sedimentary succession overlying the basement was designed (Fig. 2).Basement-sediment and the upper sediment interfaces were modeled using sigmoidal functions.Their shape was inspired from the geometry of the North American Atlantic passive margin (Watts and Thorne, 1984), which may have been a modern analogue to the pre-collisional northern Indian passive margin (Brookfield, 1993).Total length of the margin sediment succession was set to 600 km, following Brookfield (1993) and in agreement with back-stripping reconstructions of Liu and Einsele (1994) and Guillot et al. (2008).Maximum sediment thickness was set to a mean value of 4 km (uncompacted), based on recent estimations of syn-rift/post-rift Neo-Tethyan margin sediment thicknesses of Sciunnach and Garzanti (2012).
Although the lithology of the margin was variable, the proportion of carbonate sediments and organic matter may have been important (Beck et al., 1995;Liu and Einsele, 1994;Sciunnach and Garzanti, 2012).Average contents of 50 and 1 wt% were chosen for carbonate and organic carbon content, respectively.Uncompacted margin sediments were given a density of 2 t m −3 as calculated from data of Sciunnach and Garzanti (2012).

Oceanic crust and pelagic sediments
In the "carbonate subduction factory" model, CO 2 produced during oceanic subduction processes originates from deep metamorphic decarbonation of subducted crust and sediments (carbonate and organic matter), and is assumed to be released at volcanic arcs following partial melting of the subducting oceanic crust and metasomatism of the overlying mantle (Hilton et al., 2002;Gorman et al., 2006).This common statement was followed in the CSFM model (Fig. 3), therefore ignoring possible additional CO 2 sources, in particular decarbonation of the overlying crust (Lee et al., 2013).The amount of CO 2 released through arc volcanism was calculated as follows: first, for each time step, the total volume of subducting sediment and crust was computed.We assumed this volume to be similar to that encompassing metamorphic carbon loss in the sub-arc zone (∼ 120 ± 40 km depth; England and Katz, 2010) (i.e., no variation of volume during the subduction process before decarbonation).Then, the amount of CO 2 emitted at the surface was calculated from arc decarbonation efficiency, defined as the ratio of the volcanic gas CO 2 flux to the input of subducted carbon (e.g., Johnston et al., 2011) and Muttoni (2013) to perform similar calculations (10 %), based on 10 Be data in arc volcanoes of Central America (Tera et al., 1986).Finally, the time lag between decarbonation at depth and gas emission at the surface was set to 2 Ma, averaging timescales of Turner (2002) (0.4 to 4 Ma).

Continental crust and Indian margin sediments
Due to the lack of aqueous fluids in continental crust, continental subduction zones are expected to be devoid of significant syn-subduction arc volcanism in the overlying plate (Zheng, 2012).Although volcanism may have continued in Tibet after 50 Ma (Ji et al., 2009;Rohrmann et al., 2012), in the model oceanic slab-related metamorphic decarbonation and magma generation was considered to last until the arrival of the continental lithosphere at sub-arc depth (i.e., 80 km) (Fig. 3).Using preferred geometric parameters of Leech et al. (2005) for subduction of the Indian plate, this depth is reached ∼ 1.5 to 2 Ma after the initiation of continental subduction.Despite cessation of volcanic activity, subduction of continental margin sediments may have been associated with active CO 2 degassing at springs or vents as a result of efficient metamorphic sediment decarbonation at T > 300 • C (e.g., Becker et al., 2008;Evans et al., 2008).Kerrick and Caldeira (1993) suggested that limited collision-related prograde metamorphism of marly lithologies may induce a CO 2 loss of ∼ 10 wt%, equivalent to a decarbonation efficiency of ∼ 50 % for sediments with a carbonate content of 50 wt% (= 22 wt% CO 2 ).This value may represent an upper estimate as shown by thermodynamic modeling of Massonne (2010).Above-mentioned studies focus on collision rather than continental subduction, for which to our knowledge no estimations of CO 2 outgassing fluxes or decarbonation efficiency are available.To avoid overestimations of CO 2 production, we assumed that only limited margin sediment decarbonation may have occurred after the onset of continental subduction at low-grade conditions, with a 1 to 10 wt% efficiency.Time necessary for subducted margin material to reach the 300 • C isotherm after the onset of continental subduction at ∼ 55-50 Ma (corresponding to 25 km depth with a normal-subduction geothermal gradient of 15 • C km −1 ) was set to 0.5 Ma, as calculated with parameters of Leech et al. (2005).Circulation of CO 2 -rich fluids along large-scale collision-related thrust detachments has been proposed as an efficient way to promote degassing at the surface (e.g., Kerrick and Caldeira, 1993;Becker et al., 2008).Following Skelton (2011), who suggested that gas produced during low-  grade metamorphism may be rapidly released to the surface (∼ 4000 yr), we considered immediate release of CO 2 to the atmosphere (Annex 2).

Tethyan subduction rate
During the Paleocene (65 to 56 Ma), the mean subduction rate (i.e., all plates) has a constant value of ≈ 5.5 cm yr −1 (Fig. 4a).Increased rates (up to 8.3 cm yr −1 ) are computed between 56 and 53 Ma, before a gradual decrease to 3 cm yr −1 at 35 Ma.India-Asia convergence, which reaches up to 16.7 cm yr −1 at 53-52 Ma, exerts the main control on high early Cenozoic subduction rates.

Greenhouse gas production
It is important to note that decarbonation efficiencies may have strongly varied with time, depending particularly on the plate age and sediment thickness (Peacock, 2003;Connolly, 2005;Johnston et al., 2011).However, according to Johnston et al. ( 2011) the decarbonation efficiency is only roughly correlated with convergence (subduction) rate.Therefore, excess CO 2 fluxes calculated at minimum (15 %) and maximum (60 %) efficiencies correspond to extreme scenarios that very likely encompass true excess CO 2 fluxes related to Neo-Tethys closure.

With Indian continental subduction
Decarbonation of Greater Indian margin sediments, added to the last volumes of pelagic sediments at sub-arc depth, results in a peak of CO 2 production ∼ 2 Ma after the onset of continental subduction, considering a constant decar- bonation efficiency (Fig. 4d).In our model, continental subduction must start at 52.5 Ma (consistent for example with stratigraphic arguments of Najman et al., 2010 and paleomagnetic arguments of Huang et al., 2015b) for maximum CO 2 emissions to occur at ∼ 51 Ma, i.e. coeval to maximum recorded temperatures during the EECO (Zachos et al., 2008).In this case, CO 2 degassing flow rates are in the range 0.6-3.35× 10 18 molCO 2 Ma −1 (1/15-10/60 % efficiencies for margin/pelagic sediments, respectively), corresponding to 6-84 % of the modern CO 2 outgassing rate.

Modeling the impact of neothetys closure
To test the influence of calculated excess CO 2 fluxes on Paleocene/Early Eocene climate, we carried out simulations using the GEOCLIM model (Donnadieu et al., 2004;Goddéris et al., 2008).This model couples a 3-D general circulation model (GCM) called FOAM (Jacob, 1997) to a box model of geological carbon-alkalinity cycles called COM-BINE (Goddéris and Joachimski, 2004).The GCM FOAM is used in mixed-layer mode, where atmosphere is linked to a 50 m mixed-layer ocean, which parameterizes heat transport through diffusion, in order to reduce computation time (one GEOCLIM simulation needs up to 12 GCM simulations).This GCM is forced by a large range of pCO 2 (200 up to 4200 ppmv) to generate an offline catalogue of continental air temperature and continental runoff with a spatial resolution of 7.5 • long × 4.5 • lat.For each corresponding atmospheric pCO 2 value, the GEOCLIM model calculates the temperature and the runoff of each grid cell through a linear interpolation procedure from the climatic catalogue.This procedure is repeated until a steady-state is reached that corresponds to a stable atmospheric CO 2 and temperature.The model uses an ocean geometry divided into two polar oceans (including a photic zone and a deep ocean reservoir), a low-to mid-latitude ocean (including a photic zone, a thermocline and a deep ocean reservoir), two epicontinental seas (both with a photic zone and a deep epicontinental reservoir) and the atmosphere.A full description of GEOCLIM and its components COMBINE and FOAM can be found in Goddéris and Joachimski (2004) and Donnadieu et al. (2006).

GEOCLIM simulations
We first calculated the steady-state pCO 2 , assuming that the total CO 2 consumed by continental silicate rocks weather-ing equals the total solid Earth CO 2 degassing flux (Walker et al., 1981).Due to the non-consensus about the Earth degassing rate for the last 200 Ma, the degassing flux was assumed constant and fixed at a modern value of 6.8 × 10 18 molCO 2 Ma −1 , which is required in the model to balance the global consumption through the weathering of silicate lithologies (Donnadieu et al., 2006).Each terrestrial grid was prescribed a similar lithology, in which basalt weathering reaches a 30 % contribution of the total silicate weathering flux taken at present day (Dessert et al., 2003) (similar to UNI configuration of Lefebvre et al., 2013).Lefebvre et al. (2013) have shown that with this configuration steady-state pCO 2 is similar at 65, 52 and 30 Ma (320-350 ppm), despite variations in paleogeography.An Early Eocene (52 Ma) paleogeographic reconstruction was thus used in the simulation, which runs from 65 to 40 Ma.Land-ocean configuration was built from a synthesis of paleomagnetic data and geologic constraints (Besse and Courtillot, 2002;Dercourt et al., 1993).Obliquity and radiation solar constant were assumed to equal present-day values.
The main geological forcing tested in the simulation is the additional CO 2 fluxes calculated from CFSM.CO 2 fluxes of Kent and Muttoni (2013) have also been tested, using decarbonation efficiencies of 15 and 60 % in addition to the original value of 10 % proposed by the authors.Computed CO 2 outgassing rates resulting from Neo-Tethys closure were integrated to GEOCLIM, in an age step of 1 Ma.

pCO 2 evolution during the LPEE and the EECO
If minimum decarbonation efficiency (15 %) is considered, pCO 2 increase following excess CO 2 flux is negligible (Fig. 5c).For example, the addition of continental subduction, which results in higher CO 2 fluxes, allows reaching maximum pCO 2 of only 360-365 ppm at 51 Ma (i.e.close to steady state values).
If maximum decarbonation efficiency (60 %) is considered, calculated excess CO 2 fluxes lead to pCO 2 of 430-450 ppm from 65 to 54 Ma (Fig. 5c).Without continental subduction, between 54 and 51 Ma, pCO 2 increases up to 500-550 ppm.It then decreases to steady state pCO 2 values from 48 Ma.With continental subduction, pCO 2 can reach much higher values.In the preferred scenario (initiation of continental subduction at 52.5 Ma), pCO 2 strongly increases at 54 Ma up to reach a peak of ∼ 770 ppm at 51.5 Ma.It then decreases to values close to steady state at 47 Ma (Fig. 5c).If continental subduction begins at 55 Ma, a peak of similar amplitude (770 ppm) occurs at 53 Ma (Fig. 5d).In contrast, a 50 Ma age for the initiation of subduction results in two peaks of smaller amplitude, with pCO 2 values of ∼ 520 and ∼ 570 ppm at ∼ 52 and 48 Ma, respectively (Fig. 5d).

Impact of Neo-Tethys closure on
Paleocene-Eocene climate It has long been suggested that Paleocene/Eocene warming was not due to an increase of mantle degassing, calling for additional sources of atmospheric CO 2 (Engebretson et al., 1992;Kerrick and Caldeira, 1993;Hilting et al., 2008;Van der Meer et al., 2014).Kerrick and Caldeira (1993) Edmond and Huh (2003), Johnston et al. (2011) and Kent and Muttoni (2013).Values calculated from data of Kent and Muttoni (2013) (< 1.3 × 10 17 molCO 2 Ma −1 from 80 to 50 Ma) fall largely below those required by modeling of Lefebvre et al. (2013) to reach estimated pCO 2 , largely because of their choice of limited decarbonation efficiency during subduction (10 %).In contrast, estimations of Edmond and Huh (2003) and Johnston et al. (2011) vary from 0.5 to 4 × 10 18 molCO 2 Ma −1 for the entire Tethyan arc.According to results of Lefebvre et al. (2013), the higher range of these values should allow to sustain a warm climate during the Paleocene and the lower Eocene.However, excess CO 2 flux calculations of Edmond and Huh (2003) and Johnston et al. (2011) represent only average values based on simple assumptions such as a constant subduction rate for the entire Upper Cretaceous-Lower Cenozoic.Nevertheless, these estimates are generally higher than those calculated using CFSM.We rather suggest that between 65 and ∼ 55 Ma (i.e., before the oldest possible age of Indian continental subduction), Neo-Tethys closure may have released less than ∼ 10 18 molCO 2 Ma −1 , in particular owing to subduction rates lower than the one used by Edmond and Huh (2003) and Johnston et al. (2011) (∼ 5.5 vs. 8 cm yr −1 , respectively).Using the GEOCLIM model, CO 2 outgassing values obtained with maximum decarbonation efficiency (60 %) allow to reach a pCO 2 of ∼ 430 ppm, which is in agreement with proxies for the Early Paleocene (65-60 Ma) (Beerling and Royer, 2011) (Fig. 5c).Therefore, our modeling suggests that high decarbonation efficiency was a prerequisite for the "carbonate subduction factory" to have a significant impact on global climate at that time.In addition, GEOCLIM seems unable to explain the onset of Paleocene-Eocene warming at 58 Ma, coevally to an increase of at-mospheric pCO 2 (Beerling and Royer, 2011) (Fig. 5a, b,  c).Similar conclusions can be drawn from climatic simulations performed with excess CO 2 fluxes of Kent and Muttoni (2013), even using a decarbonation efficiency of 60 % www.clim-past.net/11/1751/2015/Clim.Past, 11, 1751-1767, 2015 (Fig. 5e).In that case, pCO 2 only steadily increases during the Paleocene and remains lower than values suggested by proxies (Fig. 5c).
A more significant contribution of Neo-Tethyan closure to global warming may have occurred close to the EECO (∼ 52-49 Ma), owing in particular to an increase of the Indian subduction rate (Fig. 4a).This contribution is also conditioned by maximum decarbonation efficiency, and by the onset of Indian continental subduction at ∼ 52 Ma.With these two conditions fulfilled, the model allows to reach pCO 2 values lower than, but close to proxy ones (< 850 ppm vs. ∼ 1000 ppm, respectively) (Fig. 5b).Based on calculations of climate sensitivity to atmospheric CO 2 of GEO-CLIM performed by Godderis et al. (2014) (2.4 • C for a pCO 2 doubling at 52 Ma), they may have resulted in a related global atmospheric temperature increase of ∼ 2 • C compared to the Paleocene.In contrast, if India-Asia continental subduction occurred much later (i.e., equivalent to no collision), Neo-Tethys contribution to the EECO remained negligible, even with a decarbonation efficiency of 60 % (Fig. 5c).Calculations performed with input data of Kent and Muttoni (2013) lead to the same interpretation for the Early Eocene (Fig. 5e).
Finally, our study allows to clearly moderate the impact of the Neo-Tethyan "carbonate subduction factory" on Paleocene-Eocene greenhouse, at odds with Edmond and Huh (2003), Johnston et al. (2011) and Kent and Muttoni (2008), but in accordance with recent conclusions of Kent and Muttoni (2013) and Lee et al. (2013).As a consequence, the strong decrease of CO 2 production after India-Asia collision was not a driver of pCO 2 decrease and global cooling recorded after the late Eocene (Kent and Muttoni, 2013).

Large igneous provinces
Since the role of Neo-Tethys closure on the onset of the LPEE and the EECO has likely been limited, other sources of excess greenhouse gases should be called for.These should ideally explain the decrease of marine and terrestrial δ 13 C during the LPEE, and its slight increase during the EECO (Zachos et al., 2001) (Fig. 5a).Numerous geological explanations have previously been postulated for the entire early Cenozoic greenhouse, among which a flare up in the activity of igneous provinces is the most common (e.g., Eldhom and Thomas, 1993;Reagan et al., 2013) 2013), Fig. 6), the calculated excess CO 2 flow rate falls to only ∼ 2.3 × 10 17 molCO 2 Ma −1 on average, 1 order of magnitude below fluxes necessary to reach a pCO 2 comparable to proxies using the GEOCLIM model.This value can be compared to that of Eldholm and Thomas (1993), who calculated that more than 2.3 × 10 18 molCO 2 may have been released to the atmosphere by the NAIP only, from 58 to 52 Ma (revised to 61-53 Ma by Menzies et al., 2002), corresponding to a flux of ∼ 3 × 10 17 molCO 2 Ma −1 .We infer, on the base of GEOCLIM modeling, that the effect of enhanced magmatic activity on LPEE warming and the EECO may also have been limited, unless related fluxes have been severely underestimated in the literature.As discussed later, this conclusion is consistent with previous studies of carbon cycle dynamics during the early Cenozoic (e.g., Hilting et al., 2008;Komar et al., 2013).

Metamorphic decarbonation
Lee et al. ( 2013) argued that the decarbonation of platform carbonates stored on the continental upper plate during subduction-related magmatism may have been far more efficient in driving early Cenozoic greenhouse than the activity of igneous provinces or of the "subduction factory".Lee et al. (2013) calculated that global CO 2 degassing could have reached 3.7-5.5 times the present day value from ∼ 140 to 50 Ma, making 2.2-3.7 × 10 19 molCO 2 Ma −1 for a present day value of 6.8 × 10 18 molCO 2 Ma −1 (as calculated with GEOCLIM).Cooling initiation during the late Eocene would then have resulted from a transition from a continentaldominated to an island arc-dominated world ca.52 Ma.The EECO would represent a "last spurt of CO 2 production associated with an Eocene magmatic flare-up in western North America", based on previous work of Nesbitt et al. (1995) and Kerrick and Caldeira (1998).In light of the GEOCLIM model, we argue that CO 2 fluxes of Lee et al. ( 2013) are clearly overestimated for the EECO, as pCO 2 close to proxies would be obtained for excess fluxes lower by approximately 1 order of magnitude.If they are applied to the EECO, fluxes calculated by Kerrick and Caldeira (1998) for the 60-55 Ma period (3 × 10 18 molCO 2 Ma −1 ) seem more reasonable.Similar to the decarbonation of pelagic carbonate sediments, crustal decarbonation related to magmatic or metamorphic events should lead to a positive shift of exogenic δ 13 C ( Lee et al., 2013), in agreement with proxies for the EECO.In contrast, the LPEE was characterized by a related negative shift in δ 13 C (Zachos et al., 2001), suggesting additional or alternate sources of excess isotopically light CO 2 .

Organic carbon sources
Several authors have thus proposed organic carbon (C org ) to be a significant source of excess CO 2 during the LPEE and/or the EECO, mostly based on carbon cycle models (e.g., Beck et al., 1995;Kurtz et al., 2003;Hilton et al., 2008;Kroeger and Funnel, 2012;Komar et al., 2013).The importance of organic carbon dynamics was clearly highlighted by Hilting et al. (2008), who used a carbon cycle model tuned with marine δ 13 C data to calculate Paleocene and Eocene pCO 2 values.These authors managed to reproduce pCO 2 values globally consistent with observations, even though background volcanic-metamorphic CO 2 degassing was kept constant.
According to their simulation, large changes in pCO 2 (and thus, temperature) may occur independently of the endogenic carbon cycle.Similar conclusions were drawn by Komar et al. (2013) on the basis of coupled LOSCAR-GEOCARB carbon cycle modeling.They showed that a mantellic source of excess CO 2 during the LPEE would have led to a deepening of the CCD much more important that evidenced from observations.Although based on a different approach, this conclusion is in good accordance with our suggestion of a moderate impact of LIPs on the LPEE (and the EECO).Instead, Komar et al. (2013) proposed that perturbations of the carbon cycle observed during the LPEE were likely controlled by a decrease of net organic carbon burial, either through increased C org oxidation, or through decreased C org burial.Suitable sources of organic carbon to the exogenic system include methane hydrates, which may have accumulated in marine sediments during the early Paleocene, and collapsed during the LPEE (Komar et al., 2013), terrestrial organic matter previously accumulated in swamps (Kurtz et al., 2003), or important reservoir petroleum generation (Kroeger and Funnel, 2012).However, the timing of maximum hydrocarbon production calculated by Kroeger and Funnel (2012) likely occurred during the EECO, which is hardly reconcilable with coeval increase in marine δ 13 C unless net C org burial was significantly higher than during the LPEE (e.g., Kurtz et al., 2003).Finally, Beck et al. (1995) postulated that Neo-Tethyan marine organic matter accumulated on Eurasian and Greater Indian margins may have been oxidized during India-Asia collision and subsequent exhumation, provided collision occurred no later than ∼ 60 Ma.About 1.6 × 10 18 molC Ma −1 may have been released during the first 4 Ma of the LPEE, enough to explain the concurrent negative shift in δ 13 C. Using our model, we calculate that the organic carbon contained within Greater Indian margin alone (∼ 3.8 × 10 6 km 3 ) amounts to ∼ 8 × 10 18 molC (for a sediment organic carbon content of 1 wt%), corresponding to a flux of ∼ 2 × 10 18 molC Ma −1 (i.e., close to estimates of Beck et al. (1995)) if all C org was oxidized during exhumation.This was probably not the case, and our estimate is likely overestimated.Nevertheless, it shows that oxidation of Neo-Tethyan marine C org may have contributed to the LPEE if collision occurred earlier than assumed in our model (e.g., Hu et al., 2015), to an extent that deserves to be quantified more accurately in future studies.
In contrast to the LPEE, the rise in marine δ 13 C from 52.5 to 50 Ma suggests that the EECO was characterized by increased net organic carbon burial (as proposed by Komar et al. (2013) with the methane hydrate hypothesis), or as we test in this paper, by the addition of excess CO 2 derived from one or several sources with heavier δ 13 C signatures.Both explanations need to be reconciled with recent observations that silicate weathering may have been reduced during the EECO, as discussed below.

Are modeled silicate weathering fluxes overestimated for the EECO?
Most carbon cycle models agree that during the Early Eocene volcanic degassing alone was insufficient to sustain the high pCO 2 values required by proxies, due to important weathering rates at that time (e.g., Berner, 2006;Lefebvre et al., 2013;Komar et al., 2013).For example, Berner (2006) found, based on the time evolution of seawater 87 Sr / 86 Sr that weathering was mainly controlled by increased basaltic alteration, resulting in a pCO 2 of ∼ 700 ppm at 50 Ma, i.e. lower than observations (∼ 1000 ppm).Indeed, decreasing of seawater Sr isotopic signature during the Paleocene and the Early Eocene is consistent with the alteration of igneous provinces such as the Deccan Traps or the NAIP (Hoddell et al., 2007).In detail, most paleogeographic reconstructions show that the highly weatherable Deccan traps reached the equatorial humid belt (between 5 • S and 5 • N), where weathering is maximum, at ∼ 55 Ma, with a maximum area of between ∼ 50 and ∼ 35 Ma (Dercourt et al., 1993;Besse and Courtillot, 2002;Van Hinsbergen et al., 2011a, 2012).Accordingly, taking explicitly into account the impact of paleogeography on the long-term carbon cycle as done by the GEOCLIM model has led Lefebvre et al. (2013) to suggest that the EECO was characterized by high weathering rates related to weathering of the Deccan traps.As a consequence, the model calculates a very low equilibrium pCO 2 of 340 ppm at that time.Calculations of Kent and Muttoni (2013) also pointed to an increase of silicate weathering, and thus of CO 2 consumption during the LPEE and the EECO due to the arrival of Greater India in the equatorial humid belt at that time.
Even though most models and reconstructions seem to agree with the hypothesis of increased weathering during the Early Eocene, several proxy-based observations rather suggest that silicate weathering may have been reduced specifically during the EECO.The first one is based on the estimation of the carbonate compensation depth (CCD) during the Paleogene (Hancock et al., 2007;Leon-Rodriguez and Dickens, 2010;Pälike et al., 2012;Slotnick et al., 2014).During the LPEE, deep-sea carbonate records show a progressive deepening of the CCD that can be attributed to enhanced alkalinity supply to the oceans as a result of enhanced weather-ing (Komar et al., 2013).During the EECO, high pCO 2 values (∼ 1000 ppm) similarly should have sustained high silicate weathering and thus favored a deep position of the CCD.In contrast, available records suggest its strong shoaling at that time (Pälike et al. 2012;Slotnick et al., 2014), which according to Komar et al. (2013) hardly conforms to the intense weathering deduced form carbon cycle models.However, previous GEOCLIM modeling showed that a constant silicate weathering flux does not mean a fixed pCO 2 (and thus a fixed CCD), due to the major role played by continental configuration on pCO 2 values (Donnadieu et al., 2006).The second observation is based on δ 7 Li chemistry of Paleocene and Eocene marine sediments, compiled by Misra and Froelich (2012).These authors argued that Li isotopes allow discriminating between periods of high tectonic uplift associated with important physical and chemical weathering, and periods of low alteration (Froelich and Misra, 2014).According to Froelich and Misra (2014), the LPEE and the EECO were characterized by slow weathering rates, as shown by a low and constant δ 7 Li trend.This strong discrepancy with previous interpretations is attributed to the absence of continental reliefs at that time, preventing significant weathering of uplifted, fresh silicate rocks.According to Froelich and Misra (2014), only moderate additional CO 2 may have allowed increasing pCO 2 and global temperature up to the end of the EECO.Note that for the LPEE, this interpretation is in contradiction with that of Komar et al. (2013) based on the CCD.In addition, recent modeling of Vigier and Godderis (2015) suggests that the oceanic δ 7 Li record of the Early Eocene could equally be explained by intense soil production rates (i.e, intense chemical weathering).These contradictory observations show that the intensity of silicate weathering during the LPEE and the EECO still suffers from strong uncertainties, as already highlighted by Kent and Muttoni (2013).Additional proxy-based observations are thus needed to calibrate weathering rate values obtained from models, which for example still lack the explicit integration of uplift on carbon cycle evolution (Lefebvre et al., 2013).

Conclusions
In order to test the role that Neo-Tethys closure may have exerted on warm Paleocene/Early Eocene climate through CO 2 degassing at arc volcanoes, we have calculated the volume of buried pelagic sediments and associated volcanic CO 2 release during the LPEE and the EECO, and its impact on atmospheric pCO 2 and atmospheric temperature at that time.To do so, we have applied most recently published convergence rate parameters and decarbonation efficiencies to a simplified Neo-Tethyan geometry, and integrated calculated excess CO 2 fluxes in a state-of-the-art carbon cycle model (GEOCLIM).
We show that Neo-Tethys closure was able to bury significant volumes of pelagic sediments at that time.The inset of Indian continental subduction at 55-50 Ma may have potentially given rise to important volumes of excess CO 2 , through decarbonation of thick margin sediment accumulations.However, GEOCLIM modeling demonstrates that these volumes do not generally allow reaching pCO 2 (and thus temperatures) as high as those inferred from geochemical proxies.Atmospheric CO 2 concentration may have only been able to reach significantly high values during the EECO (up to 770 ppm), but only if decarbonation efficiency was at its maximum at that time.This finding leads us to temper the impact of Neothetys closure on the LPEE and the EECO, calling for additional sources of excess CO 2 .
Among these, GEOCLIM modeling suggests that in light of available published data, the volume of CO 2 released by Large Igneous Province volcanism was 1 order of magnitude too low to have had a significant impact on climate during the Paleocene and the Early Eocene.Other recently proposed mechanisms of CO 2 release such as a decrease of net organic carbon burial may have been more efficient in driving Paleocene-Eocene warming.
Finally, an alternate explanation may be that CO 2 consumption may have been lower than suggested by carbon cycle models, calling for a better calibration of early Cenozoic weathering rates.

Figure 1 .
Figure 1.Simplified paleogeographic maps showing the positions of Africa, Arabia, India and Eurasia at 65, 52.5 and 35 Ma (3-D Globe projection; rotation poles of Müller et al. (2008), fixed Eurasian frame).The bearing and/or length of subduction trenches used in the model are represented as black lines.Flow lines (65-35 Ma) of three points representative of the central syntaxis of each plate are also reported (see text for present locations).Extension of Greater India during the Upper Cretaceous is based on the Greater India Basin hypothesis of Van Hinsbergen et al. (2012).

Figure 2 .
Figure 2. (a) Geometry of the Eastern US Atlantic coastal margin showing analogous positions for tectonic units of Indian passive margin (modified from Brookfield, 1993).(b) Geometry of Indian margin used in Carbonate Subduction Factory Model (CSFM).
(Annex 2).Decarbonation efficiency values were based on modern decarbonation efficiencies calculated recently by Johnston et al. (2011) using a mass balance approach, and those computed by Gorman et al. (2006) based on thermodynamic modeling.Decarbonation efficiencies at sub-arc depth of Johnson et al. (2011) vary from 0.1 to 70 % for 10 subduction trenches, with most values ranging between 18 and 54 %.They are quite similar to the mean values of Gorman et al. (2006) (∼ 16 and ∼ 63 %, if volcanic CO 2 is derived from decarbonation at sub-arc depth only, or both at fore-arc and sub-arc depths, respectively), which are based on 41 subduction zones.We have retained values of 15 to 60 %.Note that such values exceed the one used by Kent

Figure 3 .
Figure 3. Sketches illustrating the carbon input and outputs considered in the model during subduction.(a) General sketch for subduction of oceanic crust and pelagic sediments (Africa, Arabia and India before Indian continental subduction).Carbon deposited as carbonate and organic carbon in pelagic sediments, as well as crustal carbonate, is partly recycled at sub-arc depth and incorporated in arc magmas.(b) Similar sketch at the inset of Indian continental subduction.(c) Subduction of northern Indian margin before it reaches sub-arc depth.Carbon originating from low-grade metamorphism of Indian margin sediments is partly released to the atmosphere, making additional greenhouse gas to that released at arc volcanoes.(d) Arc volcanism stops as Indian continental crust reaches sub-arc depth.

Figure 4 .
Figure 4. (a) Calculated mean Tethyan subduction rate over the period 65-35 Ma, compared with individual subduction rates of Africa, Arabia and India beneath Eurasia.Upper and lower limits of the shaded areas are maximum and minimum velocities, respectively, corresponding to points located on the western and eastern syntaxes of each plate.Rates calculated using rotation parameters of van Hinsbergen et al. (2011a, b).(b) Amount of CO 2 produced by the subduction of the Tethys under Eurasia for the same period (green lines), using plate velocities calculated from van Hinsbergen et al. (2011a), for 15 and 60 % efficiencies.Upper and lower limits of the shaded areas are maximum and minimum gas flux rates computed for each efficiency, respectively.(c) Same as (b) but for Indian only.(d) Same as (b) but including Indian margin subduction at 55, 52.5 and 50 Ma.(e) Same as (d) but for India only.

Figure 5 .
Figure 5. (a, b) Global oceanic benthic δ 18 O (a) and δ 13 C (b) foraminiferal compilation based on data from Cramer et al. (2009).Data were smoothed using a 10-point running average.Temperatures calculated from δ 18 O values assume an ice-free world (after Komar et al., 2013).(c) GEOCLIM modeling results of atmospheric pCO 2 resulting from excess CO 2 release associated with Neo-Tethys closure using plate velocities calculated from van Hinsbergen et al. (2011a), for 15 and 60 % decarbonation efficiencies.Orange curves assume no Indian continental subduction, green lines correspond to an initiation of Indian continental subduction at 52.5 Ma.Individual (gray diamonds) and mean (black line) atmospheric pCO 2 recorded by paleoproxies are also shown (from Beerling and Royer (2011)).Black arrows and associated number refer to pCO 2 values too high to be displayed.(d) Same as (c) but for an initiation of Indian continental subduction at 55 Ma (red curves) and 50 Ma (blue curves).(e) GEOCLIM modeling results of atmospheric pCO 2 resulting from excess CO 2 release as calculated from data of Kent and Muttoni (2013) using efficiencies of 10, 15 and 60 %.
first showed, on the basis of a simple carbon cycle model that the minimal value of additional CO 2 necessary to drive climate warming during the LPEE and the EECO (≥ 1 • C) may have been close to ∼ 10 18 molCO 2 Ma −1 .More recently, Lefebvre et al. (2013) used the GEOCLIM model to calculate that a higher flux of ∼ 3.4 × 10 18 molCO 2 Ma −1 , corresponding to a 50 % increase of global CO 2 degassing rate, was needed to reach a pCO 2 value of 930 ppm consistent with geochemical proxies compiled byBeerling and Royer (2011).Estimations of CO 2 outgassing resulting from Neo-Tethys closure during the Cretaceous and the Paleogene have been previously proposed by . Reagan et al. (2013) presented a review of Late Paleocene to Early Eocene magmatism, characterized by the significant activity of at least three major igneous provinces: the North Atlantic Igneous Province, the Siletzia terrane of the northwestern United States and the Yakutat block in southern Alaska.Added to enhanced activity of Neo-Tethys and eastern Pacific subduction-related volcanism, Reagan et al. (2013) concluded for an overall excess CO 2 production of ∼ 2.3 × 10 18 molCO 2 for the late Paleocene-Early Eocene period.Even though this value may appear important, it encompasses a time duration of several million years.Assuming a 10 Ma duration for significant magmatism (seeReagan  et al. (