Research ArticleGEOLOGY

High sensitivity of Bering Sea winter sea ice to winter insolation and carbon dioxide over the last 5500 years

See allHide authors and affiliations

Science Advances  02 Sep 2020:
Vol. 6, no. 36, eaaz9588
DOI: 10.1126/sciadv.aaz9588


Anomalously low winter sea ice extent and early retreat in CE 2018 and 2019 challenge previous notions that winter sea ice in the Bering Sea has been stable over the instrumental record, although long-term records remain limited. Here, we use a record of peat cellulose oxygen isotopes from St. Matthew Island along with isotope-enabled general circulation model (IsoGSM) simulations to generate a 5500-year record of Bering Sea winter sea ice extent. Results show that over the last 5500 years, sea ice in the Bering Sea decreased in response to increasing winter insolation and atmospheric CO2, suggesting that the North Pacific is highly sensitive to small changes in radiative forcing. We find that CE 2018 sea ice conditions were the lowest of the last 5500 years, and results suggest that sea ice loss may lag changes in CO2 concentrations by several decades.


Summer sea ice in the Arctic Ocean has been shrinking in recent decades (1) in tandem with increasing CO2 emissions (2). However, winter Bering Sea sea ice extent (Fig. 1), which forms in winter and is absent in the summer under modern climate (3), has remained relatively stable and/or has increased (4) over the satellite record, suggesting that winter sea ice extent is less vulnerable to anthropogenic climate change and is more dependent on ocean-atmosphere circulation variability (5). Long-term projections predict a 34% loss in winter (February) sea ice extent for the Arctic as a whole by CE 2081–2100 using Coupled Model Intercomparison Project 5 (CMIP5) projections under representative concentration pathway (RCP) 8.5 (6). However, Bering Sea winter sea ice extent in CE 2018 and CE 2019 was 60 to 70% lower than the previous mean spring (February, March, April, and May) extent from CE 1979 to CE 2017 (1), suggesting that Bering Sea winter sea ice is diminishing more rapidly than models predict. The decline in these years was attributed to anomalous southerly atmospheric flow that also increased near-bottom water temperatures (7). How this recent warming and sea ice loss in the Bering Sea fits into the long-term context of climate change remains unresolved because of spatial gaps and low temporal resolution of regional paleoclimate and paleo–sea ice records. This is due in part to depositional limitations on the shallow Bering Shelf that underlies much of the Bering Sea, which has been more prone to erosion and low, irregular sediment accumulation during the Holocene.

Fig. 1 Map of the Arctic and selected sites discussed in the study, including the study site.

(A) Polar view of the study region, including 1. St. Matthew Island (this study), 2. Sea ice diatom record from GC-33 (14), 3. Uk’37 alkenone paleotemperature record (13), 4. Agassiz ice core record (8), and 5. IP-25 record (9). Winter (blue) and summer (red) 1981–2010 median ice edges are shown (1). Adak and Bethel, two towns with GNIP isotope data, are shown (42). (B) St. Matthew Island showing area of inset (white box) (Image: Google Earth, accessed 19 June 2019). (C) St. Matthew Island peatland inset, showing coring location (red circle) (Image: Google Earth, accessed 19 June 2019). The North Pacific is distinguished from the Bering Sea by the Aleutian Island arc, which terminates to the west of Adak. The Bering Sea is bounded in the north by the Bering Strait.

The radiative forcing from increasing anthropogenic CO2 concentrations has led to the rapid retreat of perennial summer sea ice in the Arctic Ocean basin today over the last several decades (2), reversing late Holocene cooling trends. However, rising atmospheric CO2 [~10 parts per million (ppm)] and other greenhouse gases, during the mid to Late Holocene [~6 thousand years (ka) ago to preindustrial present], coincided with cooling temperatures (8) and expanded sea ice (9) in the Arctic Ocean, suggesting that the region’s sea ice is more strongly forced by decreasing summer insolation (~25 W m−2) through ice-albedo feedbacks than the relatively small changes in preindustrial CO2 (~1 W m−2) (10). More broadly, a global proxy compilation of Holocene temperatures suggests that global cooling has occurred since the mid-Holocene (11), contrasting with warming recorded in Earth system models due to the radiative forcing of rising greenhouse gases in the atmosphere (12), suggesting that proxy reconstructions are regionally or seasonally biased. This mismatch in the proxy data and model results, referred to as the Holocene temperature conundrum, remains unresolved (12). Previous data compilations show sparse coverage over the North Pacific, where sea surface temperatures (SSTs) warm (13) and winter sea ice in the Bering Sea decreases (14) over the mid- to Late Holocene. The resulting effect is an asynchrony between the Arctic and North Atlantic Ocean basins and North Pacific proxy records.

