Accumulation reconstruction and water isotope analysis for 1736‒1997 of an ice core from the Ushkovsky volcano, Kamchatka, and their relationships to North Pacific climate

. An ice core was retrieved in June 1998 from the Gorshkov crater glacier at the top of the Ushkovsky volcano, in central Kamchatka. This ice core is one of only two recovered from Kamchatka so far, thus ﬁlling a gap in the regional instrumental climate network. Hydrogen isotope ( δ D) analyses and past accumulation reconstructions were conducted for the top 140.7 m of the core, spanning 1736–1997. Two accumulation reconstruction methods were developed and applied with the Salamatin and the Elmer/Ice ﬁrn-ice dynamics models, revealing a slightly increasing or nearly stable trend, respectively. Wavelet analysis shows that the ice core records have signiﬁcant decadal and multi-decadal variabilities at different times. Around 1880 the multi-decadal variability of δ D became lost and its average value increased by 6 ‰. The multi-decadal variability of reconstructed


Introduction
Alpine ice core records have two important roles.One is to provide data for studying local climate characteristics.The other is to fill geographic gaps in the network of glaciochemical and climatic records used to reveal how regions are affected by global climate change (Wagenbach, 1989;Cecil et al., 2004).There are several glacial regions around the North Pacific Ocean that can provide ice cores for paleoclimate reconstructions.The Ushkovsky/K2 ice core is one of two ice cores recovered from Kamchatka (Shiraiwa et al., 2001;Matoba et al., 2011).It is a paleoclimate record for this region that complements regional climate reconstructions from ice cores in Alaska and Yukon (e.g.Holdsworth et al., 1989;Wake et al., 2002;Shiraiwa et al., 2003;Fisher et al., 2004;Yasunari et al., 2007;Fukuda et al., 2011) (Fig. 1a).
Important characteristics of climate in the North Pacific region are decadal climate variabilities, the Pacific Decadal Variability (PDV) or North Pacific Decadal Variability (NPDV).One is the Pacific Decadal Oscillation (PDO), the first leading mode of North Pacific sea surface temperature variability (Mantua et al., 1997;Minobe, 1997).Another is the North Pacific Gyre Oscillation, the second mode of sea surface height anomalies (Di Lorenzo et al., 2008).These variabilities play important roles in modulating precipitation, temperature and other climatic elements around the North Pacific Ocean (e.g.Mantua and Hare, 2002;Di Lorenzo et al., 2008).
proxies.North Pacific climate variability has therefore been reconstructed using tree rings from the Northeast Pacific region (e.g.Biondi et al., 2001;D'Arrigo et al., 2001).Summer mean maximum temperature reconstructed from Kunashir tree-ring width indices (Jacoby et al., 2004) is negatively correlated with the summer PDO index (Mantua et al., 1997) and has similar multi-decadal spectral properties.Climate variability has also been detected from glaciers in Alaska and Kamchatka that changed their mass balance at the climate regime shift in 1976/1977 (Hodge et al., 1998;Shiraiwa and Yamaguchi, 2002), a timing that coincides with a major PDO shift (Hartmann and Wendler, 2005;Josberger et al., 2007).Holdsworth et al. (1989) reconstructed snow accumulation rates at Mt. Logan and found a shift in the dominant period of accumulation variation around 1850.The accumulation increased notably after the mid-19th century, paralleled by a warming over northwestern North America (Moore et al., 2002).
To provide a further climate proxy for the western North Pacific we analysed water isotope ratios and reconstructed accumulation rates from the Ushkovsky/K2 ice core (Shiraiwa et al., 2001).Ushkovsky volcano (56.04 • N, 160.28 • E, 3903 m a.s.l.) is located in the central part of the Kamchatka Peninsula (Fig. 1b).The summit of the volcano is covered by a glacier of 43 km 2 .Two craters, Gorshkov and Herz, lie near the highest part and are both filled by glacier ice.The larger, Gorshkov, is ∼ 750 m in diameter with a maximum depth of ∼ 240 m determined by ice-penetrating radar (Matsuoka et al., 1999).In June 1998, at the K2 site near to this deepest point, a 211.7 m long ice core was drilled (Shiraiwa et al., 2001).It was the first ice core drilled in the western North Pacific region.
The firn temperature was measured continuously for a year at the BH1 site ∼ 200 m south of K2.At 10 m depth the annual mean firn temperature is −15.8 • C. The monthly mean surface temperature is −5.8 • C at its maximum in August.Melting occurs only in the surface layer and the meltwater refreezes in the surface snow due to its low temperature.Melting does not, therefore, significantly disturb the annual stratigraphy preserved in the core (Shiraiwa et al., 2001).The average accumulation rate from 1969 to 1996 was determined as 0.57 m w.eq. a −1 (w.eq.: water equivalent) by shallow ice core analysis (Shiraiwa et al., 1997).
The isotopes δ 18 O and δD in the K2 core have already been analysed from the surface down to 110 m.The ice at this depth fell as snow in approximately 1823 based on dated layers of volcanic ash and counting of annual layers identified by changes in the water isotope ratio (Shiraiwa and Yamaguchi, 2002).Accumulation in Ushkovsky correlates with the Kliuchi meteorological station winter and hydrological year precipitation (Solomina et al., 2007).
We extended the record down to a depth of 140.7 m, equivalent to a date of 1735.To determine accumulation, we applied two glacier flow models to reconstruct annual layers down to the depth of the oldest dated volcanic ash layer (Murav'yev et al., 2007).We used the Salamatin (Salamatin et al., 2000) and Elmer/Ice (Zwinger et al., 2007) models.In addition, we investigated decadal and multi-decadal oscillations of the δD and accumulation records and carried out comparisons with local and large-scale climate data, demonstrating that this ice core provides climatic information for the greater region, and fills a gap in the regional instrumental climate network.

Isotope measurement and dating
To measure isotope values the K2 ice core was sliced in the low temperature room (−20 • C) at the Institute of Low Temperature Science, Hokkaido University.To ensure that annual layers could be accurately identified, slices were taken at intervals small enough to provide at least eight samples for each annual layer.The intervals varied in thickness from ∼ 100 mm near the top of the core, ∼ 50 mm near the middle of the core, and ∼ 30 mm near the bottom of the core.After slicing, each ice sample was melted in an individual sealed plastic bag.The resulting water was stored in 30-50 mL glass vials pending isotope measurements.
The water isotopes for samples from the first 110 m were measured with a Finnigan MAT Delta S mass spectrometer at the Graduate School of Environmental Science, Hokkaido University.Absolute measurement errors were 0.1 ‰ for δ 18 O and 1 ‰ for δD (both at 2σ confidence).The isotopes in slices from 110.0-140.7 m were measured at the Institute of Low Temperature Science.The δ 18 O isotope was measured using a Finnigan Delta Plus mass spectrometer (error 0.06 ‰ at 2σ confidence), and δD using an Isoprime-PyrOH mass spectrometer (error 0.3 ‰ at 2σ confidence).Two substandard water probes (secondary water standards) were used in each analysis, and results were corrected by the SMOW-SLAP scale (Hagemann et al., 1970).
The annual layers in the K2 core were dated by reference to four volcanic ash layers.The ages of these ash layers are known based on their chemistry (Murav'yev et al., 2007, Fig. 7).They are also identified by their mineral composition and stratigraphic features (Ovsyannikov et al., 2001).These layers are from the 1985 Bezymianny eruption (12.04 m below surface), 1956 Bezymianny (35.49 m), 1829 Klyuchevskoy (102.82 m), and1737 Bezymianny (138.45 m).The tephrochronological interpretation of the dates of the four ash layers is also supported by analytical calculations (Salamatin et al., 2000).Shiraiwa and Yamaguchi (2002) confirmed the ages of the 1737 Bezymyanni, 1829 Klyuchevskoi and 1956 Bezymyanni layers by counting annual layers in the ice core.They were identified by locating the minima of the seasonally varying amounts of δD.The dating error is ±2 yr at 102.82 m depth (Shiraiwa and Yamaguchi, 2002).

Accumulation-rate reconstruction
The accumulation rate, i.e. the amount of snow that falls each year, is reconstructed by combining core data with models.The core data provide the age of the annual layers and their thickness.The thinning rate, how quickly annual layers compact, can be estimated using ice flow models, measured vertical velocities, and radio echo sounding results (Paterson and Waddington, 1984;Raymond et al., 1996;Schwerzmann et al., 2006).Shiraiwa and Yamaguchi (2002) reconstructed past accumulation rates with the glacier flow model of Salamatin et al. (2000) that includes firn dynamics.We applied two accumulation reconstruction methods for glacier flow models with firn dynamics; the Lagrangian backtrajectory method for the Salamatin model and the Eulerian method for the Elmer/Ice model.Salamatin et al. (2000) devised a firn and ice flow model for a crater glacier and applied it to the Gorshkov crater glacier.It is a thermo-mechanically coupled two-dimensional glacier flow-line model that considers the effects of a flow tube of variable width.Isotropic polycrystalline ice behaves as a linear-viscous body at small stresses (Lipenkov et al., 1997;Salamatin et al., 1997).This linear rheology relating stresses and strain rates (Salamatin and Duval, 1997) is included in the model.Using a normalized ice-equivalent vertical coordinate, analytical solutions for the velocities and consequently the strain rates are obtained.The settings of the glacier model are given in Salamatin et al. (2000) and Shiraiwa et al. (2001).

Salamatin model
The trajectories of ice particles can be computed with the velocity field given by the analytical solution of the model.It is therefore possible to estimate the change in annual layer thickness by using the back-trajectory method.The origin (position at time t 0 = June 1998) of the back-trajectory of each annual layer is set to the central point of the layer at position x c (where x = (x, y, z) is the position vector, x and y span the horizontal plane and z points upward), and the end of the back-trajectory is at the glacier surface.The backtrajectory of an annual layer is given by where the integration is carried out backward in time, t, until the glacier surface is reached, and v is the analytical solution for the velocity field provided by the Salamatin model.The vertical strain rate, D zz , along each trajectory follows as The thinning function, R(x c ), is defined as the ratio of the annual layer thickness at the position x c in the glacier and at the surface, both counted in metres of ice equivalent.It is obtained by the exponential of the integral of the vertical strain rate along the back-trajectory, where t s is the time in the past when the back-trajectory ends (meets the glacier surface).The accumulation rate, a s , that

Elmer/Ice
Elmer/Ice is a thermo-mechanically coupled, threedimensional flow model (e.g.Gagliardini and Zwinger, 2008;Durand et al., 2009;Zwinger and Moore, 2009;Seddik et al., 2011Seddik et al., , 2012)).It is based on the open-source multiphysical simulation software Elmer (elmerice.elmerfem.org)which employs the finite element method.In the version we used (Zwinger et al., 2007), Elmer/Ice accounts for compressible firn (firn model by Gagliardini and Meyssonnier, 1997).It diagnostically solves the full Stokes equations for the Gorshkov crater glacier using the density profile from the K2 borehole and the present-day surface temperature.
We have improved the approach by Zwinger et al. (2007) in two major ways: -In order to avoid the necessity of describing a boundary condition at the northern outflow of the crater glacier, we have added to the model domain a 500 m long downhill slab with a constant thickness of 50 m and a slope of 0.2 (estimated from the map shown in Fig. 2 of Shiraiwa and Yamaguchi, 2002).
-In parameterizing the geothermal heat flux (Zwinger et al., 2007, Eq. 27), we optimized the parameters q min geo (minimum geothermal heat flux at the deepest point of the crater glacier) and m (exponent of the bed elevation) for best agreement with the temperature measured and the age profiles at K2. Optimizing was carried out with APPSPACK ("Asynchronous Parallel Pattern Search") software and produced values of q min geo = 17.578 mW m −2 and m = 2.3359.Applying the model provided a velocity field in the northsouth transect of the crater glacier through the K2 drillsite (Fig. 2).
In order to compute layer thinning in the Gorshkov glacier from the three-dimensional velocity field from the Elmer/Ice simulation, we applied the thinning function R(x) introduced above.It is equal to unity at the surface, and decreases with depth.In contrast to the Salamatin model, the argument of R is not limited to the positions x c of annual layers in the K2 ice core; it can now be any position x in the June 1998 configuration of the glacier.Whereas the computed velocity field accounts for compaction, this is neglected in the definition of the thinning function R(x).Therefore, an auxiliary thinning function R(x) is defined by where ρ i = 910 kg m −3 is the density of pure ice and ρ(d) the variable density of the firn at depth d.This function is governed by the equation The vertical strain rate is given by D zz = ∂v z /∂z.According to Eqs. ( 5) and ( 6) the surface boundary condition for Eq. ( 7) is where ρ s is the firn density at the glacier surface.We set ρ s = 0.45 ρ i (Zwinger et al., 2007).We implemented the numerical solution of Eq. ( 7) for R(x) in Elmer/Ice, and the original thinning function R(x) results from Eq. ( 6).For a dated annual layer at position x = x c with thickness h (in metres of ice equivalent), the accumulation rate a s is given by Eq. (4).

Isotope (δD) record and accumulation-rate reconstruction
Annual layers were more than adequately sampled.There were 1461 samples covering the 170 yr 1828-1997, equivalent to 8.55 samples per annual layer.There were 1132 samples from the 92 yr 1736-1827, equivalent to 11.44 samples per layer.Isotope measurements were hugely variable over the years (Fig. 3).The standard deviation of δD was 15.8 ‰ in the previous analysis (from the top to 110.03 m depth), and was 18.4 ‰ in the newly measured part (110.03-140.60 m depth).These values are more than ten times larger than the measurement errors (Sect.2.1).The annual mean isotope ratio was estimated using the dating results (Fig. 4).The 20 yr running mean of δD is also shown.The average annual mean isotope We also determined the deuterium excess d = δD − 8 δ 18 O because this is a good indicator of changes in atmospheric circulation, precipitation seasonality and moisture source (Dansgaard, 1964).However, the standard deviation of the annual mean deuterium excess was small, ∼ 1.3 ‰, for the entire period (although there were seasonal cycles).Because the measurement errors are large (Sect.2.1), it is difficult to analyse the annual mean deuterium excess of the ice core in a reasonable way.
Accumulation reconstruction using the two glacier flow models indicated that the ice at 140 m depth thins by a factor of 2 (water equivalent depth) in the Elmer/Ice model and by a factor of 1.5 in the Salamatin model (Fig. 5).The substantial difference between these two reconstruction methods is mainly due to the different rheologies.The constitutive flow law used in Elmer/Ice is a power law which approaches the conventional Glen flow law with a stress exponent of n = 3 in the limit of pure ice, whereas in the Salamatin model a linear rheology is assumed.The deeper ice is therefore softer in the Elmer/Ice model than in the Salamatin model.Because of this, and other differences between the models as detailed in Sect.2.2, we feel unable to judge which model is more appropriate.However, the difference between the two reconstructions provides a good estimate of the uncertainty associated with reconstructed accumulation.The trend of the accumulation rate has been nearly stable (Elmer/Ice) or slightly increasing (Salamatin) over the last 260 years (Fig. 6).The mean accumulation rates from 1736 to 1997 were 0.62 m w.eq. in the Elmer/Ice model, and 0.51 m w.eq. in the Salamatin model.There were clear multidecadal fluctuations of both reconstructed accumulation rates (20 yr running mean).They were strong in ∼ 1730-1850 and became weak after 1950.

Wavelet analysis
Morlet wavelet analysis was applied to detrended (i.e.we removed the linear trend over the entire period of the signal) annual mean δD (Fig. 7a) and detrended reconstructed accumulation rates (Fig. 8a) in order to examine in more detail decadal and multi-decadal oscillations and its variation.The global wavelet spectrums are also shown (Figs.7b and 8b).These were tested for significance at the 95 % confidence level against red noise (Torrence and Compo, 1998).Wavelet results for accumulation rate reconstructions derived by the Salamatin model and Elmer/Ice are roughly the same, since the linear detrending operation removes their greatest difference, and so only the Elmer/Ice results are plotted here.The isotope δD expresses decadal (10-20 yr band) and multi-decadal (40-60 yr band) variabilities that are statistically significant.Similar variabilities are found in the global wavelet spectrum for periods of 11-13 and 35-52 yr.The reconstructed accumulation rates also show significant decadal (10-20 yr band, before 1850) and multi-decadal (20-50 yr band) variabilities (28-35 yr in the global wavelet spectrum).For δD, the multi-decadal variabilities decrease after ∼ 1840 and fall below the confidence interval at ∼ 1880.It matches the shift in average value of δD.Both mean values and the dominant periods of δD changed between the end of the Little Ice Age and the present.The variability of the reconstructed accumulation rates also shows this change in significance.Multi-decadal significant variability (40-60 yr band) found in ∼ 1740-1830 weakens thereafter.Then the 20-40 yr domain becomes significant in ∼ 1860-1940.Some paleoclimatic proxies retrieved from the eastern North Pacific also show the change of decadal time periods of precipitation after the mid-19th century.Reconstructed accumulation rates derived from the Mt.Logan ice core (Holdsworth et al., 1989) also show that the dominant time period changed in the mid-19th century.There was a 36 yr peak before 1860 but it became less significant after that period whereas there were 4, 11 and 21 yr peaks after 1880.Tree ring records in western North America (Gedalof and Smith, 2001;Gray et al., 2003) also show the shift of multidecadal variabilities in the mid-19th century.The match of the Ushkovsky/K2 ice core results with other paleoclimatic studies in this region demonstrates that the mid-19th century climate change affected North Pacific surface precipitation conditions.However, the mechanism is not necessarily the same in both places.Water isotope records from the Mt.Logan Eclipse ice core suggest the shift between mixed and zonal flow regimes of water vapour transport occurred around 1840 (Fisher et al., 2004).However, there is a clear difference between the water isotope ratios of the Mt.Logan core and the Ushkovsky core.They decrease in the Mt.Logan core but increase in the Ushkovsky core between the Little Ice Age and the present.Increasing ratios were also found in the Altai ice core and this increase was interpreted as a warming signal (Henderson et al., 2006).

Comparison with local climate and glaciation data
We compared the annual mean δD with precipitation and air temperature records from the Kljuchi weather station (56.32 • N, 160.83 • E, approx.40 m a.s.l., WMO Global Surface Network).Monthly precipitation  and monthly mean air temperature  were used to calculate the annual mean air temperature and annual total precipitation.The annual mean δD was not significantly correlated with the annual mean air temperature (r = −0.22,n.s.) or the annual total precipitation (r = −0.14, n.s.).Here and below we lowered the effective degrees of freedom in order to reduce the effect of the autocorrelations in the two time series being correlated which would otherwise inflate the probability of obtaining significant correlations (Trenberth, 1984).The monthly-precipitation-weighted temperature (Steig et al., 1994) was derived to account for seasonal variation in precipitation.This was significantly correlated with the weighted annual mean temperature at Kljuchi station and δD (r = 0.45, p < 0.05) in 1967-1989.This suggests that δD reflects seasonal variation of precipitation rates and temperature though it is difficult to obtain a simple relationship between δD and temperature.
There are peaks at 1810-1860, 1920 and 1970 in the 20 yr mean reconstructed accumulation rates.There are troughs around 1760-1770 followed by a steady increase to large values at ∼ 1850.This again indicates that the reconstructed accumulation rate is related to the history of glacier advance and retreat in the region.The peaks of ∼ 1810-1860 and 1920 coincide with glacier advances on Kamchatka Peninsula in the mid-19th and early 20th century (Solomina et al., 1995).Increases around 1970 coincide with positive mass balances of glaciers in Kamchatka (Shiraiwa and Yamaguchi, 2002).Reconstructed, detrended annual accumulation from this ice core coincides to a degree with past winter accumulation rates reconstructed using a model of Koryto glacier, located halfway up the east coast of Kamchatka (Yamaguchi et al., 2008).Solomina et al. (2007) found that the Ushkovsky accumulation record is in phase with the regional tree ring chronology KAML for ∼ 1830-1880.Our results show that this is also the case for ∼ 1730-1830.The mass balance phase of Kozelsky glacier is generally in phase with the Ushkovsky accumulation.There is an anti-phase relationship between Ushkovsky accumulation rates and Asian PDO expression (D' Arrigo and Wilson, 2006).This is consistent with the antiphase relation between accumulation rates and PDO index (Shiraiwa and Yamaguchi, 2002).

Comparison with large-scale climate data
We also compared our δD signal and reconstructed accumulation rates with surface temperatures (ECMWF ERA-40 reanalysis, 2 m monthly mean temperature fields for the period 1958-1996, in 3 yr means;Uppala et al., 2005).For δD, there is a significant positive correlation with mid-latitude North Pacific (20-30 • N) surface temperatures, and a significant negative correlation with subpolar (40-50 • N) North Pacific surface temperatures (Fig. 9a).The correlation map of reconstructed accumulation rates and ERA-40 2 m air temperatures shows a significant negative correlation with the subpolar North Pacific (40-60 • N, 180-150 • W), and a significant positive correlation with the western coast of North America (40 • N, 125 • W and 60 • N, 145 • W) (Fig. 10a).This evidence indicates that the ice core δD records reflect extratropical North Pacific surface climate conditions.
The pattern of the correlation map of δD resembled the second leading mode of the North Pacific surface temperature (Bond et al., 2003), which is related to the North Pacific Gyre Oscillation (NPGO) (Di Lorenzo et al., 2010;Furtado et al., 2011).The correlation of annual mean δD with the annual mean NPGO index (Di Lorenzo et al., 2008) was r = 0.27 (p < 0.10) for 1950-1997 and r = 0.70 (p < 0.01) for 1979-1997.These results suggest that δD reflects the NPGO after the climate regime shift in 1976/1977, but not before that.Since then, the location of the Aleutian Low was shifted westward and the central pressure was deepened (Nitta and Yamada, 1989;Trenberth, 1990).The reconstructed accumulation rates also correlated with the NPGO index for 1950-1997 (r = 0.29, p < 0.10); however, for 1979-1997 the correlation was smaller (r = 0.21, n.s.) than that for δD.
We further compared ice core signals with a sea surface temperature (SST) record, namely the Extended Reconstructed Sea Surface Temperature (ERSST) version 3 (Smith et al., 2008), in three year means (1854 to 1995,  (Numaguti, 1999).The ocean is thus the main source of water vapour for the Kamchatka Peninsula, particularly in winter.The precipitation record at Kljuchi showed that winter precipitation is higher than during other seasons, which suggests that there might be high winter precipitation ratios in the annual layers of the K2 ice core.
Since the main source of the ice in the core is precipitation from the North Pacific, the variation of water isotope ratio and reconstructed accumulation rates reflect the conditions of the region.Numaguti (1999) differentiates between western and eastern Pacific sources, and shows that the main water vapour source is the West Pacific.The back-trajectory analysis from the Ushkovsky volcano (Kawamura et al., 2012) indicates that the air masses come from the Northwest Pacific and Siberia.Therefore, one of the main water vapour sources is the Northwest Pacific.By contrast, the Arctic Ocean is not a main source, even though it is important for the material transport.
The NPGO index is defined by the sea surface height anomaly of the California current region in the Northeast Pacific (Di Lorenzo et al., 2008).However, it is not directly reflected in the ice core records because of intervening processes.Transportation of the signal via the ocean through Rossby wave propagation takes several years at least (Ceballos et al., 2009).Other possibilities are connections via the atmosphere.There is a significant correlation between the low-passed North Pacific Oscillation (NPO) (Walker and Bliss, 1932;Rogers, 1981) and the SST anomaly around Kamchatka (Furtado et al., 2012).This represents atmospheric forcing of the NPGO (Ceballos et al., 2009;Furtado et al., 2011Furtado et al., , 2012)).Linkin and Nigam (2008) showed that this sea level pressure pattern, NPO, is linked with the upperair west Pacific connection pattern (WP pattern; Wallace and Gutzler, 1981), and that the NPO/WP is a prominent mode of winter mid-latitude variability.The correlation of the winter term (January-March) WP index with annual mean δD is significant (r = 0.50, p < 0.05) for 1977-1997, though it is not significant for 1950-1997 (r = 0.13, n.s.).This implies that the correlation between δD and the NPGO index after the climate regime shift in 1976/1977 reflects that beween δD and the WP index.
In contrast, the relationship between δD and SST for 1854-1995 differs significantly from that between δD and the ECMWF ERA-40 surface temperatures for 1958-1996 discussed above (Fig. 9b).(This difference also occurs if the ERA-40 surface temperatures are replaced by SST for 1958SST for -1996; not shown.The crucial point is thus the different time scale.)The structure shown in Fig. 9b resembles that of the PDO.There is a negative correlation between the winter (November-April) PDO index and precipitation and temperature in Kamchatka (Mantua and Hare, 2002).Winter precipitation increases and decreases cause corresponding changes in annual mean δD through changes in winter layer thickness.One of the reasons for the structure of Fig. 9b would be this precipitation-weight effect.The negative relationship between winter precipitation in Kamchatka and PDO index (Mantua and Hare, 2002) also agrees with the negative phase relationship between reconstructed accumulation and PDO index (Shiraiwa and Yamaguchi, 2002).
The water isotope ratio of precipitation depends on its sources and the condensation temperatures during transport (Dansgaard, 1964).Although variation in δD is caused by its dependence on the source conditions, this is not the cause since it also depends on the precipitation-weighted local temperature, as suggested above.Detailed water isotope modelling analysis (cf.Field et al., 2010;Porter et al., 2014) of Kamchatka and this ice core site can improve understanding of the detailed mechanism behind water isotope variation.

Summary and conclusions
We presented and analysed water isotope records for the top 140.7 m of the Ushkovsky/K2 ice core.The average level of δD changed by 6.0 ‰ from 1736-1880 to 1910-1997.This may well have been caused by an increase in average temperature between the end of the Little Ice Age and the present.It is difficult, however, to ascertain a simple relationship with annual mean temperature because annual δD is modified by seasonal precipitation change, as is suggested by the comparison with precipitation-weighted temperature.
Lagrangian and Eulerian accumulation reconstruction methods were developed for the thermo-mechanical coupled 2-D Salamatin (Salamatin et al., 2000) and 3-D Elmer/Ice (Zwinger et al., 2007) glacier models including firn dynamics.These reconstructions showed that the layer of snow deposited in a year on the surface compresses by 0.50 to 0.75 at 140 m depth at the K2 drillsite.Reconstructed accumulation rates in the late Little Ice Age were slightly less than, or almost the same, as today.One of the large uncertainties is the different firn and ice rheology employed in the two models.Glaciological observations of firn properties at the site and transient simulations using a firn layer model would improve the accumulation reconstruction.The highest reconstructed accumulation rates, at ∼ 1810-1850, 1910 and 1970, coincide with glacier advances and mass balance changes and moraine records in glaciers in Kamchatka (Solomina et al., 1995;Yamaguchi et al., 2008).Therefore, the variation in accumulation can be used as a proxy for glacier accumulation changes in central Kamchatka.It is thus likely that Siberian High and Aleutian Low activities affect the changes of Kamchatkan glaciers (Barr and Solomina, 2013).Analysis of the activity of these pressure systems, their spatial extent and elevation will improve understanding of the forcing factors for glacier mass balance in the region.
The δD record and the reconstructed accumulation rates show significant decadal and multi-decadal variabilities and changes of multi-decadal time periods in the mid to late 19th century.The time at which the multi-decadal oscillation signal in δD was lost coincided with changes in its average value between the end of the Little Ice Age and the present.The shifts of the dominant time period for the reconstructed accumulation also appears in the 1850s.
We found significant correlations between the ice core records and ERA-40 surface temperatures in the North Pacific, somewhat more pronounced for δD than for the reconstructed accumulation rates.The significant relationship between the δD record and the winter NPGO index is caused by its correlation with the winter WP index; the NPO/WP pattern is the forcing of NPGO (Ceballos et al., 2009;Furtado et al., 2011Furtado et al., , 2012)).
However, there is no significant relationship between δD and the NPGO and WP indices before the climate regime shift in 1976/1977.The correlation maps of the ice core and 140-yr SST records imply that there is a different pattern on this time scale.To determine which components are most important for the region, atmosphere-ocean interaction analysis for precipitation or temperature in this region would be required.Long-term isotope climate modelling will further enhance understanding of the relationships between δD and the PDV or other connection patterns and, more generally, of climate change between the Little Ice Age and the present.
The NPGO and WP indices are only defined since the 1950s due to a lack of observations.Past NPGO and WP reconstructions based on models would therefore be required to extend our analysis.Comparison of the water isotope records from the Ushkovsky/K2 ice core with coral records in the Northwest Pacific (Felis et al., 2009(Felis et al., , 2010) ) could also help further understanding of the connection between the source of water vapour and precipitation in Kamchatka.

T.
Sato et al.: Ice core from the Ushkovsky volcano, Kamchatka corresponds to an annual layer at position x c with thickness h (in metres of water equivalent) is

Fig. 2 .
Fig. 2. Velocity field (absolute values as texture and direction vectors) simulated with Elmer/Ice in a north-south transect of the Gorshkov crater glacier through the K2 drillsite.
Fig. 7. (a) Wavelet plot for δD: binary logarithm of the power (absolute square of the wavelet transform) in ‰ 2 .Black contours represent the 95 % confidence level against red noise (Torrence and Compo, 1998).Areas outside the cone of influence (indicated by the thick white line) are to be discarded.(b) Global wavelet spectrum for δD.The thin dashed line indicates the 95 % confidence level against red noise.
Fig.8.a) Wavelet plot for reconstructed accumulation rates (Elmer/Ice): binary logarithm of the power (absolute square of the wavelet transform) in m 2 a −2 .Black contours represent the 95 % confidence level against red noise(Torrence and Compo, 1998).Areas outside the cone of influence (indicated by the thick white line) are to be discarded.(b) Global wavelet spectrum for reconstructed accumulation rates (Elmer/Ice).The thin dashed line indicates the the 95 % confidence level against red noise.
Fig. 9. (a) Correlation map of the ice core δD with ERA-40 3 yr mean surface (2 m) temperature for 1958-1996 (39 yr).Black contours indicate the land-water margin and the 90 % confidence level.b) Correlation map of the ice core δD with ERSST 3 yr mean sea surface temperature for 1854-1995 (142 yr).Black contours indicate the land-water margin and the 90 % confidence level.