δ Earth and Planetary Science Letters

Contents lists available atScienceDirect

Earth and Planetary Science Letters

www.elsevier.com/locate/epsl

Late quaternary climate, precipitation δ ¹⁸ O, and Indian monsoon variations over the Tibetan Plateau

Jingmin Li

, Todd A. Ehlers

, Martin Werner

, Sebastian G. Mutz

, Christian Steger

, Heiko Paeth

aDept.ofGeosciences,Univ.ofTübingen,Germany bClimateScienceDivision,AlfredWegenerInstitute,Germany cInstituteofGeographyandGeology,Univ.ofWürzburg,Germany

a r t i c l e i n f o a b s t ra c t

Articlehistory:

Received19February2016

Receivedinrevisedform16September 2016

Accepted21September2016 Availableonline4November2016 Editor: A.Yin

Keywords:

Tibet

precipitationδ¹⁸O monsoonintensity mid-Holocene LastGlacialMaximum GCM

The Himalaya–Tibet orogencontains oneof thelargest moderntopographic and climate gradientson Earth.Proxydata fromthe regionprovideabasisforunderstandingTibetanPlateaupaleoclimateand paleo elevationreconstructions. Paleoclimate model comparisonsto proxy data compliment sparsely located data and can improve climate reconstructions. This study investigates temporal changes in precipitation, temperatureandprecipitationδ¹⁸O(δ¹⁸Op)overtheHimalaya–TibetfromtheLastGlacial Maximum (LGM) to present. We conduct a series of atmospheric General Circulation Model (GCM, ECHAM5-wiso) experiments at discrete time slices including a Pre-industrial (PI, Pre-1850 AD), Mid Holocene(MH,6 kaBP)andLGM(21 kaBP)simulations.Modelpredictionsarecomparedwithexisting proxyrecords.ModelresultsshowmutedclimatechangesacrosstheplateauduringtheMHandlarger changesoccurringduringtheLGM.DuringtheLGMsurfacetemperaturesare∼².0–4.0^◦C loweracross theHimalayaandTibet,and>5.0^◦CloweratthenorthwestandnortheastedgeoftheTibetanPlateau.

LGMmeanannualprecipitationis200–600 mm/yrloweroverontheTibetanPlateau.Modelandproxy datacomparisonshows agoodagreementfortheLGM,butlargedifferencesfortheMH.Largedifferences are alsopresent betweenMHproxy studiesnear eachother.The precipitationweightedannualmean δ¹⁸Op lapse rate at the Himalaya is about 0.4h/^{km larger during the MH and 0}.2h/^{km smaller} during the LGM than during the PI. Finally, rainfall associated with the continental Indian monsoon (between 70^◦E–110^◦Eand 10^◦N–30^◦N)is about44% less inthe LGMthanduring PI times.The LGM monsoon periodisabout onemonthshorterthaninPItimes. Takentogether,theseresultsdocument significantspatialand temporalchangesintemperature,precipitation,and δ¹⁸Op overthelast∼^{21 ka.}

ThesechangesarelargeenoughtoimpactinterpretationsofproxydataandtheintensityoftheIndian monsoon.

©^{2016ElsevierB.V.Allrightsreserved.}

1. Introduction

Paleoclimateandpaleoenvironmentalproxyrecordsfromthe TibetanPlateauanditssurroundingareasprovideabasisforrecon- structing Tibetan Plateau climate change and paleo-elevation re- constructions.Proxyrecordsarereconstructedfromdifferenttypes ofmaterials(e.g.,fossils;lacustrine,aeolian,andglacialsediments;

soilcarbonates;vegetationandpollen)whichcontaindifferentele- mentsandisotopes(e.g.,O,C,H,etc.).Theserecordsareoftenused to document changes in precipitation and temperature. Tibetan proxyrecordshavepreviously beenusedtoreconstructpaleo cli-

*

E-mailaddress:todd.ehlers@uni-tuebingen.de(T.A. Ehlers).

mate ranging from the Early Cenozoic (∼^{50 Ma; e.g. Graham et} al., 2005; Kent-Corson et al., 2009) to Late Quaternary (∼^{2.1 ka} BP topresent;e.g. Herzschuhetal., 2006). Agrowing numberof proxystudiescomefromLateQuaternarylacustrinesedimentsand icecores,andareusedtoinvestigaterecentclimatechangeonthe Tibetan Plateau (see data compilation insupplemental Tables S1, S2). Mostexisting studies contain recordsofHolocenevariations, andafewrecordsextendbacktotheLastGlacialMaximum.Inter- estingly, resultsfrompreviousobservationalrecordsindicatelarge spatial andtemporalvariations inclimaterelativetothemodern, but with the shortcoming that it is difficult to extrapolate be- tweenindividualstudyareasandtounderstandthenaturalspatial variabilitythatcouldbepresentacrosstheplateau.Thisstudyad- dressesthisshortcomingusinganisotope-trackingGCM.

http://dx.doi.org/10.1016/j.epsl.2016.09.031 0012-821X/©^{2016ElsevierB.V.Allrightsreserved.}

Previous paleoclimate reconstructions from observational recordson the Tibetan Plateau indicate a large degree of spatial variationinplateauclimateandprecipitationδ¹⁸O(δ¹⁸Op).Forex- ample, MH proxy records from sediment cores in Cuoe Lake in centralTibet (Wu etal., 2006) and ZabuyeLakeinwestern Tibet (Wangetal.,2002) suggestedawarmandwetterMHclimatethan today.Proxy studiesfromLake Kuhaiinthenortheastern Tibetan Plateau (Wischnewski etal., 2011) andAhungCo Lakein central Tibet(Morrilletal.,2006) suggesteddrierMHconditions.However, someMHclimatereconstructionsfromtheTibetanPlateauarein- consistent with each other. Disagreements between MH climate reconstructionsexist in recordsfrom the same,and neighboring, geographiclocations.Forexample,twonearbylakes(SumxiCoand LakeBangong) inwestern of Tibet,Gasse etal. (1991)suggested warmandwetterMHconditionsthanpresent,whereasVanCampo etal. (1996) suggesteddrier conditionsthan present. In contrast to the MH, LGM climate reconstructions from previous Tibetan Plateau studies provide a more consistent picture of the paleo precipitationandtemperaturedistribution(e.g.Kotliaetal.,2010;

Hodelletal.,1999; Herzschuhetal.,2006; Wangetal.,2002).Pre- vious studies document a colder and drier than modern climate acrossTibetexceptonestudythatsuggestscolderLGMconditions, butsimilarprecipitationastoday, inEasternNepal(Asahi,2010).