To examine controls on Bering Sea winter sea ice, we use modern simulations from an isotope-enabled global spectral model (IsoGSM) to interpret peat cellulose oxygen isotopes from a peat core on St. Matthew Island, Alaska (Fig. 1). We used this record to infer atmospheric circulation changes and Bering Sea winter sea ice extent over the last 5.5 ka.


The St. Matthew Island peat core began accumulating peat ~5.5 ka ago (fig. S1). Macrofossil assemblages were dominated by Carex spp. (sedge), Sphagnum, and brown mosses, with mosses dominating the record until ~2.8 ka ago and sedges thereafter (fig. S2). To evaluate the relationship of δ18O cellulose (δ18Oc) as a function of both plant species variability and peatland water, we collected plant and water samples across St. Matthew Island in June 2018. We found that water collected from streams, peatland pools, and ponds/lakes fell along the global meteoric water line (GMWL) (fig. S3A), suggesting that peatland water reflects unfractionated precipitation. We found no statistical difference among the δ18Oc of peatland plants (fig. S3B) and therefore did not include any species corrections to the δ18Oc record.

To interpret downcore δ18Oc changes from a peat core on St. Matthew Island, we modeled atmospheric conditions over the instrumental period (CE 1979–2018) for the St. Matthew Island grid cell using the IsoGSM model (15) (Fig. 2 and figs. S4 and S5), trends that were confirmed with the shorter LMDZiso model (16) (fig. S5), and compared δ18Oc anomalies to sea ice extent (Fig. 3). The years with δ18Oc values outside of 1 SD on either end (i.e., low and high composites or the years that were either higher or lower than 1 SD of the mean) were used to identify atmospheric patterns and sea level pressure (SLP) anomalies associated with those high and low composite anomalies (Fig. 2). These results confirm interpretations from previous studies (17, 18) that higher than average isotopic values are associated with winds dominating from the south, and those with lower than average isotopic values are associated with winds dominating from the north (Fig. 2).

Fig. 2 Isotope-enabled GCM results.

Sea level pressure (SLP) anomalies, wind anomalies, and isotopic anomalies for high-isotope composite years (left) and low-isotope composite years (right).

Fig. 3 The relationship of IsoGSM (15) modeled mean FMAM δ18O values (red) to maximum winter sea ice extent in the Bering Sea (black).

Gray bar on CE 2018 represents the standard deviation of the δ18Oc from peat moss corrected for modern water values measured that year, and the black bar represents the mean value. VSMOW, Vienna Standard Mean Ocean Water.

Peatland water δ18O reflects the atmospheric δ18O (fig. S3) and is an annual integration of the precipitation that falls in a given year. To determine the months that influence the variance in δ18O values the most, we examined the total monthly precipitation and the variance of the δ18O values for the St. Matthew Island grid cell (fig. S4). The lowest precipitation anomalies and lowest overall precipitation occurred in the summer during the instrumental record, and highest overall precipitation anomalies occurred in the autumn and winter (September to February). The interval with the highest monthly deviation in δ18O values between the “high-composite” and “low-composite” years occurs in February to May (FMAM) (see Materials and Methods and fig. S4), a period that historically coincides with maximum sea ice extent in the Bering Sea. Therefore, although the peatland water represents an annual integration of precipitation that falls in any given year, variability is most heavily influenced by conditions during the winter and spring months when sea ice is either expanding or retreating.

