Complexity of Arctic Ocean water isotope (δ¹⁸O, δ²H) spatial and temporal patterns revealed with machine learning

The Arctic is warming 4 times faster than the rest of the Earth (Rantanen et al., 2022). Warming is rapidly changing the movement of water through the Arctic, which is leading to reductions in sea ice (Hinzman et al., 2013; Bintanja and Andry, 2017; Thoman et al., 2020; Cooper et al., 2022) and an amplified Arctic water cycle in which more water is cycling through Arctic systems (Vihma et al., 2016; Shupe et al., 2022). These changes in the transport of water in the Arctic can also have implications for lower latitude hydrological cycles, which can be traced with isotopes of water (e.g., δ¹⁸O, δ²H; Puntsag et al., 2016; Boutt et al., 2019; Bailey et al., 2021; Mellat et al., 2021).

The heavy and light forms, or isotopes, of oxygen (¹⁸O and ¹⁶O) and hydrogen (²H and ¹H) respond differently to fluctuating environmental conditions (Gat, 1996). Isotopes can serve as tracers of present and past changes to the water cycle due to the preferential evaporation of the lighter forms and corresponding condensation of the heavier forms (Alley, 2000; Klein et al., 2015). Shifts in the ratio of heavy-to-light isotopes are expressed relative to a standard (e.g., δ¹⁸O and δ²H) in per mil (‰). Moreover, the parameter deuterium excess (d-excess), defined as δ²H – 8*δ¹⁸O, which is a measure of the deviation from the global meteoric water line, can reveal moisture source history and other environmental information (Gat, 1996). Therefore, a better understanding of the variability of Arctic Ocean water isotopes (both water vapor and ocean water) will greatly increase the utility of water isotopes as tracers to understand the movement of Arctic water both within and beyond the Arctic.

Early studies investigated the variation of δ¹⁸O in seawater (Epstein and Mayeda, 1953). Then field studies of simultaneous seawater and water vapor isotopes in the North Pacific, at varying relative humidity (RH) levels, resulted in the Craig and Gordon model, which explored how evaporation of seawater from many different water surfaces influenced the fractionation of stable isotope ratios of water vapor (Craig and Gordon, 1965). This model development was followed by work that invoked an evaporation mechanism and the influence of RH to explain change (enrichment) of δ¹⁸O during evaporation of seawater (Lloyd, 1966). More recent experimental work showed how surface water temperature can impact evaporative fractionation and water vapor isotope values (Cappa et al., 2003). There have also been field measurements of water vapor isotopes above the ocean in southern (Uemura et al., 2008; Uemura et al., 2010; Thurnherr et al., 2021) and northern (Kurita, 2011) polar regions, discrete seawater (Cooper et al., 1997; Cooper et al., 2022), and continuous water vapor isotopes (Klein and Welker, 2016) in northern polar regions from different expeditions, and continuous water vapor and discrete seawater isotope measurements from both polar regions (Bonne et al., 2019). Additionally, previous work used interpolation (relying on the δ¹⁸O to salinity relationship) to create a gridded data set of surface seawater δ¹⁸O values that covered Earth’s oceans (including the Arctic Ocean) at a resolution of 1° (about 87 km; LeGrande and Schmidt, 2006). This first attempt at creating gridded ocean water surface δ¹⁸O values offered the promise that future work could incorporate more parameters that impact ocean circulation, which might better reveal a more realistic depiction of the sharp gradients in surface δ¹⁸O values (e.g., changes across smaller spatial areas).

We present here the first continuous in-situ measurements of water vapor and sea surface isotopes from a 2016 western Arctic research expedition on the United States Coast Guard icebreaker Healy. These water isotope data are combined with accompanying oceanographic and meteorological measurements to initiate machine-driven learning that models water vapor and sea surface water isotope composition beyond the track of the vessel and extending across the western Arctic Ocean and into the sub-Arctic Bering Sea and Gulf of Alaska. Previously, isoscapes have been used to depict the complexity and patterns of precipitation at continental and global scales (Welker, 2012; Terzer-Wassmuth et al., 2021). The novel isoscapes presented here indicate sharp gradients in surface water δ¹⁸O values that were assumed to exist (LeGrande and Schmidt, 2006). Thus, these machine learning-driven isoscapes support the idea that oceanic water isotopes (including the Arctic) are indeed more complex than previous investigations revealed. Exploring the increased complexity of Arctic ocean water isotopes is important for understanding the modern water cycle (Klein et al., 2015; Puntsag et al., 2016; Bailey et al., 2019; Fiorella et al., 2021), resolving climate proxy records (Jouzel et al., 2005; Klein and Welker, 2016; Tabor et al., 2021), and forecasting climate (Nusbaumer et al., 2017; Dee et al., 2018).

The 3 main objectives of this study were to: (1) explore variability of seawater and water vapor isotopes on a ship track; (2) demonstrate how machine learning can use these data to create isoscapes of Arctic Ocean water isotopes; and (3) present a first step showing how machine learning can help to understand the complexity of ocean water isotopes and ocean dynamics. In time, with improvement of these methods in future iterations, new machine learning-driven understanding of the complexity of Arctic Ocean water isotopes may help to better elucidate the past, present, and future global water cycle.

Water isotope (vapor and seawater), meteorological, and oceanographic data were collected during a 34-d transect on the United States Coast Guard Cutter Healy that began July 6, 2016, in the Bering Sea, passed through the Bering Strait into the Chukchi Sea, and reached its northernmost point in the Chukchi Borderlands of the Arctic Ocean before turning south through the Bering Strait and Aleutian Islands into the Gulf of Alaska (Figure 1) on August 9, 2016. The ship left Seward, Alaska, on July 2 and returned to Seward on August 10, with data collected from July 6 to August 9.