Thus, despite the advances of previous studies, they leave some uncertaintyconcerningthespatialandtemporalvariations inLate QuaternaryTibetanclimate,andwarrantalternativeapproachesto understandclimatechangeintheregion.

The complex spatial distribution in oxygen isotopes over the Tibetan Plateau region has been also observed from direct pre- cipitationmeasurements of modern δ¹⁸Op (e.g. Tianetal., 2007;

Yaoetal., 2013)andalsowithin streamwater(Hrenetal.,2009).

Thosestudiessuggestastrongtemperatureeffectonδ¹⁸Op inthe northeastoftheTibetanPlateauandastrongprecipitationamount effectinthemiddle andsouth oftheTibetan Plateau. Besidethe δ¹⁸Opisotopedistribution,anothervariableinvestigatedhereisthe δ¹⁸O_p isotope lapserate. It is often used forclimate reconstruc- tions andfor estimates of the paleoelevation of an orogen (e.g.

ChamberlainandPoage,2000),butitvarieswithclimateandenvi- ronment.Elevationreconstructionsbasedonmodernδ¹⁸O_pisotope lapseratemaycontain climaticuncertaintiesthat overwhelmthe elevationsignal(e.g.EhlersandPoulsen,2009).

Understanding the spatial and temporal variations in climate, δ¹⁸Op, and the driving forces for these changes requires a sys- tematicinvestigation ofthe relevant processes active atdifferent time intervals. Atmospheric GCMs provide one means of investi- gatingspatialvariationsinpaleoclimateatspecifictimeslices.For example,Dallmeyerand Claussen (2011)investigateprecipitation changes from the MH to modern in the Asian monsoon region using the coupled ECHAM5/JSBACH model.They found that pre- cipitationchangesintheEastAsianmonsoon regionare notonly dependent on changes in the Indian summer monsoon circula- tion, but also on changes in the mid-latitudinal westerly winds that dominatecirculationduring thepre-monsoon season. Zheng et al. (2004) conducted paleo climate time slice experiments at 6 ka BP and 21 ka BP for the East Asian Monsoon region us- ingaregionalclimatemodel(RegCM2) withdetailedlandsurface processes. Their results show the East Asian summer monsoon strengthening, a precipitation increase, and a shift of the mon- soon rain belt westwards and northwards during the MH. They alsodemonstratedawintermonsoonstrengtheningduetolowat- mospherictemperaturesduring theLGM. Despitetheseadvances, the previous studies do not provide estimates for proxyrecords suchasδ¹⁸Opbecauseisotope-trackingGCMswerenotused.

Inthepastdecade,isotope-trackingGCMshavebecomeapow- erfultoolforunderstandingtheinteractionsofdifferentEarthsys- tems andvariations in δ¹⁸Op (e.g. Feng et al., 2013). This study

presents results from an isotope-tracking GCM (ECHAM5-wiso) for two paleo time slices including the Middle Holocene (MH, 6 ka BP) andLast Glacial Maximum(LGM, 21 ka BP). Thesetwo time slicesare combinedwith a simulation of Pre-Industrial (PI, Pre-1850AD)conditionstocapturetemperatureandprecipitation, and δ¹⁸Op variations during the most recent glacial–interglacial cycle. We complementprevious model work by (1) investigating climatechangeusingtime-specificboundaryconditionsfortheMH and LGM, and (2) evaluating changes in δ¹⁸Op in the MH and LGM.WeinvestigatethesechangesincomparisontoPre-Industrial (PI) conditions.In the following,we presentspatial andseasonal variationintemperature,precipitation,δ¹⁸Op,andtheδ¹⁸Op lapse rateforalltimeslices.Wecomparemodelpredictionswithproxy data, interpret variations in the monsoon duration and intensity, andidentify theclimate controlsδ¹⁸Op over the TibetanPlateau.

Finally, GCM predictions of δ¹⁸O_p are compared to the simpler approach of a 1D Rayleigh Distillation Model (RDM) to assess if thissimplerapproachcapturestherelevantpastandpresentfrac- tionationprocessesactive inthehigh-reliefareassurrounding the TibetanPlateau.

2. Modelsetupandmethods 2.1. Modelsetupfordifferenttimeslices

Climateandδ¹⁸O_p variationovertheTibetanPlateauareinves- tigated using an atmosphericgeneral circulationmodel equipped withisotopecalculationcapabilities(ECHAM5-wiso).Inthismodel, water isotopes (HDO, H₂¹⁶O, andH₂¹⁸O, treated as independent tracers) experience kinetic fractionation and equilibrium during phase transitions(e.g.vapor, cloud, snowetc.)in theatmosphere (Werneretal., 2011). Severalpublicationshaveevaluatedits per- formanceanddemonstratedagoodglobalandregionalagreement ofobservational data withthesimulatedisotopic fractionofpre- cipitation(e.g.Werneretal.,2011; Langebroeketal.,2011).Inthe Tibetanregion,Yaoetal. (2013)demonstratedthathigh-resolution atmosphericmodels(includingECHAM5-wiso)capturethetempo- ralandspatialdistributionofwaterisotopeinprecipitation.