In comparing the mean modeled δ18O of precipitation from the IsoGSM model for the St. Matthew Island grid cell and Bering Sea FMAM sea ice extent for the CE 1979–2018 records, we found a strong negative correlation when we lagged the sea ice record by 1 year (−0.773, P < 0.00001; Fig. 3). This suggests that the δ18Oc of the St. Matthew Island peat core records not only changes in atmospheric circulation but also sea ice extent, where lower values coincide with more northerly flow (Arctic), and an expansion of sea ice and higher values coincide with more atmospheric flow from the south (North Pacific).

The St. Matthew Island peat δ18Oc record shows lower values before 3.2 ka ago, shifting to a higher mean thereafter. A return to lower δ18Oc values occurs ~1.5 to 0.8 ka ago, before a stepwise increase toward present day. Modern (CE 2018) δ18Oc values are the highest of the previous 5.5 ka, consistent with the anomalously low sea ice extent observed from the instrumental record (1) of that year.


The relationship of Bering Sea oxygen isotopes to sea ice extent

Although stable oxygen isotopic (δ18O values) changes in precipitation across much of the Arctic are considered sensitive to temperature, the isotopic ratio of precipitation in the Bering Sea is highly sensitive to the dominant wind direction and, by extension, sea ice extent (5). This is because the precipitation over St. Matthew Island is essentially a recorder of marine precipitation, unlike the central portion of the ice sheets where most ice cores are derived and where precipitation is sensitive to orographic and continental effects. The δ18O values and their relationship to sea ice over the instrumental record in the North Pacific region more broadly, including the Bering Sea, are governed primarily by variations in the winter storm tracks (17) linked to the strength and position of the Aleutian Low (AL) (5, 17).

The modeled precipitation values over the island are negatively correlated with sea ice extent over the same period (lag +1 year, −0.773, P < 0.00001; Fig. 3), where shifts in the wind direction drive variability in sea ice extent (Fig. 2). The correlation suggests that ocean warming further drives heat into the Bering Sea via southerly atmospheric flow in the subsequent year (19), setting up a positive feedback loop toward more heat transport into the Bering Sea (and vice versa). Other factors influencing oxygen isotope fractionation of precipitation include the transport distance of moisture, where precipitation becomes more depleted in 18O with distance from the source, such as with overland or across sea ice transport, and temperature (Supplementary Materials) (18). Further amplification of these effects has previously been explained by internal forcing mechanisms, such as positive phase of the Pacific Decadal Oscillation (PDO) and enhanced El Niño–Southern Oscillation (ENSO) (5), which lead to deepening of the AL and, hence, stronger and more persistent North Pacific winds entering into the Bering Sea (5) when the low pressure center is in a more northwestern position (Fig. 2).

The IsoGSM simulations suggest that the loss of spring sea ice in the eastern Bering Sea (14) manifests as higher oxygen isotopes of precipitation (Figs. 2 and 3 and fig. S5). While the δ18Oc values reflect an annual integration of precipitation on St. Matthew Island (fig. S3), the largest deviations in the seasonal cycle in the present day occur in FMAM (fig. S4). This indicates that these months are the most crucial for interpreting past changes in both precipitation and, ultimately, maximum sea ice extent from St. Matthew Island and the Bering Sea, respectively. The modern plant samples collected from the peatland in June 2018 correspond with the highest δ18Oc values of the entire 5.5-ka peat core record, consistent with the lowest sea ice extent over the observational period (7). We therefore interpret the record to show that the Bering Sea has gradually transitioned from greater to lower sea ice extent over the last 5.5 ka, overprinted with larger excursions consistent with internal variability, such as enhanced ENSO and AL activity after ~3.2 ka (19). This interpretation agrees with a spring sea ice diatom proxy study from the Bering Sea (Fig. 4B) (14), where more perennial or extensive spring sea ice conditions are present before ~3.2 ka, after which the St. Matthew Island δ18Oc increases and the sea ice diatoms diminish, suggesting less extensive winter/spring sea ice, driven by greater tropical North Pacific influence (19) on conditions in the Bering Sea. The trend toward increasing δ18Oc in the most recent part of the record suggests both overall increased influence of North Pacific storms entering the Bering Sea and a retreat of winter sea ice. The δ18Oc data suggest that CE 2018 Bering Sea winter sea ice extent was the lowest of the last 5.5 ka.