A Picarro L2130-i analyzer measured vapor δ¹⁸O, δ²H, and water concentration. The L2130-i, a cavity ring down-spectroscopy (CRDS) analyzer (hereafter analyzer), is based on cavity-enhanced, near-infrared laser absorption spectroscopy procedures, tuned on a narrow spectral region (Crosson et al., 2002; Brand et al., 2009; Gupta et al., 2009). The instrument was set up under the bow of the ship and samples were collected through a tube connected to a PVC pipe, which extended 1 m away from the ship and about 11 m above the ocean surface, analogous to previous similar field studies (Klein and Welker, 2016). The water vapor isotope ratios were measured approximately every second by the analyzer. Twice a day (approximately every 12 h), standard waters (United States Geological Survey 45, –2.24‰ δ¹⁸O and –10.30‰ for δ²H; and 46, –29.8‰ δ¹⁸O and –235.8‰ for δ²H) were injected for calibration (Text S1) of values for each 12-h period (Figure S1 and Table S1).

The end of the intake pipe was fitted with an approximately 0.5 cm cotton fiber membrane, which hinders particles and liquid water from entering the inlet line and reaching the cavity of the CRDS during water vapor isotope collection, but does not inhibit air flow. The air flow ran continuously through Teflon-based tubing that was heated and had a flow rate of 5 L s⁻¹, which helped to minimize vapor condensation in the inlet line. The analyzer’s water vapor concentration was calibrated with known humidity values from meteorological instrumentation, which resulted in specific humidity values about 19% lower than the water concentration and was similar to other estimates of Picarro water vapor concentration corrections (Dennis and Jacobson, 2014). Due to logistical difficulties, a water vapor humidity response experiment was not performed during the 2016 campaign. However, a water vapor isotope humidity response experiment was done in the field at a later time on the same Picarro analyzer (Akers et al., 2020). Based on the experiment by Akers et al. (2020), no varying effects of humidity are expected above about 4,000 ppm (Figure S1). As the humidity range in this study never dropped below 6,400 ppm, isotope-humidity dependency corrections have not been applied.

Sea surface water samples were collected from the uncontaminated science seawater system (USSS) science sample station, which continuously pumps water from about 8 m below the sea surface at a continual flow rate of about 16 L min⁻¹ (Amos and Wickham-Rowe, 2010). An intake valve linked the USSS science sample station to a Picarro continuous water sampler (CWS), which is a module that connects to a Picarro L2130-i analyzer (different from the instrument collecting water vapor isotopes) and allows for continuous analysis of water isotope values. The liquid water flow rate is approximately 2 mL min⁻¹, δ¹⁸O precision is 0.1‰, and δ²H precision is 0.2‰ for the CWS. Water isotope measurements were made about every 5 s and then incorporated into the 5 min mean seawater isotope values used in this study.

The Healy USSS science sample station has a sink with a continuous flow of seawater passing through it. Water flows through the line at about 2.5 m s⁻¹, and the flow path is about 42 m (Amos and Wickham-Rowe, 2010). Therefore, water flows through the line for about 16 s before reaching the sink. A bucket (about 2 L) was placed at the bottom of this sink and filled with this continuous flow water. The CWS sample intake line was placed at the bottom of this bucket and held in position with a weight. Based on the flow rate of the seawater (about 16 L min⁻¹), the residence time of water in the bucket was about 8 s, which is less than the 5 min mean seawater isotope values used in this study. Water temperature at the USSS science sample station matches that of the ambient sea surface water, suggesting that waters from the USSS represent surface waters. The USSS intake is not blocked by ice while the vessel is in an ice pack. Further, the supply line is not infused with air bubbles, as this step would disrupt the thrusters and dynamic positioning system on the ship (Amos and Wickham-Rowe, 2010). Therefore, there was likely minimal to no physical change, or isotope fractionation, during the movement of water from near the sea surface to where it was sampled.

The CWS allows for the connection and use of standard waters for episodic calibration, but due to problems unserviceable while at sea, these standard waters were not available for use. Instead, periodic samples (n = 14) were collected from the Healy science sample station in opaque 25 mL Nalgene sample bottles (Table S2). Standard water isotope sample collection procedures were followed during collection of these samples (e.g., capped without headspace, stored in a cool and dark room). The relationship between the CWS analyzer and lab analyzer values of these seawater samples, with known times and dates, were then used to calibrate the CWS values and apply a humidity correction (Text S2 and Figures S2–S5). Similar to the water vapor isotope ratios, the 5 min means of seawater isotope ratios were used for analysis.

For water isotope analysis of these periodic samples, 2 mL of each sample were placed in septa-capped glass vials in the lab at the University of Alaska Anchorage (UAA; Klein et al., 2016; Bailey et al., 2019). Isotopic analyses of seawater samples were conducted on a liquid water isotope analyzer (e.g., Picarro L2130-i) at the UAA Stable Isotope Lab about 2 months after collection. Results of stable isotope analyses are presented using the delta notation reported relative to the Vienna-Standard Mean Ocean Water (V-SMOW) standard, with precision of ±0.2‰ for δ¹⁸O and ±2‰ for δ²H. Each of the samples was injected 6 times and each liquid injection was measured for 5 min. To account for memory effect between samples, the values from the first 3 injections were discarded. The mean value of the remaining 3 injections was retained. Reanalysis of the sample occurred if the standard deviation of the 6 replicates was greater than 0.3‰ for δ¹⁸O and/or 3‰ for δ²H, or if the internal standard for the run differed from the accepted value by greater than ±0.2‰ or 2‰ for δ¹⁸O and δ²H, respectively.

