Relaxation models applied to paleoclimate dynamics: southern ocean mechanisms controlling glacial-interglacial cycles

CARMEN HERRERO CARMEN HERRERO RELAXATION MODELS APPLIED TO PALEOCLIMATE DYNAMICS: SOUTHERN OCEAN MECHANISMS CONTROLLING GLACIALINTERGLACIAL CYCLES RELAXATION MODELS APPLIED TO PALEOCLIMATE DYNAMICS: SOUTHERN OCEAN MECHANISMS CONTROLLING GLACIALINTERGLACIAL CYCLES Barcelona October 2015 UNIVERSIDAD DE LAS PALMAS DE GRAN CANARIA Facultad de Ciencias del Mar ANEXO I Dª MARÍA ISABEL PADILLA LEÓN, SECRETARIA DE LA FACULTAD DE CIENCIAS DEL MAR, ÓRGANO RESPONSABLE DEL PROGRAMA DE DOCTORADO EN OCEANOGRAFÍA, DE LA UNIVERSIDAD DE LAS PALMAS DE GRAN CANARIA. CERTIFICA Que el Consejo de Doctores del Programa de Doctorado en Oceanografía, en su sesión de fecha 28 de octubre de 2015, tomó el acuerdo de dar el consentimiento para la tramitación, de la tesis doctoral titulada: “Relaxation Models Applied To Paleoclimate Dynamics: Southern Ocean Mechanisms Controlling Glacial-Interglacial Cycles ” presentada por la doctoranda: Dª Carmen Herrero Navarro dirigida por el Doctor D. Antonio García-Olivares Rodríguez Y para que así conste, a efectos de lo previsto en el Artº 6 del Reglamento para la elaboración, tribunal defensa y evaluación de tesis doctorales de la Universidad de Las Palmas de Gran Canaria, firmo el presente en Las Palmas de Gran Canaria, a veintiocho de octubre de dos mil quince. PÁGINA 1 / 1 ID. DOCUMENTO XwwN32p7awxT6rz34GMFFg$$ FIRMADO POR FECHA FIRMA ID. FIRMA 43646105V ISABEL PADILLA LEÓN 29/10/2015 13:05:48 NTI0NTY= Documento firmado digitalmente. Para verificar la validez de la firma copie el ID del documento y acceda a / Digitally signed document. To verify the validity of the signature copy the document ID and access to https://sede.ulpgc.es:8443/VerificadorFirmas/ulpgc/VerificacionAction.action Tesis doctoral presentada por Carmen Herrero Navarro dirigida por el Dr. Antonio García-Olivares Rodríguez para obtener el grado de Doctora por la Universidad de Las Palmas de Gran Canaria, Departamento de Física, DOCTORADO EN OCEANOGRAFÍA Bienio 2012-2014 En Barcelona, a octubre de 2015 RELAXATION MODELS APPLIED TO PALEOCLIMATE DYNAMICS: SOUTHERN OCEAN MECHANISMS CONTROLLING GLACIAL-INTERGLACIAL CYCLES MODELOS DE RELAJACIÓN APLICADOS A LA DINÁMICA PALEOCLIMÁTICA: MECANISMOS DEL OCÉANO AUSTRAL CONTROLANDO LOS CICLOS GLACIALES-INTERGLACIALES El director La doctoranda A la meva família Contents 12 ............. Preface 14 ............. Acknowledgements 22 ............. Glossary 23 ........ INTRODUCTION 30 ....... Chapter 1. RELAXATION MODELS APPLIED TO LATE PLEISTOCENE CLIMATIC OSCILLATIONS 33 .................1.1. INTRODUCTION 34 .................1.2. PP04-DERIVED MODELS 37 ................ 1.2.1. Biological export production model 43 ................ 1.2.2. Two response times for CO2 44 ................ 1.2.3. Two response times for ice volume 46 ................ 1.2.4. Oceanic pulse exponentially dependent on stratification 46 ................ 1.2.5. 3τ model 47 ................ 1.2.6. Local stratification model 48 .................1.3. MODELS PERFORMANCE 50.................1.4. DISCUSSION 56 ........ Chapter 2. ON THE PHYSICAL MECHANISMS BEHIND GLACIAL-INTERGLACIAL DYNAMICS 59 ................. 2.1. INTRODUCTION 60................. 2.2. RELAXATION MODELS 63 ................. 2.3. NON-LINEAR ANALYSIS 63 ................ 2.3.1. Fourier transform 16 Al setembre del 2009, ja amb la carrera pràcticament acabada, vaig prendre la decisió de lluitar per un futur en el món de la recerca. Va ser aleshores quan vaig presentar-me al despatx d’un dels meus professors de la carrera per demanar si podia fer pràctiques amb ell. Aquell professor no era ni més ni menys que l’Antonio García-Olivares, director d’aquesta tesi, tot i que ni ell ni jo no n’érem conscients del llarg camí que ens esperava al davant. No sabria dir si les oportunitats apareixen soles o són fruit de buscar-les, però al cap de molt pocs dies de començar a fer les pràctiques, em van trucar dient que hi havia una vacant per la beca que havia demanat per anar uns mesos a la Universitat de les Illes Balears, i que inicialment m’havien denegat. Vaig estar sis mesos a Esporles, aprenent de grans investigadors com en Damià Gomis, en Biel Jordà, en Francisco Calafat i sobretot, na Marta Marcos, amb qui vaig tenir el privilegi de treballar activament. Vaig conèixer gent meravellosa amb els que encara, tot i que esporàdicament, compartim trucades i anècdotes: na Chus Alonso, en JuanJo Ensenyat, en GianPa Candini, en Miquel Gomila... Va ser una experiència genial i un bon aprenentatge. Però tenia clar a on volia estar i aquell lloc era, clarament, l’ICM. Va ser llavors, al setembre del 2010, quan vaig decidir embarcar-me en aquesta gran aventura. Nada de esto hubiera sido posible sin ti. Gracias Antonio por confiar siempre en mi y ofrecerme la oportunidad de aprender, de crecer, de mejorar y de conocer nuevos aspectos de la vida; gracias por enseñarme y por transmitirme esa gran pasión por la ciencia. Sense cap mena de dubte l’altre gran suport d’aquesta tesi ha estat en Josep Lluís Pelegrí. Gràcies per donar-me l’oportunitat de treballar amb vosaltres, per ajudar-me dia rere dia a què madurés tant com a investigadora com a persona. Gràcies també per confiar sempre en el meu criteri i per oferir-me grans oportunitats. Treballar a l’ICM ha estat un grandíssim plaer. No només per què es tracta d’un dels centres de recerca més potents de tot el Mediterrani, si no perquè la gent que el forma es meravellosa. He compartit despatx amb grans persones com l’Evan Mason, 17 Agraïments / Agradecimientos / Acknowledgements la Patricia De La Fuente, en Miquel Rossell, en Suso Peña, en Sergio Ramírez... tots ells grans investigadors i uns companys excel·lents. També hi ha en Marc Gasser amb el que, tot i no haver coincidit explícitament a un despatx, hi ha una connexió especial; gràcies Marc per totes les converses apassionants que hem tingut i per estar sempre disposat a ajudar o a fer una coca-cola. L’oficialment extingit DOF (Departament que Organitza més Festes, Claret et al. 2011) ha donat pas al DOFT (Departament Organitzador de Festes i Trasllats) però per sort la gent que el forma segueix sent fantàstica. Començant pel «nostre» passadís, amb els dos Jordis (Font i Salat), i l’Emilio García, seguint transversalment amb en Kintxo Salvador, en Pere Fernández, la Maribel Lloret i en Jose Pozo, continuant amb l’Álvaro Víudez, el mejor compañero de despacho que alguien pueda imaginar (mil gracias por aguantarme estos últimos meses), i més endavant amb en Mikhail Emelianov i en Jordi Isern. Una planta més amunt trobem en Jose Antonio Jiménez, en Jordi Solé, en Jaume Piera, la Carine Simon, en Miguel Ángel Rodríguez i tot el grup SMOS amb l’Antonio Turiel, en Quim Ballabrera, la Caro Gabarró, en Marcos Portabella, la Maria Piles, la Vero González, la María Belmonte, l’Estrella Olmedo, la Marta Ramírez, la Nina Houreau, la Marta Umbert... i per descomptat, en Fernando Pérez i en Justino Martínez, sempre amb la solució als nostres problemes. Alguns dels que ja han «volat del niu» són en Pedro Llanillo o la Rocio Rodríguez, grandíssimes persones. Hi ha quatre persones, però, que es mereixen una menció especial, per ordre d’aparició: Paola, Mariona, Maria Rosa i Dorleta. Paola, ha sido un auténtico placer estar a tus espaldas; tienes un gran corazón y una energía aclaparadora, estoy segura que llegaras allí donde te propongas. Gracias, por supuesto, por todos los coffees (aunque yo nunca he tomado café) y por siempre encontrar un hueco para escucharme y darme un buen consejo. La Mariona, ai, Mariona... aquests dies em balla pel cap l’època del B36 i el meu famós «Disfruta, que això només es fa un cop a la vida»; deixam dir-te que se t’ha trobat a faltar per «aquestos lares», però et desitjo que allà on et porti el vent, ja sigui a una banda de l’Atlàntic o a una altra, siguis molt i molt feliç. Maria Rosa, ¿què puc dir-te que no sàpigues?, ha estat un immens plaer poder compartir aquests últims anys al teu costat, parlant de viatges i de noves tec- 18 nologies; ha estat un gran honor conèixer-te una miqueta més a fons. Per cert, ets la meva millor alumna, sense cap mena de dubte! Gràcies per escoltar-me i aconsellar-me amb la teva immensa saviesa. Dorle, ¿quién nos iba a decir que acabaríamos como dos sin papeles en Salvador de Bahía? Ha sido genial el haberte conocido, eres fantástica y estoy segura que llegarás muy lejos. Gracias por siempre tener una palabra de ánimo a mano y por afrontar la vida con una sonrisa. ¡Viva la biodramina! No, espera, permíteme rememorar... ¡Ojo, que volcamos!. Per sort, l’ICM no es redueix a un únic departament, hi ha vida social més enllà dels físics i majoritàriament la trobem a Biologia. Als ICM-food, per tots els dinars (amb o sense tuper), coffees, pels cines, per les 4h d’espera a la Crêperie, per les barbacoes i activitats variades... gràcies! Mireia, Elisa, María D. L. F., Xavi, Ramiro, Rachelle, Isabel, Denisse, Carol, Yaiza, Estela, Sdena, Rosana, Norma, AnaMari, Fran C., Laia, Daffne, Cat y por supuesto, Fran A. Habéis hecho que ir a trabajar sea absolutamente fenomenal. Estimado Sr. Don. Permítame que me dirija a usted en particular para transmitirle mi más sentido agradecimiento. En este tiempo hemos compartido grandes partidas de DS, visto todos los capítulos de Futurama (al menos yo, a usted le ha faltado alguno), comentado grandes películas de la saga Jurassic Park (las buenas y las no tan buenas), ha cambiado usted unas 13 veces de teléfono (todas ellas justificadas), pero lo más importante es que siempre ha estado ahí para animarme cuando más falta ha hecho. Muchos ánimos en este último sprint, tú puedes con esto y más!!! Sin más, me despido. Suya, siempre, una Nobel. Gràcies també a tot el personal d’administració i serveis, amb la Núria Angosto, l’Eva López, en César García, en Jordi Estaña, en José Fortuño, en Jose María Anguita, en Xavi Leal, la Rosa Cabanillas, la Mariví Martínez de Albéniz i especialment, la Conchita Borruel, qui és capaç de resoldre qualsevol problema. Mil gràcies a tots i a totes per ajudar-nos amb la burocràcia i per aconseguir que l’ICM funcioni «a las mil maravillas». A l’ Eli Broglio, gràcies per ajudar-nos amb la divulgació i per organitzar les activitats més divertides i interessants que algú pugui imaginar. A Anita Bonilla de Pelopantón, gracias por ilustrar la ciencia en color y traer un poco de diseño al aburrido mundo científico, y por supuesto, por ser capaz de capturar en 19 Agraïments / Agradecimientos / Acknowledgements papel mis ideas. És més que possible que m’estigui oblidant més d’un i més d’una, així que mil gràcies a tots i a totes per aquests grandíssims 5 anys. Evidentment, no tot es redueix a l’ICM. A todo el grupo PalMA, de la Universidad Complutense de Madrid. Gracias Marisa por contestar al correo de una estudiante sin mucha idea diciendo que le gustaba esto de los modelos. Gracias por acogerme, por ofrecerme la oportunidad de trabajar con vosotros y por enseñarme más de cerca el universo CLIMBER. Gracias Jorge y Alex por todas las discusiones apasionantes y por ayudarme a entender un poquito más esto de los modelos, incluso los de hielo. A Rubén, por no ser solo un compañero sino un gran amigo con el que poder comentar la jugada, ¡incluso la del mountain bike! A Etor, por todas las conversaciones apasionantes ya sean de física o de fotografía; a Laura, por su entusiasmo y energía positiva; a Fidel, a Núria, a Edmundo, a Elena, a Ángela, a Jorge N., a Jorge C., y seguro que me estoy dejando a algun@ por el camino. A todo PalMA, gracias, de verdad, por hacerme un huequito semanal esporádicamente en vuestras vidas y por todas esas palmeras de chocolate que hemos compartido (las mejores del mundo mundial, doy fe) y por las que nos quedan por hacer. Thank you so much to the UCSB for an amazing oportunity. My most sincere gratitude to Lorraine Lisiecki, for welcoming me and offering me the chance to work with her. A special thanks to Vivian Stopple and the administration group, for all the help with the documentation and their warm welcome. Thanks to Rachel Spratt for all the interesting talks, for showing me around and for all the lunches together; to my office-mates, Deborah Khider and Cedric Twardzik, for their help and kindness; to Will Gray, for all the table tennis games (I want my revenge). Also thanks to Frank Kinnaman, for being a great friend in Goleta and for sharing my passion for photography. A huge thanks to the Gans Family, for all the amazing times and for making us feel like home. On the other hand, I’m deeply grateful to all the Intense Family: Steberelli, Dettmers, Herrick, Niermeijer, Bible, Palmer... Thanks so much for being absofrackinglutely rad #intenseforlife. Canviem radicalment d’òrbita. Gràcies als amics que he anat fent en cadascuna de les etapes de la meva vida. A la gent del Roser, per compartir amb mi la infància i la joventut; especialment a la Marta, amb qui he compartit grans moments (com 20 una pizza amb gust de castanya). A les nenes del batxillerat: l’Ester, la Mar, la Cris, l’Anna i la Elena; gràcies per acollir-me quan vaig arribar a una escola desconeguda i fer-me un lloc en les vostres vides. A tota la gent que em va acompanyar al llarg de la vida universitària, especialment als Freaksics: Fran, Laly, Tofu, Tinoco, Guillem, Palau, Pablo, Sebas, Sergi, Alberto, Jordi C., Jordi R., Mario, Alba, Laura, Albert... i segur que em deixo a algú; a tots, gràcies. No m’oblido d’aquells que van compartir amb mi les aules per primer cop: en Daniel (junt amb la Julia i l’Èrika), qui ha passat de ser una unitat individual a una molècula de 3 àtoms, amb qui seguim compartint la vida encara que ens separin 19300 km; sembla mentida que amb els anys que han passat i la de coses que han canviat hi hagi certs aspectes de la vida que segueixen igual, com que hi ha converses que malauradament només tú i jo trobem fascinants. And of course, la Sílvia, qui fa massa anys que és lluny però amb qui sempre ens quedarà una màquina de cafè i un «sense got». Als òptics de la UAB i derivats, sempre disposats a fer un sopar, un vídeo o un club de fans: en Menxi (qui té ulleres de putilla) i la Lina, amics d’aquells que sempre ve de gust veure; en Mompi, un dels millors professors que he tingut mai; en David, sempre aventurer i somrient; la Sònia, amb la seva energia positiva i de qui sóc presidenta en superposició del seu club de fans; i per descomptat, dues persones clau: en Bensi i la Car. Albert, ets un crack; deixa’m dir-te que has estat el meu model a seguir a la tesi i hem viscut mil i una aventures, però crec que em quedo amb INVERNAAAAALIAAAA i la superposició presidentística (amb el moment estel·lar a l’Abacus, explicant que havíem comprat 100 stands de Quàntic Love i no ens els havien enviat). Car, CarOlid, Carpeta, Cartró i derivats similars... quantes coses hem viscut eh! Car, si algú ens hagués dit en aquella pràctica de col·lisions relativistes tot el que viuríem juntes, no ens l’haguéssim cregut. Gràcies per tots els bons moments, totes les idees esbojarrades, les pelis de Spiderman (i Kill Bill) que hem compartit... i per tots els bons moments que encara estan per venir. A tots els amics bikers, especialment a la Vane i en Miki, qui sempre tenen una estona per ajudar en tot el que faci falta (per quan un sushi?); gràcies per estar sempre al nostre costat, sou molt grans! Evidentment no hi ha prou espai per anomenar a tothom (tot i que m’he esforçat en mirar de fer-ho), així que gràcies a totes aquelles persones que m’heu acompanyat al llarg de la meva petita història. 21 Agraïments / Agradecimientos / Acknowledgements Un lugar especial lo reservo a mis padres, Lucía y Emilio, porque sin ellos absolutamente nada de esto hubiera sido posible. Gracias por creer siempre en mi, por enseñarme, por educarme, por darme vuestro apoyo incondicional y también, por malcriarme (y cuánto he disfrutado con ello). Gracias por hacer que me cuestionara las grandes decisiones para estar segura de que era lo que realmente quería (aunque siempre fuera nada más levantarme) y gracias por darme las herramientas necesarias para ser capaz de asumir las consecuencias derivadas. Nada de lo que pueda escribir será suficiente para agradeceros todo lo que habéis hecho por mi, pero permitidme insistir: gracias por estar ahí, siempre. Gracias, también, a la familia que, pese a estar lejos, siempre están al otro lado del teléfono dispuestos a contarme las últimas novedades o dispuestos a compartir los grandes momentos. A Carmina, a los Pepes, a los tíos, a los primos y a los chicos (Victoria, Pablo y Ana), gracias. Y, por supuesto, gracias por todas esas velitas que habéis encendido y que me han iluminado el camino. A la familia Guàrdia - Pascual, per acollir-me a casa seva i fer-me part del seu dia a dia. A l’Antonia i en Josep, a en Dídac i la Mireia, a la Montse i en David, a les iaies Rosita i María, als tiets i als cosins: gràcies. Last but not least, a en Bernat. Tot va començar a les aules de física, ara ja fa uns quants anys, amb un «nosaltres anem a seure allà, vens?» i des d’aleshores la Terra gira a ritme de pedal. Gràcies per compartir la vida amb mi i fer que cada dia sigui més especial que l’anterior. Gràcies per escoltar-me, per ajudar-me, per animar-me, per fer-me riure cada matí i per fer-me lluitar pels meus somnis. Gràcies per haver cregut sempre en mi i per fer-me ser cada dia millor persona. I don’t know how to put this, but you’re kind of a big deal. Barcelona, 21 d’octubre de 2015 AABW Antarctic Bottom Water ACC Antarctic Circumpolar Current AMOC Atlantic Meridional Overturning Circulation AP After Present BP Before Present CDW Circumpolar Deep Water CO2 Carbon dioxide CRP Cross Recurrence Plot DIC Dissolved Inorganic Carbon EMIC Earth System Model of Interme diate Complexity HS1 Heinrich event Glossary ITCZ Intertropical Convergence Zone LCDW Lower Circumpolar Deep Water LGM Last Glacial Maximum MOC Meridional Overturning Circulation NADW North Atlantic Deep Water NH Northern Hemisphere NMOC Northern Meridional Overturning Circulation ppm Parts per million SH Southern Hemisphere SMOC Southern Meridional Overturning Circulation SO Southern Ocean YDS Younger Dryas 23 Let me start with a fundamental question: why should we study climate? Climate affects our daily life in many ways: on the food we eat, on the houses we live in, on our work, on how we travel... Even affects our culture, our spare time or our health. But all of those are only locally important. What we need to see here is the big picture: climate affects the way all living species have adapted to the biosphere, it is key for our survival, but we are drastically influencing it and we cannot predict the consequences of our actions. It is, then, fundamental to broaden our knowledge on the climate system of Earth, the only home we have ever known, to study the past climatic variability on various time scales to obtain clues that will help society face future climate change. As Carl Sagan foreseen on Cosmos, one of his most known books (Sagan [1980]): Our intelligence and our technology have given us the power to affect the climate. How will we use this power? Are we willing to tolerate ignorance and complacency in matters that affect the entire human family? Do we value short-term advantages above the welfare of the Earth? Or will we think on longer time scales, with concern for our children and our grandchildren, to understand and protect the complex life-support systems of our planet? The Earth is a tiny and fragile world. It needs to be cherished. Introduction What is a scientist after all? It is a curious man looking through a keyhole, the keyhole of nature, trying to know what’s going on. — Jacques-Yves Cousteau 24 Paleoclimatology is the science that studies those changes in climate on a really long time scale; as long as the entire history of Earth, 4.54 billion years. It uses data previously preserved on a wide representation of environments (i.e. rocks, sediments, ice sheets, corals...) to reconstruct the past state of the Earth climate and its corresponding variability. To find detailed data as old as Earth is almost an impossible task, but we can certainly use a variety of proxy methods to obtain reliable time-series with a time scale of millions of years. In this work, we have used proxy records with a maximum time domain of 5 million years, and our results are reduced to what we may consider the most recent time, only 800,000 years, or in paleoclimate units, 800 kyr. To understand climate in detail, we should look into all relevant processes of the system, although this might be slightly terrifying, as the number of these processes which must be understood is astonishing. To make it more interesting, the amount of true knowledge that we have is limited and partial. Should we cry on despair? Certainly not. Luckily, paleoclimate, as all sciences, is based on the scientific method, which help us on developing an understanding of reality using terms as mechanism, model and theory. We may describe a theory as a detailed mathematical explanation of phenomena that has sufficient relevance to make predictions from fundamental principles, or in other words, as stated by Truesdell and Toupin [1960], a theory is no more than a mathematical model for nature. Truth is, there is no clear path to obtain a reasonable theory of climate in the near future. Rather, we all hope to develop models of climate, in which general conservation principles and assumptions on local mechanisms of feedback allow us to derive mathematical algorithms and predictions that need to be compared with experimental observations. As long as the predictions and the experiments match, we may consider the model valid to describe a reality, a complex object called climate system. Historically, many inroads have been made to model the structure of long-term climate variability. In the 19th century, Louis Agassiz (Agassiz [1838]) proposed for the first time that climate on Earth could have been much colder in older periods of time, at least in the Northern Hemisphere (NH), but under a catastrophist perspective. Later on, Joseph Adhémar (Adhémar [1842]) formulated the first astronomical RELAXATION MODELS APPLIED TO PALEOCLIMATE DYNAMICS: SOUTHERN OCEAN MECHANISMS CONTROLLING GLACIAL-INTERGLACIAL CYCLES 25 Introduction theory of glacial ages, suggesting that astronomical parameters should modify climate, and explained the position of ice sheets according to the precession of equinoxes, and the subsequent change of location of the perihelion (nearest point to the Sun on Earth orbit) in relation to seasons. Nevertheless, Adhémar formulation was largely criticized, mostly with unfounded arguments, although it had indeed problems. James Croll (Croll [1867]) worked on a more elaborated version of this theory, considering the role of eccentricity as a modulation of the precessional forcing, and following Adhémar ideas, insisted on the role of snow accumulation during winters as lead of glaciations. It was on the 20th century when Milutin Milankovitch (Milankovitch [1941]) formulated a theory that is still valid nowadays. Croll’s main problem was to consider winter as the critical season for the ice sheet evolution, while Milankovitch demonstrated that summer melting is far more important for the ice mass balance, considering summer insolation a key parameter on his theory. Furthermore, Milankovitch considered a third astronomical parameter, the obliquity of the Earth axis (axial tilt), setting the basis of the present astronomical theory of climate. On the other hand, other hypothesis were developed following a completely different point of view, the geochemical theory. Svante Arrhenius (Arrhenius [1896]), inspired by the previous work of Joseph Fourier or John Tyndall, stated the role of CO2 on climate. He correctly computed how a variation of atmospheric CO2 levels affects the surface temperature through the greenhouse effect, considering the atmosphere as a carbon reservoir able to control ice ages. Under this perspective, CO2 and global climate are driving the ice sheet changes, while Milankovitch point of view states that ice sheets controlled by summer northern insolation are the ones driving global climatic changes. Both theories have limitations and problems, and as the reader may suspect, both are essentially valid and not exclusive. Hays et al. [1976] demonstrated that variations of glacial-interglacials states during the Quaternary may be ultimately caused by changes in parameters of the Earth’s orbit (Milankovitch theory), as the astronomical periodicities are found in paleoclimatic records. They also showed that at the obliquity (41 kyr) and precessional bands (23 and 19 kyr) exist a direct climatic response to high latitude summer 33 .................1.1. INTRODUCTION 34 .................1.2. PP04-DERIVED MODELS 37 ................ 1.2.1. Biological export production model 43 ................ 1.2.2. Two response times for CO2 44 ................ 1.2.3. Two response times for ice volume 46 ................ 1.2.4. Oceanic pulse exponentially dependent on stratification 46 ................ 1.2.5. 3τ model 47 ................ 1.2.6. Local stratification model 48 .................1.3. MODELS PERFORMANCE 50.................1.4. DISCUSSION Contents 33 1.1 INTRODUCTION Some simple relaxation models have been proposed to explain the glacial-interglacial cycles. Particulary, Paillard and Parrenin [2004] proposed a model (hereafter PP04) that incorporates very simple parameterizations attempting to represent dense water formation in the SO. Though the model is simple, the results are encouraging because they correctly reproduce the pace and termination times of all the glacial cycles observed. Here, the PP04 model has been generalized and calibrated to the δ18O and CO2 time series available for the last 800 kyr Before Present (BP) (Petit et al. [1999]; Monnin et al. [2001]; Pepin et al. [2001]; Lisiecki and Raymo [2005]; Siegenthaler et al. [2005]; Luthi et al. [2008]). The objectives of this chapter are: Chapter 1 RELAXATION MODELS APPLIED TO LATE PLEISTOCENE CLIMATIC OSCILLATIONS We embarked on our journey to the stars with a question first framed in the childhood of our species and in each generation asked anew with undiminished wonder: What are the stars? Exploration is in our nature. We began as wanderers, and we are wanderers still. We have lingered long enough on the shores of the cosmic ocean. We are ready at last to set sail for the stars. ‒ Carl Sagan, COSMOS RELAXATION MODELS APPLIED TO LATE PLEISTOCENE CLIMATIC OSCILLATIONS 34 The next section describes the structure of the initial PP04 model and its best calibrations to the δ18O and CO2 data for the last 800 kyr BP. A biological carbon export rate is next introduced instead of the original mechanism, and the modified model is calibrated and analyzed. An alternative model with two different response times for the CO2 emission and absorption is also considered. The model is then modified by two different response times for accumulation and ablation of ice, and by adding an exponential oceanic pulse. Next, we compare the patterns found when using different model calibrations, and analyze different mechanisms controlling the oceanic pulse that triggers the deglaciations. Finally, the results and conclusions are summarized. 1.2. PP04 DERIVED MODELS The set of equations of the original PP04 model are the following: (i) to obtain the set of parameters that best fit the late Pleistocene time series in PP04’s original model; (ii) to investigate several generalizations of Paillard’s model by introducing some simple sub-models for the export production, the strength of deep stratification and other mechanisms; and (iii) to analyze how the different sub-models affect the data fit. 1.1 dV = (Vr ‒ V ) , τV dt 1.3 dC = (Cr ‒ C), τC dt 1.2 dA = (V ‒ A) , τA dt Chapter 1 35 where V and A are dimensionless indexes for non-Antarctic ice volume and extent of Antarctic ice sheet respectively,C is a dimensionless index for atmospheric CO2, I65 is the daily insolation at 65° N on 21 June, P(F) is the deep ocean contribution to the reference CO2, which in the original model is P(F) = H(F), where H is the Heaviside function (H = 1 if F < 0; H = 0 otherwise), and F is the salty bottom water formation efficiency parameter. F increases with ice volume V and decreases when continental shelf areas are reduced (through A) and when I60 (daily insolation at 60 S on 21 February) increases. The main variables V, A, and C tend to decay exponentially to a reference state Vr, V, and Cr in characteristic times τV , τA and τC , respectively. Reference volume (Vr) decreases when C increases, and when I65 is high. The CO2 of reference (Cr) increases when I65 is high through parameter α and when the oceanic pulse of CO2 is on (through parameter γ), and decreases when the ice volume increases. The β parameter represents the feedback between ice volumeV and CO2. Daily insolations of the last 800 kyr BP are obtained from the Berger [1978a] and Berger and Loutre [1991] software (Appendix A). We have used negative and positive values when respectively referring to times before present (BP) and after present (AP), where year zero is taken as 1950 AD. Insolations have been normalized with their standard deviations to obtain I65 and I60. No lower bound is imposed on the V, C and A variables, which may take positive or negative values because they represent relative and not absolute variations. Note that all results of the model have been normalized by their standard deviation before plotting. V r = ‒ xC ‒ yI65 + z,1.4 Cr = αI65 ‒ βV + γP ( ‒F ) + δ,1.5 F = aV ‒ bA ‒ cI60 + d, 1.6 RELAXATION MODELS APPLIED TO LATE PLEISTOCENE CLIMATIC OSCILLATIONS 36 In Paillard and Parrenin [2004] no precise tuning of the parameters was performed but a table with central value and range for the parameters was presented. In original PP04, γ parameter had a value of 0.5 due to a mistake. Noticed by Michael Crucifix and Didier Paillard in a later revision, the correct value is 0.7 (Table 1.1). Three different integration methods were used for the differential equations (Mathematica internal algorithms, Matlab internal algorithms1 and Predictor-Corrector explicit method with time step of 50 years) and the results were coincident. Given the strong nonlinearity in the variable C dynamics, at least 16000 steps are needed in the latter algorithm to grant that results are independent of the number of steps. A systematic search of the set of parameters that lead to maximum correlation was implemented around this central case. Correlations of the V and C time series with the (respectively) δ18O and CO2 experimental series available for the last 800 kyr BP were calculated to quantify the maximum observational explained variance, giving us a first approximation of the similarity between simulated and observational time series. The expression for the statistical correlation is 1Mathematica is a trademark of Wolfram Research and Matlab is a trademark of The MathWorks. where cor[s1, s2] is the correlation between the time series s1 and s2, E[x] is the expected value (mean) of x, the dot symbol denotes scalar product, σ[x] is the standard deviation of x, and N is the number of components of arrays s1 and s2. The correlations obtained are shown in Table 1.1. As can be observed in column 2 of the table, the correlations of the central case of PP04 can be improved from 0.59 and 0.63 to 0.67 and 0.71 for CO2 and ice volume, respectively, when the α parameter takes the zero value (model PB). This parameter represents a direct forcing cor[s1, s2] = (s1 ‒ E [s1]) • (s2 ‒ E [s2]) , (N ‒ 1)σ[s1] σ[s2] 1.7 Chapter 1 37 of I65 on atmospheric CO2 and its physical interpretation is not clear, given that the CO2 response to ice volume is already included in the model. Thus, we have eliminated this parameter from the model hereafter. The second and third panels of Figures 1.1. and 1.2. show the ice volume and CO2 time series simulated with these two parameter sets. The main difference between the two cases appears in the CO2 time series, where the oscillations of 41 and 21 kyr in the case with α = 0 present a lower amplitude than in the case with α = 0.15. 1.2.1. Biological export production model Martínez-García et al. [2009] observed that the export productivity in the subantarctic Atlantic during glacial stages grew exponentially. This suggests that the export production in the SO could be driven by changes in the supply of iron by dust, which is known to Vl mechanism, is able to generate fits of similar correlation to the ones obtained with the original model (PP04). With this purpose, the oceanic pulse in expression Cr is replaced by an expression that represents the biological rate of carbon exportation of the SO as where Cr = αI65 ‒ βV + B + δ,1.8 B = ɛ(ekV ‒ 1), V > 0 . 0, V < 0 1.9 RELAXATION MODELS APPLIED TO LATE PLEISTOCENE CLIMATIC OSCILLATIONS 38 PARAMETERS PP04 PB BIO 2τ4τEP 3τLS τV15000 15000 3667 12000 17006 4600 16585 11325 τV2- - - - 3797 1100 3105.5 2325 τC5000 5000 100 800 6796 3040 13505 2793 τC2- - - 4500 17667 7600 -8414 τA12000 10500 -11333 8089 12000 9004 10266 x1.3 1.32 0.85 1.29 0.767 1.525 0.905 0.669 y0.5 0.45 10.5 0.442 0.2 0.489 0.527 z0.8 0.8 0.92 0.85 1.033 0.96 0.946 0.761 α0.15 - - - - - - 0.237 β0.5 0.496 0.76 0.496 0.406 0.476 0.336 0.793 γ0.7 0.506 0.49 0.513 1.642 0.46 2.044 1.955 δ0.4 0.4 0.4 0.434 0.407 0.445 0.228 0.146 a0.3 0.3 -0.3 0.395 0.3 0.54 - b0.7 0.71 -0.71 0.8 0.71 1.205 0.936 c0.01 0.005 -0.005 ---0.533 d0.27 0.27 -0.27 0.27 0.255 0.483 0.069 ɛ- - 0.8 - - - - - k- - 0.47 - - - - - Vm- - 0.5 - - - - - λ- - - - - 37 - - RV0.63 0.71 0.58 0.82 0.89 0.85 0.88 0.87 RC0.59 0.67 0.45 0.75 0.79 0.77 0.79 0.76 Table 1.1 Central values of the models studied. Columns 2 to 9 correspond to: (2) original Paillard and Parrenin [2004] model with γ = 0.7 (PP04); (3) best fit of previous model (PB); (4) model with biological export production (BIO); (5) model with two relaxation times for C (2τ); (6) model with two relaxation times for C and two relaxation times for V (4τ); (7) model with oceanic pulse exponentially dependent on stratification (EP); (8) model with one relaxation time for C and two relaxation times for V (3τ) and (9) model with local stratification parameters (LS). RV and RC represent the correlation between proxy and modeled data for global ice volume, V, and atmospheric CO2 concentration, C, respectively. Chapter 1 39 MODELS DESCRIPTION NUMBER OF PARAMETERS PP04 Original Paillard and Parrenin [2004] model with γ = 0.7 14 PB Best fit of Paillard and Parrenin [2004] model 13 BIO Biological rate of exportation included at the CO2 reference state 11 2τ Two different response times for the accumulation and absorption of CO2 τC = τC1 for the periods where C < Cr and τC = τC2 when C > Cr 14 4τ Two different response times for the accumulation and absorption of CO2 and two different response times for accumulation and ablation of ice: τV = τV1 for periods when V < Vr and τV = τV2 when V > Vr 14 EP Oceanic pulse exponentially dependent on stratification 14 3τModel with one relaxation time for C and two relaxation times for V13 LS Model with local stratification parameters 15 Description of the different models. Table 1.2 RELAXATION MODELS APPLIED TO LATE PLEISTOCENE CLIMATIC OSCILLATIONS 40 Figure 1.1 Best fits obtained for the ice volume time series. From top to bottom: (1) original Paillard and Parrenin [2004] model with γ = 0.7 (PP04); (2) best fit of previous model (PB); (3) model with biological export production (BIO); (4) model with two relaxation times for C (2τ); (5) model with two relaxation times for C and two relaxation times for V (4τ); (6) model with oceanic pulse exponentially dependent on stratification (EP); (7) model with one relaxation time for C and two relaxation times for V (3τ ) and (8) model with local stratification parameters (LS). Proxy records of δ18O time series (Lisiecki and Raymo [2005]) are superimposed in every panel. Blue bands represent interglacial periods considering benthic δ18O below 3.8 per mil. PP04 I -2.5 0 2.5 PB -2.5 0 2.5 ExP -2.5 0 2.5 2τ -2.5 0 2.5 4τ -2.5 0 2.5 EP -2.5 0 2.5 3τ -2.5 0 2.5 LS Time (kyr) -800 -700 -600 -500 -400 -300 -200 -100 0 -2.5 0 2.5 II IIIIV VVI VIIVIIIIX Chapter 1 41 Figure 1.2 Best fits obtained for the CO2 time series. From top to bottom: (1) original Paillard and Parrenin [2004] model with γ = 0.7 (PP04); (2) best fit of previous model (PB); (3) model with biological export production (BIO); (4) model with two relaxation times for C (2τ); (5) model with two relaxation times for C and two relaxation times for V (4τ); (6) model with oceanic pulse exponentially dependent on stratification (EP); (7) model with one relaxation time for C and two relaxation times for V (3τ ) and (8) model with local stratification parameters (LS). Records of experimental CO2 time series (Petit et al. [1999]; Monnin et al. [2001]; Pepin et al. [2001]; Siegenthaler et al. [2005]; Luthi et al. [2008]) are superimposed in every panel. Blue bands represent interglacial periods considering experimental CO2 concentrations above 250 ppm. PP04 -2 0 3 PB -2 0 3 ExP -2 0 3 2τ -2 0 3 4τ -2 0 3 EP -2 0 3 3τ -2 0 3 LS Time (kyr) -800 -700 -600 -500 -400 -300 -200 -100 0 -2 0 3 III IIIIV VVI VIIVIIIIX RELAXATION MODELS APPLIED TO LATE PLEISTOCENE CLIMATIC OSCILLATIONS 48 obtained with a = 0. For this reason, we used the final expression (Eq. 1.18) that uses one parameter less. Stratification is thus considered dependent only on regional variables (the extent of Antarctic ice sheet A and SO temperature, related to C). In parallel with it, it was assumed that SO temperature, as represented by variable C, led also the evolution of ice sheet A (Eq. 1.15). This change has no effect on the correlation obtained for V but it improves slightly the correlation obtained for C. The final equations are: 1.3. MODELS PERFORMANCE The analysis of Table 1.1 gives some hints on the conditions that lead to the best data fits. As can be observed in Figure 1.1, all models have problems to accurately simulate the minimum of V at –200 kyr (cycle II), the timing of the V minimum at –490 kyr (Termination VI) and the maximum ice volume at –750 kyr (except, to some extent, model BIO). dA = (‒C ‒ A) , τA dt 1.15 Vr = ‒ xC ‒ yI65 + z,1.16 Cr = αI65 ‒ βV + γH ( ‒F ) + δ,1.17 F = ‒ bA ‒cC + d. 1.18 1.13 dV = (Vr ‒ V ) , τV dt 1.3 dC = (Cr ‒ C), τC dt Chapter 1 49 It may be observed that all models have problems simulating the glacial cycle between –600 and –500 kyr, especially the timing of its termination. Inaccurate fitting of the timing of terminations probably derives from limitations of both model and δ18O records. A common way of estimating time in paleoclimate records is to stretch, squeeze and shift a record’s chronology in order to align it with a template indicative of changes in the Earth’s orbital and rotational configuration, a process generally referred to as orbital tuning (Huybers [2011]). For this reason, it is not clear that the time scale of the data is better or worse than the time scale of the model (within, let’s say, a quarter period of the fastest forcing, the precession, i.e. ±6 kyr). This is particularly true for the timing of terminations, which apparently are not linearly correlated with the astronomical forcing intensity (see e.g. Parrenin and Paillard [2003]). For this reason, a “better fit” of the timing of terminations of the δ18O data may not be an improvement, because the data are not necessarily very accurate in their estimation of the precise dates when terminations take place. Time series for the fit error were obtained by subtracting the simulated and experimental time series resulting from dividing every time series by its standard deviation. The correlations obtained with the EP model do not improve upon the ones obtained with the 4τ model (Table 1.1) so, we can conclude that the results obtained do not justify the introduction of an additional parameter. However, the extra contribution to glacial CO2 produced in this model makes it possible to fit the experimental CO2 series with a lower value of the y parameter (sensitivity to I65 forcing), which produces a V time series with the high-frequency oscillations more damped. The “texture” of this series, as can be observed in Figures 1.1 and 1.2, resembles that observed in the experimental data. The performance of 4τ and 3τ models in the simulation of V is quite similar, except for Termination V and the glacial cycle II, where 4τ performs slightly better. For these reasons, 4τ correlation for V is slightly better than 3τ (0.89 instead of 0.88). The use of a different relaxation time for emission and absorption of CO2 allows 4τ to match CO2 in glacial cycle V better than 3τ. However, the extreme simplicity of the carbon model shows up in both models though the lack RELAXATION MODELS APPLIED TO LATE PLEISTOCENE CLIMATIC OSCILLATIONS 50 of high-frequency oscillations in the evolution of C. The LS model has a content of high-frequency variance larger than that of 3τ because it permits a direct forcing of I65 to C but its correlation is not better because this variance is frequently out of phase with the observational one. 1.4. DISCUSSION Starting from the Paillard and Parrenin [2004] model, several box models incorporating simple parameterizations of the oceanic CO2 pumping and response times of carbon and ice volume have been developed. The models’ parameters were calibrated to provide the best fit to the δ18O and CO2 experimental time series available for the last 800 kyr BP. The PP04 model is insensitive to the α parameter (direct forcing between I65 insolation and CO2) and the Sun’s effect seems to emerge only through the filter of global ice volume. The fit of this model to observational data may be improved if different response times are assumed both for absorption/ emission of CO2 and for ablation/accumulation of ice. Correlations between simulated and experimental time series increase from 0.59 and 0.63 to 0.79 and 0.89 for CO2 and V, respectively. Several modifications of the PP04 model that lead to the right timing of the last nine terminations were tested. In particular, the qualitative behavior of the last eight glacial-interglacial cycles may be roughly reproduced with an export production model with export dependent on ice volume (BIO). However, in this model, the best fit corresponded to a dependence between CO2 export and V that was not exponential as proposed by Martínez-García et al. [2009], but a square function. In our formulation, biological export alone was not able to simulate the experimental data as accurately as the other models, even though some kind of biological export of CO2 is very probably acting in synergy with physical processes in the glacial-interglacial dynamics. One limitation of the methodology used in this study, common in paleoclimate models, is that there is no distinction between calibration data and validation data. For Chapter 1 51 this reason, it may be useful to know whether parameter optimizations obtained on half of the data interval lead to nearly the same best estimates. Our models seem to fit better the sawtooth oscillations that climate has shown in the last –400 kyr than the more irregular oscillations occurring between –800 kyr and –400 kyr. In order to include both kinds of dynamics, we used the interval –600 to –200 kyr as calibration data in the LS model. This interval led to a best fit that is coincident with the one shown in the last column of Table 1.1, and the correlations obtained were 0.872 and 0.763, which are slightly better than those obtained for the complete interval (0.866 and 0.757). We have rounded the latter results to two decimal digits in the last column of Table 1.1. When a second and more exhaustive searching is implemented around the values obtained by the genetic algorithm, it is possible to find a different parameter combination with even better correlations (0.884 and 0.774). However, the case is then so tightly fitted to the calibration interval that it fails to reproduce the pace of the whole interval (and correlations drop to 0.43 and 0.22). On the other hand, the experimental data used (δ18O: Lisiecki and Raymo [2005]; CO2: Monnin et al. [2001]; Petit et al. [1999]; Siegenthaler et al. [2005]; Luthi et al. [2008]) contain much high-frequency variance, representing (among other things) weather noise, measurement error, and annual-to-centennial climate variability. PP04-derived models do not contain mechanisms able to produce significant variance in these time scales. For this reason, in one of our runs we smoothed the proxy data with a moving average of 2 kyr, and then calibrated the LS model with these time series and the complete interval (–800 kyr, 0 kyr). The correlations obtained were 0.869 and 0.761, respectively, whereas the result with unfiltered data was 0.866 and 0.756, respectively. In conclusion, the smoothing of the experimental time series provides an improvement of less than 1% in the fits of the two variables. Four models (4τ, EP, 3τ and LS) show a good performance in the simulation of the ice volume observed in the last eight glacial-interglacial cycles and even in many details of that signal belonging to the scale from 23 to 41 kyr. The good fits obtained suggest that density of deep water and its rate of formation may be important factors controlling the oceanic pulse that triggers the deglaciations. Schmittner [2007] has confirmed, with a global climate model, the sensitivity of atmospheric CO2 to processes that affect stratification in SO waters. In our analysis, oceanic pulses that RELAXATION MODELS APPLIED TO LATE PLEISTOCENE CLIMATIC OSCILLATIONS 52 strongly depend, in a non-linear way, on the deep ocean stratification are necessary to trigger deglaciations. The pulses obtained in the best models are always very close to the Heaviside function proposed by the original work of PP04. Oceanic CO2 pulses with a duration of between 10 and 20 kyr were found at the beginning of the nine last deglaciations according to the best models analyzed. In addition, different response times for emission and absorption of CO2 and for accumulation and ablation of ice are necessary to obtain the best fit to the available data. Accumulation and ablation of ice are different physical mechanisms that could have different characteristic times. However, the two response times for emission and absorption of CO2 are more difficult to explain, as we have said before. To decide whether these two response times correspond to real physical mechanisms, a deeper understanding of the long-term transfer rates in the carbon cycle should be achieved. Björkström [1979] identified the two slowest response times for the emission of CO2 to the atmosphere to be related to the remineralization of organic carbon of dead material in soils, and the long travel time of the CO2 from the deep ocean compartments to the surface. These authors used 1000 year as an order of magnitude for both parameters. Brovkin et al. [2002] used the range 400 - 1000 years to model the slow soil response with the CLIMBER-2 paleoclimatic model. The characteristic time of CO2 absorption in the long-term scales is related to the CO2 advection from the surface into the deep ocean according to Björkström [1979] and it should be, again, of the order of 1000 yr. Our result τC2 > τC could be meaningful if the characteristic time of deep water circulation was longer than the long-term emission time of soils. Montenegro et al. [2007] found that 25% of an instantaneous CO2 release to the atmosphere remains there after 5000 years. Archer et al. [2009b] reviewed the literature on the carbon cycle, which agrees that 20% to 35% of an instantaneous CO2 release remains in the atmosphere in 200 to 2000 years, after equilibration with oceans. Due to oceanic acidification, subsequent dissolution of CaCO3 increases the uptake capacity of oceans in the scale of 3 to 7 kyr. These results are apparently coherent with the order of magnitude of the τC2 parameter that we have found. However, more investigation is needed to have a precise understanding of the long-term carbon cycle. Chapter 1 53 A different problem is deciding what confidence we can give to those models with higher correlations, given that the record used for the past ice volume probably suffers different biases. Indeed, the record is a stack of δ18O data from benthic foraminifera, and the ratio O18/O16 is known to depend on both the isotopic composition and the temperature of the water where the foraminifera develops. Siddall et al. [2010] found that ice volume becomes increasingly sensitive to temperature change at low temperatures. Waelbroeck et al. [2002] found that the relationship between δ18O and ice volume is not linear, since δ18O decreased faster than the increase in ice volume at the beginning of the last glaciations, and then progressively more slowly until the ice sheets reached their maximum size. This may produce uncertainties greater than 20% for the ice volume estimations in some periods. Therefore, the models with the largest correlations between predicted ice volume and δ18O are not necessarily better than other models with slightly lower correlation because of the uncertainty in the observational data. This is a complex and important question requiring an additional mathematical analysis that we leave for a future work. Fortunately, there is no such problem for CO2, which is directly measured in ice cores. It may be argued that it is not surprising to get good fits (Table 1.1) using models with many parameters. This is particularly true when experimental data are fitted to generic mathematical functions that are dependent on an arbitrary number of parameters. However, it is not so easy to get the same good fit with mathematical expressions that imitate geophysical mechanisms such as feedbacks and pumping rates. An additional outcome from relaxation models of this kind is that they point to specific physical mechanisms that are potentially drivers of the observed changes. Thus, when good fits are obtained with such models, the specific physical mechanisms modeled should be investigated in more detail to confirm them (and not an alternative mechanism producing a similar behavior) as inducers of the observed dynamics. To sum up, we have obtained eight different fits with ice volume correlations between 0.58 and 0.89 which show a good quantitative and qualitative agreement with the empirical time series, especially 4τ, EP, 3τ and LS; although the warm event of –500 kyr is not properly reproduced by any model. The 4τ model improves the PP04 correlations (0.59 and 0.63) to 0.79 and 0.89 using the same number of parameters RELAXATION MODELS APPLIED TO LATE PLEISTOCENE CLIMATIC OSCILLATIONS 54 (namely, 14). The EP model uses 15 parameters but the additional parameter, which was aimed at improving the form of the oceanic pulse function, did not improve the correlations and can be considered useless. The 3τ model obtains almost the same correlations as 4τ using only 13 parameters. The LS model (15 parameters) does not improve the correlations of 4τ (which are slightly higher) but incorporates a functional form for stratification F that seems more consistent with both the mechanisms suggested by Paillard and Parrenin [2004]. If our aim were to select the model with the best explained variance per parameter, the choice would be 3τ . However, a second objective of this work is to determine whether models with an explained variance similar to the one obtained by 3τ contain parameterizations that can be related as realistically as possible to observed mechanisms. LS can offer a valuable insight into the real meaning of the good performance of PP04-derived models. In the next chapter, the dynamics and mechanisms incorporated in 3τ and LS models will be analyzed and compared. Chapter 2 ON THE PHYSICAL MECHANISMS BEHIND GLACIAL-INTERGLACIAL DYNAMICS Chapter based on the following published articles: García-Olivares and Herrero [2013], Simulation of glacial-interglacial cycles by simple relaxation models: consistency with observational results. Climate Dynamics 41:1307-1331 Herrero and García-Olivares [2014], Non-linear analysis methods applied to observational and simulated climatic time series. Proceedings International work-conference On Time Series. Volume 2, p10681080. I.S.B.N: 978-84-15814-97-4 Herrero and García-Olivares [2015], Non-linear analysis methods applied to observational and simulated climatic time series. Boletín Geológico Minero, accepted for an ITISE special volume. 64 ON THE PHYSICAL MECHANISMS BEHIND GLACIAL-INTERGLACIAL DYNAMICS As shown in Figure 2.3 top and Figure 2.4, the 100-kyr band is the dominant period in both the δ18O and V time series of the 3τ and LS models, followed by the 41-kyr band, and power in both frequencies is similarly distributed in observational and simulated series. The δ18O time series contains periodicities in the range 1-10 kyr; in contrast, our models do not generate any significant periodicity under 10 kyr. The 23-kyr period is weakly present in both experimental and 3τ time series. The LS and 3τ models show the same patterns in their wavelet diagrams, which are almost indistinguishable from each other. For CO2 (Figure 2.3 bottom and Figure 2.5), the power distribution in the simulated and observational time series are similar for the 100-kyr band. Power in the 41-kyr band is somewhat more homogeneously distributed in the observational time series than in the simulated time series. The 23-kyr contribution is weaker in both series and is distributed in a different way in the modeled and observational time series. The difference is greater in the first three cycles, which are relatively difficult to model. The complex cross-wavelet transform of two time series (Figures 2.6 and 2.7) can be interpreted as the shared power in a given periodicity band (absolute value) and the phase difference between the two series in time frequency space (Grinsted et al. [2004]). Lighter purple color indicates greater shared power in that time and periodicity band, and the arrow angle indicates the phase between the observational and simulated series. As can be observed in Figure 2.6, observational and simulated ice volume share a great common power in the 100-kyr and 41-kyr band. In the 100-kyr band, the two periodicities are in phase, but in the 41-kyr band the modeled periodicity tends to lead the observational one at some moments, with a phase between 0 and π/4. In the 41 and 23-kyr band, the two periodicities tend to be in phase when the common power is high and out of phase when the shared power is low. The two models studied (3τ and LS) show very close patterns in their cross-wavelet transforms for V. Regarding CO2 (Figure 2.7), the largest shared power is in the 100-kyr band, especially after ‒500 kyr, when the fit between simulated and observational cycles is somewhat better (see Figure 1.2 and 2.1). In the 41-kyr band the shared power is lower and has relative maxima around ‒600 kyr, between ‒400 and ‒200 kyr and between ‒150 and -50 kyr. These intervals of good coincidence can be also Chapter 2 65 0 20 40 60 80 100 120 0 500 1000 1500 2000 Spectral Power V 0 20 40 60 80 100 120 0 500 1000 1500 2000 Spectral Power C Period (Kyr) At top, spectral power of normalized ice volume as predicted with 3τ model (blue), LS model (red) and proxy δ18O records from Lisiecki and Raymo [2005] (grey line); at bottom, normalized CO2 as predicted with the 3τ (blue line) and LS (red line) models and proxy data (Petit et al. [1999]; Indermuhle et al. [2000]; Monnin et al. [2001]; Siegenthaler et al. [2005]; Luthi et al. [2008]) (grey line). Figure 2.2 observed in Figures 1.1, 1.2 and 2.1 for the LS model. Low shared power is observed in the 23-kyr band, with an out-of-phase pattern in some intervals. Wavelet coherence of two series can be interpreted as a localized correlation coefficient in time frequency space. It is useful to find locally phase-locked behavior, that is, 66 ON THE PHYSICAL MECHANISMS BEHIND GLACIAL-INTERGLACIAL DYNAMICS 1/64 1/32 1/16 1/8 1/4 1/2 1 2 4 8 16 32 64 Time (Kyr) Frequency (Kyr) −700 −600 −500 −400 −300 −200 −100 0 0.20 0.39 0.78 1.56 3.13 6.25 12.5 25 50 100 200 Wavelet diagram of proxy δ18O records from Lisiecki and Raymo [2005] (top) and proxy CO2 data from Petit et al. [1999]; Indermuhle et al. [2000]; Monnin et al. [2001]; Siegenthaler et al. [2005]; Luthi et al. [2008] (bottom). Figure 2.3 1/64 1/32 1/16 1/8 1/4 1/2 1 2 4 8 16 32 64 Time (Kyr) Frequency (Kyr) −700 −600 −500 −400 −300 −200 −100 0 0.20 0.39 0.78 1.56 3.13 6.25 12.5 25 50 100 200 Chapter 2 67 1/64 1/32 1/16 1/8 1/4 1/2 1 2 4 8 16 32 64 Time (Kyr) Frequency (Kyr) −700 −600 −500 −400 −300 −200 −100 0 0.20 0.39 0.78 1.56 3.13 6.25 12.5 25 50 100 200 Wavelet diagram of the simulated V time series with 3τ model (top) and LS model (bottom). Figure 2.4 1/64 1/32 1/16 1/8 1/4 1/2 1 2 4 8 16 32 64 Time (Kyr) Frequency (Kyr) −700 −600 −500 −400 −300 −200 −100 0 0.20 0.39 0.78 1.56 3.13 6.25 12.5 25 50 100 200 68 ON THE PHYSICAL MECHANISMS BEHIND GLACIAL-INTERGLACIAL DYNAMICS 1/64 1/32 1/16 1/8 1/4 1/2 1 2 4 8 16 32 64 Time (Kyr) Frequency (Kyr) −700 −600 −500 −400 −300 −200 −100 0 0.20 0.39 0.78 1.56 3.13 6.25 12.5 25 50 100 200 Wavelet diagram of the simulated C time series with 3τ model (top) and LS model (bottom). Figure 2.5 1/64 1/32 1/16 1/8 1/4 1/2 1 2 4 8 16 32 64 Time (Kyr) Frequency (Kyr) −700 −600 −500 −400 −300 −200 −100 0 0.20 0.39 0.78 1.56 3.13 6.25 12.5 25 50 100 200 Chapter 2 69 1/64 1/32 1/16 1/8 1/4 1/2 1 2 4 8 16 32 64 Time (Kyr) Frequency (Kyr) −700 −600 −500 −400 −300 −200 −100 0 0.20 0.39 0.78 1.56 3.13 6.25 12.5 25 50 100 200 Cross-wavelet transform between proxy δ18O records from Lisiecki and Raymo [2005] and simulated V time series with 3τ model (top) and LS model (bottom) . Figure 2.6 1/64 1/32 1/16 1/8 1/4 1/2 1 2 4 8 16 32 64 Time (Kyr) Frequency (Kyr) −700 −600 −500 −400 −300 −200 −100 0 0.20 0.39 0.78 1.56 3.13 6.25 12.5 25 50 100 200 70 ON THE PHYSICAL MECHANISMS BEHIND GLACIAL-INTERGLACIAL DYNAMICS 1/64 1/32 1/16 1/8 1/4 1/2 1 2 4 8 16 32 64 Time (Kyr) Frequency (Kyr) −700 −600 −500 −400 −300 −200 −100 0 0.20 0.39 0.78 1.56 3.13 6.25 12.5 25 50 100 200 Cross-wavelet transform between proxy CO2 data from Petit et al. [1999]; Indermuhle et al. [2000]; Monnin et al. [2001]; Siegenthaler et al. [2005]; Luthi et al. [2008] and simulated C time series with 3τ model (top) and LS model (bottom). Figure 2.7 1/64 1/32 1/16 1/8 1/4 1/2 1 2 4 8 16 32 64 Time (Kyr) Frequency (Kyr) −700 −600 −500 −400 −300 −200 −100 0 0.20 0.39 0.78 1.56 3.13 6.25 12.5 25 50 100 200 Chapter 2 71 0 0.1 0.2 0.3 0.4 0.5 0.6 0.7 0.8 0.9 1 Time (Kyr) Frequency (Kyr) −700 −600 −500 −400 −300 −200 −100 0 0.20 0.39 0.78 1.56 3.13 6.25 12.5 25 50 100 200 Wavelet Coherence between proxy δ18O records from Lisiecki and Raymo [2005] and simulated V time series with 3τ model (top) and LS model (bottom). Figure 2.8 0 0.1 0.2 0.3 0.4 0.5 0.6 0.7 0.8 0.9 1 Time (Kyr) Frequency (Kyr) −700 −600 −500 −400 −300 −200 −100 0 0.20 0.39 0.78 1.56 3.13 6.25 12.5 25 50 100 200 72 ON THE PHYSICAL MECHANISMS BEHIND GLACIAL-INTERGLACIAL DYNAMICS 0 0.1 0.2 0.3 0.4 0.5 0.6 0.7 0.8 0.9 1 Time (Kyr) Frequency (Kyr) −700 −600 −500 −400 −300 −200 −100 0 0.20 0.39 0.78 1.56 3.13 6.25 12.5 25 50 100 200 Wavelet Coherence between proxy CO2 data from Petit et al. [1999]; Indermuhle et al. [2000]; Monnin et al. [2001]; Siegenthaler et al. [2005]; Luthi et al. [2008] and simulated C time series with 3τ model (top) and LS model (bottom). Figure 2.9 0 0.1 0.2 0.3 0.4 0.5 0.6 0.7 0.8 0.9 1 Time (Kyr) Frequency (Kyr) −700 −600 −500 −400 −300 −200 −100 0 0.20 0.39 0.78 1.56 3.13 6.25 12.5 25 50 100 200 Chapter 2 73 moments in which both series oscillate with the same frequency and a given phase difference. Figure 2.8 shows the wavelet coherence of the δ18O and the simulated V series for 3τ and LS and Figure 2.9 shows the wavelet coherence of the observational CO2 and the simulated C series for the same models. It can be observed that the correlations of V for the two models are always in phase, are larger than 0.9 in the 100-kyr band for all times and are especially high in the four last cycles (‒400 kyr to 0). In the 41-kyr band the simulated and observational series oscillate together with a correlation that exceeds 0.9 around ‒600 kyr, between ‒400 and ‒350 kyr, and is always greater than 0.7 between ‒130 and ‒90 kyr. The coherence is not significant between ‒750 and -630 kyr and between ‒540 and -380 kyr, corresponding to intervals in which the models are worse at matching the observational series (see Figures 1.1 and 2.1). Regarding CO2 (Figure 2.9), in the two models the coherence at the 100-kyr band is in phase and high after ‒500 kyr and weaker and out of phase before ‒500 kyr coinciding with the first glacial cycles (see also Figure 2.2), which are more difficult to simulate. In the 41-kyr band the coherence is high in roughly the same intervals as in the V series, and in the 23-kyr band it is practically inexistent. This can be attributed to the combined effects of high-frequency damping produced by parameter τC (and/or τC2 in LS) and poor representation of the carbon dynamics by PP04-derived models. 2.3.3. Phase space portraits Phase space portraits and embedding attractor techniques can also be useful for quantifying the performance of a simulated time series in matching the dynamical properties of an observational time series. Time series of 16000 linearly interpolated, equally spaced data are used to reconstruct the attractors. To avoid autocorrelated effects we use the method of mutual information (Marwan et al. [2007]), a well-established measure to detect nonlinear dependencies within a time series, which estimates an optimal time lag of 1433 for V and 1371 for C. The first step 80 ON THE PHYSICAL MECHANISMS BEHIND GLACIAL-INTERGLACIAL DYNAMICS Embedded attractor of the simulated C time series with LS model. Figure 2.15 Chapter 2 81 2.3.4. Cross recurrence plots We then obtained the cross recurrence plot (CRP) of observational and simulated trajectories which, despite its name, does not represent recurrences but rather the conjunctures of states of the two systems. The CRP reveals all the times when the phase space trajectory of the first system visits roughly the same region in the phase space where the phase space trajectory of the second system is (Marwan et al. [2007]). Figures 2.16 to 2.19 show the CRP obtained with the above parameters for the simulated series of V and C obtained by the models. The Matlab tool created by Marwan et al. [2007] was used. Last glacial cycle does not show up because this period data is smaller than D times l, where D is the dimension (D = 3) and l is the lag (l = 1433 for V and l = 1371 for C). We can see that the persistent intervals of good coincidence between the V and δ18O series (Figure 2.1) coincide with intense black patches in Figures 2.16 and 2.17, especially spots following Termination IX and Termination VII, 5000 < N < 6000, 7000 < N < 9000, 9200 < N < 9500, 10500 < N < 11000 and 11500 < N < 12500 (Terminations are shown in Figures 1.1 and 1.2). Many patterns occurring in the main diagonal appear deformed in lines and columns, indicating that both series are cyclical and the way in which V(t) coincides with δ18O(t) at time t shows similarity to the way in which V(t) coincides with δ18O(t‒T), where T is the 100-kyr period. Roughly 85% of the time the predicted and observational V trajectories are neighbors at a distance of 0.46 times the standard deviation of the phase space. The C trajectories are much more poorly simulated, with the predicted and observational trajectories coinciding less than 50% of the time (Figures 2.18 and 2.19). However, the dynamics in the intervals when the CO2 is maximum is normally well simulated, which seems to be sufficient to allow the V dynamics to obtain a good similarity to the observational dynamics. 2.4. ALTERNATIVE PROXY FOR ICE VOLUME As pointed out in Chapter 1, the Lisiecki and Raymo [2005] record is a stack of δ18O data from benthic foraminifera, and the O18/O16 ratio is known to depend on 82 ON THE PHYSICAL MECHANISMS BEHIND GLACIAL-INTERGLACIAL DYNAMICS Cross-recurrence plots obtained for smoothed proxy δ18O records from Lisiecki and Raymo [2005] and simulated V time series with 3τ model. Figure 2.16 Chapter 2 83 Cross-recurrence plots obtained for smoothed proxy δ18O records from Lisiecki and Raymo [2005] and simulated V time series with LS model. Figure 2.17 84 ON THE PHYSICAL MECHANISMS BEHIND GLACIAL-INTERGLACIAL DYNAMICS Cross-recurrence plots obtained for smoothed proxy CO2 data from Petit et al. [1999]; Indermuhle et al. [2000]; Monnin et al. [2001]; Siegenthaler et al. [2005]; Luthi et al. [2008] and simulated C time series with 3τ model. Figure 2.18 Chapter 2 85 Cross-recurrence plots obtained for smoothed proxy CO2 data from Petit et al. [1999]; Indermuhle et al. [2000]; Monnin et al. [2001]; Siegenthaler et al. [2005]; Luthi et al. [2008] and simulated time series with LS model. Figure 2.19 86 ON THE PHYSICAL MECHANISMS BEHIND GLACIAL-INTERGLACIAL DYNAMICS both the isotopic composition and the temperature of the water where foraminifera developed. Waelbroeck et al. [2002] found that the relationship between δ18O and ice volume is not linear, since δ18O decreases faster than the increase in ice volume at the beginning of the last glaciations, and then progressively more slowly until the ice sheets reached their maximum size. This may produce uncertainties of 20% for the ice volume estimations at the beginning of glaciations. Consequently, there is some risk of over-interpretation when model parameters are adjusted to fit that isotopic record tightly. If that error was assumed to take place during 20 kyr at the beginning of each glacial period, it would cause a 7% loss of explained variance or a 4% loss of correlation. This is slightly above the correlation differences in the models that we have studied, so their performance must be considered equivalent within the observational uncertainty. Bintanja et al. [2005] tried to eliminate the temperature effects from the δ18O series by using a model for air and deep ocean temperatures. The resulting sea-level time series ca be used as an alternative to δ18O. When we calculate the correlations between this time series and the V predicted by our models, we obtain the surprising result that they increase from 0.88 to 0.90 for 3τ and from 0.87 to 0.89 for LS. This shows that the calibration of our selected models to an isotopic series with experimental biases has produced a robust result in spite of the observational uncertainty, although a non-negligible component is that Bintanja’s time series is slightly smoother than δ18O. In Figure 2.20 the V predicted by the 3τ and LS models is compared with sea level time series (Bintanja et al. [2005]) and δ18O records (Lisiecki and Raymo [2005]). In the bottom panel, the power spectrum of all series is shown. As can be observed, the behavior of the 3τ and LS models compared with Bintanja et al. [2005] sea level is good, reinforcing the robustness of the fits implemented. Although the correlation improves, the results do not offer any new dynamical insights. 2.5. DYNAMICS OF THE MODELS A common feature of PP04-derived models is their lack of sensitivity to the I60 forc- Chapter 2 87 ing. This result could be interpreted in three alternative ways: (i) Our function F (V, A) is an approximate representation of stratification but I60 is not a good proxy for Southern Ocean temperature. (ii) Stratification and/or CO2 release in the SO is not controlled by regional temperatures. (iii) Release of CO2 in the SO is not controlled by stratification and the function F (V, A) in our models represents a different process. A possible interpretation belonging to group (iii) is that A might represent sea ice instead of the ice-sheet area and that F (V, A) could be a proxy for a rate of emission of CO2 that is controlled directly by sea ice area, A, and SO deep temperature. SO deep temperature could be controlled by NH temperature (through V) more strongly than by I60. One possible mechanism able to produce this effect was proposed by Gildor and Tziperman [2001]: stratification in the SO is composed of cold, fresh water above warm, salty water. Glacial conditions in the NH cool the North Atlantic Deep Water (NADW) and, consequently, via the southward flow and upwelling of NADW, lower the deep temperature in the SO. Because of the permanent ice cover over Antarctica, surface ocean temperature in the SO near Antarctica is close to freezing point during the entire glacial cycle, so it cannot cool very much even during glacial conditions. Glacial conditions therefore increase the density of deep SO water but not of surface SO water and this strengthens the vertical stratification there. As a result, vertical mixing in the SO is expected to be reduced. This mechanism implies that temperature changes in the NH (supposedly correlated with global ice volume V) lead to the CO2 changes in the SH. This mechanism is plausible dA = (V ‒ A ‒ cI60) . τA dt 2.7 88 ON THE PHYSICAL MECHANISMS BEHIND GLACIAL-INTERGLACIAL DYNAMICS −800 −700 −600 −500 −400 −300 −200 −100 0 −2 −1 0 1 2 Time (Kyr) Global Ice Volume 0 20 40 60 80 100 120 0 500 1000 1500 2000 Period (Kyr) Spectral Power Top panel: normalized ice volume as predicted with the 3τ (blue line), LS (red line) models, proxy δ18O records from Lisiecki and Raymo [2005] (grey line) and Bintanja et al. [2005] sea-level time series (green line). Bottom panel: power spectrum of normalized ice volume as predicted with 3τ model (blue), LS model (red), proxy δ18O records from Lisiecki and Raymo [2005] (grey line) and Bintanja et al. [2005 sea-level time series (green line). Figure 2.20 Chapter 2 89 but interpreting F as a model for such mechanism is uncertain, because correlation between V and North Atlantic temperature time series (Bard [2003]) is low, only 0.64 for the last 800 kyr. On the other hand, to determine whether A is sensitive to I60, Equation 2.2 of the model 3τ was modified as However, the best fits obtained after this change did not improve upon those obtained by our 3τ model, showing that Equation 2.2 is also insensitive to I60. These results are consistent with the following conclusions: a variable A that follows the evolution of V with a delay of 5 to 7 kyr seems necessary to obtain good fits. This variable can be interpreted, as Paillard and Parrenin [2004] do, as the ice sheet extent on the Antarctic shelf, but it can also be interpreted as the whole SO ice extent since the capping effect (Keeling and Stephens [2001]) and the control of A on the residual circulation act in the same direction as the deep stratification. CO2 release to the atmosphere in the SO may be controlled by either stratification or sea ice extent, or by both. However, neither of the two corresponding variables in our models (F and A) seems to be sensitive to I60 insolation. The conclusion is that either some of these mechanisms are not sensitive to SO temperature, which is difficult to conceive, or I60 is not a good proxy for SO temperature. Huybers and Denton [2008] pointed out an important fact that had gone unperceived and can bring new light to these conclusions: Antarctic summer duration is highly correlated with I65 insolation. In addition, it is plausible that the Antarctic climate remains in near radiative equilibrium with local heat accumulation, which is controlled by summertime duration. In contrast, northern changes are mediated through the response of the northern ice sheets, which are much more sensitive to insolation at the solstice (I65) than to summertime duration. Thus, SO temperature could be a crucial variable that has influence on the southern sea ice or the density of deep water, or both; but this southern temperature does not need to be controlled via teleconnection with the north, it may respond to some regional astronomical forcing different to I60. 96 ON THE PHYSICAL MECHANISMS BEHIND GLACIAL-INTERGLACIAL DYNAMICS the westerlies belt (Fischer et al. [2010]). However, observations of wind and rain patterns in HS1 and YDS suggest that the Inter-tropical Convergence Zone (ITCZ) and trade winds shifted southward (Peterson et al. [2000]; Leduc et al. [2009]; Saikku et al. [2009]; Wang et al. [2004]; Grffiths et al. [2009]). Two warming pulses in the SH coincided with HS1 and YDS, respectively, suggesting a north-south connection mechanism such as the bipolar seesaw suggested by Broecker [1998]. The 231Pa/230Th ratios from core GGC5 off the Bermuda Rise (McManus et al. [2004]) reflect a strong reduction in Atlantic overturning circulation during HS1 and a moderate reduction during YDS. In contrast, biogenic opal flux in the SO, interpreted as a proxy for changes in upwelling south of the Antarctic Polar Front (Anderson et al. [2009]) shows an increase in upwelling during HS1, a decrease during the Bolling-Allerod period and a further increase in overturning during YDS. This suggests that the shutdown in the Atlantic MOC (AMOC) removed a source of dense water to the ocean interior and that this “density vacuum” (Broecker [1998]) precipitated an increase in AABW formation to fill that vacuum (Sigman et al. [2010]). These pieces of evidence have been interpreted similarly by different authors (see Cheng et al. [2009]; Denton et al. [2010]; Shakun et al. [2012]). The sequence of events would start with a rising boreal summer insolation controlled mainly by inclination and precession cycles (Huybers [2011]). This insolation produces a gradual northern warming in the presence of a massive, isostatically depressed NH ice sheet which has growth gradually for 80-100 kyr. In these conditions, an ice instability mechanism makes the Laurentide, Greenland and European sheets ablate and collapse locally, producing a strong meltwater pulse into the North Atlantic, which initiates a Heinrich event. This slows the AMOC and lowers the south-north Atlantic heat flux, with stadial conditions in the north and southern warming and southward movement of climatic zones in the SH. The above three studies agree that this sequence is the deglacial trigger. A number of mechanisms are potentially able to connect this trigger with an increased release of CO2 in the SO and subsequent atmospheric CO2 rise: (i) Southward movement of climatic zones could include a southward shift in the westerlies, resulting in enhanced wind-driven Chapter 2 97 upwelling on the Antarctic divergence (Marchitto et al. [2007]). This would promote ventilation of CDW and the observed productivity peak (Anderson et al. [2009]), would erode the Antarctic salinity-driven stratification and would thus facilitate the formation of deep water in Antarctica (Toggweiler et al. [2006]). The importance of this mechanism deserve further investigation, as recent studies have shown an oscillatory behavior of AMOC between weak and strong states due to variations on the density gradient in the Atlantic Ocean, derived from variations of CO2 and strength of SO winds (Banderas et al. [2015]). (ii) Warming of the SO could melt sea ice, leading to an increase in ventilation as a result of enhanced residual northward transport and also to degassing effects. As suggested by Fischer et al. [2010], the increase and decrease in sea ice coverage could have modulated the annual net heat gain of the SO surface, and the reduced buoyancy gain would have weakened the MOC during glacial periods through its control of the residual circulation. This would have limited the upwelling of CO2-enriched deep water, acting in the same direction as the capping effect that winter ice extent is assumed to have on deep-water outgassing of CO2 (Stephens and Keeling [2000]). (iii) Warming associated with southerly shifts could reduce Patagonian glaciation, lowering the flux of dust and iron from Patagonia to the SO. This would reduce the fficiency of the biological pump (Martin and Fitzwater [1988]). (iv) Increased CO2 release from Antarctic divergence would decrease the concentration of dissolved inorganic carbon (DIC) and thus increase the carbonate ion concentration in the deep ocean. This would increase the burial rate of calcium carbonate, forcing a decrease in whole-ocean alkalinity, decreasing solubility and putting an additional fraction of CO2 into the 98 ON THE PHYSICAL MECHANISMS BEHIND GLACIAL-INTERGLACIAL DYNAMICS atmosphere. According to Sigman et al. [2010], this carbonate compensation in cooperation with switch on/off of deep water formation or a degassing effect is able to account for 40 ppm of pCO2 difference between interglacial and glacial conditions. It is probable that neither Antarctic overturning nor gas exchange were completely shut off/on. However, more complete nutrient consumption and lower export of organic matter have been observed in the Antarctic during glacial periods, suggesting a reduced supply of nutrients and water from the deep to the surface ocean (Sigman et al. [2010]). Thus, overturning decreased in the glacial SO while productivity did not decrease as much as the nutrient supply from below did, making pCO2 decrease as a result. This mechanism would have complemented the previous ones, generating the atmospheric pCO2 difference of 40 ppm. (v) In the sub-Antarctic zone, the sinking flux of organic matter was greater in glacial periods than today, with some evidence for more complete nutrient consumption (Sigman et al. [2010]). The effect of this zone on alkalinity and availability of nutrients at low-latitudes is great, reaching 40 ppm for the pCO2 glacial-interglacial difference (Hain et al. [2010]). As observed by Sigman et al. [2010] the declines in atmospheric pCO2 to their minima during peak ice ages tend to occur over tens of thousands of years and/or in steps. This may result from the progressive activation of different pCO2 reducing processes (Hain et al. [2010]). In particular, the Antarctic cooling early in the last ice age, 115 kyr BP, has been interpreted as reduced Antarctic overturning or increased sea-ice suppression of gas exchange (Peacock et al. [2006]) with an additional fraction of CO2 uptake due to ocean cooling. In contrast, the second major decline in atmospheric pCO2, 70 kyr BP, coincides with a major dust flux increase in Antarctic and sub-Antarctic zones (Martínez-García et al. [2009]), which may have caused iron fertilization of the SO (Watson et al. [2000]). Additionally, this period may have suffered the sharpest Chapter 2 99 transition from NADW to GNAIW formation (Hodell et al. [2003]), which may have improved the ability of the deep SO to lower atmospheric pCO2 (Sigman et al. [2010]). 2.7. DISCUSSION To what extent are PP04-derived models consistent with the observational dynamics summarized in last section? First, the sequence of activation of processes that release CO2 during the terminations, suggested by the triggering sequence mentioned above, is very different to the sequence of steps that reduce CO2 during glacial times, and this difference generates the typical saw-tooth shape of the CO2 glacial-interglacial oscillations. The lack of a complete model for the biological and carbonate pumps in the PP04-derived models makes it impossible to precisely reproduce the series of events that characterize glacial periods, and using different relaxation times for moments of increasing and decreasing CO2 is the easiest way to take into account this asymmetry and to reproduce the double slope cycle. When a single relaxation time τC is used (τC2 = 0) in the LS calibration, the correlation obtained for V decreases by less than 1% (0.85 instead of 0.86) but the correlation for C decreases by 8% (0.68 instead of 0.74) because of the far more abrupt ups and downs of the simulated C during the glacial periods. The 3τ model is much less sensitive to this τC2 parameter, and maintains a very high correlation when it is eliminated (0.88 and 0.79, for V and C, respectively, in the 3τ model). We thus obtain a model with almost the same correlations as 4τ and with one parameter less than PP04. The good performance in relation to PP04 must be attributed only to the introduction of a different time for accumulation (τV) and for ablation of ice (τV2). Though the glacial decay is only roughly simulated by PP04-derived models, the precise fitting of pace and intensity of interglacial CO2 that these models generates is sufficient to obtain a good match of the ice volume evolution. The term cI60 has been removed from Equation 2.6 in 3τ model as one conclusion of this study is that PP04-derived models are not sensitive to I60 southern insolation. In these models, C seems to be a much better proxy of Antarctic temperature than I60 100 ON THE PHYSICAL MECHANISMS BEHIND GLACIAL-INTERGLACIAL DYNAMICS is. This indirectly reinforces the claim of Huybers and Denton [2008] that Antarctic temperature is not sensitive to insolation during a particular summer day and that it could be sensitive to another kind of variable, such as austral summer duration. The stratification function, F, is a crucial variable of PP04-derived models, which must be interpreted as stratification produced by density difference between AABW and NADW with a possible minor contribution from stratification produced by density difference between NADW and surface water. F controls the upwelling of CO2 in a nonlinear way: an abrupt oceanic release of CO2 similar to a rectangular pulse of 10 to 20 kyr, very sensitive to F, is necessary in these models to generate the terminations. Mixing rates may decrease as density differences increase (Watson and Garabato [2006]) if we use the internal wave parameterization for the vertical mixing coefficient suggested by Gargett [1984]. If the density difference between AABW and NADW is small, the blending process of the two water masses may be more efficient and a larger fraction of the AABW wells up with the mixed Lower Circumpolar Deep Water (LCDW) into the two loops of the meridional overturning circulation (NMOC and SMOC). Bouttes et al. [2012] have shown the plausibility of this stratification-dependent mechanism with a model of intermediate complexity. The changes in ocean circulation associated with westerlies, deep ocean stratification, expansion of sea ice, variation of AMOC efficiency, together with variations of biological uptake and carbonate compensation seems to be the main drivers for the CO2 atmospheric difference between glacial and interglacial states. Ferrari et al. [2014] has recently proposed a new mechanism able to link Antarctica sea ice expansion with the temperature drop, the rearrangement of deep water masses and the change of circulation in LGM, proving that they are not independent mechanisms, but feedbacks. The precise mechanism is lead by the isopycnals of the SO, which control the upwelling and mixing of AABW and NADW, separating the NMOC and SMOC cells. The slope and position of this isopycnal change according with the expansion of summer sea ice, mixing actively AABW and NADW in the interglacial period, but confining the mixing-driven upwelling of abyssal waters at the LGM, closing the SMOC and preventing the mixing of NADW and AABW. As a result, AABW filled oceans basins up to 2 km depth, instead of the present 4 km depth. Chapter 2 101 This closed abyssal overturning cell would act as a ocean storage of carbon. This approach, considers the expansion of quasi-permanent sea ice in the SH as the driver of the termination. This suggest that sea ice extent may have an effect on the CO2 release intensity which works in the same direction than stratification during glacial and interglacial periods; therefore, we cannot discard a possible parallel contribution of sea ice area in the intensity of the CO2 release. PP04-derived models need a quick response of regional SO temperature (represented by C in LS and by –V in 3τ) to I65 rising, as quick as τV2 = 3 kyr in 3τ and τC = 3 kyr in LS. The real climate system could produce a quick response of this kind during deglaciations by the bipolar seesaw mechanism (Broecker [1998]). A is another crucial variable for this kind of models and follows the evolution of V and C with a delay of about 5-7 kyr. This dependence shows that A is sensitive to Antarctic temperature, and the delay suggests that it represents the extent of ice sheet on Antarctic shelves and not on SO sea ice extent. Polynyas are the places where most of the brine is produced and are located between the ice sheet edges and sea ice. Thus, the position of icesheet boundaries is crucial for the rate and density of the deep water formed. According to Anderson et al. [2002], the Antarctic icesheet edge was close to the continental slope in the Last Glacial Maximum (LGM) and it started to retreat onto the continent some thousand years after the LGM on most of the shoreline. This delay seems to be coherent with the delay obtained for the A variable in our calibrations. Therefore, PP04-derived models reinforce an interpretation of the oceanic CO2 release during deglaciations as controlled by the deep stratification rather than by a capping effect. As can be observed in Figure 2.21, A takes about 80 to 100 kyr to reach a level that makes stratification become critical (F < 0). Therefore, the 100-kyr period (main frequency of glacial cycles) must be considered as internally generated by the climate system. More specifically, it must be interpreted as the characteristic time needed for the Antarctic ice sheet to reach the continental slope. The observations of Anderson et al. [2002] on the advance and retreat of ice shelves are not in disagreement with this time scale, but additional observations are needed to contrast this prediction of PP04-derived models. 102 ON THE PHYSICAL MECHANISMS BEHIND GLACIAL-INTERGLACIAL DYNAMICS According to Tziperman et al. [2006], many of the models that best fit glacial oscillations have a “phase locking” between a Milankovitch frequency ωm and an internal cyclicity ωi. Frequently, a good fit can be obtained when the quotient ωm/ωi is a rational number. However, this good fit does not imply that the mechanisms represented by the model are the correct ones. Phase locking is a necessary condition but not a sufficient one. However, Tziperman et al. [2006] conclude that “the actual glacial cycles may also be similarly phase locked to the Milankovitch forcing”. That is, this property is probably also one of the climate system itself. If this is so, showing nonlinear phase locking is a good ‒though not sufficient‒ feature for a correct paleoclimate model. What should be sufficient? There is no clear mathematical criterion for deciding. However, if the right mechanisms were to be incorporated in a model with nonlinear locking, we would expect the model to display the following features: (i) Its predictions should have high correlations with the experimental time series. (ii) It should reproduce the ups and downs of the experimental time series in the scales between thousand to million years better than other models of similar complexity. (iii) It should generate time series with frequency composition and embedding dynamics similar to the observational ones. However, ultimately, only the experimental confirmation that the proposed mechanisms are really acting would be a definitive sufficient condition. The visual inspection of Figure 2.1 as well as the good correlations obtained show that our models offer great improvements in criteria (i) and (ii) in comparison with other simple models studied in the literature. The wavelet analysis and phase space portraits show that the models that we have discussed here satisfy the third criteria, at least for variable V. On other hand, the pace of orbital timescale is well simulated Chapter 2 103 for V and C, but in the sub-orbital timescales the phasing of C is poorly matched in many periods that include the deglaciations (Figure 2.9). In addition, the cross-recurrence analysis shows that the short term coherence between simulated and observational CO2 is only sporadic, indicating that both time series do not follow the same dynamical behavior (Figures 2.18 and 2.19). However, in the deglacial periods the two carbon series become dynamically close. Finally, neither simulated C nor V match the observations on timescales 1-10 kyr, which are the timescales where the first events of the deglacial trigger take place. The most important conclusion of this study is that a detailed dynamics of CO2 is of little relevance to obtain very good agreement between simulated ice volume and δ18O. All that is needed is a nonlinear instability that is reached after 80-100 kyr and that enables that even a modest sub-orbital variation may trigger a sudden release of atmospheric CO2. This nonlinear instability seems to be controlled by V or C (here interpreted as Antarctic temperature) and by A (interpreted as Antarctic ice sheet) in PP04-derived models. More specifically, Antarctic cooling would press toward greater stratification because a cooler regional temperature would precondition the surface layer temperature making easier the water column on the shelf to get the ‒1.9 to ‒2°C needed for freezing and brine formation. On the other hand, ice sheet extent would tend to block the active polynyas, producing an opposite effect. With the progress of glaciation the latter effect would become dominant and would place the system close to instability, due to the decrease in stratification. A precise simulation of the detailed set of events constituting the deglacial trigger would require adding some model for the North Atlantic ice sheet instability to our models as well as more complex carbon models. PP04-derived models produce good fits in spite of their simplicity in the modeling of oceanic CO2 release, which is made to depend only on stratification. The other contributions (carbonate compensation, shoaling or deepening of NADW formation, change of productivity in SO and slowing down of SO overturning) are not included in the model, but should be included in future improvements to obtain a closer representation of all the physical processes that are apparently involved. Heinrich events seem to belong to a cyclic process of calving with a recurrence peri- 104 ON THE PHYSICAL MECHANISMS BEHIND GLACIAL-INTERGLACIAL DYNAMICS od much shorter than deglaciations (Bond and Lotti [1995]). However, according to the interpretations of Cheng et al. [2009]; Denton et al. [2010]; Shakun et al. [2012], Heinrich events are able to produce large iceberg discharges and long stadials only when ice sheets become large enough. Long stadials are needed to produce a strong CO2 oceanic release able to start the termination according to Denton et al. [2010]. However, large oceanic CO2 emissions can be produced only if stratification has weakened at the end of the glacial period (Burke and Robinson [2012]; Martinez-Boti et al. [2015]) for the lack of brine formation (Paillard and Parrenin [2004]; Bouttes et al. [2012]) or may be weakened by the expected SO upwelling activation. Ultimately, stratification change plausibly controls a large fraction of the oceanic CO2 release and is synchronous with other mechanisms which contribute an additional fraction in warm periods, such as biological and carbonate pumps in the sub-Antarctic zone and sea ice extent which controls the residual circulation (Hain et al. [2010]; Sigman et al. [2010]) and the depth of upwelled water (Ferrari et al. [2014]). In this regard, stratification may play the same role as an “order parameter” (Haken [1987]) plays in a synergetic nonlinear change. IMPACT OF ANTHROPOGENIC CO2 ON THE NEXT GLACIAL CYCLE 112 The Intergovernmental Panel on Climate Change (IPCC) published a Special Report on Emission Scenarios (SRES) that contains 40 scenarios for future fossil fuel production to assess future climate change (SRES [2000]). For example, the SRES A1 set uses an average accumulated carbon emission between 1750 and 2100 AD of approximately 1995 Gtce, and the A1F1 extreme scenario implicitly assumes 4100 Gtce of URR, according to the Hubbert linearization done by Berg and Boland [2013]. However, these scenarios are based on simple extrapolations of present emission rates with no consideration on actual URR values. Höök et al. [2010] and Berg and Boland [2013] show that this report is based on unreasonably optimistic expectations about future fossil fuel production which do not take into account realistic estimates for proven and probable fossil fuels reserves. In order to obtain a realistic estimate for U we have used the analysis of both Laherrere [2006] and the Energy Watch Group [2007] on the URR for the three main fossil fuels (coal, oil and gas). The Energy Watch Group [2007] estimates 479 Gt for bituminous coal and anthracite, 272 Gt for sub-bituminous coal, and 158 Gt for lignite. Assuming 92%, 40% and 66% of carbon respectively for these three minerals, we obtain approximately 558 Gtce and the peak of coal-derived CO2 emissions on year 2042 AD. The analysis of Laherrere [2006] gives an estimate of 3000 Gb for all kind of oils (Gb = Gigabarrel, one oil barrel is very close to 159 liters of oil), which amounts to 353 Gtce if we use the conversion factor recommended by EPA2: 0.118 tons of carbon per barrel. The International Energy Agency (IEA [2010, 2012]) recognizes that conventional oil production is at or very close to a maximum since 2006 AD and that this production will hardly increase much more, which amounts to an implicit recognition of the oil peak arrival (Laherrere [2012a]). Regarding gas, Laherrere [2006] predicts an URR of 9100 Tcf, which amounts to about 120-132 Gt of carbon. Recently, Laherrere [2012b] has updated the URR coal estimate to 750 gigatons of oil equivalent (Gtoe). Using a carbon content of 25.7 tce/TJ for average coal, in the upper range of reported values (Gassan-zade [2004]), this amounts to 808 Gtce. Given these considerations, the total fossil fuels URR amounts to approximately 1300 Gtce, being this our best estimate for the integrated anthropogenic emission of carbon. 2 ht t p : //www . e p a . go v/ g r ee np owe r/ pub s/ c a lc m e t h.h t m Chapter 3 113 Figure 3.1 shows the historical data on fossil fuels emissions during the 1800 to 2010 AD period and the corresponding fit using a Lorentz function (Equation 3.1), obtained setting U = 1300 Gt, which corresponds to a peak emissions on year 2037 AD. This fit has been done such that the integrated carbon emission at 2010 AD has exactly the same area than the historical time-series until 2010 AD (dashed lines in Figures 3.1 and 3.2). This modified Lorentz function is the one we use to produce the full time series for the accumulated anthropogenic emission of fossil fuels (Figure 3.2). The peak of emissions occurs in year 2037 AD, coincident with the inflection point in the accumulated curve, and tends to zero near year 2324 AD, when the a ccu m u l a t ed emissions reach t h ei r m a x i m u m v a l u es . 1800 1850 1900 1950 2000 0 2 4 6 8 10 Time !Year AD" Gt#yr Historical data of fossil fuels emissions during the 1800-2010 period (continuous line) and slightly modified Lorentz function such that the area until 2010 equals the accumulated historical emissions (dashed line). Figure 3.1 IMPACT OF ANTHROPOGENIC CO2 ON THE NEXT GLACIAL CYCLE 114 3.3.2. CO2 atmospheric response The relatively fast, at the 10 to 100 yr scale, anthropogenic CO2 emission will lead to a rise in atmospheric concentration but it will also produce an increased storage rate into the ocean, ground vegetation, trees, detritus and soil (Kheshgi and Jain [2003]). Our models do not include this type of relatively rapid feedback climatic mechanisms. To take them into account we require a transfer function of accumulated carbon emissions (Ac, in Gt) to atmospheric CO2 concentration (pCO2, in ppmv). With this objective we turn to the global compartmental model ISAM (Jain et al. [1994]; Kheshgi and Jain [2003]), which may be interactively run online3. 1800 1900 2000 2100 2200 2300 0 2 4 6 8 10 12 1800 1900 2000 2100 2200 2300 0 2 4 6 8 10 12 Time !Year AD" Gt#yr 10 2 Gt Modified Lorentz function illustrating the anthropogenic carbon release per year (dashed line) and the accumulated carbon release (solid line). Figure 3.2 Chapter 3 115 In this model, the rates of transfer from the atmosphere to the above five compartments have been calibrated to match the average projections (by six dynamic global-vegetation models and 10 coupled ocean-atmosphere models) for CO2 uptake between years 2000 and 2200 AD. In order to obtain the Ac to pCO2 transfer function, we use the ISAM model with the A1T scenario parameters, as proposed by the IPCC (SRES [2000]). We run the ISAM model for values of URR between 0 and 3000 Gt, at increasing URR intervals of 100 Gt, replacing the fossil fuels input of A1T scenario with the modified Lorentz function. For each URR value we obtain the corresponding atmospheric response; these final response points are then linearly interpolated in order to have a continuous curve (Figure 3.3). The corresponding transfer function should be adequate to model the short-term (at decadal and century time scales) response in the range from 0 to 3000 Gt; in particular, it may be used to model the results with our URR best estimate, U = 1300 Gt (Figures 3.1 and 3.2). Note that contributions from future deforestation are considered by the ISAM model, therefore they are reflected in the transfer function of Figure 3.3. We may now finally estimate the time evolution of the anthropogenic atmospheric CO2 pulse. The historical data of atmospheric concentrations is used for the period between 1750 and 2010 AD. Thereafter, and until the end of the anthropogenic emission in year 2324 AD, the atmospheric CO2 is estimated using the predicted accumulated emissions from the modified Lorentz function (Figures 3.1 and 3.2) and the transfer function as deduced from the ISAM model (Figure 3.3). In this way the accumulated carbon emission is 361 Gt in year 2010 AD, 642 Gt in year 2037 AD, and reaches a maximum value of 1300 Gt by year 2324 AD, and the atmospheric concentrations at those dates can be obtained from the corresponding ordinates in Figure 3.3. 3http://climate.atmos.uiuc.edu/isam2 IMPACT OF ANTHROPOGENIC CO2 ON THE NEXT GLACIAL CYCLE 116 3.4. PROJECTIONS FOR THE NEXT GLACIAL CYCLE Before making a projection for the next 300 kyr we must take into consideration if there are other plausible mechanisms which are not incorporated in our models. In this work we consider two potential contributions: the very well established weathering compensation mechanism and the less known emission of methane from clathrates. We have chosen these two mechanisms not only because they are indeed relevant effects but also because they are illustrative of the sort of modifications that feedback mechanisms may cause in the Pleistocene dynamics. According to Archer and Ganopolski [2005], the silicate weathering cycle will cause that 7% of the anthropogenic CO2 will still remain in the atmosphere 100 kyr after an anthropogenic perturbation. This implies that, after ending the anthropogenic emission, about 9% of the CO2 concentration will decay with a time constant of approximately 400 kyr. We take into consideration this non-linear buffer (ocean carbon chemistry) effect by simply assuming that 91% of the anthropogenic perturbation will follow the dynamics given by Equations 2.1 to 2.6, while the remaining 9% will have a long-term decay, with a 400 kyr time decay constant. A similar approach was used by Paillard [2006]. Archer et al. [2009a] have studied the greenhouse effect expected to occur in the next 10 kyr as the rising temperatures lead to the emission of a fraction of the methane clathrates from continents and shelves. Their projection has a large uncertainty related to determining the critical fraction of methane bubbles that are able to reach the ocean surface. Taking this fraction as 2.5% and considering a scenario of 1000 Pg of carbon emission, suffciently close to the one used in our models, he predicts an increase of 0.4 – 0.5°C related to the escape of methane during some 10 kyr. Taking a climate sensitivity of 3°C for a doubling of CO2, the effect of the predicted methane emission, Cmet , is equivalent to 40.5 ppmv of CO2 concentration during a period of about 10 kyr. We have introduced this contribution in our projections by assuming that Cmet (additional equivalent CO2 concentration derived from methane) grows linearly from 0 to 40.5 ppmv between 0 and 2 kyr, remains constant between 2 and 10 kyr, and decreases linearly from 40.5 ppmv to 0 between 10 and 20 kyr. Chapter 3 117 Following these considerations, we are now ready to start the projection. The model equations are integrated forward, with the anthropogenic carbon emission scenario in Figure 3.2 as the initial condition and properly forced by the insolation at 65°N (Berger [1978b]), to predict the evolution of the interglacial-glacial transitions during the next 300 kyr. In order to obtain a CO2 projection in dimensional units (ppmv), rather than the model non-dimensional units, a concentration of 280 ppmv is assumed for year 1750 AD and a dimensional scale is obtained by considering the range of atmospheric concentrations in the Taylor Dome Ice Core between the last glacial maximum (181 ppmv) and year 1750 AD (280 ppmv) (Indermuhle et al. [2000]). Similarly, in order to transform from model ice volume to benthic δ18O (per mil), the conversion used is V = (Vnormalized + std(δ18O)) + mean(δ18O),3.2 where Vnormalized is the modeled ice volume normalized with its own standard deviation and mean. All other variables remains non-dimensional as there is no dimensional scale available for their transformation. The projection for the next 300 kyr is computed for both unperturbed and perturbed conditions using either the 3τ or LS models; Figure 3.4 illustrates, as an example, the predicted evolution of atmospheric CO2 following the anthropogenic CO2 pulse. Both models display an exponential decay in atmospheric CO2, with characteristic time decay somewhat longer for 3τ than for LS. The two models, despite having different parameterizations for ocean stratification and ice sheet evolution, do produce a very similar outcome for the following glacial cycle. Hereafter, we will focus on the future projection obtained with the 3τ model, a comparison of both models shows analogous results. Figure 3.5 shows the CO2 evolution predicted by model 3τ following the anthropogenic short-term pulse, after adding the long-term weathering compensation mechanism and further including methane emissions from clathrates. It may be observed that in all cases the timing for the next interglacials does not substantially change. The main effect is caused by the weathering compensation which causes somewhat IMPACT OF ANTHROPOGENIC CO2 ON THE NEXT GLACIAL CYCLE 118 higher concentrations and a slightly longer duration of the present interglacial. The effect of the methane emission is much reduced except for a relatively short plateau in CO2 values after ending the anthropogenic pulse. The introduction of the anthropogenic CO2 pulse clearly perturbs the natural cycle for all four model variables during the forthcoming 300 kyr (Figure 3.6). Notice that the values for variables I65, V, A and F in Figure 3.6 are relative variations, only C has been translated to physical units (ppmv); further note that a minimum bound for V has been set to an additional 10 meters sea level rise with respect to year zero, which corresponds to a melting about 30% higher than the Greenland ice sheet (IPCC [2001]) and, for this reason, V doesn’t take values under ‒ 0.93/ ‒1.94 for 3τ/LS models. 0500 1000 1500 2000 2500 3000 0 100 200 300 400 500 600 700 Ac !Gt" pCO 2 !ppmv" Relation between the accumulated anthropogenic carbon emissions (Ac) and the atmospheric CO2 concentration (pCO2). The data is obtained from an aggregate of A1T scenarios (SRES [2000]) as predicted by model ISAM (Jain et al. [1994]), available at http://climate.atmos.uiuc.edu/isam2/, for every possible value of URR. Figure 3.3 Chapter 3 119 A comparison between the different panels in Figure 3.6 illustrates the influence of I65 on the sequence of events, as discussed with detail in Chapter 2 and García-Olivares and Herrero [2013]. In both the perturbed and unperturbed simulations, a V minimum (maximum) is produced a few thousand years after every positive (negative) peak in I65. During a glaciation, the long-term increase in V produces an increase in A which tends to decrease the F value (due to the inverse dependence of F on A, Equation 2.6). The transit to the interglacial relies heavily on F turning negative (through the Heaviside function in Equation 2.5): F reaches negative values only during periods of advanced glaciation because at these times A is also large. In these situations, peaks in I65 lead to relative minima in V and a subsequent decrease in both A and F, which may induce negative F values that trigger the oceanic CO2 emission. The condition to generate a deglaciation is, thus, the coincidence of a maximum of I65 in a period when V is large, as we have seen previously in Chapter 2. The model predicts a maximum CO2 concentration of 519 ppmv in 2300 AD followed by an exponential decay. To define interglacial conditions we use two simple criteria (Crucifix [2011]): CO2 concentrations above 250 ppm and benthic δ18O below 3.8 per mil; the CO2 criterion has been used also by IPCC [2007]. According to the CO2 criterion the current interglacial would end 40 kyr AP (35 kyr AP) for 3τ (LS), and according to the δ18O criterion it would end 46 kyr AP (40 kyr AP) for 3τ (LS). Therefore, the length of this interglacial would be 54 ‒ 47 kyr (according to the CO2 values) and 58 ‒ 50 kyr (according to the δ18O) as predicted by the 3τ and LS models. This interglacial would be followed by a glacial period lasting until 113 kyr AP. The model also predicts the disappearance of the next interglacial unperturbed age, which should start at 64 kyr AP. This is due to the abnormally reduced ice volume and ice sheet area predicted for the present interglacial, which takes long to recover. The ice sheet, A, follows the evolution of V with a delay of about 7 kyr (Equation 2.2 and τA value in Table 1.1). The F function, however, turns negative only when A is large enough (due to the high value b = 1.205). At 64 kyr AP, V decreases (C rises) in response to the positive forcing of the I65 insolation but A is not high enough for F to become negative. If V and A were large, IMPACT OF ANTHROPOGENIC CO2 ON THE NEXT GLACIAL CYCLE 120 as occurs in the unperturbed case, F would be close to the threshold value (F = 0) and the decrease of V (rise of C ) at that time would drive F to negative values, triggering the oceanic pulse. This delayed recovery of the ice volume V (and A) is caused by the anthropogenic pulse introduced in the model, which produces 20 kyr of abnormally high greenhouse effect. Under these conditions V needs a longer time to reach the glacial maximum and the next glacial cycle moves 44 kyr forward in time. At 235 kyr AP perturbed and unperturbed interglacials coincide again and the further evolution of all variables remain in phase, suggesting the recovery of the natural periodicity. 0 50 100 150 200 250 300 200 250 300 350 400 450 500 Time (kyr) C (ppmv) Prediction of the dimensional atmospheric CO2 concentration for the next 300 kyr using, as an initial condition, the anthropogenic carbon release scenario shown in Figure 3.2. The blue and red lines correspond to the projection by the 3τ and LS models, respectively. Figure 3.4 Chapter 3 121 We have also explored how sensitive are the model results, in particular the starting time for the next interglacial cycle, to different URR values. One main result is that a progressive increase in the anthropogenic pulse leads to smaller, with lower CO2 values, interglacial cycles. A quite relevant result is that the timing for the next interglacial will change discretely as U exceeds different threshold values (inset in Figure 3.5). For example U = 475 Gt represents the URR threshold beyond which the next interglacial will experience a delay of 44 kyr; the following thresholds will occur at 1725 and 2100 Gt, with corresponding additional delays of 42 and 72 kyr (inset in Figure 3.5). 0 50 100 150 200 250 300 200 250 300 350 400 450 500 550 Time (kyr) C (ppmv) 0 500 1000 1500 2000 2500 3000 50 100 150 200 250 300 Total emission (Gt C) Time of next interglacial (kyr) Dimensional atmospheric CO2 concentration for the next 300 kyr as predicted by model 3τ following the anthropogenic short-term pulse (red line), after adding the long-term weathering compensation (green line), and after further incorporating the methane emissions from clathrates (blue line). The inset shows how the start of the next interglacial shifts depending on the total emissions of carbon. Figure 3.5 Chapter 4 INSIGHT TO MARINE ISOTOPIC STAGE 13 USING LATE PLEISTOCENE RELAXATION MODELS AND SEA LEVEL STACK Ch a p t er based on the following published a r t icl es : Herrero, Lisiecki and García-Olivares [2015], Insight to Marine Isotopic Stage 13 using late Pleistocene relaxation models and sea level stack. In preparation. 131 .............. 4.1. INTRODUCTION 132 .............. 4.2. RELAXATION MODELS 134 .............. 4.3. NON-LINEAR ANALYSIS OF BOTH SIMULATED AND PROXY TIME SERIES 137 .............. 4.4. DYNAMICS OF SEA LEVEL RE-CALIBRATED MODELS 137 .............. 4.5. DISCUSSION Contents 131 4.1. INTRODUCTION As we have seen in previous chapter, glacial-interglacial oscillations of late Pleistocene climate (last 800 kyr) reveal a characteristic 100-kyr ice-age cycle, assumed to be mainly derived from orbital parameters and from internal feedbacks of the climate system (Hays et al. [1976]; Archer et al. [2000]; Paillard [2010, 2015]). Following Paillard and Parrenin [2004], several relaxation models, based on simple parameterizations of deep ocean stratification, have been developed (García-Olivares and Herrero [2012, 2013]). Two of these models, 3τ and LS, have good skills reproducing the conditions during the last eight glacial cycles (Chapter 1). The 100-kyr glacial-interglacial periodicity is internally generated through three coupled variables: atmospheric CO2 concentration, global ice volume and the extension of the Antarctic ice shelf (for a detailed dynamics, see Chapter 2) Chapter 4 INSIGHT TO MARINE ISOTOPIC STAGE 13 USING LATE PLEISTOCENE RELAXATION MODELS AND SEA LEVEL STACK It seems to me that the natural world is the greatest source of excitement; the greatest source of visual beauty; the greatest source of intellectual interest. It is the greatest source of so much in life that makes life worth living. - David Attenborough 131 .............. 4.1. INTRODUCTION 132 .............. 4.2. RELAXATION MODELS 134 .............. 4.3. NON-LINEAR ANALYSIS OF BOTH SIMULATED AND PROXY TIME SERIES 137 .............. 4.4. DYNAMICS OF SEA LEVEL RE-CALIBRATED MODELS 137 .............. 4.5. DISCUSSION INSIGHT TO MARINE ISOTOPIC STAGE 13 USING LATE PLEISTOCENE RELAXATION MODELS AND SEA LEVEL STACK 132 δ18O data from benthic foraminifera (Lisiecki and Raymo [2005]) may be used as a proxy for ice volume (Shackleton et al. [2000]; Shakun et al. [2015]; Waelbroeck et al. [2002]), although the O18/O16 ratio is known to depend on both the isotopic composition and temperature of the water where foraminifera develop. Waelbroeck et al. [2002] found that the relationship between δ18O and ice volume is not linear, causing some uncertainties in the ice volume variations. Global ice volume changes derived from growth and retreat of continental ice sheets may be associated with sea level variations in the glacial-interglacial cycles (Chappell and Shackleton [1986]; Waelbroeck et al. [2002]; Lambeck et al. [2014]), a proxy with less uncertainties. Some reconstructions of sea level from ocean sediment core data have been performed using different proxies and models, each of them limited by measurement error, local variations in salinity and temperature, and assumptions particular to each technique. Spratt and Lisiecki [2015] have compiled a wide representation of these reconstructions, developing a sea level stack, which represents the eustatic sea level record more accurately than each of the individual reconstructions. In this work, we use Spratt and Lisiecki [2015] sea level stack to analyze and identify differences between δ18O and sea level as proxy data for the ice volume. In next section, the models are presented and some non-linear methods are applied in Section 4.3 to outline the differences between δ18O and sea level stack. Section 4.4 shows the recalibration and optimization of 3τ and LS models to the new set of sea level data, and the discussion of the results, as well as some brief conclusions, are presented in last section. 4.2. RELAXATION MODELS There are several differences between models 3τ and LS (Figure 4.1, bottom panel). One difference is the value used for the reference Antarctic ice sheet, either ‒C (representing the inverse effect of Antarctic temperature in Antarctic ice sheet extent) for model LS or V for model 3τ. Another minor difference is the inclusion of I65 in model LS when specifying the reference atmospheric CO2 concentration value, Cr. Chapter 4 133 However, the main difference between both models, is their parameterization of the ocean state, with F = F (V, A) in model 3τ and F = F (C, A) in model LS. Both V and C are good proxies for the Southern Ocean (SO) regional temperature and either F = F (V, A) or F = F (C, A) are plausibly ways to model the local formation of brines. The fact that the dependence F = F (V, A) performs somewhat better than F = F (C, A) may be due to the non-negligible role that V has on stratification through teleconnections (García-Olivares and Herrero [2013]; Gildor and Tziperman [2001]) or to a possible larger effect of sea level than Antarctic temperature on brine formation (Chapter 2). −3 −2 −1 0 1 2 3 Proxy −800 −700 −600 −500 −400 −300 −200 −100 0 −3 −2 −1 0 1 Time (Kyr) Models 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 Top panel: proxy δ18O records from Lisiecki and Raymo [2005] (grey line) and sea level stack from Spratt and Lisiecki [2015] (green line); numbers represent Marine Isotopic Stage labels for the early Pliocene from Lisiecki and Raymo [2005]. Bottom panel: normalized ice volume as predicted with the 3τ (blue line) and LS (red line models. Figure 4.1 INSIGHT TO MARINE ISOTOPIC STAGE 13 USING LATE PLEISTOCENE RELAXATION MODELS AND SEA LEVEL STACK 134 The best-fit correlations are 0.87/0.88 (between δ18O and modeled V ) and 0.76/0.79 (between the reconstructed atmospheric CO2 concentration and modeled C ) for models LS/3τ. The parameter values that lead to the best data fit for either model are shown in Table 4.1. The two different proxies, δ18O records from Lisiecki and Raymo [2005] and sea level stack from Spratt and Lisiecki [2015], have similar dynamics but the timing of deglaciations is slightly shifted in all cases, as well as depth of glacial cycles, specially in the first four cycles (from ‒800 to ‒400 kyr) (Figure 4.1, top). The difference is, somehow, more evident around ‒500 kyr where the sea level data reconstruct the cycle with less variance. A comparison between sea level stack and benthic δ18O has been performed in Spratt and Lisiecki [2015], showing that the relationship between benthic δ18O and sea level is well-described by a linear function in the first four cycles (from ‒800 to ‒400 kyr) and a quadratic function in the last cycles (from ‒400 to 0 kyr). 4.3. NON-LINEAR ANALYSIS OF BOTH SIMULATED AND PROXY TIME SERIES Some non-linear analysis have been performed to analyze internal frequencies of each time series. Spectral power of observational and simulated time series are obtained with Fast Fourier Transform (Figure 4.2). The sea level stack has more power than any other series in the 100-kyr band, but the power in the 41-kyr band is visibly lower than any other time series. As it seems, the sea level stack tends to over predict the power for long wave frequencies and to under predict it for short wave frequencies when compared with δ18O data. The complex cross-wavelet transform of two time series can be interpreted as the shared power in a given periodicity band (absolute value) and the phase between the two series in time frequency space (Grinsted et al. [2004]). Lighter purple color indicates greater shared power in that time and periodicity band, and the arrow angle indicates the phase between the observational and simulated series. It is in- Chapter 4 135 teresting to see that in the main 100-kyr band the shared power is maximum (as expected) and the phase is really similar for all three cases (Figure 4.3), but the most interesting features are the blue spots in the 41-kyr band at the sea level and δ18O case. Apparently, there are some significant differences between both proxies around ‒500 and ‒250 kyr, in terms of obliquity content, which made both proxies dynamically different at those times. 0 20 40 60 80 100 120 0 500 1000 1500 2000 2500 Spectral Power Period (Kyr) Spectral power of normalized ice volume as predicted with 3τ model (blue), LS model (red), proxy δ18O records from Lisiecki and Raymo [2005] (grey line) and Spratt and Lisiecki [2015] sea-level stack (green line). Figure 4.2 INSIGHT TO MARINE ISOTOPIC STAGE 13 USING LATE PLEISTOCENE RELAXATION MODELS AND SEA LEVEL STACK 136 1/64 1/32 1/16 1/8 1/4 1/2 1 2 4 8 16 32 64 Time (Kyr) Frequency (Kyr) −700 −600 −500 −400 −300 −200 −100 0 0.20 0.39 0.78 1.56 3.13 6.25 12.5 25 50 100 200 1/64 1/32 1/16 1/8 1/4 1/2 1 2 4 8 16 32 64 Time (Kyr) Frequency (Kyr) −700 −600 −500 −400 −300 −200 −100 0 0.20 0.39 0.78 1.56 3.13 6.25 12.5 25 50 100 200 1/64 1/32 1/16 1/8 1/4 1/2 1 2 4 8 16 32 64 Time (Kyr) Frequency (Kyr) −700 −600 −500 −400 −300 −200 −100 0 0.20 0.39 0.78 1.56 3.13 6.25 12.5 25 50 100 200 Cross-wavelet transform between proxy sea level stack from Spratt and Lisiecki [2015] and δ18O records from Lisiecki and Raymo [2005] (top), simulated V time series with 3τ model (middle) and LS model (bottom). Figure 4.3 Chapter 4 137 4.4. DYNAMICS OF SEA LEVEL RE-CALIBRATED MODELS Considering the existing dynamic differences between the sea level stack and the benthic δ18O, a new calibration has been done for both 3τ and LS models using the sea level (SL) stack from Spratt and Lisiecki [2015]. This has derived in new 3τ and LS models, called 3τSL and LSSL hereafter. The comparison between parameter values of all models is shown in Table 4.1. The correlations between the original models and the new recalibrate models do not change much; for V, the correlation has increased from 0.88 (3τ) to 0.89 (3τSL), just 1%, while there is no change for the LS/LSSL case. For C, however, the correlation has slightly decreased in both cases from 0.79 (3τ) to 0.77 (3τSL) and from 0.76 (LS) to 0.73 (LSSL), a 2% and a 3% respectively. This change in the correlation should not be important, as the uncertainties in the δ18O may produce an over interpretation of the data, meaning that these correlation differences of less than 4% are included under the observational uncertainty. Nevertheless, the most remarkable feature of these new optimizations is a double peak shown in the 3τSL ice volume around ‒500 kyr. This double peak was not present in the original 3τ models, and it is not shown in any LS model, suggesting that there are a significant dynamic difference between both proxies and/or both models. 4.5. DISCUSSION Marine Isotopic Stage 13 (MIS-13) interglacial occurred approximately 500,000 years ago (‒500 kyr). It is of particular interest, as it suffered severe summer monsoons simultaneously with increasing marine oxygen isotope and decreasing Antarctic ice core records in temperature compared to other interglacials (Yin and Guo [2008]; Lang and Wolff [2011]; Muri et al. [2012, 2013]). All of these anomalies indicate a warm Northern Hemisphere (NH) and a cool Southern Hemisphere (SH), and consequently a strong asymmetry of hemispheric climates during MIS-13 (Guo et al. [2009]). 145 I n this final chapter I will briefly revise the results presented in this thesis. In Chapter 1 several relaxation models incorporating a wide representation of physical mechanisms, like the oceanic CO2 pumping or the response times of carbon and ice volume, have been developed. The models’ parameters have been calibrated to provide the best fit to the δ18O and CO2 experimental time series available for the last 800 kyr BP. We have described eight different sub-models, all derived from the original Paillard and Parrenin [2004], to evaluate how those different parameterizations affect the data fit. Some of those mechanisms like the biological exportation production or an exponential CO2 oceanic pulse does not offer new insights. On the other hand, the sub-models with different response times for accumulation and ablation of ice, and/or emission and absorption of CO2, show a good quantitative and qualitative agreement with the empirical time series, especially 4τ , EP, 3τ and LS, suggesting that different relaxation times are a possible way to reproduce the asymmetry in the glacial cycles and that the mechanisms that these models incorporate may be important factors controlling glacial-interglacial oscillations. Conclusions We on Earth have just awakened to the great oceans of space and time from which we have emerged. We are the legacy of 15 billion years of cosmic evolution. We have a choice: We can enhance life and come to know the universe that made us, or we can squander our 15 billion-year heritage in meaningless self-destruction. What happens in the first second of the next cosmic year depends on what we do, here and now, with our intelligence and our knowledge of the cosmos. ‒ Carl Sagan, COSMOS 146 The 4τ model improves the PP04 correlations (0.59 and 0.63) to 0.79 and 0.89 using 14 parameters. The EP model uses 15 parameters but the additional parameter, which was aimed at improving the form of the oceanic pulse function, did not improve the correlations and can be considered useless. The 3τ model obtains almost the same correlations as 4τ using only 13 parameters. The LS model (15 parameters) does not improve the correlations of 4τ (which are slightly higher) but incorporates a functional form for stratification F that seems more consistent with the mechanisms suggested by Paillard and Parrenin [2004]. 3τ is, then, the model with the best explained variance per parameter, and LS, having a explained similar variance, contain parameterizations that can be related to the observed mechanisms. Both models offer a valuable insight into the real meaning of the good performance of PP04-derived models, and deserve further analysis. In Chapter 2, the dynamics and mechanisms incorporated in 3τ and LS models have been analyzed and compared, to explore the most plausible physical interpretations of the mathematical expressions. First, we have evaluated their non-linear dynamics using wavelet transform, cross-wavelet transform, wavelet coherence, Fourier analysis, attractors and cross-recurrence plots. We have, then, identified the mechanisms and dynamics that lead to good data fit and we have compared the results with the observational dynamics that supposedly cause glacial-interglacial oscillations. We have demonstrated that, in fact, relaxation models are a useful tool to analyze and study the physical mechanisms of the climate system, as well as to identify feedbacks and the involved variables. We have shown that the models are not sensitive to the southern insolation; rather, they respond to a different variable, C, which seems to be a much better proxy of Antarctic temperature. One important conclusion if that the detailed dynamics of CO2 are not important to obtain a good match between the models and the proxy data. What is really important is a nonlinear instability, reached after 80 - 100 kyr (allowing V to reach a maximum level), allowing any modest increase of the insolation to cause a release of atmospheric CO2. Variable F is crucial in the models. We have pointed out that an abrupt oceanic release of CO2, similar to a rectangular pulse of 10 to 20 kyr, is necessary to trigger the termination, and F is the non-linear control of this upwelling of CO2. It can be Conclusions 147 described as a deep ocean control parameter. We believe that F can be physically considered the stratification synchronously combined with the biological and carbonate pumps in the sub-Antarctic zone, and with the sea ice extent which controls the residual circulation and the depth of upwelled water. A precise simulation of the detailed set of events constituting the deglacial trigger would require adding some parameterization for the North Atlantic ice sheet instability to our models as well as more complex carbon dynamics. PP04-derived models produce good results in spite of their simplicity in the modeling of oceanic CO2 release, and to obtain a closer representation of all the physical processes involved in the climate system, more mechanisms should be included in future versions. In Chapter 3 we have illustrated an application of the relaxation models, useful too to answer more direct dynamic questions. We have used 3τ and LS models to predict the future evolution of global Earth variables during the forthcoming 300 kyr, with and without the atmospheric CO2 perturbation caused by anthropogenic fossil fuels emissions. The anthropogenic CO2 pulse produces 20 kyr of abnormally high greenhouse effect, involving a delay in the future advance of the ice sheet over the Antarctic shelf. As a result, the corresponding peak of northern insolation, causing a termination in the unperturbed scenario, will have no effect on the stability of the developing glaciation. However, the following insolation peak will take place in an appropriate state of the climate system and will be sufficient to induce the new deglaciation, moving, accordingly, the next glacial cycle 44 kyr forward in time. After three cycles, perturbed and unperturbed interglacials coincide again and the further evolution of all variables remain in phase, suggesting the recovery of the natural periodicity. On the other hand, a progressive increase in the anthropogenic pulse leads to smaller, with lower CO2 values, interglacial cycles, and the timing for the next interglacial will change discretely as U exceeds different threshold values. In Chapter 4 the models have been re-calibrated using a new sea level data stack and are used to understand the dynamics of a very particular interglacial event, the MIS-13. We have compared 3τ and LS models with Spratt and Lisiecki [2015] sea 148 level stack to analyze and identify differences between δ18O (Lisiecki and Raymo [2005]) and sea level as proxy data for the ice volume. The comparison of the deep ocean parameter, F , with a proxy of deep stratification (benthic δ13C data, Lisiecki [2010]) shows a mismatch in the first cycles (from ‒800 to ‒400 kyr), suggesting that F might be representing not only the deep water formation but also other variables, like the Antarctic ice sheet. On the last cycles (from ‒400 to 0 kyr), F seems to be in good pace with the stratification, although slightly shifted, suggesting a representation of a mixed state of the deep ocean. This reinforces the results of Chapter 2, where F represents not only the stratification but also the sea ice extent. The appearance (or not appearance) of a double CO2 peak at MIS-13 in our models, suggest that the parametrization F = F (V, A), affected by both NH and SH dynamics, is more appropriate to simulate the climate. We may then assess that relaxation models can be extremely useful tools to characterize the complex climate system, and helpful to examine specific questions as the future evolution of climate under different scenarios of anthropogenic fossil fuels emissions. Relaxation models may contribute to broaden the knowledge of the mechanisms controlling the variability of the late Pleistocene glacial cycles; nevertheless, their results should be compared with those obtained from an intermediate complexity model, which should include conservation equations, realistic geometry and long-term global carbon dynamics. The analysis of the present results under the complementary perspective given by these complex models would be a natural continuation of this thesis. Appendix 150 ............ Appendix A. On insolation forcing 152 ........... Appendix B. Genetic algorithms for optimization 155 ............ Appendix C. Resumen en castellano 150 Appendix A We have used Berger [1978a] and Berger and Loutre [1991] sofware to calculate insolation time series at different latitudes, covering a time domain from last 800 kyr to the next 300 kyr in time steps of 100 years. We have used negative and positive values when respectively referring to times before present (BP) and after present (AP), where year zero is taken as 1950 AD. In Chapters 1, 2 and 4 we have taken into account only the past time (last 800 kyr) but in Chapter 3 we have made some projections, using the next 300 kyr domain. Two sets of insolation data has been developed for the past time. On one hand, insolation regarding one single day of the year at one specified latitude, this is Northern Hemisphere summer insolation (65°N) on 21st June, named I65, and on the other hand, late Austral summer insolation (60°S) on 21st February, named I60 (Figure A.1). For Chapter 3, a new set of Northern Hemisphere summer insolation (65°N) on 21st June has been calculated, this time considering the whole time domain (from ‒800 to +300 kyr), shown in Figure A.2. A. ON INSOLATION FORCING 151 Appendix A Figure A.1 Northern Hemisphere summer insolation (top) and late Austral summer insolation (bottom) for the past 800 kyr −800 −700 −600 −500 −400 −300 −200 −100 0 −3 −2 −1 0 1 2 3 i 65 −800 −700 −600 −500 −400 −300 −200 −100 0 −3 −2 −1 0 1 2 i 60 Time (kyr) Figure A.2 Northern Hemisphere summer insolation (65°N) from the last 800 kyr to the next 300 kyr. −800 −700 −600 −500 −400 −300 −200 −100 0 100 200 300 −3 −2 −1 0 1 2 3 i65 Time (kyr) 152 To perform optimizations of mathematical expressions, we have implemented a multi-objective genetic algorithm with the aim of maximize correlation between experimental and modeled data. This has been done with Global Optimization Toolbox from Matlab1. First of all, we must answer a key question: what exactly is a genetic algorithm? According to Charbonneau [2002], genetic algorithms are, fundamentally, a class of search techniques that use simplified forms of the biological processes of selection/ inheritance/variation. This optimizers are based on natural selection, which states that individuals better adapted to their environment, i.e., for whatever reason better at obtaining food, avoiding becoming lunch, and finding/attracting/competing for mates, will, on average, leave behind more offspring than their less apt colleagues. B. GENETIC ALGORITHMS FOR OPTIMIZATION 1Matlab is a trademark of The MathWorks. Appendix B 153 For natural selection to lead to evolution, two more essential ingredients are required: 1inheritance: offspring must retain at least some of the features that made their parents fitter than average, otherwise evolution is effectively reset at every generation. 2variability: at any given time individuals of varying fitnesses must coexist in the population, otherwise natural selection has nothing to operate on. To better understand the mechanism, here we show a genetic optimization problem based in our model. One is given a model that depends on a set of parameters u (like 3τ), and a functional relation f (u) that returns a measure of quality, or fitness, associated with the corresponding model; in our case, this function is the correlation between the modeled output and the experimental data. The optimization task usually consists in finding the “point” u* in parameter space corresponding to the model that maximizes the fitness function f (u); in our case, we search in a parameter space of 11 to 15 coordinates, depending on the model. We now define a population as a set of Np realizations of the parameters u. A top-level view of a basic genetic algorithm is then as follows: 1Randomly initialize population and evaluate fitness of its members (parents). 2Breed selected members of current population to produce offspring (child) population (selection based on fitness). 3Replace current population by offspring population. 4Evaluate fitness of new population members. 5Repeat steps (2) through (4) until the fittest member of the current population is deemed fit enough. Appendix B