Fig. 4 Composite summary figure.

(A) St. Matthew Island δ18Oc (gray) and smoothed (maroon) with the cellulose values corrected to reflect precipitation, Gulf of Alaska Uk′37 alkenone paleotemperature record (13), and the winter insolation curve at 60°N, (B) eastern Bering Sea diatom record (14), (C) Agassiz ice cap δ18O record (8) and the summer insolation curve at 90°N (20), (D) IP-25 sea ice proxy record from the Canadian archipelago (9), and (E) composite atmospheric CO2 concentrations, spanning instrumental (Mauna Loa) and Antarctic Ice core records [European Project for Ice Coring in Antarctica (EPICA) Dome C and West Antarctic Ice Sheet (WAIS) Divide Ice core] (21). Shading indicates cooler intervals in the Bering Sea associated with lower CO2 concentrations and higher sea ice extent.

Drivers of Holocene Bering Sea sea ice changes

The long-term increase in δ18Oc from St. Matthew Island toward warmer Bering Sea conditions and decreasing winter sea ice extent, similar to other regional records of sea ice extent (14) and SSTs (13), follows mid- to Late Holocene increases in both winter insolation and atmospheric CO2 concentrations (Figs. 4 and 5). Despite decreasing mean annual solar insolation during the Holocene, winter and spring insolation increase (20), which can explain decreasing winter sea ice in the Bering Sea over the duration of the record (Figs. 4 and 5), along with the long-term increase in atmospheric CO2 concentrations. However, on shorter multidecadal to centennial time scales, shifts in δ18Oc correspond to perturbations in atmospheric CO2 concentrations, which suggests that Bering Sea ice extent is sensitive to shorter-term changes in atmospheric CO2 (Fig. 5, B and C). The strength and persistence of directional winds can drive greater shifts in δ18Oc values at St. Matthew Island (Fig. 5C) (5, 17) than the perturbations in CO2 given the nonlinearity of the relationship between Bering Sea winter sea ice and CO2.

Fig. 5 The relationship of St. Matthew δ18O inferred sea ice extent to atmospheric CO2.

(A) Preindustrial and industrial atmospheric CO2 (21) and preindustrial radiative forcing (10) versus inferred sea ice extent. (B) Lag correlation relationship between St. Matthew Island and the Bering Sea diatom record (gray) (14), Gulf of Alaska Uk′37 alkenone paleotemperature record (GoA, red) (13), the atmospheric CO2 record (purple) (21), (14), including age model uncertainty for all records. (C) St. Matthew Island peat cellulose isotopes plotted as inferred winter sea ice extent (cyan) versus atmospheric CO2 concentrations (purple) (21), with the St. Matthew record shifted by 100 years.

In this study, the relationship between preindustrial atmospheric CO2 and winter sea ice extent (Fig. 5A) suggests that CO2 plays a sensitive role in the variability of winter sea ice in the Bering Sea, either directly, through radiative forcing, or indirectly, through dynamical changes in ocean-atmosphere circulation. While this study is not equipped to resolve or explore the influence of dynamical changes in ocean-atmosphere circulation on sea ice extent in the Bering Sea, it is possible to invoke the radiative impact of increasing CO2. The coherence between radiative forcing of preindustrial CO2 (21) and St. Matthew Island δ18Oc values suggests that even small (~0.5 W m −2) changes in the radiative effect of increased CO2 results in a decrease in winter sea ice in the Bering Sea. We estimate that the preindustrial radiative CO2 forcing of ~1 W m−2, combined with the increase in winter insolation of ~4 W m−2, has led to a ~42% loss of winter sea ice in the Bering Sea over the last 5.5 ka (Fig. 5C). In the Arctic Ocean basin, where perennial (summer) sea ice was increasing over the Late Holocene (Fig. 4D), the CO2 forcing was small compared with the radiative cooling of ~25 W m−2 resulting from decreasing summer insolation over the same period (Fig. 4, C and D) (20). Over the observational record, the rapid (decades to century) increase in anthropogenic CO2 forcing of ~2 W m−2 was found directly responsible for the reduction in Arctic Ocean basin summer sea ice extent (2), as it far outweighs the change in insolation forcing (Fig. 4C) over the same time period.