Meteorological measurements taken aboard the Healy (from about 20 m above the water) and incorporated into analyses included: location (Kongsberg Differential Position Sensor GPS); air temperature and wind speed and direction (RM Young); sea surface temperature (Seabird SBE3S); barometric pressure (Paroscientific, Inc. MET-3A); and RH (Viasala HMP 110). Samples were collected for water isotope analyses from different conductivity, temperature, and depth (CTD) stations (Figure 1). Sea ice cover was estimated at 5 min intervals from photographs taken by a camera mounted to the aloft tower on the Healy, which resulted in repeat images from the same location and direction (in front of the bow of the ship). All the photographs were loaded into R programming language and converted to raster format. A script created 350 equally spaced points (35 columns and 10 rows) across an image around the bow of the ship (approximately 10 km² and 270° to 90° view angle) and determined if each point was open water or ice by interpretation of pixel color. The total value for each class was divided by the number of points (n = 350) to arrive at a sea ice cover value, ranging from 0% to 100% ice cover for 6,605 images. Sea ice data were grouped, and then isotopic, meteorological, and oceanographic data were analyzed relative to percent sea ice cover. This high spatial and temporal resolution sea ice estimation only included sea ice cover and did not account for vertical depth and density. Further, this technique likely underestimated ice coverage over the whole cruise because the icebreaker was often preferentially navigated to open water instead of entering sea ice.

To better understand variability of ocean water isotope values (both seawater and vapor), we developed sets of exploratory and predictive models of the spatial patterns of isotopic variability, similar to what has been previously done with large scale precipitation isotope data (Dutton et al., 2005; Vachon et al., 2010). The exploratory models were trained on the observed water isotope and oceanographic data, whereas the predictive isoscape models were trained using measured isotope values from the southbound transects, attributed with a subset of predictors (i.e., those with continuous geographic information system [GIS] datasets across the study extent), and validated with observed isotope values from the northbound transect.

Gradient boosting machines (GBM) are a type of machine learning algorithm that use sequential, binary, recursive decision trees to classify data. GBM are nonparametric and capable of accounting for nonlinear, multivariate interactions to analyze large, sometimes inconsistent datasets. They are effective for describing and predicting the complexity of ecological systems (Li et al., 2011; Wiersma et al., 2011; Humphries et al., 2018; Baltensperger et al., 2020). Predicted results are data-driven and not fit to an a priori model as in more commonly used frequentist methods (Breiman, 2001; Friedman, 2002), and so GBM are effective tools for deciphering highly variable patterns among isotopes in the marine environment. Binary decision trees classify data into terminal nodes that minimize within-node variances and sequential trees are run to explain the variance remaining after previous trees runs. Because GBM are based on random samples of data and subsets of predictors (“bagging”) and run iteratively to resolve remaining variance (“boosting”), correlative structure within the dataset is disassociated. The resulting models, grown with a slow learning rate and assessed internally using 10-fold cross validation, should be robust to overfitting. In addition to using environmental predictors to develop models, gbmidwcv (Li et al., 2019) also incorporates an inverse distance weighting function to tune predictions to appropriate numeric scales. This element becomes especially important when correlative patterns between environmental predictors and training points become less informative as the number of model training points decreases (Li et al., 2011).

To identify covariates that could serve as powerful predictors in isoscape models, we trained preliminary, exploratory GBM models on measured δ²H, δ¹⁸O, and d-excess values, attributed with 25 predictors observed at sampling locations. We then explored the relative influence of these predictors on isotope values using the spm package in R (gbmidwcv). Observed predictors used in the exploratory model included: air temperature, barometric pressure, chlorophyll-a, CO₂, dewpoint, water depth, water oxygen, radiation, RH, salinity, sea-surface temperature, ship speed, wet bulb temperature, wind direction, and wind speed. We grew models to 5,000 trees, using a learning rate of 0.01, 3 observations per node, and inverse distance power of 2. Default settings were used for all other model parameters. The gbmidwcv function calculates model accuracy based on an internal, 10-fold cross-validation on withheld subsamples of training data versus model predictions. We identified the most accurate models as those with the lowest root-mean square error (RMSE). From the most important predictors we identified those for which spatially continuous GIS datasets were available for the study area and timeframe, to be used in subsequent predictive isoscape models.

From the exploratory models, we selected 4 of the most important predictors for which continuous, spatial datasets were available for our study area and timeframe (Table 1). As isotopic values can change over short timeframes and increased variation may confound model predictions, we segregated the full training dataset into northbound and southbound datasets (Figure 1), limiting the timespan represented by each model to 15 d and 20 d, respectively. This approach resulted in large improvements in model accuracy for models trained on single transects versus those trained using the roundtrip dataset.

We proceeded to develop 2 sets (northbound and southbound) of 9 spatially predictive isoscape models across the study extent using gbmidwcv. To train spatial models, we attributed sets of observed isotope values (δ²H water, δ¹⁸O water, and d-excess water; δ²H vapor, δ¹⁸O vapor, and d-excess vapor) with 4 remotely sensed and spatially modeled environmental predictors (bathymetry, sea surface temperature, salinity, and chlorophyll-a concentrations; Table 1). We used the midpoint date of the cruise (July 20, 2016) to select the environmental predictor variables (e.g., sea surface temperature) from the available daily datasets that were used for modeling.

