Mass Balance of Novaya Zemlya Archipelago, Russian High Arctic, using Time-Variable Gravity

25 Pages·2017·12.61 MB·English

Preprints (www.preprints.org) | NOT PEER-REVIEWED | Posted: 3 October 2018 doi:10.20944/preprints201810.0056.v1
Peer-reviewed version available at Remote Sens. 2018, 10, 1817; doi:10.3390/rs10111817
Article
Mass Balance of Novaya Zemlya Archipelago,
Russian High Arctic, using Time-Variable Gravity
from GRACE and Altimetry Data from ICESat and
Cryosat-2
EnricoCiracì1,†,‡ ID,IsabellaVelicogna1,2‡ ID,TylerClarkSutterley3‡ ID
1 DepartmentofEarthSystemScience,UniversityofCalifornia,Irvine,CA92697,USA;
2 JetPropulsionLaboratory,CaliforniaInstituteofTechnology,Pasadena,CA91109,USA;
3 UniversitiesSpaceResearchAssociation/NASAGoddardSpaceFlightCenter;
[email protected](I.V.);[email protected](E.C.);
* Correspondence:[email protected];Tel.:[phone]
Abstract: WeexaminethemassbalanceoftheglaciersintheNovayaZemlyaArchipelago,located
intheRussianHighArcticusingtimeseriesoftime-variablegravityfromtheNASA/DLRGravity
RecoveryandClimateExperiment(GRACE)mission,laseraltimetrydatafromtheNASAIceCloud
and land Elevation Satellite (ICESat) mission, and radar altimetry data from the ESA CryoSat-2
mission. Wepresentanewalgorithmfordetectingchangesinglacierelevationfromthesesatellite
altimetrydataandevaluateitsperformanceinthecaseNovayaZemlyabycomparingtheresultswith
GRACE.WefindthatthemasslossofNovayaZemlyaincreasedfrom10±5Gt/yrover2003-2009to
14±4Gt/yrover2010-2016,withabriefperiodofnearmassbalancebetween2009and2011. The
resultsareconsistentacrossthegravimetricandaltimetricmethods. Furthermore,theanalysisof
elevationchangefromCryoSat-2indicatesthat60%ofthemasslossoccursatlowelevation,where
thinning rates are highest. We also find that marine-terminating glaciers in Novaya Zemlya are
thinningsignificantlyfasterthanland-terminatingglaciers,whichindicatesanimportantroleofice
dynamicsofmarine-terminatingglaciers. Wepositthattheglacierchangeshavebeencausedby
changesinatmosphericandoceantemperatures. Wefindthattheincreaseinmasslossafter2010is
associatedwithawarminginairtemperatures,whichincreasedthesurfacemeltrates. Thereisno
enoughinformationontheoceantemperatureatthefrontoftheglacierstoconcludeontheroleof
theocean,butwepositthatthetemperatureofsubsurfaceoceanwatersmusthaveincreasedduring
theobservationperiod.
Keywords: NovayaZemlya;Altimetry;Gravity;RussianHighArctic;GlaciersandIceCaps;Remote
Sensing;ClimateChange.
1. Introduction
Inthetwentiethcentury,withamarkedincreaseinglobalatmospherictemperaturesattheendof
theLittleIceAge,themassbalanceofmountainglaciersandicecaps(GIC)hasbeenpredominantly
negative[1–4]. Availableobservationssuggestthattheshrinkageprocesscontinuedduringthefirst
decadesofthe21stcentury[4,5]. Recently,Reageretal.[6]estimatedthatglaciersoutsideGreenland
andAntarcticahavebeenlosingmassatarateof253±35Gt/yrbetween2002and2014. Thisresult
placesGICasamajorcontributortopresentandnear-termeustaticsealevelrise(SLR).Morethan70%
ofthetotalglaciermasslosstakesplaceintheArctic[7],wherealbedofeedbackcausedbydecreasing
seaiceandsnowcoveryieldswarmingrates2-3timeslargerthantheglobalaverage[8–10]. This
studyfocusesontheglaciersoftheNovayaZemlyaArchipelago(NZEM),intheRussianHighArctic
(RHA).TheRHAhaslosticemassatarate9.9±8Gt/yrbetween2002and2014, withmorethan
80% of the mass wastage attributed to the glaciers of Novaya Zemlya [6,7,11,12]. Recent studies
© 2018 by the author(s). Distributed under a Creative Commons CC BY license.
Preprints (www.preprints.org) | NOT PEER-REVIEWED | Posted: 3 October 2018 doi:10.20944/preprints201810.0056.v1
Peer-reviewed version available at Remote Sens. 2018, 10, 1817; doi:10.3390/rs10111817
2of25
suggestthatbothnegativesurfacemassbalance(SMB)andacceleratingicedischarge(D)increased
theglaciermasslossafteryear2000[13–15]. TheseresultsdifferfromotherArcticregions,e.g. the
Canadian Arctic Archipelago [16–19] or Alaska [20], where SMB has been identified as the main
processdrivingthemassloss. Yet,adefinitiveassessmentofthepartitioningofthemasslossbetween
DandSMBishoweverlackingfortheNZEM.Atmospherictemperatureswilllikelycontinuerisingin
thenextcentury[5]. ModelprojectionsindicatethatRHAwilllosealmost50%ofitstotalvolume,
correspondingtoanequivalent24±8mmofSLR,by2100[21,22],whichsuggeststhattheregional
masslosswillcontinuetoincreaseinthecomingdecades. Modelprojections, however, arebased
on surface mass balance models that do not include ice dynamics and hence may underestimates
the mass loss. Additional observations will be necessary to determine the processes dictating the
present-daymasschangeandtobetterconstrainprojectionsoftheevolutionoftheRHAglaciers.
Here, weevaluatethemassbalanceofNZEMusingtimeseriesoftime-variablegravityfrom
the NASA/DLR GRACE mission between April 2002 and August 2016. We compare the GRACE
results with independent estimates obtained from satellite altimetry data acquired by the NASA
ICESatmissionbetweenOctober2003andOctober2009andtheESACryoSat-2missionbetweenJuly
2010andJuly2016. Wepresentanovelalgorithmformappingiceelevationchangeusingaltimetry
data designed to optimize the spatial coverage of measurements over the region of interest. We
calculateregionalmeanelevationchangesusingCryoSat-2andcomparethemwithGRACEduring
theoverlappingtimeperiodforevaluation. Wecomparethinningratesoverland-terminatingversus
marine-terminatingglacierstodeterminethemaincontrolonthemassloss. Wediscusstheinterannual
variabilityinmasschangeanditsrelationshipwithclimate,mostlyprecipitationandairtemperature,
togaininsightintotheprocessesdrivingtheglaciermassloss. Weconcludeontheevolutionofthe
regionalmasslossofNZEMinthepastdecades.
1.1. Studyregion
TheRHAincludes51,592km2 ofglaciersandicecaps,whichis13%ofthetotalglaciatedarea
outside of Greenland and Antarctica. The glaciated areas of the RHA are distributed across four
regions: NovayaZemlya(22,379km2),SevernayaZemlya(16,382km2),FranzJosefLand(12,756km2),
andUshakonIsland(359km2). NZEM,liesbetweentheBarentsandKaraSeasandconsistsoftwo
majorislands,SevernyIslandandYuzhnyIsland,alongwithseveralsmallerislands. Severnyand
YuzhnyIslandsareseparatedbytheMatochkinStrait,a600meterwidechannelthatstretches100km
andconnectstheBarentsandKaraSeas. Anaxialmountainrangeextendsovertheentirelengthofthe
NZEMwithamaximumelevationof1,340monYuzhnyIslandand1,596monSevernyIsland[22].
92%ofthetotalglacierarea(20,784.4km2)islocatedonSevernyIslandandconcentratedinasingle
largeicecap(NorthernIcefield)thatcovers45%oftheislandandconstitutesthelargestbodyoficein
theEurasianContinent[23,24](Fig. 1).
TheNorthernIcefieldismorethan400kmlong,withameaneast-westwidthof80km[24]. The
ice cap has an average elevation of 800 m and an average thickness of 400-450 m [25]. 71% of the
NorthernIcefieldarea(16,064km2)isdrainedbymarine-terminatingglacierswhilemostoftheglaciers
separatedfromthemainicecapterminateonland. Glacierswereclassifiedascold-based[26];but
recentstudiessuggestawiderrangeofthermalregimeandatransitiontowarm-based/polythermal
glaciers[22,27]. Theclimateintheregionisdeterminedbyacomplexinteractionbetweenmultiple
oceanicandatmosphericforcings. OntheNorth-West,mildtemperaturesarefavoredbytheadvection
of warm and salty water transported by the North Cape current from the Atlantic Ocean [28–30].
Abundantannualprecipitation(400mm/yr)isfavoredbytheinfluxofairmassesrichinmoistureby
theAtlanticCyclone[28,30]. Atmosphericconditionsbecomegraduallydrierandcoldertowardthe
Southeast. Thecentralmountainchainprovidesanorographicbarrierfortheeastwardpenetration
oftheAtlanticcyclonicsystem[28]andprecipitationgraduallydecreasestowardtheSouth. Onthe
KaraSeacoast,annualprecipitationamountsto250mm/yr. Here,loweratmospherictemperatures,
Preprints (www.preprints.org) | NOT PEER-REVIEWED | Posted: 3 October 2018 doi:10.20944/preprints201810.0056.v1
Peer-reviewed version available at Remote Sens. 2018, 10, 1817; doi:10.3390/rs10111817
3of25
especiallyinWinterandSpring,arecontrolledbysemi-permanentseaiceandcoldArcticwater[28].
JanuaryandFebruaryarethecoldestmonthsoftheyearwithminimumairtemperaturesdecreasing
fromtheNorth-West(-15◦C)totheSouth-East(-21◦C).Maximumtemperature(+2◦C)isrecordedin
August,withnosignificantdifferencebetweenthetwocoasts. Maximumprecipitationisbetween
SeptemberandOctober. AprilandMayarethedriestmonthsoftheyear.
Figure1.(a)NovayaZemlyaglaciers.Marine-terminatingglaciersarered(KaraSea)andblue(Barents
Sea).Landterminatingglaciersaregreen(KaraSea)andyellow(BarentsSea).Minorglaciersoutside
theNorthernIceFieldareorange.GlaciersoutlinesarefromRastneretal.[24].(b)LocationofNovaya
Zemlya.
1.2. Priorstudies
NZEM’s present-day average temperatures are 3◦C higher than in the 19th century [31].
Atmosphericwarmingthroughoutthefirstdecadesofthe20thcenturycontributedtoglaciershrinkage
[32],duringwhichglaciersintheregionlostmasswithanaverageratebetween0.16and0.30meters
ofwaterequivalentperyear[3,13,22,33,34]. Theprocesswasnotuniformintime,butconsistedofa
warmfirstphasebetweenthe1920’sandthe1960’swithextremelyfastretreatrates(greaterthan300
m/yr),followedbyasecondphasebetweenthe1960’sandthe1990’scharacterizedbylessnegative
Preprints (www.preprints.org) | NOT PEER-REVIEWED | Posted: 3 October 2018 doi:10.20944/preprints201810.0056.v1
Peer-reviewed version available at Remote Sens. 2018, 10, 1817; doi:10.3390/rs10111817
4of25
mass balance and lower retreat rates (50 –150 m/yr). [22,32]. Zeeberg and Forman [32] observed
thatthedecadalvariabilityinmassbalanceonNZEMislinkedtolong-termshiftingofatmospheric
andoceaniccirculationpatternscausedbytheNorthAtlanticOscillation(NAO).Duringapositive
phaseoftheNAO,sea-levelpressureathighlatitudesisbelownormal,whileabovenormalconditions
areregisteredatmid-latitudes. AnincreasedpressuregradientovertheNorthAtlanticcontributes
to an intensified and more zonally oriented jet-stream. When the NAO is in a negative phase, a
decreasedpressuregradientbetweenhighandmid-latitudesweakensthejet-streamwhichbecomes
moremeridionallyoriented[35]. WhentheNAOisinaprolongedpositivephase,thestraightening
ofthepolarvortexcontributestoanincreasedfluxofairmassesfromtheAtlanticOcean,resulting
inelevatedwinterprecipitationoverNZEM.AstrongerNorthAtlanticpressuregradientintensifies
theadvectionofwarmwaterfromtheAtlanticintotheBarentsSea[32]. TheAtlanticWater(AW)
releasesheatintotheatmosphereandyieldswarmerthannormalwinterandsummertemperatures.
DuringextendednegativeNAOphases,thejetstreamweakensandatmosphericconditionsinthe
regionbecomedrierandcolder. Betweenthe1920’sand1950’s,theNAOwaspredominantlypositive
and higherthan normal atmospheric andocean temperatures favored fast glacial recession. After
1960,theNAOturnedtonegative,i.e. colderatmosphericconditionsandloweroceantemperatures,
which slowed down the glacial retreat [32]. Several studies have shown that the NZEM has been
experiencingastateofprolongednegativemassbalance[6,11,22,36]. Moholdtetal.[11]usedsatellite
laseraltimetryfromICESatandgravitydatafromGRACEbetweenOctober2003andOctober2009
toestimateaglacierlossof7.6±1.2to5.8±3.0Gt/yr,withrapidglacierthinningatelevationsbelow
500 m (0.9 m/yr). The analysis of elevation change, however, revealed no significant difference
betweenfrontalthinningofmarine-andland-terminatingglaciers(0.94mvs0.89m),butsignificant
differencesbetweenglaciersontheNorth-West(BarentsSea)(0.45m/yr)vsSouth-East(KaraSea)(0.25
m/yr),consistentwithSMBdominanceobservedinotherArcticregions[16–19]. Zhaoetal.[37]used
multipleactiveandpassivemicrowavedatabetween1995and2011toreportaprogressiveincrease
inmeltseasonlength,consistentwithincreasedatmospherictemperatures. Increasedsnowmeltand
summerrunoffcontributedtothehighicethinningratesobservedbyMoholdtetal.[11]. Yet,new
studies suggest that increased ice discharge from marine terminating glaciers may be significant
[13,14,36,38]. Carr et al. [14] documented NZEM retreat between 1973 and 2015 for both marine-
andland-terminatingglaciers. Marine-terminatingglaciersretreated3.5fasterthanland-terminating
glaciers (46.9 vs 13.8 m/yr) and faster in the Barents Sea vs Kara Sea (55.9 vs 37.2 m/yr). Glacier
recessionacceleratedafteryear2000ontheBarentsSeaandafter2003ontheKaraSea[36,38]. Marine
glaciersretreatratespeakedat85.4m/yrfortheBarentsSeaand64.8m/yrfortheKaraSea. These
retreatscoincidewithwarmeratmosphericandoceantemperaturesandlowerseaiceconcentrations
comparedtothepreviousdecades. TheseresultswereevaluatedbyMelkonianetal.[13]whomapped
icevelocitybetween1999and2014andelevationchangebetween1952and2013/2014tofindthe
highesticespeedsandthinningratesatthefrontofmarine-terminatingglaciersfacingtheBarents
Sea. Melkonianetal.[13]documentedanaccelerationinicevelocityandthinningratesafteryear2000.
Steadyincreasesinmarine-terminatingglaciervelocitywereobservedbySunetal.[39]between2002
and2015andbyStrozzietal.[15]between2008and2016. Allthesefindingsindicateanincreasingly
negativeregionalglaciermassbalanceaftertheyear2000.
2. Dataandmethods
2.1. Time-VariableGravity
TheGravityRecoveryandClimateExperiment(GRACE)isajointoperationbyNASAandthe
German Aerospace Center (DLR) to measure changes in the Earth gravity field at a monthly time
scale[40]. Here,weuse156monthlyGRACERelease-5(RL05)gravityfieldsolutionsprovidedby
theCenterforSpaceResearchattheUniversityofTexas(CSR)fortheperiodApril2002–August2016
[40,41]. MonthlygravityfieldanomaliesfromCSRaredistributedintheformoffully-normalized
Preprints (www.preprints.org) | NOT PEER-REVIEWED | Posted: 3 October 2018 doi:10.20944/preprints201810.0056.v1
Peer-reviewed version available at Remote Sens. 2018, 10, 1817; doi:10.3390/rs10111817
5of25
sphericalharmonicsuptodegreeandorder60. WereplacetheC coefficientswithmoreaccurate
20
estimatesfromthesatellitelaserrangingmissions(SLR)[42]. Weaccountforthevariationofdegree-1
using coefficients calculated from a combination of GRACE coefficients and ocean model outputs
followingSwensonetal.[43]. Weincludeanadditionalpoletidecorrectiontoremovethelong-period
poletidesignalsnotincludedinthestandardcorrections[44]. WeremovetheeffectsfromGlacial
IsostaticAdjustment(GIA),i.e. theviscoelasticresponsebythesolidEarthfromchangesiniceloading
sincethelastglacialmaximum(LGM),usingcoefficientsfromAetal.[45]. Toisolategravitychanges
relatedtoicemasschange, wesubtractthesignalcontributiongeneratedbychangesinterrestrial
hydrologyfromthemonthlyterrestrialwaterstorageestimatesfromGRACE(i.e. groundwater,soil
moistureandsnowcover). Weestimatethesecomponentsusingtheaverageoutputsfromtwoland
surfacemodels: 1)theCommunityLandSurfaceModel(CLM)version4.5[46]and2)theGlobalLand
DataAssimilationSystem2(GLDAS-2)model,versionNOAH-3.3[47].
We calculate time series of glacier mass change using the least squares fit mascon approach
presentedinJacobetal.[12]. Wecoverallregionswithaglacierarea > 100km2 withoneormore
mascons. For each mascon we calculate a set of Stokes coefficients associated to a unitary mass
equivalent to 1 cm of water uniformly distributed over the mascon surface Tiwari et al. [48]. We
calculatethemonthlymassanomalyassociatedtoeachmasconbysimultaneouslyfittingthemascons
StokescoefficientstotheGRACEmonthlycoefficientscorrectedforGIAandtheTerrestrialHydrology
[12].
Wefittheregionaltimeseriestoanannual,semiannual,lineartrend,andconstantcomponentsin
ordertoinvestigatethepresenceoflong-termtrendsintheglaciermassbalance. Weaddaquadratic
term to assess potential acceleration in mass loss. The total uncertainty in mass loss rate is the
sum in quadrature the GRACE measurement errors, error in GIA correction, error in hydrology
correction, leakage error caused by the assumption of uniform mass change within each mascon,
ocean/atmosphericmasscorrectionerrors,andthestatisticaluncertaintyofthemodelfit. Tocompare
thequadraticandlinearmodels,weuseavariantoftheAkaikeInformationCriterion, AIC ,with
c
smallsamplesizedatasets[49]thataccountsforthegoodnessofthefitandnumberofparameters,
andidentifieswhichmodelbestfitsthedata. WecalculateGRACEmeasurementerrorsfollowing
Wahretal.[50]. GIAandhydrologyuncertaintiesarecalculatedasinJacobetal.[12]andGardneretal.
[51]. Toevaluatetheleakageerroronthetrendintroducedbyassumingauniformmassdistribution
withineachmascon,weuseaMonteCarloapproach. Wegenerate10,000syntheticconfigurationsby
distributingthemonthlymasschanges,(fromthemasconinversiondescribedabove),non-uniformly
withineachmascon. The10,000randomconfigurationsaregeneratedusingpseudo-randomweights
extractedfromaGaussiandistributionwithstandarddeviationequalto1andzeromean. Ateach
iteration,weapplythemasconfittothesyntheticdataandcalculatetheregionalleakageerrorasthe
differencebetweenthetrendobtainedapplyingtheinversiontoGRACEandthetrendobtainedusing
therandomlygeneratedsyntheticdata. Ourfinalestimateoftheuncertaintyassociatedtotheleakage
errorissetequalto2timesthestandarddeviationofallerrorsmeasuredinthe10,000iterations. The
leakageerrorassociatedwithourmasconconfigurationistakenatthe2σerror. Trendsuncertainties
duetotheoceanmasscorrectionarecalculatedfollowingVelicognaandWahr[52].
Weinvestigatetheinterannualvariabilityoftheglaciermassanomalybysmoothingtheregional
timeserieswitha13-monthmovingaveragefilter. Wedefineaglaciologicalyearasthetimeperiod
separatingtwoconsecutiveablationseasons,i.e. the1stSeptemberinyearnand31stAugustinyear
n+1. Wecalculatetheregionalglacierannualmassbalanceasthedifferencebetweenvaluesfromthe
smoothedtimeseriesintwoconsecutiveSeptembermonths.
2.2. ICESataltimetry
TheNASAIce,Cloud,andLandElevationSatellite(ICESat)wasanEarthorbitinglaseraltimeter
designedtomeasureelevationchangesoficesheetsandarcticicecapsmissionoperativebetween
February2003andOctober2009. ThemaininstrumentcarriedbythesatellitewastheGeoscience
Preprints (www.preprints.org) | NOT PEER-REVIEWED | Posted: 3 October 2018 doi:10.20944/preprints201810.0056.v1
Peer-reviewed version available at Remote Sens. 2018, 10, 1817; doi:10.3390/rs10111817
6of25
LaserAltimeterSystem(GLAS),asingle-beamlaseraltimeteroperatingatawavelengthof1064nm
[53]. Thesatelliteflewatanelevationof600kilometerswithanorbitinclinationof94degrees. GLAS
tookmeasurementsat40-Hz,thereforesamplingtheEarth’ssurfaceelevationevery172metersalong
thesatelliteflightdirectionwithafootprintsizeofalmost70meters[54]. Here,weuselevel-2Global
LandSurfaceAltimetryHDF5(GLAH14)release34,whichincludequalityattributesandelevation
correctionsforeachfootprint. Thedataqualityattributesincludeawavesaturationflagtoindicate
saturationofthesensorwhenrecordingthereturnedpulseandacorrectionofapotentialbiasinthe
extracted elevation [55,56]. We apply the saturation and the inter-campaign bias corrections. The
Gaussian-Centroid correction has been already applied to this data release [57,58]. We use glacier
outlinesfromthe6threleaseoftheRandolphGlacierinventory[23]toisolateelevationmeasurements
over the NZEM. Following Moholdt et al. [11] and Gardner et al. [16], we do not apply the data
parameter-based cloud filtering [59], as it tends to be too strict and to discard a high number of
elevationmeasurementsthatcouldbeotherwiseusedtomeasureelevationchanges. OverNovaya
Zemlya,theapplicationofthedata-cullingprocedurereducesthenumberofavailablemeasurements
by65%.
Weevaluateglacierelevationchangesemployinganalong-trackplane-fittechniquefollowing
Moholdtetal.[60].Thisapproachevaluatesiceelevationchangeratesandterrainslopesimultaneously
byleastsquaresfittingatimevariableplanefunctiontoasetofmeasurementswithinaspatialrange
of700m. Theplanefunctionispresentedinequation(1)wherex,y,trepresenteasting,northingand
time,respectively.
dz
z(x,y,t) = α +α x+α y+ t (1)
0 1 2
dt
Theplanefitisappliedatmultiplelocationsconsideringplaneswithcentroidsequallyspacedalong
each satellite track. The spacing is chosen equal to one half the size of the plane radius (350 m),
consecutiveplanesarethereforepartiallyoverlapping[11]. Toreducetheeffectoferroneouselevation
measurements,theleastsquaresinversionisappliediterativelyateachcentroidlocationdiscardingall
measurementswitharesidualvaluelargerthan3timesthestandarddeviationofallresiduals(3-σ
Clip). Thisprocedureisrepeatediterativelyuntilallobservationshavearesidualcontainedwithinthe
3-σinterval[61,62]. Inaddition,weconsideronlymeasurementscontainingatleast15observations,
frommorethan4differentsub-tracks,anddistributedonatemporalintervalofatleast2years. We
rejectelevationchangesderivedfromplaneswithanestimatedslopegreaterthan10◦.
2.3. Cryosat-2
TheESACryoSat-2satellitesystemwaslaunchedinApril2010tostudychangesaffectingthe
arcticregions[63]. CryoSat-2payloadconsistsofahighlysophisticatedradaraltimetersuitableto
study land and sea ice. The use of radar technology allows the altimeter to work in all weather
conditionswithnolimitationsrelatedtocloudcoverage,whichispersistentintheArcticandhighly
affectslaseraltimeters. CryoSat-2fliesatanelevationof717kmwithanorbitinclinationequalto92◦
thatallowsanalmostunprecedentedcoverageofpolarregions. ThepayloadoftheCryoSat-2mission
consistsintheSAR/InterferometricRadarAltimeter(SIRAL).SIRALisclassifiedasapulselimited
delay-doppleraltimeter(workingintheKuband-13.573GHz,wavelength∼2.2cm)withafootprint
sizeof1.65x1.65km2. Insteadoftransmittingasinglepulsepertime,SIRALtransmitsaburstof64
pulsesevery50microsecond. Thereturningechoesareaffectedbysmallshiftinphaseduetothe
movementofthesatelliteduringtheintervalbetweenthebursttransmissionandacquisition(delay
doppler effect) [64]. The phase shift between the different echoes is used to increase the altimeter
resolutionalongthedirectionofflightofthesatellite(azimuthdirection)uptoamaximumof400
metersonaflatterrain[65]. Thesmallsizeoftheradarfootprintmakespossibleforthefirsttimethe
applicationofsatelliteradaraltimetryinregionsliketheicesheetmarginsandsmallarcticicecaps
Preprints (www.preprints.org) | NOT PEER-REVIEWED | Posted: 3 October 2018 doi:10.20944/preprints201810.0056.v1
Peer-reviewed version available at Remote Sens. 2018, 10, 1817; doi:10.3390/rs10111817
7of25
wherethepresenceofsteepslopesanddiscontinuoussurfaceslimitedtheapplicationofprevious
microwavealtimeters[63].
Over NZEM, SIRAL operates in SAR Interferometry mode (SARIn). In this configuration, a
secondary antenna that receives each backscattered signal from a slightly different position with
respecttotheprimaryoneisactivated. Thesignalsreceivedbythetwoantennasarecharacterized
by a shift in phase due to the different paths linking them to the ground scatterer. Analyzing the
differenceinphasedelaybetweenthetwosignals,SIRALdetectsnotonlytherangedistancebetween
thesatelliteandthegroundscattererbutalsoitslocationintheacross-trackplane[63,66,67]. Major
limitationstotheapplicationofCryoSat-2toestimateelevationchangesoverGICareassociatedtoa
radarfootprintoneorderofmagnitudelargerthanforalidarsystem. Radarperformancesignificantly
deterioratesoverterrainswithslopeshigherthan1◦. However,Wangetal.[68]demonstratedthat
evenintheseoperativeconditions,theSIRALaccuracyisstilloneorderofmagnitudehigherthanthe
accuracyofpriorradaraltimeterse.g. ERS-1/2andEnvisat. AsecondlimitationofCryoSat-2isrelated
tothetemporalvariabilityinscatteringpropertiesoficeandsnow. Asreportedin[69–72]changes
insnowpackreflectionintroduceepisodicartifactsthatinfluencetheestimationofelevationchange
especiallyovershorttimeperiods. Elevationmeasurementsarenotalignedalongtheflightdirection,
butscatteredaroundit. Aftermultiplepassagesoverthesameregion,theelevationmeasurementsare
distributedoveralargearea,complicatingtheapplicationoftheclassiccross-overorco-linearanalysis
usedwithotheraltimeters. Thenon-uniformdistributionofthereturningechoesintroducesbiasesin
thederivationoficeelevationchanges. Elevationmeasurementsaremainlyathigherelevations,while
depressionsandlowelevationregionsareundersampled[69].
Inthisstudy,weemployLevel-2Baseline-CSARinmodeelevationdataprovidedbytheEuropean
Space Agency [73] and freely available at (ftp://science-pds.cryosat.esa.int). We locate elevation
measurementsovericeusingtherecentlyreleasedglacierinventorybyRastneretal.[24]. Wemap
elevation changes over NZEM using the plane fit employing all available elevation measurement
locationsasplanecentroidsasinWoutersetal.[74]. Givenasinglecentroid,welocateallelevation
datawithin1,000mfromtheplanecentroid. Weleastsquaresfitthetimevariableplanefunctionof
theelevationobservations. Werepeatthefititerativelyuntilallelevationmeasurementsarewithinthe
3-σinterval. SimilartoourICESatanalysis,weleast-squaresfitthetime-variableplanefunctiontothe
elevationobservationsfromCryoSat-2. GiventhelargerglacierareaconsideredintheCryoSat-2case,
werepeattheinversionusingtwomoremodelfunctions:
dz
z(x,y,t) = α +α x+α y+α x2+α y2+ t (2)
0 1 2 3 4
dt
dz
z(x,y,t) = α +α x+α y+α x2+α y2+α xy+ t (3)
0 1 2 3 4 5
dt
Weobtainourfinalestimateselectingtheelevationtrendvalueassociatedtothemodelproviding
thebestfittotheconsideredelevationmeasurements. WeusetheR2 asamodelselectioncriterion
adj
(i.e. wechosethedh/dtvalueassociatedwiththemodelfunctioncharacterizedbythelargestR2 ).
adj
Weusetheplane-fitresidualtocalculateregionalmonthlytimeseriesoficeelevationchangefollowing
Woutersetal.[74]andNoëletal.[75](Formoredetails,seeAppendixA).
2.3.1. AveragingProcedure
WeaveragetheelevationchangemeasurementsobtainedwithICESatandinCryoSat-2ona1
kmgrid. ThegridisdefinedonthestandardNSIDC/NorthPolarStereographicProjection(EPSG:
3413)(https://nsidc.org/data/polar-stereo/ps_grids.html). Theaveragingprocedureisappliedasin
Woutersetal.[74]followingthestepspresentedbelow:
• Iftheconsideredgridcellcontainsasinglemeasurement,weassignthisvaluetotheentiregrid
cell;
Preprints (www.preprints.org) | NOT PEER-REVIEWED | Posted: 3 October 2018 doi:10.20944/preprints201810.0056.v1
Peer-reviewed version available at Remote Sens. 2018, 10, 1817; doi:10.3390/rs10111817
8of25
• If the considered grid cell contains multiple measurements, we set its elevation change rate
equal to the measurements mean. We use a 3-σ iterative filter in order to discard erroneous
observationsthatcouldinfluenceourfinalelevationchangeestimate.
2.4. SpatialExtrapolation
Theevaluationofthetotalregionalicevolumechangerequirestheextrapolationofelevation
change rates at the locations where measurements are not available. In this study, we update the
approachpresentedinMoholdtetal.[11]bydesigninganextrapolationschemethatconsidersnot
onlytherelationbetweenelevationchangeandabsoluteelevation,butalsoitsspatialvariability. We
extrapolateelevationchangeratesforgridcellswithnoassignedelevationchangevaluebyapplying
thefollowingsteps. Givenagridcellwithnovalue:
• wefindallthemeasurementsavailablewithinadistancesmallerthanorequaltoaselectedsearch
radiusandweevaluatethequadraticpolynomialfunctionthatprovidesthebestparametrization
oftherelationbetweentheconsideredelevationchangeestimatesandtheirmeanelevation.
• Weassigntheconsideredgridcellanelevationchangevaluecalculatedusingtheparametersfrom
thequadraticpolynomialfitandtheaverageelevationofthegridcellprovidedbyArcticDEM.
ArcticDEMisahigh-qualityDEMusingopticalstereoimageryfromtheDigitalGlobeconstellationand
highperformancesoftwarefromtheNationalGeospatial-IntelligenceAgency(NGA)andtheNational
ScienceFoundation(NSF)[76],similartowhatwasproducedforGreenland[77]butextendedtothe
entireArctic.
In order to define a realistic distance to be used in extrapolation procedures, we evaluate the
semivariogram of the elevation changes measured using the plane fit. We find an autocorrelation
distance (lag-correlation) equal to 50 kilometers. We therefore assume that ice elevation changes
remaincorrelatedfordistancessmallerthan50kmandwesetthesearchradiusequaltoonehalfof
theestimatedautocorrelationdistance. Furtherdetailsregardingthespatialextrapolationalgorithm
areprovidedinAppendixB.
2.5. FromVolumeToMassChange
Theicevolumechangeforeachgridpointiscalculatedbymultiplyingtherelativemeanelevation
change value with the grid cell area. The total ice volume change is calculated by summing the
contributionsfromallgridcells. Thelackofdataregardingpossiblechangesiniceelevationdueto
snow/firncompactioncomplicatestheconversionoficevolumechangetomassvolumechange[78].
FollowingMoholdtetal.[11],wecalculatethetotalicemasschangebymultiplyingthetotalvolume
changewiththedensityofice(0.917g/cm3)andconsidertheeffectofthisassumptioninourerror
budget. Weprovideourfinalestimatesconsideringuncertaintytermsrelatedto: elevationchange
measurement and extrapolation error; error in the considered glacier area; error associated to the
non-uniformdistributionoftheelevationchangemeasurementsontheglaciersurface;errorassociated
withthevolumetomassconversion. Adetaileddescriptionofthecalculationofthedifferenterror
componentsisprovidedinAppendixC.
2.6. AtmosphericTemperaturesandTotalPrecipitation
WeanalyzesurfacetemperaturesusingERA-Interimreanalysisdata[79]availableat(http://
www.ecmwf.int/en/research/climate-reanalysis/era-interim). Weemploy"monthlymeansofdaily
means"of2-metertemperature. ERA-Interimhasaresolutionof0.75◦ thatiscoarseversusthesizeof
theregionofinterest. WethereforecomparetheERA-Interimlinearlyinterpolateddatatomonthly
temperaturemeasurementsobtainedfromtheHydrometeorologicalInformation-WorldDataCentre
BaselineClimatologicalDataSets(http://meteo.ru/english/climate/temp.php)attwometeorological
stationsavailableintheregion: MalyeKaramakuly(59E,70.4N-WMOID:20744)andE.K.Fedorova
Preprints (www.preprints.org) | NOT PEER-REVIEWED | Posted: 3 October 2018 doi:10.20944/preprints201810.0056.v1
Peer-reviewed version available at Remote Sens. 2018, 10, 1817; doi:10.3390/rs10111817
9of25
(52.7E,72.3N-WMOId:20946)(SeeFigure1).Thecomparisonshowsanexcellentagreementbetween
thetwodatasets.Inbothcases,wefindarootmeansquareddifferencebelow0.7◦Candacorrelation>
0.98.WethereforepursuetheregionalanalysisusingonlythemonthlytemperaturesfromERA-Interim.
Weevaluatetheseasonalmeantemperaturesandtemporalanomaliesinreferencethetimeperiod
1979–2002inordertoinvestigateclimatevariabilityintheregion.WealsouseSynopticMonthlyMeans
ofTotalPrecipitationtoevaluateaccumulationvariabilityduringtheperiodunderanalysis. Wedo
notemploydatafrommeteorologicalstationinthiscasesincethesedataareconsideredbiasedand
underestimatesolidprecipitation[11,80].
3. Results
ThemasschangetimeseriesforNZEMbetweenApril2002andAugust2016areshowninFigure
2. BetweenApril2002–August2016,theNZEMglacierslostmassatanaveragerateof8±4Gt/yr.
Theregionalmassbalancewasconsistentlynegativebetween2002and2016withtheexceptionofthe
periodbetween2008and2011. After2011theglaciermassbalancewasagainstronglynegative. The
yearswiththelargestmasslosswere2004,2013,2015and2016(SeeTable1). Theglaciermasschange
timeseriesexhibitsapronouncedinterannualvariabilitycharacterizedbythethreemainphases: a
firstphaseofstablemassloss(2002–2008),asecondphaseofpositivemassbalance(2008–2010),anda
thirdandfinalphasewithincreasinglynegativemassloss(2010–2016). TheCryoSat-2andGRACE
timeseriescompareverywellandexhibitasignificantcorrelationof0.92.
Figure 2. In black, NZEM glacier mass change time series derived from GRACE using a mascon
inversionforthetimeperiodbetweenApril2002andAugust2016(monthlymeasurementerrorisgrey
shade). Inred,regionalglaciermasschangeusingICESatbetweenOctober2003andOctober2009
(rangeofvariabilityinmasschangeisreddash).Ingreen,regionalmasschangetimeseriescalculated
usingCryoSat-2betweenJuly2010andJuly2016withverticalerrorbars.
Preprints (www.preprints.org) | NOT PEER-REVIEWED | Posted: 3 October 2018 doi:10.20944/preprints201810.0056.v1
Peer-reviewed version available at Remote Sens. 2018, 10, 1817; doi:10.3390/rs10111817
10of25
Table1.AnnualMassBalanceof
NovayaZemlya
expressedinGigatonnes[Gt].
Year AnnualMassBalance
2003 -0.1
2004 -25.9
2005 -6.0
2006 -11.9
2007 -11.8
2008 -9.5
2009 4.5
2010 9.8
2011 -9.2
2012 -9.9
2013 -26.1
2014 -4.8
2015 -23.1
Total -124
Table2showstheicemasschangetrend(anditscomponents)measuredbyGRACEfortheentire
mission (2002–2016), the periods overlapping with ICESat (2003–2009) and CryoSat-2 (2010–2016).
Thiscomparisonshowsthatafter2010,theNZEMglacierslostmassatarateof14.3±4Gt/yr,hence
significantlyhigherthanfortheentireperiod. However,whenanalyzingthemasschangetimeseries
overtheentireperiodunderobservation,the AIC indicatesthataquadraticfitisnotasignificantly
c
betterfittotheregionaltimeseriesthanalinearmodel. Thismeansthatdespitetheincreasedglacier
masslossafter2010,theentiretimeseriesdoesnotexhibitthecharacteristicsofastatisticallysignificant
acceleratedmassloss. Thisresultsisconfirmedbythe R2 criterioncalculatedforthetwomodels.
adj
We find the same value 0.8 for both the linear and quadratic models. Between 2003 and 2009, we
calculate a mass loss of 9.9±5 Gt/yr with a significant positive acceleration equal to +3±1 Gt/yr
peryear,indicatingaslowdowninmasslossafter2008. Theregionalglaciermassbalancebecomes
positive for two consecutive years. For the years 2009 and 2010, we calculate a total annual mass
balanceof+4.5and+9.8Gt/yr(SeeTable1). After2010,NZEMglacierslosemassatanincreasingrate.
Between2002and2017,glaciermassbalancewasnegativefor11ofthe13consideredglaciological
years,contributingtoacumulativeicelossof124Gtequivalentto0.32mmofSLR.
Table2.AnnualMassBalanceof
NovayaZemlya
expressedinGiga-tons[Gt].
TimePeriod TotalMassBalance GRACE GIA Hydrology
04/2002-07/2016 −8±4 −4.5±2.9 −3.2±1.6 -0.2±0.2
10/2003-10/2009 −10±5 −6.2±2.6 −3.2±1.6 −0.3±0.1
07/2010-07/2016 −14.3±4 −10.9±2.4 −3.2±1.6 −0.2±0.1
Our GRACE estimates agree with Reager et al. [6] for the time period 2002–2014 and Matsuo
andHeki[81]forthetimeperiod2004–2012,bothobtainedusingGRACEdataproductsfromthe5th
release. Wefindalargedisagreement,however,withSunetal.[39]whoreportedamasslossrateof
1.04±0.25Gt/yrbetween2003and2014(versus7.5±4Gt/yrhere). Weattributethediscrepancy