The response of Bering Sea winter sea ice to increasing winter insolation and atmospheric CO2 highlights the seasonal and spatial bias in existing Holocene paleoclimate compilations that exhibit a cooling trend over the Holocene (11). The warming that is modeled over this time period in response to rising CO2 concentrations (12) falls in line with the results of this study, suggesting that, by including more winter proxy records, the Holocene temperature conundrum could be resolved. Furthermore, the limited spatial coverage of the North Pacific in these studies suggests that processes related to the interocean seesaw (22), where warming in the North Atlantic is coupled with cooling in the North Pacific leading to climate stability, could be missed.

We found that Bering Sea ice extent was related to preindustrial CO2 concentrations over the Holocene (Fig. 5). Assuming an increase in CO2 translates to a loss of sea ice extent, the relationship would suggest a complete loss of sea ice by 285.9 ± 5 ppm by volume (ppmv) under the preindustrial relationship (Fig. 5A), equivalent to concentrations ~CE 1870. However, this relationship would suggest a lagged response to Bering Sea ice extent loss to the radiative warming of CO2, indicating that current winter sea ice conditions are out of equilibrium with modern climate. To evaluate whether a lag relationship exists between this study and Bering Sea diatom sea ice proxy (14), North Pacific alkenone SST (13), and the Holocene CO2 record (21), we performed a lag correlation analysis where we evaluated the relationship of the proxy records that were offset from one another by up to 400 years including age model uncertainties for all records (see Materials and Methods). While a lag correlation cannot be ruled out, the response time of interest falls within the range of age model uncertainty, but correlations remain high for all records (r > 0.55 and P < 0.05) for a lag of up to 100 years (Fig. 5B). At worst, a lag of 100 years would suggest the Bering Sea is committed to complete winter sea ice loss. At best, supposing that Bering Sea winter sea ice is in equilibrium with modern climate and following the postindustrial CO2–sea ice relationship (Fig. 5A), Bering Sea winter sea ice would be lost by the time atmospheric concentrations reach 685 ppmv, which under RCP 8.5 would be reached in ~CE 2070 and CE 2090 under RCP 6.0. It also suggests that a complete loss of Bering Sea winter sea ice could be mitigated by adopting CO2 reduction strategies in line with RCP 4.5 or lower (6).

A lagged response in the Bering Sea could explain the lack of a steep rise in the St. Matthew Island δ18Oc record when compared with the inferred abrupt warming from the Agassiz ice core record (8) that occurs in step with the exponential rise in anthropogenic CO2 (Fig. 4E). This difference can be explained by seasonal representation (summer versus winter) and the response to the respective forcing. A potential explanation for the lack of a steep rise in δ18Oc values at St. Matthew Island is that this small island is influenced more heavily by the slower thermal inertia of the ocean. While the radiative forcing on the atmosphere occurs in step with the rise in CO2, the heat exchange between the atmosphere and the surface layers of the ocean, as well as the deep ocean, i.e., the thermal inertia of the ocean (23), can lag a rise in CO2 by decades to a century (2325). Furthermore, a dynamical lag resulting from ocean circulation change in response to multidecadal modes of climate variability (e.g., PDO) (26) and ocean heat export (27, 28) can create an ocean temperature lag responsible for the possible lag in winter sea ice loss in the Bering Sea.

There are numerous sources of uncertainty in this projection. Nonetheless, even with an increase or decrease in the range for ice-free conditions by several decades, our data suggest that the potential ocean-atmosphere heating required to generate a complete loss of Bering Sea winter sea ice may already exist. Furthermore, the window for trajectories of winter sea ice loss in the Bering Sea is highly sensitive to internal forcing mechanisms, such as the phase of the PDO and AL, which can alter dominant wind directions (5, 17) to amplify sea ice loss or gain on interannual to decadal time scales (26). However, given the relationship of arctic amplification to sea ice loss in the Arctic, the change in temperature and heat absorption could increase disproportionally with ever-decreasing sea ice (29), resulting in a substantial shortening of the winter ice loss timeline. Under future warming scenarios, AL intensity is expected to increase and be driven farther north (30), both of which would amplify winter sea ice loss in the Bering Sea by increasing southerly flow (as in the high composite scenario in Fig. 2).