For each model, we varied the number of trees grown (1,000–10,000), learning rate (0.0001–0.01), minimum number of observations per node (2–10), and the inverse distance power (1–2) in order to identify hyper-parameters which produced the most accurate model for each isotope. Default settings were used for all other model parameters. We selected the most accurate model RMSE. Models were then scored (i.e., applied) to a set of regularly spaced grids of points (2 km resolution, which was the highest resolution possible given computing constraints) attributed with the same 4 spatial environmental predictors in order to make continuous predictions across the study area. Predicted values (‰) at points in the 2 km grid were smoothed using the inverse-distance weighting function and clipped to the study area extent in ArcGIS 10.8 (ESRI, Inc., Redlands, CA) to generate spatially continuous isoscape raster maps of δ²H water, δ¹⁸O water, d-excess water, δ²H vapor, δ¹⁸O vapor, d-excess vapor, δ²H water to vapor difference, and δ¹⁸O water to vapor difference.

To test the validity of the values in the isoscapes, the modeled sea surface isotope values were compared with observed values from a different expedition that collected data from the Chukchi Sea in a similar vicinity and time (i.e., mid- to late July) as our modeled values (Cooper et al., 2017). From these data, 167 δ¹⁸O water values were selected (ranging from –2.37‰ to –0.52‰ δ¹⁸O), which included all values less than 80 m deep, as values above this depth are generally part of the near surface mixed layer with similar characteristics (Figure S6). This comparative study collected multiple samples from each depth, so the mean δ¹⁸O value from each depth was used for comparison with our data (n = 28). This comparative analysis of modeled (predicted) isoscape δ¹⁸O water values with those from the observed data showed that some of the predicted values are similar to the observed values with a low calculated mean squared error (MSE) of 0.59 (Figure S7). While this comparison has values off the 1:1 line and outside the 0.25‰ range, it also shows that many are within this range. Moreover, the only observed samples used for comparison (Cooper et al., 2017) are from various depths, none of which are actually at the depth from which we collected samples for the continuous seawater isotopes (about 8 m). Additionally, while observations were from a similar time and location, they are not exactly the same temporally and spatially as the model predictions. Therefore, some isotopic variability (e.g., Figure S6) could be expected when they are compared with our modeled values. Even though this comparison has limits, it shows that the modeled values are close, especially considering the spatial, temporal, and depth differences between them.

To provide additional model validation, we compared southbound-derived model predictions with observed water isotopes along the northbound transect for both seawater and vapor isotopes. The RMSE of the difference between the predicted and observed values was used to measure the strength of the model (Table 2). The RMSE for the δ¹⁸O seawater model is the lowest (1.6), which indicates an excellent fit to the data (0 would be a perfect fit). While not quite as low as seawater, the δ¹⁸O water vapor model (4.8) is still a good fit to the data. There is more variability in the RMSE values for δ²H, with seawater having a much lower value than water vapor. The RMSE values for the d-excess models are similar for both water vapor and seawater.

However, these RMSE values include some large spatial variations in squared error values. Most of the individual squared error values at sample locations are low, but there is a small number that are very large, which then pulls the RMSE values upward (Figure 2). Upon closer examination, for most models there are 2 primary regions where these larger squared error values are found across δ¹⁸O, δ²H, and d-excess models: (1) between St. Lawrence Island and the Seward Peninsula, and (2) near the edge of the Canada Basin (Figures 2 and S8–S12). The main exception to this trend is the seawater d-excess model (Figure S10), which has larger squared error values in more locations than the other models, but still has a similar overall RMSE value to d-excess vapor. All environmental predictors are present in the areas identified with larger RMSE values, so the larger errors are not a result of missing data. While model limitations and possibilities for improvement are discussed later in this article, this validation exercise shows that the models are accurate for the vast majority of the study area, despite temporal variation in predictors and observed isotope values between north and south transects.

Water sources for collected samples were estimated using δ¹⁸O and salinity values in a stable isotope-based Bayesian mixing model in the SIMMR package for R (Parnell and Inger, 2016). The 3 end members in this mixing model were sea ice melt (SIM), river runoff (RR), and Arctic Ocean surface water, which was our saline end member (Figure S13, Table 3). The Arctic Ocean end member data for salinity and δ¹⁸O were selected from data during the expedition when the ship was not close to the shore nor to sea ice, so as to minimize the influence of freshwater. Additionally, these values agree with previously published studies of Arctic Ocean water end members: Timmermans and Marshall (2020) for salinity and Yamamoto-Kawai et al. (2008) for δ¹⁸O.

The data collected were partitioned into 4 geographic regions: (1) Gulf of Alaska, (2) Bering Sea, (3) Chukchi Shelf, and (4) Chukchi Borderlands (Figure 3). Unimak Pass is used to separate the Gulf of Alaska and Bering Sea, while the Bering Strait marks the start of the Chukchi Shelf. The shallower Chukchi Shelf is separated from the deeper Chukchi Borderlands at a latitude of 73° (Sharma, 1979).

Water vapor isotope ratios (δ¹⁸O and δ²H) were generally more depleted in the Arctic Ocean in the Chukchi Borderlands and Shelf than at lower latitudes in the Bering Sea and Gulf of Alaska, but there is still variability within the Arctic Ocean (Figure 3). For example, the δ¹⁸O and δ²H values after day of year (DOY) 200 in the Chukchi Borderlands are generally lower than those before. Additionally, the decoupling of the trends in δ¹⁸O and δ²H (e.g., separation between relative changes in δ¹⁸O and δ²H) increases north of the Bering Strait. Similar to δ¹⁸O and δ²H, there is variability in the Arctic Ocean d-excess water vapor values, but they were generally lowest in the Chukchi Borderlands, which had the greatest sea ice cover (Figure 3).