In this study, simulations are conducted at a resolution of T106L31 (spatial resolution of 1.1^◦×¹.1^◦, and 31 vertical lev- els). The simulations are set up with LGM, MH andPI boundary conditions. The simulations are run for >13 model years. Since isotopes arenot tracked insubsurfacewaters,simulations havea short(<2 yr)spin-uptimeandthelast>10 yrareusedforanal- ysis. Pre-industrialoceanboundaryconditionsare takenfromthe AMIP2 observational climatologic data. MH andLGM sea surface temperatures (SST) and sea ice concentrations (SIC) are derived from the coupledatmosphere–ocean model ECHO-G (Lorenz and Lohmann,2004) asdescribedinDietrichetal. (2013).Thevegeta- tionboundaryfortheMHandLGMisconstructedfromthePale- oclimateModellingIntercomparisonProject(PMIP)PhaseIIobser- vationaldata(http://pmip2.lsce.ipsl.fr; Braconnotetal., 2007) and its gapsare filledwithvegetationmodel predictionsfromArnold etal. (2009). FortheLGMclimate,the land–seadistribution,and icesheetextentandthickness,followtheguidelinesofPMIPPhase III (https://pmip3.lsce.ipsl.fr). Thegreen houseconcentrationsand orbitalparameterduringthePI,MHandLGMarebasedontable 1 inDietrichetal. (2013)andfromOtto-Bliesneretal.(2006).

2.2. Rayleigh DistillationModel(RDM)

TheRDMmodeliswidelyusedintheobservationalstudiesfor paleo elevation reconstructionsdue to its simplicity (e.g.Rowley and Garzione, 2007). The modelis based on the conservation of moist static energy. It assumes the condensation of a single air parcel during adiabatic cooling. When an unsaturated air parcel

ascends,itfirstexperiencescoolingatthedryadiabaticlapseun- tilits temperaturedropstothedewpoint.Followingthis,itcools at the moist adiabatic lapse rate. Under the assumption that all condensedvaporprecipitates,theremainingvaporfractionandits deltavalue can be calculated based on the liquid–vapor equilib- riumandtheestimatedtemperatureandelevationcurves.Thede- tailedalgorithmusedcanbefoundinRowleyandGarzione (2007) andFengetal. (2013).Inthisstudy,RDMδ¹⁸Opiscalculatedstart- ingwithaGCM-derivedmoisturesourceatlowelevation(foreland oftheHimalaya).TheGCMsimulatedδ¹⁸O_p changewithelevation isthencomparedtotheRDMpredictiontoevaluatethedifference betweenthe approaches over the large elevation gradient ofthe Himalaya.

The GCM predicted δ¹⁸Op considers more complex physical processes than the RDM.This is because an isotopic counterpart (H₂¹⁶O, H₂¹⁸O and HDO) for all components of the hydrological cycle(e.g. vapor, cloud, snow anddrainage) aretreated indepen- dently in the GCM prediction and undergo kinetic fractionation andequilibriumduringtheirphasechangesintheatmosphere.The isotopeeffects neglected bythe RDMinclude: e.g. kineticeffects, partialsub-cloudre-evaporationofraindroplets.

2.3. Trajectoryanalysisforvapor sourcechanges

Lagrangian trajectory analysis is a method for quantifying air masstrajectoriesthat controlthetransportofδ¹⁸OinaGCM.For thetrajectoryofaspecificinfinitesimally smallairparcel,thepar- celpositionatacertain time isdependentonthe velocityvector of that time andits previous location. Using the GCM predicted 3Dwindfields,thebackwardstrajectoriesofvaporcanbetracked.

A comprehensivedescriptionoftheLagrangianbacktrajectorycal- culationusedherecanbefoundinFengetal. (2013).

Thetrajectoryanalysisisappliedtothe modelsimulations us- ing6-hourlymodel outputincludingthewind velocityin apres- surelevel system. The wind velocity is linearly interpolatedto a 20-mintime interval. The target location (i.e. the location where thebacktrajectoryiscalculatedfrom)usedinthisstudyisinthe southeastern Himalaya.Thislocation was selected becauseit isa potentialvapor sourcefortheTibetanregionduringtheMonsoon season. The windfield isaveraged from 100–850 hPa wherethe watervapor isconcentrated.The windvectorsatthetarget loca- tioninsidetheT106resolutiongridboxwereinterpolatedusinga bilinearinterpolationandtrackedbackwardsfortendays.

3. Results

3.1. Simulatedanomaliesoftemperature,precipitationandδ¹⁸Op

We compare the 10-year mean annual surface temperature, precipitationandδ¹⁸Opforthethreetimeslices(PI,MHandLGM), andshowtheresultsasthedifference invaluescomparedtothe PI(e.g.MH-PI,LGM-PI).Hereafter,werefertothesedifferencesas theanomalies (relativeto PI). Astatisticalt-test isappliedto as- sess the significance of the difference between MH, LGMand PI (Figs. 1–3).Areasthatarecoloredshowdifferencesthataresignif- icantatthe95%confidencelevel,thewhiteregions indicateareas wherethedifferencebetweensimulationsisresolvedatthe<95%

confidencelevelandnotinterpreted.

SimulatedPItemperatureandMH-LGMtemperatureanomalies are shownin Fig. 1.Simulated PI temperatures (Fig. 1A) indicate that in the eastern Tibetan Plateau (at longitudes ∼>90^◦E), the meanannualtemperatureis∼−³.0 to 0.5^◦C. Coldertemperature arepredictedforthesouthwestplateau(theeasternHimalaya,Na, Fig. 1) and northwest (the Kunlun mountains, KU, Fig. 1) with mean annual temperatures between ∼−^{10 to} −⁷.0^◦C. The MH temperaturedeviateslittlefromthePI temperatureconditionson

Fig. 1.Simulated annual meantemperatureanomalies. (A):simulatedPI annual meantemperature.(B)and(C):simulatedMHandLGMannualmeantemperature anomalies.Thedifferencesthataresignificantbelowthe95%levelareshadedin white.Theredlinecirclestheregionwhereelevationexceeds1500 m.Theloca- tionsmarkedhereareDu:Dushanbe,Nn:Nangaparbat,Ku:KunlunMountains,Ka:

Kathandu,Lh:Lhasa,Na:NamchabarwaandCh:Chengdu.(Forinterpretationofthe referencestocolorinthisfigurelegend,thereaderisreferredtothewebversionof thisarticle.)

theplateau.Morespecifically,duringtheMHsurfacetemperatures are ∼0.5^◦C colder on thecentral Tibetan Plateau, ∼1.0^◦C colder attheHimalayanfront,and∼^{0.5◦CwarmerthanPI}temperatures over the northern part of the Tibetan Plateau. In contrast, LGM temperaturesaretypicallycoolerovertheregion.DuringtheLGM surface temperaturesare ∼².0–6.0^◦C cooler across the Himalaya and Tibet, and >6.0^◦C cooler in the northwest and southeast edgesofthePlateau.

Simulated PI precipitation (Fig. 2A) shows higher valuesfrom Dushanbe (>2.0 m/yr, Du Fig. 2A) located west of the Tibetan PlateauandcontinuingalongtheHimalayatotheEast.Thesouth- eastPlateaureceiveslargeamountsofprecipitation(1.0–2.5 m/yr), whereasthenorthwestPlateauisdrierwithprecipitationamounts of<0.5 m/yr.PredictedmeanannualMH precipitationanomalies show a generaltrend of wetter conditions in the MH and drier intheLGMcomparedtoPI conditions(Figs. 2B and C).Meanan- nual MH precipitationis ∼<⁰.1 m/yr wetter onthe plateauand 0.3–0.5 m/yrwetteracrosstheHimalayacomparedtothePI.Dur- ingtheLGM,meanannualprecipitationis0.2–0.6 m/yrdrierover the Tibetan Plateau relative to PI conditions(Fig. 2C). Acrossthe westernHimalaya,LGMmeanannualprecipitationis0.4–0.8 m/yr drierthanthePI.AlongtheeasternHimalaya,thepredictedmean annualprecipitationis wetter(>0.8 m/yr)closeto Kathanduand

Fig. 2.Simulated annualmeanprecipitationanomalies. (A):simulatedPIannual meanprecipitation.(B)and (C):simulatedMHand LGMannualmeanprecipita- tionanomalies.Thedifferencesthataresignificantbelowthe95%levelareshaded inwhite.LabelsandtheredlinearethesameasinFig. 1.(Forinterpretationofthe referencestocolorinthisfigurelegend,thereaderisreferredtothewebversionof thisarticle.)

Namcha,butshowsnosignificantchangeinsurroundingareasrel- ativetoPIvalues(whitecolorsinFig. 2C).

PIδ¹⁸O_p(Fig. 3A)ischaracterizedbyadepletionofδ¹⁸O_p from thenorthwest ofthe Tibetan Plateau (∼−¹⁰h^{) to the southeast} oftheTibetanPlateau(∼−¹⁹h^{).Locallowsin}δ¹⁸Op valuesexist inthe Pamir region (−¹⁶h ^to −¹⁹h^{) northwestof the Tibetan} PlateauandnorthofKathmandutoLhasa(−¹⁹h ^to−²²h^{) and} in the southeast plateau. In comparison to PI conditions, simu- latedMH δ¹⁸Op is 1.0h ^{depleted on theplateau and 1}.0h ^en- richedacross,andsouthof,theHimalaya(Fig. 3B).LGM δ¹⁸Op is 2.0–6.0h^{enrichedrelativeto}pre-industrialpredictionsacrossthe plateau,and continuing eastward to the Pacific Ocean. However, northwestof the plateau, LGM δ¹⁸Op values go the opposite di- rectionandaremoredepleted (Fig. 3C).The causesfortheδ¹⁸Op changes will be discussed later in the context of variable vapor source, temperature, precipitation amount, and evaporative recy- cling.

3.2.Seasonalandspatialvariationsofδ¹⁸O_p,temperatureand precipitation

Thespatialandseasonal(winterdefinedasthemonthsDJF,and summerdefinedasJJA)variations in δ¹⁸O_p andclimatewere in- vestigated for the MH and LGM. To simplifythe presentationof

Fig. 3.Simulated annualmeanδ¹⁸Op anomalies.(A): simulatedPI annualmean δ¹⁸Op.(B) and(C):simulatedMHandLGMannualmeanδ¹⁸Opanomalies.Thedif- ferencesthataresignificantbelowthe95%levelareshadedinwhite.Labelsand redlinearethesameasinFig. 1.(Forinterpretationofthereferencestocolorin thisfigurelegend,thereaderisreferredtothewebversionofthisarticle.)

Fig. 4.(1)South,north,westandeastregionsontheTibetanPlateauforinvestigat- ingthevariationofδ¹⁸Op andtheclimateduringtheMHandLGM(blackboxes:

∼³⁵⁰×^{350 km).(2) 4 mountainslopes atthe edgeofthe TibetanPlateaufor} isotopelapseratecalculation(redboxes),regionsincludethesouthslopeofthe Himalaya,northslopeoftheTibetanPlateau,westslopeofthePamir,eastslope oftheTibetanPlateau.(Forinterpretationofthereferencestocolorinthisfigure legend,thereaderisreferredtothewebversionofthisarticle.)

resultswe focuson fourdifferentdomains neartheedges ofthe TibetanPlateau(black boxesinFig. 4).Thesizeofeachdomainis 3×^{3 modelgridcellstoprovidearegional}representation.Thedif- ferentdomainsshowninFig. 4aretheninvestigatedforvariations inδ¹⁸Op,temperature,andprecipitationforeachofthethreetime

Fig. 5.Boxplotsforδ¹⁸Op,temperatureandprecipitationof2seasons(DJFandJJA),4regions(South,north,westandeast)and3timeslices(LGM,MHand PI).Eachbox plotrepresentsthemaximum,75pthpercentile,median,25thpercentileandminimumvalueofthevariables.Thestatisticalcalculationsarebasedonanumberof>270 samples(9 gridboxesineachregion,3months,and>10simulationyears).Winter(DJF)andsummer(JJA)resultsarerepresentedwithredandblueboxesrespectively.The numbersoftheprecipitationamountareaddedtoFig. 5I–L.(Forinterpretationofthereferencestocolorinthisfigurelegend,thereaderisreferredtothewebversionof thisarticle.)

slices(PI, MHandLGM). The resultsprovide thefollowing infor- mation:(1)Seasonalvariations,(2)spatialvariationsbetweeneach region, and(3)paleo to modern variation of δ¹⁸Op, precipitation and temperature. The correlation coefficients of δ¹⁸Op, tempera- ture, andprecipitationare notshownhere, butwillbepresented inthenextSection4.3.

The seasonal andspatial variation ofδ¹⁸O_p, temperature, and precipitationduring the LGM, MHandPI are presentedinFig. 5.

Each bar-and-whisker plot represents the maximum value, 75th percentile,median,25percentile andminimumvalue ofthevari- ablesplotted.Thestatisticalcalculationsarebasedon>270values (9 gridboxesineachregion,3 months,and>10simulationyears).

Winter (DJF) andsummer (JJA) resultsare presented asred and blueboxesrespectively.

Theδ¹⁸Op variationsinFigs. 5A–Dillustrateseveralkeypoints.

First, for all time slices(LGM, MH andPI), the median δ¹⁸Op is

∼².0–5.0h ^{more enriched during winter than summer for the} southandeastregion(Figs. 5Aand D),andis∼^10to−¹³h^more depleted during the winter than summer for the northern and westernregion(Figs. 5Band C).Thesevariationsinδ¹⁸Opindicate that the three δ¹⁸Op distributions zones (southeast zone, transi- tionalzone,andnorthwestzones)identifiedbasedonthemodern dayδ¹⁸O_pvariability(Yaoetal.,2013) alsoexistedduringtheMH andLGM as well. Second, LGMand MH climate caused a larger (∼¹.0–5.0h^{)differenceinmedian}δ¹⁸Op betweenthewinterand summeracrosstheTibetan Plateau thanduring PIclimate.Inthe southregion,forexample,theseasonaldifferencesare ∼⁶.0h ⁱⁿ the LGM, ∼¹⁰h ^{in the MH, and} ∼⁵.0h ^{for the PI} simulations (Fig. 5A). Finally (third),theLGMandMH δ¹⁸Op anomalies show bothregionalandseasonalvariations.Inthesouthernandeastern regions, LGMδ¹⁸Op ismoreenriched than PIvalues forboththe winterandsummer.Thelargestenrichment(∼⁴.0h^)issimulated for the eastern region during the winter (Fig. 5D). In the west- ern andnorthern regions (Figs. 5B and C),the LGMmean δ¹⁸O_p is ∼<0.5h ^{moredepleted inthe winterbut}∼<0.5h ^moreen-

richedinsummer.MHδ¹⁸Op anomaliesshowa complextrendas well.Inthesouthandeastern regions(Figs. 5Aand D),MHδ¹⁸O_p is approximatelythe sameasPI δ¹⁸Op in thewinter,butheavily depleted inthesummer(∼⁴.0–6.0h^{). Inthenorthernandwest-} ern regions,MHmedianδ¹⁸Op is∼².0h ^{moredepletedthanthe} PI medianin winter,and ∼⁰.5–4.0h ^{more enriched than thePI} medianinthesummer.

Seasonaltemperaturevariations(Figs. 5E–H)alsovarybetween the different time slices and regions. The summer–winter tem- perature amplitudeishigherduring the LGMandMH(18–25^◦C) than in the PI (14–16^◦C). For example, in the southern regions thewinter–summertemperaturedifferenceis∼^{20◦CfortheLGM,}

∼^{18◦C forthe MHand} ∼^{14◦C forthe PI (Fig. 5E).}Furthermore, the LGMandMH temperatureanomalies (relativeto thePI) also show more pronounced seasonal variation. The LGM is colder (8.0–10^◦C)in winterbutLGM summertemperatureis aboutthe sameasduringthePItimes.MHtemperatureishigher(2.0–7.0^◦C) in summer but colder (3.0–5.0^◦C) in winter than during the PI times.ThewesternregionischaracterizedbythelargestMHtem- perature differences, i.e.a summerwarming of 7.0^◦C andwinter coolingof5.0^◦C(Fig. 5G).

Seasonal precipitationvariations (Figs. 5I–L)demonstratethat:

first,astronginfluenceofthemonsooninthesouthernandeast- ern regions for the LGM and MH simulations. For example, the southernregionLGMprecipitationincreasesfrom18.1 mm/month in the winter to 288 mm/monthin the summer (Fig. 5I). In the eastern region,precipitationincreasesfrom6.5 mm/monthinthe winter to 180 mm/monthin summer (Fig. 5L). Second, the LGM andMH precipitation anomaliesalso show differentregionaland seasonal variation. The LGM winter is drier than the PI winter across the Tibetan Plateau exceptin thesouth region. Forexam- ple,LGMwinterprecipitationis6.0 mm/monthhigherthanduring PItimes(18.1 vs12.1 mm/month)(Fig. 5I).LGMsummersarewet- ter in the southern (LGM: 288 mm/month, PI: 243 mm/month) andnorthernregions(LGM:58.8 mm/month,PI:28.0 mm/month)

Fig. 6.ModelsimulatedisotopelapserateduringthePI(toprow),MH(middlerow)andLGM(bottomrow).Theanalysishasbeenconductedat4mountainslopes:thesouth slopeoftheHimalaya(25^◦N–30^◦N,84^◦E–92^◦E)(redcycle),westslopeofthePamir(36^◦N–38^◦N,65^◦E–70^◦E)(orangesquare),northslopeoftheTibetanPlateau(36^◦N–39^◦N, 85^◦E–89^◦E)(bluetriangle)andtheeastslopeoftheTibetanPlateau(32^◦N–36^◦N,100^◦E–110^◦E)(greendiamond).Theprecipitationweightedtemporaldifferencesare presentedfortheannualmean(leftcolumn),winterseason(DJF)(middlecolumn)andmonsoonseason(JJA)(rightcolumn).Theuncertaintyofslopesandtheirsignificant levelarelistedinTable S3inthesupplementary material.(Forinterpretationofthereferencestocolorinthisfigurelegend,thereaderisreferredtothewebversionofthis article.)

(Figs. 5I and J). MH summer precipitation is lower than the PI precipitationin the western region (MH: 28.4 mm/month, vsPI:

59.0 mm/month),andMHwinterprecipitationis22.3 mm/month butPIprecipitationis25.5 mm/month(Fig. 5K).

Taken together, the above results suggest: (1) LGM and MH bothhavehigherseasonal(winter–summer) amplitudesinδ¹⁸O_p, temperatureandprecipitationvaluesthanduringPI times.Ofthe four regions analyzed, these seasonal differences are most pro- nounced in the western region. (2) Based on the climatological meananomaliesdiscussedintheprevioussection,theLGMisgen- erallycolder anddrierandthe MHiscolder andwetter thanthe PIclimate.(3) Thethreeδ¹⁸Op distributionszonesobservedinthe modern(Yaoetal.,2013) alsoexistedduringtheMHandLGMas well.

3.3.SpatialandtemporalvariationsinisotopelapseratesduringthePI, MHandLGM

Inthisstudy,spatialandtemporalvariationsofthemodelsim- ulatedisotopic lapse rate havebeen analyzed at four highrelief zoneson the flanks of theplateau (red boxes in Fig. 4). The lo- cations include the southern slope of the Himalaya(25^◦N–30^◦N, 84^◦E–92^◦E),western slopeofthePamir(36^◦N–38^◦N,65^◦E–70^◦E), northern slope of the Tibetan Plateau (36^◦N–39^◦N, 85^◦E–89^◦E), and the eastern slope of the Tibetan Plateau (32^◦N–36^◦N, 100^◦E–110^◦E). Theanalysis hasbeen done forthe winterseason (DJF), monsoon season (JJA), and precipitation weighted annual mean (Fig. 6). The statistical analysis of the estimated isotopic

lapseratesincludinguncertaintyandthesignificantlevelisshown inTable S3insupplementary material.

Significant temporalchanges inthe isotopiclapserateare ev- ident fromthe model results (Fig. 6). The primary trends in the model results include: (1) some regional variations in the mean annualandseasonalisotopelapserates.Forexample,thesouthern Himalayanfronthasasimilarlapserate(∼³.0h/km)astheeast- ern TibetanPlateau (redlinesoverlayorare paralleltothe green lines), but the northern slope of the Plateau (blue lines) has a lowerisotopelapseratethantheslopesinthesouthernandeast- ernregions(e.g.Figs. 6A,C,D,G, I).ThewesternPamirregionhas a higher isotope lapse ratein the winter than the other regions (Figs. 6B,E, H).Forexample,theMHwinterwestPamirlapserate is −⁴.2h/km, whereas the southern andeastern lapse ratesare

∼−³.2h/km andthenorthernlapserateis−⁰.7h/km (Fig. 6E).

(2) Different regions havedifferent seasonal variationin the iso- topelapserate. Thesouthern andeastern slopes havethelowest seasonalvariation(∼⁰.2h/^km,max0.4h/^{km)(compareFigs. 6B,} E,HwithFigs. 6C,F, I).The summerisotopelapseratecalculated for the west slope of the Tibetan Plateau (Figs. 6C, F, I) are not significant at the 95% significance level, which is due to model limitationsfor calculatinglapserate whenprecipitation ratesare verylow.Finally,(3) theprecipitationweightedannualmeaniso- topelapseratesinthethreetimeslices(PI,MHandLGM)indicate the smallest differences in the lapserate (∼⁰.1h/^{km) occur on} the eastern slope betweeneach time slice. The largestdifference (∼¹.4h/^{km)isfoundinthewesternslope wherethelargestun-} certainty is (Table S3). The northern slope has a mean annual

difference of ∼⁰.4h/km between the LGM, MH and PI whereas forthe southern slope, the isotope lapserateis about 0.4h/km largerduringtheMHand0.2h/km smallerduringtheLGMthan duringthePI(Figs. 6A,Dand G).

4. Discussion

4.1. Comparisonsbetweenthesimulatedanomaliesandthe proxy-basedreconstructionsduringtheMHandLGM

The simulated anomalies have been compared with climate changereconstructionsbasedonproxydata(seesupplementalTa- bles S1, S2 for dataused). Inmost cases climatereconstructions fromproxydataindicatea changerelative toa knowncondition, such as the PI times, and do not provide an absolute value or magnitudeofthechanges.Asaresult,thechangescanonlybede- scribedhereinrelativeterms,i.e.as‘wetterthan’,‘similarto’and

‘drier than’ for precipitation, and ‘warmer than’, ‘similar to’ and

‘colder than’for temperature. The PI climate serves asour refer- enceclimateforthesedescriptions.Theproxydatabasedchanges are plottedover the modelsimulated precipitationandtempera- tureanomaliesfortheMH(Fig. 7A)andLGM(Figs. 7Band C).

There is both an agreement and disagreement between the modelandobservationsduringtheMH(Fig. 7A).Modelandobser- vationsagreeattheHimalayanfrontthattheMHclimateiswetter thanthePI climate.OnthenortheastpartoftheTibetanPlateau, MH proxy data show a large disagreement between neighboring locations, e.g.atlocations2–6,9–10,12,18–19, 22–23and25in Fig. 7A, whereas the model shows no significant change at the 95% confidencelevel (represented by the whitecolorin Fig. 7A).

Themodel resultsfromGCMssuggest that theclimatechange in thisregion ismuted and the variations in the proxy datamight resultfromdifferentmethods(Table S1)andlocalizedclimatesig- nals.Theaboveresultsagreewithpreviouswork(Dallmeyeretal., 2013) that conducteda similar comparisonbetweentheMH and PI climates.Theirresultsshow a wetterMH climateatthesouth andmiddleTibetanPlateau,andsuggestedaregionaldissimilarity intheannualprecipitationsignalbetweentheMHandPIclimate aswell.

Thereisbetteragreementbetweenmodelandobservationsfor the LGM than for the MH (Figs. 7B and C). Predicted tempera- turechanges show good agreement withthe proxy datain both thesignandmagnitudeofthepredictedchange(Fig. 7B).Schmidt etal. (2011) derived an LGMmax T of 3.0to 4.0^◦C frombio- geographicalandphylogeneticdataofbiomarkersatlocation6in Fig. 8AwherethemodelsimulatedLGMtemperaturedifferenceis

−³.7^◦C.Herzschuhetal. (2006)deriveda∼⁴.0–7.0^◦C coldertem- peratureduringtheLGMfromthevegetationofasedimentcoreat theQilianMountainswherethemodelsimulatedresultis−⁴.7^◦C.

The predicted precipitation differences also show a good agree- mentwiththeproxydata(Fig. 7C).ThesignofthesimulatedLGM precipitationchangesagreeswiththeexistingproxydatanearthe Himalaya(locations 2,3inFig. 7C),attheZabuyesaltlakeonthe centralTibetan Plateau(location 5),atQilianMt(location 7),and atXinyunlakeintheSouthChina(locations9, 10).Modelpredic- tions alsoindicateboth wetter,drier, andsimilarto presentLGM precipitation condition for the central Himalaya (red, blue, and whitecolorsinFig. 7C),andoneproxystudiesinthisregion(loca- tion 1)suggestprecipitationconditionsthesameastoday(Asahi, 2010).

Unfortunately,a moredetailedcomparisonbetweenGCMpre- dictionsandproxydatacannotbecompleted.Thisismainlydueto thefactthatmanyoftheavailableproxydataindicateonlyarela- tivechangeinconditions,suchthatotherparameterslikechanges in the maximum andminimum predictedand observed temper- ature or precipitation cannot be compared. No better or worse

Fig. 7.(A)MHmodelandproxydatacomparison.Locationsandreferencescitedon theplot:1. Namco,Muegleretal. (2010).2. LakeKuhai,Wischnewskietal. (2011).

3. LakeKoucha,Mischke etal. (2008). 4. LakeCuoe,Wu etal. (2006).5. Ahung co, Morrillet al. (2006).6. LakeZabuyeSalt,Wangetal. (2002). 7. KaiduRiver, Wuennemannetal. (2006).8. SelinCo,Zhangetal. (2011).9. LakeQinghai,Shen et al. (2005).10. LakeQinghai,Liuet al. (2014).11. Mt.Qilian,Herzschuhetal.

(2005).12. LakeXimencuo,ZhangandMischke (2009).13. PumoyumCo,Nishimura etal. (2014).14. QaidamBasin,YuandLai (2014).15. SumxiCo,Gasseetal. (1991).

16. LakeBangong,VanCampoetal. (1996).17. LakeTsoKar,Demskeetal. (2009).

18. LakeNaleng,Krameretal. (2010),Wischnewskietal. (2011).19. LakeZigetang, Herzschuhetal. (2006).20. LakeQiliu,Hodellet al. (1999),Zhanget al. (2011).

21. LakeXingYun,Hodelletal. (1999),Zhangetal. (2011).22. LakeYidun,Shen et al. (2006), Wischnewskietal. (2011).23. RenCo,Tangetal. (2000).24. Lake Hidden, Tangetal. (2000).25. Hongyun,Zhouet al. (2002), Wischnewskiet al.

(2011).26. DundeIceCore,Liuetal. (1998),Wischnewskietal. (2011).LGMtem- perature(B)andprecipitation(C)modelandproxydatacomparison.Locationsand referencescitedontheplot:1. EasternNepalHimalaya,Asahi (2010).2. Kumaun Himalaya,Kotliaetal. (2010).3. KumaunHimalaya,Kotliaetal. (2000).5. Zabuye Saltlake,Wangetal. (2002).6. YarlungZhangbocatchment,Schmidtetal. (2011).

7. West Qilian Mt., Huet al. (2014). 8. Qilian Mt.,Herzschuh et al. (2006). 9.

10. XinyunLake,Hodelletal. (1999).(Forinterpretationofthereferencestocolor inthisfigurelegend,thereaderisreferredtothewebversionofthisarticle.)

agreement betweencertain typesofproxydataandthe GCMre- sults can be determined withthe available data.Agreement and disagreementcanbeseenamongthesametypeofproxydatawith theGCMresultsfortheMH(Table S1).Wenote howeverthatthe agreement betweenthemodelanddataforthe LGMisgenerally good, despitethe fact that the eightLGM proxystudies usedin-

Fig. 8.EIMR(Goswamietal.,1999) andWYI(WebsterandYang,1992) monsoon indexcalculatedfromthesimulatedlong-termmeandataduringthePD,PI,MH andLGM.ThemeanandstandardderivationofsimulatedyearlyEIMRandWYIare listedinTable S4insupplementmaterial.

cludedsix differentproxymethods (Table S2). From ouranalysis wecanhoweverconcludethat:(1) MHclimatechangerelativeto today is muted and the climate anomalies might lie in the un- certainty of both predicted GCM values and proxy data. Finally, (2) paleoclimatechangefortheLGMisconsistentlyreconstructed byboththeGCMandproxystudies.

4.2.LengthandintensityofIndianmonsoonduringthePI,MHandLGM

PrecipitationovertheTibetanPlateauisstronglyinfluencedby the Indian monsoon. The intensity and the period of the mon- soonsystemcanbeidentifiedwithvariousmethods.Inthisstudy, weused the ExtendedIndian Monsoonrainfall index(EIMR) and Wesbter and Yang Index (WYI) to identify the variations of the monsoonintensityduring the MHandLGM (Fig. 8). The EMIRis definedasthesummer(JJAS)rainfallinthecontinentalregionbe- tween 70^◦E–110^◦E and10^◦N–30^◦N (Goswami et al., 1999). WYI

is definedasthe difference inzonal wind speed at850 hPa and 200 hPa for the region 40^◦E–110^◦E and 0^◦–20^◦N (Webster and Yang, 1992). Using theGCMpredictedwind field andabulk tra- jectoryanalysismethod,thesourceofthevaporandthelengthof the monsoonperiod during the LGM, MHand PI aredetermined inthisstudy(Fig. 9).Thebulktrajectoryanalysisisthebacktrack- ingofthemoistureonlyintheboundarylayer(1000to850 hPa) wherethevaporisconcentrated(Fengetal.,2013).

Both the EIMR and WYI monsoon index suggest significantly reduced monsoon intensity during the LGM (Fig. 8). Relative to thePI climate,theLGMhada44.5%decreaseinEIMR(decreased from304 mm/monthduring thePIto 168 mm/monthduring the LGM),anda 27.6%decreaseinWYI(decrease from25.4 m/sdur- ing the PI to the 18.4 m/s during the LGM) in comparisonto PI conditions.The changes ofthe MH monsoonintensity relativeto the PI are not consistent betweenboth calculation methods but thosechangesareminor(Fig. 8).Forexample,theEIMRincreases from303 mm/monthduring PIto314 mm/monthduring theMH (∼^{3.6% increase),while the WYIdecreases from25.4 m/s during} thePI to24.2 m/sduring theMH(∼^{4.7%decrease).Thestandard} deviation(

σ

σ

σ

Theinter-annualbulktrajectoryanalysisshowsthatthelength oftheLGMmonsoonisaboutonemonthshorterthanduring the PI,andthattherearenoobviousdifferencesinthemonsoonlength between the MH and PI (Fig. 9). During the winter (NDJFMA) (Figs. 9A,F andKshow anexample forJanuary),the dominantly westerly winds demonstrate that no monsoon exits during this time. The onset of the monsoon is marked by the beginning of themoistureoriginatingfromtheIndianOcean,whichhappensin Mayforallthreetimeslices(Figs. 9B,Gand L).Thedeclineofthe monsoon initiated inOctoberduring thePI (Fig. 9E)andthe MH (Fig. 9J),butoccursamonthearlier,inSeptember,duringtheLGM (Fig. 9N).

Fig. 9.Interannualbulktrajectoryanalysisat25^◦N,87.5^◦E.Ittracksthebulkmoistureintheboundarylayer(1000to850 hPa)wherethevaporisconcentrated(Fengetal., 2013).ThelocationliesatthesouthoftheHimalayanfront,whichisapotentialvaporsourceofprecipitationontheTibetanPlateauduringthemonsoonseason.

Fig. 10.RDMandGCMcomparisonandrelatedatmosphericprocessesforthemonsoonseason(JJA)duringthePI,MHandLGMatHimalayafront(zonalaveragebetween 84^◦E–92^◦E).(A–C)RDMandGCMδ¹⁸Opcomparison.Modelδ¹⁸Oparepresentedassoliddotsandblacksolidline(trendlinefromalinearfitting).Theredtrianglerepresents theRDMestimationinitiatedwithGCMderivedmoisturesource.Theinitialvaporsourcearepresentedasthreevariables:Ts,qsandδ¹⁸Ov(listedinthedashbox),where Tsistheinitialvaportemperature,qsistheinitialvaporspecifichumidityandδ¹⁸Ovistheδ¹⁸Ointheinitialwatervapor.(D–F)Modelsimulatedspecifichumidityatcross section10^◦N–40^◦N.Colorscontourrepresentqorqanomalyvalues.Streamlinesrepresenttheverticalwind(omega)andthemeridionalwind(v)field.TheTibetanPlateau isshadedingray.(H–J)Estimatedδ¹⁸Ovmixingrateatcrosssection10^◦N–40^◦N,dataarezonallyaveragedbetween84^◦E–92^◦EduringJJA.Themixingrateisdefinedas thechangeofδ¹⁸Opduetothemixingeffectofwatervaporwithitsenvironment(units:h/^h).Streamlinesrepresentverticalwind(omega)andmeridionalwind(v)field.

TheTibetanPlateauisshadedingray.(K–M)Changesinδ¹⁸Opifoneprecipitationtypeisabsent.Purplerepresentsthechangesinδ¹⁸Opiftheconvectiveprecipitationis absent,orangerepresentschangesinδ¹⁸Opifthesnowtypeisabsent.TheelevationoftheTibetanPlateauisplottedasblacksolidline.(Forinterpretationofthereferences tocolorinthisfigurelegend,thereaderisreferredtothewebversionofthisarticle.)

4.3. Correlationsbetweentemperatureandprecipitationwithδ¹⁸O_p

Temperature,topographyandprecipitationarethemostimpor- tantclimatecontrolsonδ¹⁸Op.Inthisstudy,we correlatemodel- predictedδ¹⁸Op withtemperatureandprecipitation,andcorrelate temperature with precipitation (Fig. S1). The seasonal cycle has beenremovedfromthemodeloutputtofilterouttheseasonalsig- nalintheanalysis(toavoidhighcorrelationcoefficientsproduced by season-dependent co-variability of the investigated variables).

The correlation analyses are conducted foreach grid box witha populationofmorethan120datapoints.Onlythecorrelationco- efficients that are above the 95% significance level are shown in theplot (Fig. S1). The correlation analyses indicates(1) a signifi- cantpositive correlation ofδ¹⁸O_p andtemperature inthe region

>∼^{30◦N,andanegative}correlationbetweenδ¹⁸O_p andprecipita- tionintheregion<40^◦N.(2) Thedominantcontroloftemperature onδ¹⁸Op atthe highaltitudesandprecipitation onδ¹⁸Op atlow altitudes maybe explainedby an overshadowing of theeffectof temperature in latitudes governed by higher precipitation rates (theSouth).(3) OntheTibetanPlateau,thetemperaturecontrolis insignificant,possiblybecausethevariabilityintemperaturethatis relevanttoδ¹⁸Opvariabilityliesmostlyintheseasonalcyclewhich hasbeenremoved.(4) Thepatternofcorrelationcoefficientsdoes not varymuch betweenthe differenttime slices, suggestingthat

thedifferentclimatesduringtheLGM,MHandPIhavelittleorno influenceonthespatialdependenceofthecontroloftemperature andprecipitationonδ¹⁸Opinthisregion.

4.4. GCMandRDMapproachesforestimatingδ¹⁸OpduringthePI,MH andLGM

The comparison between the more simple one-dimensional RDM and more comprehensive GCM predicted δ¹⁸Op was con- ducted on the south slope of the Himalaya (84^◦E–92^◦E, 23^◦N–30^◦N) during the monsoon season JJA. The initial vapor moistureconditions(representedasTs,qs andδ¹⁸Ov)usedforthe RDMcalculationarederivedfromtheGCMoutput,whereTsisthe initial vaportemperature, qs isthe initialvapor specifichumidity andδ¹⁸O_vistheδ¹⁸Ointheinitialwatervapor.

Results show that the RDM agrees with the GCM during the PI (Fig. 10A)andLGMtimes(Fig. 10C). DuringtheMH,the RDM predicts more depleted δ¹⁸Op at ∼^{1 km elevation and} ∼².0h more enriched δ¹⁸O_p than the GCM prediction at ∼^{5 km eleva-} tion (Fig. 10B). The differencesin theRDM andGCM predictions duringtheMHcanbeexplainedinpartbywarmerconditionsand more vaporcontent intheinitial moisturesource during theMH than that during the PI for the RDM calculations(dash boxes in Figs. 10Aand B,andFigs. 10Dand E).