Implications for winter sea ice in the Bering Sea and the Arctic

The substantial rate of anthropogenic CO2 inputs into the atmosphere over industrialization suggests that a loss in Bering Sea sea ice extent is accelerating or is already committed to complete sea ice loss as a result of delayed response to anthropogenic forcing. Low winter sea ice anomalies in CE 2018 and CE 2019 indicate future conditions that favor an ice-free Bering Sea. Widespread effects of Bering Sea winter sea ice loss are expected to occur. Ecosystem responses to low sea ice in CE 2018 included altered food webs that led to sea bird die-offs and may represent a harbinger of future low sea ice extent (31). Further intensification of observed North Pacific influence in the Bering Sea leading to a reduction in sea ice can further affect heat transport to the Arctic Ocean basin. Although the Bering Strait throughflow may be relatively small (<1 Sv; 1 Sv = 106 m3 s−1), it can have a disproportionate influence on heatflux into the Arctic Ocean basin, and recent increases have been linked to weakening northerly winds (32), signifying enhanced winds originating from the North Pacific could amplify Arctic Ocean sea ice decline via increasing winds from the south. Simultaneously, the increased frequency and duration of winter cyclones in the Arctic have led to the large reductions in freezing degree days in Arctic Ocean winters (33, 34). A loss of sea ice can also increase coastal erosion and increase land temperatures that result in permafrost thaw (35), further amplifying warming (36).


St. Matthew Island (60.4°N, 172.7°W) is a small (357 km2) uninhabited island in the central Bering Sea (Fig. 1). It lies within the CE 1979–2018 median winter seasonal sea ice edge but can be surrounded by open water in some winters, as was the case in winter of CE 2018 when the Bering Sea sea ice edge was ~350 km north of its modern average position and completely absent by April when it historically reaches its maximum winter extent. Island temperatures remain cold during the winters, while summers remain cloudy and cool, fostering a vegetation community of forb tundra.

The 1.45-m St. Matthew Island peat core was collected in 2012 from a small fen near the southern margin of North Lake, on the northwestern end of the island using a Russian-style peat corer to minimize compaction. Surface vegetation at the coring location included Sphagnum fimbriatum, Carex spp., and Drepanocladus spp. in wetter depressions. The core was subsampled into 1-cm intervals, and samples were analyzed for plant macrofossils, loss on ignition, and oxygen isotopes from the analysis of bulk peat cellulose. Cellulose was extracted from the peat core at every centimeter from bulk peat using the University of Waterloo Environmental Isotope Laboratory (UWEIL) method (37).

Macrofossils from 10 samples were radiocarbon dated at Lawrence Livermore National Laboratory (4), Beta Analytic (5), and National Ocean Sciences Accelerator Mass Spectrometry (1) (table S1), and an age model was generated using Bacon (fig. S1A) (38). Two ages were rejected by the age model because one was too young for its stratigraphic position, likely having incorporated root material from the sedge leaves dated, and the other was anomalously old. Peat accumulation rates were calculated by dividing the amount of time in each centimeter interval (fig. S1B). Macrofossil abundances were analyzed at every centimeter by sieving peat through a 250-μm screen and tallied based on total relative abundance of herbaceous (i.e., sedge) peat versus bryophytic (i.e., moss) peat, using semiquantitative methods (fig. S2) (39). To calibrate the water on St. Matthew Island to the living plants, water samples were collected in June 2018 from a variety of settings (lakes, streams, bog, and fen) and compared to cellulose oxygen isotope values from a range of plant species from each respective location (table S2, A and B). We used the UWEIL method of cellulose extraction at the University of Alaska Fairbanks (37) for the peat core and a subset of modern CE 2018 samples, and we used the cuprammonium solution (CUAM) method (40) for cellulose extraction at the United States Geological Survey in Reston, VA, for the modern samples.