Seawater δ²H values changed during the transition from the shallower Chukchi Shelf to the deeper Chukchi Borderlands (Figure 3). There was also a gradual increase in seawater δ¹⁸O values beginning around DOY 195. Values for δ¹⁸O and δ²H display different trends, especially in the Chukchi Borderlands between DOY 195 and 203 when δ¹⁸O became enriched and δ²H became more depleted.

Seawater d-excess decreased gradually after the first transition to the Chukchi Borderlands and then decreased sharply around DOY 203, remaining low before increasing after moving back near the Chukchi Shelf. It then decreased again after going back into the Bering Sea.

Air temperatures were highest in the Gulf of Alaska and Bering Sea and became lower in the Chukchi Shelf before decreasing more in the Chukchi Borderlands. Air temperatures decreased further around DOY 203 and remained at their lowest values before increasing after moving back onto the Chukchi Shelf. Sea surface temperatures followed a similar overall pattern between the regions. The relationship between air and sea surface temperatures was dynamic in the Chukchi Borderlands, as air temperatures decreased below sea surface temperatures around DOY 203. Both increased toward the Bering Sea, where air temperatures were generally above sea surface temperatures, before becoming similar in the Gulf of Alaska.

Salinity decreased when transitioning from the Chukchi Shelf to the Borderlands. The salinity dropped further after DOY 195 and remained relatively stable before rising again after moving back onto the Chukchi Shelf. The water is deeper overall in the Chukchi Borderlands than the other regions, but is also variable, as might be expected in this area with seafloor canyons. There are some gaps in the bottom depth data where there was interference with sea ice, which was more extensive but still variable in the higher latitude Chukchi Borderlands. The most pronounced change in sea ice cover was the transition from sustained higher levels to lower cover after DOY 195.

Relative humidity was highest, and least variable, in the Chukchi Borderlands, but there were sustained decreases between DOY 195 and 200 (Figure 4). Wind speed was variable, with high and low values occurring throughout the different regions. Wind direction was also variable, but was most often from the north-northwest, especially in the Chukchi Borderlands. Barometric pressure ranged from about 985 mb to 1025 mb and was variable throughout the expedition.

The water vapor isotopes measured (at about 11 m) can include contributions from both local evaporation as well as more distant air parcels. Accordingly, this work assumes that water vapor is composed of contributions between these local and distant sources, which are influenced by factors such as RH and wind speed. However, even with these variable water vapor sources, the difference between ocean water and vapor isotopes may help to better understand the interaction of some oceanic and atmospheric processes. For example, the difference during periods that might favor local evaporation (e.g., lower RH and winds) can help to understand more about possible fractionation patterns. If winds are stronger (implying mixing of more distant moisture sources) and RH is higher (suggesting less local evaporation), then the differences between isotopes can help to understand the contributions of these different sources (e.g., the similarity of difference values to expected fractionation patterns). These difference values may eventually help to understand evaporation and fractionation processes that could be better constrained in the future with additional data like water vapor flux. Overall, the difference between surface water and vapor δ¹⁸O and δ²H values are variable, but show the most change at the higher latitude Chukchi Borderlands region with a large shift (δ¹⁸O becoming less negative relative to δ²H) around DOY 195, then switching back around DOY 204 (Figure 4). The relative difference changes in δ¹⁸O and δ²H values then remain similar before widening around DOY 211 in the Chukchi Borderlands, narrowing after reaching the Chukchi Shelf, and then widening again in the Gulf of Alaska.

The isoscape models focus on the southbound transect as it had the largest training dataset for the models. Some general patterns and highlights of the isoscapes (δ¹⁸O and d-excess) are mentioned, but not all areas and details are described due to the high resolution and complexity they present. δ²H isoscapes can be found in Figures S14 and S15.

The southbound δ¹⁸O seawater isoscape shows a range of values from –0.1‰ to –2.9‰, with the Chukchi Sea and Arctic Ocean generally more negative in δ¹⁸O values than the Gulf of Alaska and southern Bering Sea (Figure 5). More negative δ¹⁸O water values are predicted near the shorelines and deltas of large Arctic rivers like the Colville in the western Beaufort and Mackenzie in the eastern Beaufort Sea. There is also more negative δ¹⁸O value water predicted in the Bering Sea southwest of St. Lawrence Island. Salinity had the most influence on spatially modeled δ¹⁸O seawater values, and the RMSE between the modeled and withheld values was 0.06 (Table 4).

The isoscape of seawater d-excess reveals a range of values from –13‰ to 13‰ (Figure 6). The Arctic Ocean shows the greatest variability with the highest values (of the whole Arctic and sub-Arctic region) in the Chukchi Sea northwest of Alaska, while the Beaufort Sea generally has lower values. Similar to the δ¹⁸O isoscape, there is some high d-excess value water southwest of St. Lawrence Island. Water depth had the greatest influence on spatially modeled d-excess seawater values, and the RMSE was 0.39.