Modern water, plant cellulose, and cellulose from peat core samples were analyzed at the Alaska Stable Isotope Facility on a continuous flow isotope ratio mass spectrometer (CF-IRMS) using a High Temperature Conversion Elemental Analyzer (TCEA) attached via a Conflo IV to a Thermo Delta V+ IRMS. Analytical precision associated with the stable oxygen isotope analyses was <0.6 per mil (‰) and is expressed as 1 SD from the mean based on the results from multiple (n = 10) analyses of a laboratory standard (EMA-P1 from Elemental Microanalysis, part no. B2203, certificate no. BN/132358) conducted during the run of samples. A suite of international standards (NBS N-1, NBS-18, and NBS-19) were analyzed with the run (measured versus expected, r2 = 0.99) to allow the calibration of the data, which is expressed relative to Vienna Standard Mean Ocean Water (VSMOW).

Isotopic values of all water samples collected fell along the GMWL, suggesting that evaporation effects are minimal (fig. S3A). Stream water values were lower for both δ18O and δ2H values than from bog ponds and North Lake, likely because streams reflect a dominant snowmelt source, whereas ponds and lakes represent an annual integration of the water isotopic signature. To verify the relationship of δ18Oc values to surface water, we collected modern peat and water samples from St. Matthew Island and found a relationship of plant to water offsets of −33.0 ± 1.1‰ using the CUAM method (table S2). The modern δ18Oc values show a slightly lower offset was determined from the modern samples using the UWEIL method (−31.3 ± 2.0‰), and because this was the method used for the core extraction, we offset the values by this amount to attain the paleo-water δ18Oc values. These offsets are larger than the typically assumed −27.0 ± 0.3 ‰ offsets but in line with an analysis of latitudinal biochemical offsets determined in (41). We also measured the δ18Oc values between the dominant peat species: moss (nonvascular) and sedge (vascular) (fig. S3B). As no significant (P < 0.0001) difference between sedges and mosses were recorded, no species adjustments were made for the core δ18Oc values.

The closest sites in the Global Network of Isotopes in Precipitation (GNIP) to St. Matthew Island are in Adak, Alaska, in the Aleutian Islands to the south, which show little seasonal variation in the oxygen isotopes of precipitation (δ18Op) and are higher (~−9 ± 0.9‰) than stations to the north (Utqiaġvik: −18.5 + 3.7‰). Bethel, Alaska, located due east on the Yukon-Kuskokwim Delta displays oxygen isotope seasonality, with winters averaging −14‰ and summers averaging ~−10‰ and an annual average of −13‰ (42).

In the absence of recorded weather or climatological data from St. Matthew Island, we relied on global circulation models with isotope tracers to draw inferences about the controls on isotopic variability in this region. The analysis relied on two models, IsoGSM (15) and LMDZiso (16). Each of these models is an independent three-dimensional atmospheric circulation model that runs with prescribed SST and sea ice conditions and includes all relevant isotope tracer physics including fractionation associated with phase changes. The underlying physics for the IsoGSM model is the Experimental Climate Prediction Center (ECPC) Global Spectral Model, whereas the LMDZiso model is based on the Laboratoire Métérologie Dynamique (LMD) climate model. Both of these models were run using a “nudging” routine where the model is corrected to atmospheric conditions from reanalysis models. In this way, outputs from these simulations are sometimes referred to as isotope reanalysis, which means that model data from, for example, CE 2018 reflects isotope variability associated with the actual atmospheric conditions during that year. The IsoGSM model is nudged with the National Center for Environmental Prediction (NCEP) Reanalysis II product and LMDZiso to the European Centre for Medium Range Weather Forecasts (ECMWF) reanalysis product. The nudging procedure is in contrast to a “free run” simulation, where the atmosphere responds to external forcing but is not explicitly corrected to the atmospheric circulation patterns for a given window of time. For both models, we calculated the annually averaged precipitation weighted for the grid cell that includes St. Matthew Island. The LMDZiso data were accessed from the SWING2 database and data spanned 1994–2010, whereas the IsoGSM spans the entire period from CE 1979 to CE 2018. We analyzed the annual time series data for LMDZiso and IsoGSM models together to highlight that the interannual anomalies are well reproduced by both models, despite being nudged by different reanalysis products. Using the IsoGSM model, where we have a longer time series (CE 1979–2018), we identify a series of years that had anomalously high and low isotopic ratios of precipitation and use these to explore the cause of interannual variability. High composite and low composite analyses were performed on those years that fell above or below 1 SD of the mean, respectively. We first explore whether these interannual anomalies were caused by changes specific to a season. We thus calculated the average seasonal cycle from 1979 to 2018 and tested how this compared to the seasonal cycles during the years that went into the “low” and “high” composites. We found that isotopic anomalies were most pronounced during the FMAM period of the year (fig. S4) and thus look at atmospheric circulation patterns during this period of the year to understand the cause of interannual isotopic anomalies. For the anomaly maps, we first calculated the average FMAM conditions over the 1979–2018 period and subtracted this value from the average FMAM for the years that went into the composite. We show the SLP and wind vector anomalies for those years and also determined the δ18Op anomalies for those high- and low-composite years. We used the FMAM average oxygen isotope values to compare to sea ice extent over the same time period and used a lag + 1 correlation analysis to determine the statistical relationship of Bering Sea sea ice extent to the FMAM oxygen isotopes of each year from CE 1979 to CE 2018 (fig. S4).