The isoscape of δ¹⁸O water vapor ranged from –25‰ to –12‰, with more negative values in the Arctic (especially in the Canada Basin) and relatively more positive values in the Gulf of Alaska (Figure 7). Near the deltas of large rivers (e.g., Yukon, Kuskokwim, Mackenzie), the δ¹⁸O vapor values are more negative. Additionally, the east Bering Sea region generally has more positive (predicted or expected) values than the west Bering Sea, except for the vapor values southwest of St. Lawrence Island. Salinity had the most influence on spatially modeled δ¹⁸O water vapor values, and the RMSE was 0.29.

The d-excess isoscape of water vapor shows broad variability, with the highest values in the Gulf of Alaska (about 18‰–24‰) and the lowest values (based on the southern transit data) in the Chukchi Borderlands (about –12‰; Figure 8). Additionally, there is substantial variability in the Arctic with higher values in parts of the Canada Basin (about 14‰–22‰) than the Chukchi Sea. Sea surface temperature had the most influence on spatially modeled d-excess water vapor values, and the RMSE was 1.23.

Near the Atlantic entrance to the Arctic (e.g., Fram Strait) and at depth, the Arctic Ocean is becoming warmer and more saline, which makes it thermally and geochemically more like the Atlantic and Pacific Oceans (Polyakov et al., 2017; Thoman et al., 2020; Tesi et al., 2021). However, in the western Arctic region of this research, freshwater contributions are increasing (Cooper et al., 2022). Moreover, freshwater fluxes through the Bering Strait are rising. For example, Bering Strait mooring data from 1990 to 2015 show an increase of about 0.01 Sv yr⁻¹ in annual mean transport of Pacific waters into the Arctic, which also includes an increased freshwater flux of about 30 km³ yr⁻¹ (Woodgate, 2018). Additionally, freshwater inventories in the Canada Basin have also increased (Rosenblum et al., 2022). Some of this freshwater increase can be attributed to greater RR, but it is also influenced by SIM (Serreze et al., 2006). Reduced Arctic Ocean sea ice coupled with warmer air and water temperatures will increase the amount of Arctic-produced evaporate (i.e., moisture), which amplifies the Arctic water cycle and can be associated with transport of Arctic-sourced moisture throughout the Arctic and also into the sub-Arctic and temperate regions (Bintanja and Selten, 2014; Puntsag et al., 2016; Bailey et al., 2021; Mellat et al., 2021).

However, these Arctic Ocean changes are spatially and temporally complex and driven by interactions between the ocean water and surrounding boundary layer air. The data presented here show the water isotope signature of this complexity in an unprecedented way using water isotope maps, isoscapes, of western Arctic and sub-Arctic seawater and water vapor isotopes. These Arctic water isoscapes show new complexity of Arctic Ocean water isotopes and could lead to a transformational depiction of the Arctic’s water cycle behavior, which can drive novel discoveries of the present, past, and future Arctic and global water cycles.

For the isoscape modeling, depth was the spatial environmental predictor that most influenced d-excess of seawater (Table 4). However, given that stratification of oceanic waters generally prevents waters at depth from influencing the isotopic composition of water at the surface, there is not a direct relationship between depth and the d-excess of surface waters. The d-excess of surface waters is influenced by the SIM contribution to the total freshwater amount, which does vary by depth (Tables 5 and S3). The deeper Chukchi Borderlands surface waters have a higher sea ice cover, which may lead to a greater contribution from SIM and lower values for both salinity and d-excess. Leading up to DOY 195, sea ice cover is about 100% and sea surface temperature is relatively constant, but starting at DOY 195 there is an increase in sea surface temperature, which is followed by lower sea ice cover and a decrease in salinity (Figure 3). These changes suggest an increased SIM contribution to ocean waters. This apparent SIM influences isotopic fractionation and leads to waters with less negative δ¹⁸O values (with less relative change to δ²H), which decreases the d-excess values of the seawater (Bonne et al., 2019); less variability, however, is observed within the deeper Chukchi Borderlands waters.

Seawater d-excess values decrease steadily with more SIM until around DOY 203, when there is a sharp drop in δ²H values, which drives down the d-excess values of seawater to below zero (Figure 3). The d-excess values increase sharply when transitioning from the Chukchi Borderlands to the Shelf as the influence of SIM decreases (which is especially apparent in the less negative δ²H values) and salinity increases. The d-excess values then decrease after passing through the Bering Strait into the warmer waters of the eastern Bering Sea and Gulf of Alaska.

Additionally, there is a plume of relatively high d-excess seawater southwest of St. Lawrence Island (Figure 6). This plume appears to be influenced by Anadyr Current (AC) water in the western Bering Sea. This AC water is colder and saltier than the Alaska Coastal Current (ACC) in the eastern Bering Sea. St. Lawrence Island separates the AC and ACC in the Bering Sea, while the Diomede Islands (middle of the Bering Strait) separate them to the north in the Bering Strait (Kawaguchi et al., 2020). Further, the data from the same time period used for the predictive modeling to generate the isoscapes (Table 1) show these sea surface temperature, salinity, and chlorophyll-a differences to be related to the AC and ACC (Figures S16–S18). For example, the colder and saltier AC waters are apparent on the west side of the Bering Sea and Bering Strait. As all water was flowing north through the Bering Strait into the Arctic during the time of the modeled isoscapes (Zhuk and Kubryakov, 2021), the higher d-excess waters southwest of St. Lawrence Island appear to be the result of the colder and saltier AC waters flowing around St. Lawrence Island (Figures S16 and S17). The temperature and salinity differences between these AC and ACC waters could lead the higher d-excess water to form an eddy southwest of St. Lawrence Island, which is in between both currents. The AC waters then concentrate in the northwestern Bering Sea and on the western side of the Bering Strait before flowing north into the Chukchi and East Siberian Seas and mixing more with the ACC. Additionally, lower d-excess seawater samples (likely from the ACC) were observed (not modeled) passing north through the Bering Strait east of the Diomede Islands (Figure 6), which would be expected given that the islands traditionally separate the colder and saltier AC from the fresher and warmer ACC (Figures S16 and S17). After passing through the Strait, the d-excess values become higher and more representative of the colder AC. Additionally, the variability of the d-excess seawater values (e.g., southwest of St. Lawrence Island) also supports the utility of water isotopes to help understand distinct currents and movements of oceanic water masses (Benetti et al., 2017) beyond what might be possible with traditional temperature and salinity measurements (e.g., Figures S16 and S17) in both modern and historical contexts.

The water isotope parameter d-excess is important for understanding moisture sources, especially in the Arctic (Kurita, 2011; Klein et al., 2015; Klein et al., 2016; Puntsag et al., 2016). However, prior to this study we lacked a clear spatial understanding of the complexities of d-excess in water vapor in the western Arctic Ocean around Alaska. Some previous studies of terrestrial isotopes (e.g., precipitation and ice cores) have linked Arctic moisture sources with higher d-excess values than sub-Arctic and lower latitude water sources (Klein et al., 2015; Klein and Welker, 2016; Puntsag et al., 2016), but this study reveals more complexity to the pattern than previously recognized. For example, Figure 8 shows that the d-excess of Arctic Ocean vapor north of Alaska can range from around –12‰ to 18‰. The higher values in the Canada Basin are potentially influenced by the Beaufort Gyre, which captures Arctic Ocean waters from the Transpolar Current (Serreze et al., 2006; Morison et al., 2012). The Gulf of Alaska also displays a large gradient of water vapor d-excess values (about 10‰ to 24‰). The highest of these d-excess water vapor values are likely influenced by warmer waters (Figure S17); interestingly, some of these values are as high as the those generally associated with Arctic moisture sources (Klein et al., 2015). Additionally, among the lowest d-excess vapor isoscape values in the Arctic are in the East Siberian Sea. These values are potentially influenced by relatively warm and fresh water discharge from Russia (Figures S16 and S17).

The lowest measured d-excess vapor values (about –10‰ to –12‰) occur after DOY 195 (Figure 3). These values could be influenced by a shift to wind from the south (Figure 4), which could transport the relatively enriched δ¹⁸O vapor values from just north of the Bering Strait (Figure 7) that drive down the measured d-excess values. While the water vapor isotope values may be influenced by local evaporation, if local evaporation were a large contributor, the d-excess values could be expected to increase due to greater kinetic fractionation (e.g., more depleted δ¹⁸O vapor values), which is not observed (Figure 3). These low d-excess vapor values persist until a shift to wind from the west (Figure 4), which could bring more depleted δ¹⁸O values (Figure 7) that influence the sharp transition to relatively higher d-excess values. This atmospheric change is also evident in lower barometric pressure and air temperature (Figure 4). Additionally, wind speeds increase and RH decreases around DOY 197 (Figure 4), which potentially increases the local evaporation contribution that would also result in higher d-excess vapor values (Steen-Larsen et al., 2013; Klein et al., 2015).

Water vapor isotope values can be composed of 2 main sources: local and distant. The local sources are the result of evaporation of seawater, while the distant sources are from ambient air parcels that traveled from another location. As the water vapor intake was about 11 m above the ocean surface, we assume that water vapor includes both local evaporation as well as moisture from more distant sources. The proportion of local versus distant water vapor is influenced by factors such as RH and wind speed (Bonne et al., 2019), but for this work we are not attempting to separate these different water vapor sources. Some studies have started to investigate the influence of vertical location on distant moisture sources and local evaporation of ocean water on water vapor isotopes (Thurnherr et al., 2020; Seidl et al., 2023). Future work should collect water vapor isotope data at different heights (including right above the surface), in conjunction with sea surface isotopes, to further investigate these ocean-vapor water isotope dynamics.

While the water vapor isotopes can be influenced by both local and distant water sources, some interesting patterns can still be observed in the difference between surface water and vapor isotopes. In the measured data, a noticeable shift in the differences in δ¹⁸O and δ²H occurs around DOY 203, with the δ¹⁸O difference showing more negative change than δ²H (Figure 4). This difference could show a shift to more local evaporative vapor sources, as the air temperatures drop below sea surface temperatures, which could increase the vapor pressure deficit and fuel more evaporation.

The isoscapes presented here appear to reveal new complexities of Arctic water isotopes, but they do have limitations: they are predictive correlative models, not causative mechanistic models. Possible explanations of mechanisms influencing some components of the isoscapes are presented, but given the detailed complexity indicated by these new machine learning-driven techniques, not all aspects of potential variability in the isoscapes are discussed. Additionally, some of the newly revealed complexities contained in the isoscapes might not have a mechanistic explanation, as these machine learning-based models are fundamentally inductive using correlation, not causation. Or perhaps some of the complexities suggested are erroneous, which is possible as this work presents the first iteration of creating ocean isoscapes. Even though this initial work can be viewed as experimental, the models underpinning the isoscapes were checked internally using cross-validation and spatial predictions of the models were validated using an independent dataset. The models were also validated using the southbound-derived models to predict to the northbound observational data. The validation revealed some consistent spatial clustering of higher errors across the models, but it also showed that the models were accurate for the vast majority of the study area despite temporal variation in the predictor variables.