Sea ice extent data were retrieved from the National Snow and Ice Data Center (1). We used data from CE 1979–2018. Relationships between St. Matthew Island δ18Oc record, the Bering Sea diatom record (14), Gulf of Alaska SST (15), and atmospheric CO2 concentrations (21) were investigated, incorporating age model uncertainty. For each of the time series (St. Matthew, Bering Sea diatom, Gulf of Alaska SST, and CO2), we generated 1000 plausible age models based on assumed uncertainty. For St. Matthew Island δ18Oc record, we used the uncertainty intervals produced by the Bacon age model. For the other published time series, we applied a 10% error, which, by definition, increases with time (i.e., 10% of 5 ka is greater than 10% of 1 ka). We then interpolated each record from their native resolution to whichever was the lower-resolution time series being used for the correlation analysis. We then calculated the Pearson correlation coefficient for each of the 1000 time series and for all possible lag times between +50 and −400 years. The median correlation coefficients (based on the 1-ka age models) were plotted as dots at each lag time, and we calculated the 97.5 and 2.5% (error bars) (Fig. 5B).


Supplementary material for this article is available at

This is an open-access article distributed under the terms of the Creative Commons Attribution-NonCommercial license, which permits use, distribution, and reproduction in any medium, so long as the resultant use is not for commercial advantage and provided the original work is properly cited.


Acknowledgments: We thank D. Klein for initiating this work and D. Klein and R. Kleinfelder for collecting the 2012 core. We thank M. Romano for coordinating field research efforts, T. DeGange for collecting plant and water samples in June 2018, and the Alaska Maritime National Wildlife Refuge for making the trip to St. Matthew Island possible. Stable isotope compositions of samples were analyzed courtesy of the Alaska Sable Isotope Facility. Reviews by T. Cronin, S. Praetorius, and B. Gaglioti improved earlier versions of this manuscript. We thank two anonymous reviewers, T. Roland, and the editor for their constructive suggestions for changes that improved our work. Any use of trade, firm, or product names is for descriptive purposes only and does not imply endorsement by the U.S. Government. Funding: Funding support came from USGS Climate and Land Use Change R&D (to M.C.J.) and NSF grant no. 1502776 (to M.B.). Author contributions: M.C.J. coordinated and directed the field sampling, conceived of the study, performed data analysis, and wrote the paper. M.J.W. and K.J.K. extracted cellulose. M.J.W. measured δ18O. K.J.K. assisted with lab and data, statistical analyses, and interpretation. M.B. ran and analyzed the isotope modeling. K.Y. provided the 2018 IsoGSM data. All authors contributed to the writing of the manuscript. Competing interests: The authors declare that they have no competing interests. Data and materials availability: All data needed to evaluate the conclusions in the paper are present in the paper and/or the Supplementary Materials. Additional data related to this paper may be requested from the authors.

Stay Connected to Science Advances

Navigate This Article