Similar to prior efforts to understand ocean water isotopes at a larger scale, the machine learning-driven isoscapes presented here can be improved in the future. For example, although we used the best available datasets, if additional spatially continuous variables become available to describe oceanic and atmospheric conditions, future machine learning isoscapes could incorporate more than the 4 variables used here, which would likely improve the results. Our preliminary analyses of water isotopes and ship-based data (covariates) revealed that barometric pressure, wind speed, and wind direction were important predictors for water isotopes, but these variables were not available across a spatially continuous area for the time period of interest, and thus were not used in the models. The availability of more spatially continuous data for these variable could also help to improve the accuracy of the models in locations where they showed consistently greater errors (i.e., between St. Lawrence Island and the Seward Peninsula and near the edge of the Canada Basin). More specifically, the model validation exercise revealed that the d-excess seawater model had the largest spatial range in errors of any isotope model (Figure S10), even though the seawater δ¹⁸O and δ²H models themselves had low errors. Therefore, seawater d-excess values appear to be disproportionately influenced (relative to δ¹⁸O and δ²H) by factors not captured as well by the 4 variables in the current model and might benefit more than the other isotope models from the availability of more spatially continuous data. Lastly, and perhaps most importantly, as machine learning could help to reveal novel connections between water isotopes and different oceanic and atmospheric variables, future work could help to evaluate if the data and relationships presented here are representative of various conditions found across different Arctic locations, times, oceanic depths, and atmospheric conditions and elevations.

These simultaneous and continuous coupled water vapor and seawater analyses and machine learning-driven isoscape models indicate new complexities of the water isotope cycle in the western Arctic and sub-Arctic. Although not previously observed, this complexity is not entirely unexpected, as prior efforts to map ocean water isotopes suggested that future work integrating components of ocean circulation models (e.g., salinity, sea surface temperature, water depth), as we have done here, would create more realistic gradients in ocean water isotopes (LeGrande and Schmidt, 2006). Additionally, the apparent complexity in these high temporal resolution data provide further support to a recent study that showed seasonal and regional variability of periodic Chukchi and Beaufort Sea δ¹⁸O seawater values over a longer time, from 1987 to 2020 (Cooper et al., 2022). Our analyses provide a novel framework to estimate the spatial and temporal complexity of water isotope processes at their oceanic sources and identify different oceanic water masses. Furthermore, with application of our algorithms to 4 oceanic parameters (salinity, water depth, sea surface temperature, and chlorophyll-a), water isotope properties can be estimated with relatively low error ranges (Table 4).

Understanding the spatial and temporal variability of water isotopes at oceanic moisture sources is critical for assessing how liquid water and vapor isotopes are changing along different atmospheric trajectories before deposition (Dütsch et al., 2018). These patterns are important for interpreting past hydroclimate changes with water isotopes from paleoclimate archives (Jones et al., 2012; Klein et al., 2016), as well as using water isotope-enabled models (Fiorella et al., 2021) to estimate the spatial and temporal variation of future shifts across a changing Arctic (Tesi et al., 2021). The water isotope values and processes presented here are important for tracing water in an increasingly rainier Arctic (Bintanja and Andry, 2017; McCrystall et al., 2021) with less sea ice, more local moisture sources (Bintanja and Selten, 2014), and increased discharge of rivers (Box et al., 2019) into a warmer Arctic Ocean (Timmermans and Marshall, 2020).

Data and scripts are available from the NSF Arctic Data Center: Simultaneous ocean water and boundary layer water vapor isotopes in western Arctic, 2016. https://doi.org/10.18739/A2BZ61952

Geospatial analyses performed using ESRI ArcMap software: https://www.esri.com/en-us/arcgis/products/arcgis-desktop/resources

The supplemental files for this article can be found as follows:

Texts S1–S2

Figures S1–S18

Tables S1–S3

This material is based upon work supported by the National Science Foundation Division of Polar Programs under AON (Arctic Observing Network) grant 1504141 and MRI (Major Research Instrumentation) grant 0923571 awarded to JMW. In addition, this research was supported by the U.S. Department of Homeland Security grant DHS-14-ST-061-COE-001A-02 awarded to the University of Alaska Anchorage, and support was provided by JMW’s UArctic Research Chairship. The authors thank the United States Coast Guard Healy crew for help with field logistics. The views and conclusions contained in this document are those of the authors and should not be interpreted as necessarily representing the official policies, either expressed or implied, of the U.S. Department of Homeland Security. This manuscript was improved through the comments of 3 anonymous reviewers and the journal editors, Drs. Miller and Deming.

The authors have no competing interests to declare.

Conceptualization: ESK, JW.

Data curation: ESK.

Formal analysis: ESK, APB.

Funding acquisition: JW, ESK.

Investigation: ESK, JW.

Methodology: ESK, JW, APB.

Project administration: ESK, JW.

Visualization: ESK, AB.

Writing—original draft: ESK.

Writing—review and editing: ESK, JW, APB.

Approved the submission of this manuscript: All authors.

How to cite this article: Klein, ES, Baltensperger, AP, Welker, JM. 2024. Complexity of Arctic Ocean water isotope (δ¹⁸O, δ²H) spatial and temporal patterns revealed with machine learning. Elementa: Science of the Anthropocene 12(1). DOI: https://doi.org/10.1525/elementa.2022.00127

Domain Editor-in-Chief: Jody W. Deming, University of Washington, Seattle, WA, USA

Associate Editor: Lisa A. Miller, Institute of Ocean Sciences, Fisheries and Oceans Canada, Sidney, BC, Canada

Knowledge Domain: Ocean Science