Impact Factor 4.912 | CiteScore 5.0
More on impact ›


Front. Mar. Sci., 09 January 2019 |

Modeled Effect of Coastal Biogeochemical Processes, Climate Variability, and Ocean Acidification on Aragonite Saturation State in the Bering Sea

Darren J. Pilcher1,2*, Danielle M. Naiman2,3, Jessica N. Cross2, Albert J. Hermann1,2, Samantha A. Siedlecki4, Georgina A. Gibson5 and Jeremy T. Mathis6
  • 1Joint Institute for the Study of the Atmosphere and Ocean, University of Washington, Seattle, WA, United States
  • 2NOAA Pacific Marine Environmental Laboratory, NOAA, Seattle, WA, United States
  • 3Jacobs School of Engineering, University of California, San Diego, San Diego, CA, United States
  • 4Department of Marine Sciences, University of Connecticut, Groton, CT, United States
  • 5International Arctic Research Center, University of Alaska Fairbanks, Fairbanks, AK, United States
  • 6Georgetown University, Washington, DC, United States

The Bering Sea is highly vulnerable to ocean acidification (OA) due to naturally cold, poorly buffered waters and ocean mixing processes. Harsh weather conditions within this rapidly changing, geographically remote environment have limited the quantity of carbon chemistry data, thereby hampering efforts to understand underlying spatial-temporal variability and detect long-term trends. We add carbonate chemistry to a regional biogeochemical model of the Bering Sea to explore the underlying mechanisms driving carbon dynamics over a decadal hindcast (2003–2012). The results illustrate that coastal processes generate considerable spatial variability in the biogeochemistry and vulnerability of Bering Sea shelf water to OA. Substantial seasonal biological productivity maintains high supersaturation of aragonite on the outer shelf, whereas riverine freshwater runoff loaded with allochthonous carbon decreases aragonite saturation states (ΩArag) to values below 1 on the inner shelf. Over the entire 2003–2012 model hindcast, annual surface ΩArag decreases by 0.025 – 0.04 units/year due to positive trends in the partial pressure of carbon dioxide (pCO2) in surface waters and dissolved inorganic carbon (DIC). Variability in this trend is driven by an increase in fall phytoplankton productivity and shelf carbon uptake, occurring during a transition from a relatively warm (2003–2005) to cold (2010–2012) temperature regime. Our results illustrate how local biogeochemical processes and climate variability can modify projected rates of OA within a coastal shelf system.


The global ocean has absorbed ∼41% of all anthropogenic carbon dioxide (CO2) emissions resulting from the combustion of fossil fuels and industrial processes (Khatiwala et al., 2009), thereby mitigating some of the warming associated with climate change. When CO2 dissolves in seawater, it can react with water molecules to form a weak acid, carbonic acid (H2CO3). The dissociation of this acid shifts the marine carbonate system, generating a decline in surface pH and reduced carbonate saturation states (Ω), a process referred to as ocean acidification (OA). This process is expected to negatively impact the growth and survival of a wide-range of marine organisms, particularly marine calcifiers (Doney et al., 2009; Kroeker et al., 2013). OA has already generated a pH decline of 0.1 units globally, equivalent to a ∼30% increase in acidity since the beginning of the Industrial Revolution (Orr et al., 2005). High-latitude regions are particularly vulnerable to this expected change as these waters are already naturally low in carbonate ion concentrations, due to ocean circulation patterns and the effect of colder water temperatures increasing CO2 gas solubility (Fabry et al., 2009).

The Bering Sea is one of the most productive marine ecosystems in the world and supports a $3 billion annual U.S. fishery that constitutes approximately 40% of all U.S. commercial fish landings (Wiese et al., 2012). Many of the marine organisms that comprise this ecosystem are vulnerable to OA, posing a large risk to Alaskan commercial and subsistence fisheries (Mathis et al., 2015a). Bivalves such as the Pacific Oyster (Crassostrea gigas) display reduced growth and survival rates under conditions of low pH and Ω (Barton et al., 2012; Waldbusser et al., 2015). Survival rates of larval and juvenile red king crab (Paralithodes camtschaticus) and tanner crab (Chionoecetes bairdi) are also negatively impacted at lower pH (Long et al., 2013; Long et al., 2016), which may lead to significant reductions in catch and profits for these fisheries (Punt et al., 2016). Pelagic fish appear to have relatively greater resilience to direct impacts of OA. For instance, no significant negative impacts from OA have been found for walleye Pollock (Gadus chalcogrammus) larval and juvenile growth rates (Hurst et al., 2012; Hurst et al., 2013). However, pelagic fish may still be vulnerable due to indirect effects of OA on the food web, particularly by reducing prey abundance and behavior. Pteropods are an important prey source for many subarctic and arctic fish species and are already exhibiting shell dissolution due to OA in the Southern Ocean (Orr et al., 2005; Bednarsek et al., 2012).

Many of these biological experiments testing organismal response to OA consist of monitoring the organism in a controlled setting of constant low pH/high CO2 water. However, the natural environment displays substantial spatial-temporal variability, particularly in dynamic coastal settings where the underlying biogeochemistry is impacted by a range of land, ocean, and atmosphere processes. For example, freshwater runoff from terrestrial (Semiletov et al., 2016) and glacial (Evans et al., 2014; Siedlecki et al., 2017) sources, in addition to wind-driven upwelling (Feely et al., 2008) produce localized regions of corrosive water conditions on seasonal timeframes. Thus, while long-term trends in open ocean pH are expected to emerge within a 10–15 years timeframe (Rodgers et al., 2015; Sutton et al., 2017), trends in coastal pH have been substantially more difficult to detect (Duarte et al., 2013). Nonetheless, these coastal regions are a variable but increasing sink of atmospheric CO2 (Laruelle et al., 2018), suggesting that OA is occurring. Furthermore, the interplay of these coastal processes with OA can produce accelerated rates of pH decline compared to the open oceans (Chavez et al., 2017; Carstensen et al., 2018).

The remote geographic location, harsh winters, and large spatial extent of the Bering Sea shelf are additional challenges that observational field campaigns in this region face. Nonetheless, previous work stemming from the 2008 to 2010 Bering Sea Ecosystem Study (BEST) and Bering Sea Integrated Research Project (BSIERP) cruises illustrated distinctive variability, even at low spatial-temporal resolution, and greatly expanded our understanding of the underlying carbon cycle and vulnerability to OA (Wiese et al., 2012; Mathis et al., 2015b). Similar to other high-latitude regions, the Bering Sea shelf has a highly seasonal carbon cycle, with periods of carbon efflux to the atmosphere in winter and carbon influx to the ocean in summer and fall (Cross et al., 2014). Net carbon uptake on annual timeframes results from substantial summer-fall biological productivity, especially along the outer shelf “green belt” region where mixing processes supply off shelf nitrate and coastally derived iron (Mathis et al., 2010; Bates et al., 2011; Cross et al., 2012; Lomas et al., 2012). This organic carbon eventually sinks and is remineralized by bacteria at depth, producing conditions of aragonite undersaturation (ΩArag < 1) in bottom waters for at least several months starting in June (Mathis et al., 2014). Periods of aragonite undersaturation have also been observed in surface waters along the inner Bering Sea Shelf, near sources of freshwater runoff from the Yukon and Kuskokwim rivers (Mathis et al., 2011). This river runoff is supersaturated with terrestrial-derived carbon and supports significant respiration and net heterotrophic conditions in nearshore coastal waters (Striegl et al., 2007; Cross et al., 2012). The Bering Sea shelf also displays substantial interannual variability in sea ice extent and water temperature due to natural climate variability. This variability has recently coalesced into 4–5 years regimes of persistent cold, high sea ice (e.g., 2007–2012) or warm, low sea ice (e.g., 2001–2005) conditions, which have had significant effects on the Bering Sea ecosystem (Stabeno et al., 2012; Heintz et al., 2013; Eisner et al., 2014). This physical variability likely also produces substantial variability in the biogeochemistry and carbonate system, however, the available data are insufficient to resolve this effect (Mathis et al., 2010; Cross et al., 2012).

Computational regional models grounded in observable data have emerged as a useful tool to help elucidate spatial-temporal variability due to their much finer spatial and temporal resolution (Xue et al., 2016; Siedlecki et al., 2017). A regional physical-biological model was developed as part of BEST-BSIERP and utilized in previous work to understand biophysical variability and project changes to the Bering Sea ecosystem under future climate change scenarios (Gibson and Spitz, 2011; Hermann et al., 2016) Here, we add carbonate chemistry to this Bering Sea model and run a 10-year hindcast simulation from 2003 to 2012. The 2003–2012 timeframe was selected as it: (1) provides a full decadal timeframe, (2) encompasses the 2008–2010 period when shipboard measurements of TA and DIC were collected throughout the Bering Sea shelf, and (3) also encompasses both a warm (2003–2005) and a cold (2010–2012) temperature regime. The purpose of this hindcast simulation is to elucidate (1) the underlying mechanisms of spatial variability in ΩArag, (2) the impact of climate variability on the shelf carbon cycle, and (3) how this climate variability impacts the expected rate of OA.

Materials and Methods

Physical Model Description

We use the Regional Ocean Modeling System (ROMS; Shchepetkin and McWilliams, 2005; Haidvogel et al., 2008) configured to the Bering Sea by Hermann et al. (2013) following an implementation for the Northeast Pacific (NEP-5) described by Danielson et al. (2011). This configuration (known as the “Bering10K” model) has ∼10km spatial horizontal resolution with 10 vertical layers (Figure 1). The use of 10 sigma-coordinate vertical layers significantly decreases model run time, while still maintaining many of the physical and biological dynamics produced by a 60-vertical layer implementation (Hermann et al., 2013). The Bering10K model simulates both sea ice (Budgell, 2005) and tidal mixing. The version used here includes recent corrections to the ice thermodynamics terms, which improved the seasonal timing of ice formation and retreat. The Climate Forecast System Reanalysis (CFSR; Saha et al., 2010) and subsequent CFSv2 operational analysis (CFSv2; Saha et al., 2014) provide the atmospheric forcing of air temperature, winds, specific humidity, rainfall, and shortwave and longwave radiation.


Figure 1. Model domain with contoured color shading for model depth. The contour lines for the 50, 100, and 200 m isobaths are also labeled as these represent the generalized boundaries for the inner, middle, and outer Bering Sea shelf domains. White dots signify station locations from the 2008–2010 BEST-BSIERP cruises and the yellow star signifies the locating of the M2 mooring.

Biogeochemical Model Description

We use the biological model developed as part of the BEST program by Gibson and Spitz (2011) (BEST-NPZD). This model contains 12 total pelagic state variables (Figure 2) – large and small phytoplankton, microzooplankton, small and large copepods, euphausiids, jellyfish, iron, nitrate, ammonium, and slow (1 m/day) and fast (10 m/day) sinking detritus. Coupled to this pelagic model is a 3-component ice-biology submodel which simulates ice algae, nitrate, and ammonium within the bottom 2 cm of the ice. Benthic processes are simulated in a 2-component benthic submodel consisting of benthic fauna and benthic detritus. Phytoplankton and detritus that sinks out of the bottom water layer enters the benthic detritus component. Of this flux, 20% is removed to simulate burial and off-shelf transport, based on observations at the M2 mooring (Figure 1; Gibson and Spitz, 2011). An additional 1% is removed due to de-nitrification, which is relatively low in the Bering Sea (Gibson and Spitz, 2011).


Figure 2. Schematic diagram of the model ecosystem structure, reprinted from Gibson and Spitz (2011) with permission from Elsevier.

We add dissolved inorganic carbon (DIC) and total alkalinity (TA) cycling to this existing biological model, following the framework of Fennel et al. (2008) for DIC and Siedlecki et al. (2017) for TA. The model is coded in units of carbon and is converted to nitrogen when necessary (e.g., for nitrate, ammonium, TA) following a constant C:N ratio of 106:16 (Sarmiento and Gruber, 2006). The time evolution of DIC is as follows:


where remin is the detrital (D) remineralization rate, P is the phytoplankton group, μ is the phytoplankton growth rate which is a function of light and nutrients (N), Resp is the respiration rate, Z is the zooplankton group, Graz is the zooplankton grazing rate, wD is the sinking rate of 1 m/day for slow sinking detritus and 10 m/day for fast sinking detritus, [CO2]sol is the CO2 solubility calculated from Weiss (1974), VCO2 is the gas transfer velocity, u is the wind speed, Δz is the vertical thickness of the model grid cell, z = bottom is the bottom water model grid cell, and Sc is the Schmidt number from Wanninkhof (1992). The atmospheric CO2 concentration ([CO2]air) is set to a spatially mean value at monthly time resolution from 2002 to 2012, calculated for the Bering Sea from the NOAA Earth System Research Laboratory’s CarbonTracker 2016 product (CT2016; Peters et al., 2007). Atmospheric CO2 increases from 370 μatm to 393 μatm over the 10-year model hindcast.

Alkalinity cycling follows many of the same processes as DIC, except with the addition of nitrification and the lack of air-sea exchange:


where Nitrif is nitrification. Phytoplankton uptake of nitrate (NO3-) increases TA while the uptake of ammonium (NH4+) decreases TA (Brewer and Goldman, 1976; Goldman and Brewer, 1980; Zeebe and Wolf-Gladrow, 2001; Siedlecki et al., 2017). Nitrification decreases TA due to the production of H+ from the conversion of ammonium to nitrate (Brewer and Goldman, 1976).

DIC and TA are not explicitly simulated within the ice sub-model, however, surface layer DIC and TA concentrations are directly modified by ice algae growth, in order to ensure mass conservation. For example, growth of ice algae decreases surface water DIC, which is then transferred back to the water column following ice algae senescence, grazing, or sinking. A caveat of this simplified approach is that it ignores any additional dynamics between sea ice and TA/DIC, which may have significant local effects in Arctic environments. For example, the formation of ikaite crystals in sea ice can have a greater impact on surface water pCO2 concentrations than ice algae production (Rysgaard et al., 2012; Moreau et al., 2014). To our knowledge, this effect has not been quantified in a 3D biogeochemical model, so we leave this to future work. The presence of sea ice also acts to inhibit the magnitude of the air-sea CO2 flux according to the percentage of the grid cell covered by ice. For example, a grid cell that is 40% covered with sea ice will have the flux magnitude reduced by 40%.

Boundary Conditions

Oceanic horizontal open boundary conditions are interpolated from the CFSR/CFSv2 global model reanalysis/analysis (for all physical state variables) and climatological estimates (for BEST-NPZ state variables). Values for TA and DIC are derived from empirical relationships with salinity, calculated from shipboard data collected as part of the BEST-BSIERP program from 2008 to 2010. Sampling typically occurred between April and July, with an additional late September to early October cruise in 2009. The values are calculated using the following robust linear regressions:


This method was previously employed by Siedlecki et al. (2017) for the Gulf of Alaska, which yielded similar fits, suggesting that the fits are robust throughout the region. All of the fits display strong R2 values (0.75–0.95), except for the DIC fit at salinity values less than 32.6 psu (R2 = 0.21). This fit can be moderately improved by increasing the value of the salinity cutoff between the two fits, however, this leads to a discontinuity where the two regressions intersect. This discontinuity can generate substantial differences in DIC values between two very similar salinities and is therefore less suitable for our purpose of simulating the model boundary conditions.

Freshwater terrestrial runoff reduces salinity in nearshore gridpoints, with monthly climatological magnitude calculated from Dai et al. (2009). The freshwater runoff also has a seasonally varying concentration of DIC and TA following data collected from the mouth of the Yukon River (Table 1). The Yukon River supplies a substantial portion of the total freshwater runoff to the Bering Sea, and is therefore a reasonable representation of total runoff flux of DIC and TA to the Bering Sea. Model water column iron, nitrate, and ammonium are relaxed to climatological values on a 1-year timescale in order to reduce any long-term drift issues that may develop on decadal model timeframes (Hermann et al., 2013). Climatology values are based on observational data collected in the Bering Sea (e.g., BEST-BSIERP), with some additional tuning based on model ecosystem development and multidecadal hindcast analysis (Gibson and Spitz, 2011; Hermann et al., 2013). This climatological relaxation was found to have no discernable impact on 4-year timeframes, but is essentially for multidecadal runs (Hermann et al., 2013). Model surface salinity values are relaxed to monthly values from the Polar Science Center Hydrographic Climatology (PHC; Steele et al., 2001), using the same 1-year relaxation timescales as was used for the nutrients. Model initial conditions for all prognostic variables except TA and DIC are interpolated from a long-term spinup of the Bering10K model, started in 1970 and fully described in Hermann et al. (2013). Initial conditions for TA and DIC are calculated from the salinity regressions.


Table 1. Seasonal concentrations of DIC and TA in freshwater runoff.

Model Skill Assessment

We use observational shipboard data collected during the 2008–2010 BEST-BSIERP cruises for model validation. These data were collected during the spring (April/May) and summer (June/July) seasons, with an additional fall (September/October) sampling in 2009. The sampling covered three cross-shelf transects and one line along the 70 m isobath (Figure 1). Model output is compared to co-located discrete measurements of temperature, salinity, nitrate, DIC, and TA. Values of pCO2 and ΩArag were calculated from TA and DIC by Cross et al. (2012, 2013) using CO2SYS version 1.05 (Lewis and Wallace, 1998). Model comparisons to observed data are summarized with target diagrams (Jolliff et al., 2009). Target diagrams allow for easy visualization of a model’s bias and variance compared to the observed data. Location on the y-axis indicates either a positive (Y > 0) or negative (Y < 0) model bias, while location on the x-axis (calculated as the normalized unbiased RMSD between model and data, multiplied by the sign of modeled minus observed standard deviation) indicates whether the model has a larger (X > 0) or smaller (X < 0) standard deviation compared to the observed data. The radial distance from the origin is also related to the modeling efficiency metric (MEF; Stow et al., 2009). Model fields that lie within the circle have a MEF value of less than 1, indicating better than average model performance as compared to using the mean of the observations. A “perfect” fit would result in all points at X = 0, Y = 0.

Model Hindcast Setup

The model hindcast was run from 2002 to 2012. The year 2002 serves as the model spin-up year for the DIC and TA tracers, therefore we only display results from 2003 to 2012. To ensure sufficient model spin-up, we also tested using two consecutive years of 2002 forcing, but found no substantial subsequent model differences compared to the single year of 2002 forcing. Upon model completion, the model output is re-gridded from the curvilinear, sigma-coordinate grid to a regular latitude-longitude-depth grid to facilitate plotting and comparison to observational data.


Model Performance

The Bering10K model with the BEST-NPZD ecosystem model has been previously validated and utilized in a variety of studies exploring biophysical variability in the Bering Sea in both hindcast and decadal projection simulations (Gibson and Spitz, 2011; Hermann et al., 2013, 2016; Ortiz et al., 2016). The model version presented in this manuscript incorporates updated algorithms for the vertical absorption of shortwave radiation in the water column and sea ice growth described in Hermann et al. (2016). These updated algorithms significantly improved model-data comparisons of water temperature and sea ice coverage. The model does contain a few previously described biases, including a late sea ice retreat and colder water temperatures in the northern regions (Hermann et al., 2016). This late sea ice retreat also delays the timing of the spring bloom and tends to generate somewhat lower total phytoplankton biomass in spring and greater biomass in fall (Ortiz et al., 2016). The model version used in the present study also includes recent corrections to the ice thermodynamics terms, which reduce but do not completely eliminate these biases (K. Kearney, personal communication).

Most model fields compare well with the observations and fall within the MEF < 1 circle, suggesting that the model is skillful at reproducing the observed system (Figure 3). The lack of any significant mean biases suggest that the model is especially useful for simulating mean-state conditions, particularly near the surface. Only surface TA and NO3- values fall outside of the circle, underestimating and overestimating the observed variability respectively. This bias in variability for surface NO3- suggests that the model may be producing more variability in nutrient transport or primary productivity. The biases described above by Ortiz et al. (2016) support the latter mechanism. However, if primary productivity variability is overestimated, we may also expect to see an overestimate in DIC variability. On the contrary, DIC variability is underestimated. Thus, either this mechanism is not present, is not of sufficient magnitude to produce a similar bias in DIC, or is compensated by additional factors governing DIC variability (e.g., TA, air-sea CO2 flux, solubility).


Figure 3. Target diagrams for model variables compared to observed variables from the 2008–2010 BEST-BSIERP cruises for (left) the full water column and (right) surface layer. X-axis is the normalized unbiased RMSD between model and data, multiplied by the sign of the difference between modeled standard deviation and observed standard deviation. Y-axis is the normalized mean bias between model and data. The symbol locations thus illustrate whether the model contains a positive (Y > 0) or negative (Y < 0) bias, whether the model standard deviation is larger (X > 0) or smaller (X < 0) compared to the observations, and the normalized unbiased RMSD between model and data (absolute value of X).

Seasonal and Spatial Patterns

The direction of the air-sea flux of CO2 is determined by the gradient between the pCO2 of the surface ocean (pCO2ocn) of the atmosphere (pCO2atm). Here we define this gradient as the ΔpCO2 (pCO2ocnpCO2atm), where a negative value indicates a flux direction from the atmosphere into the ocean. Thus, regions of negative ΔpCO2 are associated with regions of net ocean carbon uptake.

Values of ΔpCO2 are generally neutral on the Bering Sea shelf throughout most of the winter and spring, but become largely negative in summer and fall (Figure 4). In particular, June–October ΔpCO2 values are strongly undersaturated with respect to the atmosphere, typically ranging from -100 to -200 μatm. However, much of the coastal ocean near western Alaska and Nunivak Island contains positive ΔpCO2 values, particularly near the Yukon River delta. In the deep-water basin off of the Bering Sea shelf, ΔpCO2 remains slightly positive or neutral for nearly the entire year. CO2 flux values display a similar seasonal pattern but with differences in relative magnitude due to the influence of wind speed and sea ice (Figure 5). For example, sea ice inhibits air-sea exchange on the inner and middle coastal shelf during winter and spring months, leading to relatively low CO2 flux values. Conversely, frequent storms and strong wind speeds produce substantial positive CO2 flux values of 20–30 mmolC/m2/day in the deep Bering Sea basin despite relatively low ΔpCO2 values (∼ 25–50 μatm).


Figure 4. (A–D) Seasonal plots of model ΔpCO2 averaged over the 2003–2012 timeframe and (E) time series of area-weighted spatial mean ΔpCO2 for the three shelf domains, denoted by the 50, 100, and 200 m isobath lines shown in Figure 1.


Figure 5. (A–D) Seasonal plots of model air-sea CO2 flux averaged over the 2003–2012 timeframe and (E) time series of area-weighted spatial mean CO2 flux for the three shelf domains. A positive value signifies a flux of carbon from the ocean to the atmosphere.

Values of surface ΩArag are relatively low throughout winter and spring, with values in coastal regions along the western Alaskan shoreline consistently below 1 (Figure 6). ΩArag values increase to 2–3 during summer and fall for most of the Bering Sea shelf, except for coastal regions near Nunivak Island and the Yukon River delta where ΩArag remains less than 2 and less than 1 in a few isolated locations. Bottom water ΩArag values are lower compared to the surface and display a more subdued seasonal cycle compared to surface ΩArag (Figure 7). Coastal waters along the western Alaskan shoreline are persistently undersaturated in aragonite through winter and spring similar to surface waters, however, coastal regions near Siberia are also undersaturated for large portions of the year. Regions along the middle and outer Bering Sea shelf also display an inverse seasonal cycle compared to surface ΩArag. That is, bottom water ΩArag values decrease during the winter-summer transition, as opposed to surface ΩArag values which increase.


Figure 6. (A–D) Seasonal plots of model surface ΩArag averaged over the 2003–2012 timeframe and (E) time series of area-weighted spatial mean ΩArag for the three shelf domains The red dots signify locations where ΩArag < 1.


Figure 7. (A–D) Seasonal plots of model bottom water ΩArag averaged over the 2003–2012 timeframe and (E) time series of area-weighted spatial mean bottom ΩArag for the three shelf domains. The red dots signify locations where ΩArag < 1. Regions deeper than the Bering Sea coastal shelf break (>1500 m) are whited out.

The spatial differences between the inner and outer Bering Sea shelf also manifest on annual timescales, with important implications to the underlying carbon chemistry (Figure 8). Salinity, DIC, and TA are all substantially lower along the inner shelf compared to the rest of the model domain, and particularly within Norton Sound, indicative of outflow from the neighboring Yukon River delta. Plotting the TA/DIC ratio (Figure 8E) illustrates the relative buffering capacity of the water, with lower values indicating a reduced buffering capacity. Values along the inner shelf are relatively lower compared to the outer shelf, indicating this reduced buffering capacity. The spatial pattern of the TA/DIC ratio correlates very strongly with the spatial pattern in ΩArag (r = 0.996, Figure 8F), illustrating the utility of this ratio in identifying regions of corrosive water conditions. Relatively low phytoplankton productivity on the inner shelf and high productivity on the outer shelf (Figure 8B) further contributes to the spatial pattern of this ratio through changes in DIC.


Figure 8. Annual average plots of model (A) surface salinity, (B) primary production integrated over the top 200 m of the water column, (C) surface DIC, (D) surface TA, (E) surface TA/DIC ratio, and (F) surface ΩArag. All variables are averaged over the 2003–2012 timeframe.

Warm and Cold Regimes

The model hindcast timeframe covers two distinct observed temperature regimes: a period of warm, relatively low sea ice extent years (2003–2005) and a period of cold, relatively high sea ice extent years (2010–2012). Surface water temperatures during the cold regime are generally 1–2°C colder throughout most of the shelf region, and up to 4°C colder along the southeastern middle and outer shelf domains (Figures 9A–C). Similarly, sea ice extent is 20–40% greater during the cold regime (Figures 9D–F). Primary production increases for the outer shelf and the northwestern shelf, but remains largely unchanged or slightly decreases through most of the inner shelf and for Bristol Bay (Figures 9G–I). Surface salinity slightly decreases for the outer and middle shelf, but slightly increases for the inner shelf (Figures 9J–L). During the transition from a warm to a cold period, surface ΩArag decreases by ∼ 0.2 units throughout most of the Bering Sea shelf, with locations south of Nunivak Island and in Bristol Bay decreasing by up to 0.3 units (Figures 10A–C). Surface ΩArag also decreases in the Bering Sea Basin, though to a lesser extent (∼0.1 units) compared to the shelf. Surface DIC increases by 20 mmol/m3 for most of the shelf, with a smaller increase of 5–10 mmol/m3 along the shelf break region (Figures 10D–F). However, the area around Nunivak Island remains largely unchanged. Surface TA remains relatively unchanged for the middle and outer coastal shelf, but decreases by ∼20 meq/m3 for large portions of the inner shelf near Nunivak Island (Figures 10G–I). The change in air-sea CO2 flux is more spatially heterogeneous than the change in ΩArag and DIC (Figures 10J–L). Following a similar pattern to primary production (Figures 9G–I), CO2 flux becomes more negative (i.e., increased ocean carbon uptake) over the outer and northwestern shelf regions, but more positive over the inner shelf. In particular, the outer shelf region near the Aleutian Island chain transitions from a neutral or weak carbon source to a significant carbon sink. Carbon efflux from the deep Bering Sea basin also decreases, though the region remains a significant carbon source on annual timeframes.


Figure 9. (left) Model annual mean variables averaged over the warm and (middle) cold temperature regimes, and (right) the difference generated by the transition from warm to cold conditions. (A–C) Surface temperature, (D–F) fraction of grid cell covered by sea ice, (G–I) production integrated over the top 200 m of the water column, and (J–L) surface salinity.


Figure 10. (left) Model annual mean variables averaged over the warm and (middle) cold temperature regimes, and (right) the difference generated by the transition from warm to cold conditions. (A–C) Surface ΩArag, (D–F) surface DIC, (G–I) surface TA, and (J–L) air-sea CO2 flux (where a positive value indicates a flux out of the ocean).

Decadal Trends

Statistically significant trends (p-value < 0.05) of decreasing ΩArag are evident throughout most of the Bering Sea shelf, with magnitudes ranging from 0.025 – 0.04 units/year (Figure 11). These regions also correspond to significant positive trends in DIC of 2.5 – 4.0 mmol/m3/year. Localized regions of negative trends in TA of -2.5 to -5.0 mmol/m3/year are also present in the southeastern inner shelf. Trends in pCO2 and ΔpCO2 are largely positive along the inner and northwestern middle Bering Sea shelf, but large portions of the outer and southeastern middle shelf display a negative trend, particularly for ΔpCO2. These areas of significant negative ΔpCO2 trends correspond to regions of significant negative trends in CO2 flux, indicating an increase in ocean carbon uptake.


Figure 11. Linear trends in annual average surface model (A) ΩArag, (B) DIC, (C) TA, (D) pCO2, (E) ΔpCO2, and (F) air-sea CO2 flux. Linear trends are calculated using the 2003–2012 timeframe. Only trend values that are significant at the 95% confidence level are displayed.

The negative trends in CO2 flux drive positive trends in surface pCO2 and DIC, along with a substantial increase in the total carbon uptake of the entire Bering Sea shelf (Figure 12). Surface ocean pCO2 largely follows the estimated, expected increase in ocean pCO2 based on the air-sea pCO2 disequilibrium calculated in 2002 and the increase in atmospheric CO2 over the model hindcast timeframe (Figure 12A). We then recalculate DIC using CO2SYS and the ocean pCO2 values calculated in Figure 12A (dashed line), along with model TA, salinity, and temperature (Figure 12B). This allows for comparison between the model increase in surface DIC and the expected increase in DIC, calculated from the increase in ocean pCO2 associated with the increase in atmospheric CO2. The goal is to understand the extent to which the change in model pCO2 and DIC is driven by the increase in atmospheric CO2. Figure 12B illustrates that the two values track closely, similar to Figure 12A. Shelf annual total carbon uptake increases from 15.9 TgC/year in 2003 to 26.3 TgC/year in 2012, and total benthic carbon increases slightly from 2.61 TgC in 2003 to 2.78 TgC in 2012 (Figures 12D,E). Surface ΩArag decreases from 1.7 to 1.5, while bottom ΩArag decreases from 1.2 to 1.1. Thus, the decrease at the surface is approximately twice as great as the decrease at the bottom.


Figure 12. Time series plots of annual average model (A) surface ocean pCO2 and (B) surface DIC. The dashed line in (A) represents the expected increase in surface ocean pCO2 based on the model air-sea pCO2 disequilibrium in 2003, combined with the yearly increase in atmospheric CO2. The dashed line in (B) represents the excepted increase in surface ocean DIC, calculated from the expected surface pCO2 values shown in the dashed line of (A), and the model values of total alkalinity, salinity, and temperature. (C) Total ocean carbon uptake and (D) total benthic carbon over the entire Bering Sea shelf. (E) Annual average values of surface ΩArag (solid line) and bottom water ΩArag (dashed line).


The model simulates significant spatial variability in carbonate chemistry on the Bering Sea shelf, resulting from local biogeochemical processes. The inner shelf region contains low values of ΩArag, with some areas near the Yukon Delta and Norton Sound experiencing nearly year-round ΩArag < 1 conditions. These corrosive water conditions result from river inputs of terrestrial carbon from the Yukon River combined with low coastal productivity, which both drive a reduced TA/DIC ratio. Previous studies using shipboard data have also noted relatively low productivity and ΩArag waters along the inner shelf (Mathis et al., 2011; Cross et al., 2014), though the corrosive conditions in Norton Sound are difficult to verify due to a lack of data. Norton Sound contains a relatively small red king crab fishery, but the life cycle of these crab is largely unknown (Norton Sound Economic Development Corporation [NSEDC], 2016). The sound may therefore be a natural analog to future OA conditions and an ideal location for experimental testing on the impacts of OA on red king crab growth, survival, and adaptability. However, observational data are first needed to verify these corrosive water conditions.

In contrast to the inner shelf, the middle and outer shelf regions have greater values of ΩArag, particularly during the summer and fall, due to substantial productivity and minimal influence of freshwater runoff. While this productivity boosts ΩArag in surface waters, carbon remineralization at depth generates bottom water ΩArag values close to 1. This result is supported by observational evidence (Mathis et al., 2014) and is noteworthy because the region is home to substantial tanner and red king crab fisheries, which are negatively impacted at ΩArag values less than 1 (Long et al., 2013; Mathis et al., 2015a; Long et al., 2016). Furthermore, while summer and fall surface ΩArag is highly saturated due to productivity, winter and early spring values are much lower at 1.0–1.5. Thus, while productivity may buffer ΩArag against OA during the summer and fall, biological organisms are still highly vulnerable to winter and early spring conditions as ΩArag values are already near 1.

Modeled total annual carbon uptake for the Bering Sea shelf is 15–25 TgC/year, which is within the 2–67 TgC/year range of previous estimates, but is significantly greater than the most recent estimate of 6.8 TgC/year (Cross et al., 2014). This difference between the latter estimate and the model results from substantially greater modeled carbon uptake in fall. Cross et al. (2014) report average CO2 flux values of -5.0 and -0.48 mmolC/m2/day in September and October respectively, whereas model average flux values are between -10 and -15 mmolC/m2/day for these months. These higher magnitude model flux values are caused by a fall phytoplankton bloom, generated by vertical mixing of nutrients under still-optimal light and temperature conditions (Figure 13). Similar fall blooms are observed on the eastern Bering Sea shelf (Sigler et al., 2014) and reproduced in other biogeochemical models (Cheng et al., 2016), but their prevalence has received less attention compared to the spring bloom. This fall biological drawdown is not apparent to the same extent in the Cross et al. (2014) analysis, though they still calculate a pCO2 drawdown of 50–100 μatm in September and October due to biological productivity. The data are also sparser in September and October compared to the spring and summer months, primarily consisting of a north-south transect line from the Aleutian Islands to the Bering Strait. In this transect, relatively larger ΔpCO2 of ∼150 μatm are located on the outer and middle shelf region south of Nunivak Island, similar to the region of enhanced model primary productivity. It is possible that the finer spatial and temporal resolution of the model is capturing this fall productivity signal that is largely absent in the observed pCO2 data, though the model is also likely biased considering the stronger magnitude fall bloom (Ortiz et al., 2016).


Figure 13. Time series plots of area-weighted spatial mean (A) 30 m column integrated phytoplankton biomass, (B) surface NO3-, and (C) mixed layer depth for model year 2003 and 2012. Phytoplankton biomass is divided into the large phytoplankton (black line), small phytoplankton (blue line), and ice algae (green line) functional groups.

The strengthening of this fall phytoplankton bloom is also the primary mechanism behind the ∼10 TgC/year increase in total carbon uptake from 2003 to 2012. This is a substantial increase in total carbon uptake over only a 10-year timeframe, though we note that this change equates to approximately one-sixth of the total variability between the various observationally based estimates for Bering Sea shelf carbon uptake. Comparing Figures 9, 10 illustrates that the increase in carbon uptake is co-located with the increase in primary production in outer shelf regions. Colder water temperatures will also decrease pCO2 values which can generate enhanced ocean carbon uptake by enhancing negative ΔpCO2 values. Indeed, the regions of greatest increasing carbon uptake (i.e., the southeastern outer shelf) are located where both mechanisms (i.e., increased productivity and colder temperatures) are present. Though there is a slight decreasing trend in annual wind speed (∼0.05 m/s/year) through large portions of the coastal shelf, the trend is not statistically significant. The substantial increase in ice cover during the cold period can also impact shelf carbon uptake by further inhibiting air-sea exchange and by meltwater dilution of surface DIC and TA concentrations. The latter mechanism is evident by the freshening of the outer shelf waters during the cold period (Figure 9L), though the effect is not as apparent with DIC and TA. Implementing mechanistic connections between sea ice formation/melt and DIC and TA concentrations (e.g., Rysgaard et al., 2012) that are currently not included in our model may help resolve this connection for future work.

Previous studies have noted significant variability in temperature and sea ice extent between the 2003–2005 and 2010–2012 timeframes, which likely impacted the marine ecosystem (Stabeno et al., 2012). There is considerable debate regarding whether warmer water temperatures and reduced sea ice extent will increase or decrease primary productivity (Lomas et al., 2012; Brown and Arrigo, 2013; Eisner et al., 2014; Banas et al., 2016; Liu et al., 2016). In the Bering Sea, warmer years have been tied to reduced large copepod and euphausiid populations and subsequent weaker summer/fall lipid content and survival of age-0 walleye Pollock and Pacific cod recruits, which consume the copepods and euphausiids (Hunt et al., 2011; Mueter et al., 2011; Stabeno et al., 2012). However, this response is thought to be tied specifically to the strength of the spring diatom bloom supporting spring zooplankton populations (e.g., copepods and euphausiids) as opposed to annual net productivity. Our model results suggest that warmer years are less productive overall, though with slightly more productivity on the inner shelf. A similar response with zooplankton (i.e., warmer temperatures generating less zooplankton overall, but slightly more on inner shelf) was noted by Hermann et al. (2016) in climate change projections out to 2040 using the Bering10K model. Thus, we expect that the modeled increase in primary productivity from 2003 to 2012 is temporary and will revert to lower values with the return of warm conditions and future warming.

For this study, our primary interest in productivity stems from the biogeochemical connection to the carbonate system. An increase in productivity and carbon export can decrease surface DIC concentrations, and thereby mask the expected increase in surface DIC and corresponding decrease in ΩArag due to OA. However, shelf DIC increases at a rate similar to that expected due to increasing atmospheric CO2 (Figure 12). This increase in DIC accounts for the rapid decrease in ΩArag, but it is perhaps less intuitive as to why the DIC is increasing at a rate comparable to the expected rate from increasing atmospheric CO2, despite the increase in productivity. This discrepancy results from the composition (i.e., small phytoplankton) and fall timing of the increase in productivity. For instance, the cold period drives enhanced negative ΔpCO2 values in fall due to colder water temperature and increased productivity. However, this productivity is in the form of small phytoplankton, which have a significantly slower sinking rate compared to large phytoplankton (0.05 m/day vs. 1.0 m/day, respectively). Thus, this organic carbon is not permanently buried (Figure 12D), but rather remineralized back to DIC during the winter and retained in the water column due to ice inhibition of air-sea exchange. This process combined with the increase in atmospheric CO2 results in a net annual increase in DIC, despite the increased productivity.

It is possible that the model is overestimating the trend in productivity and thereby shelf carbon uptake, which would exacerbate the previously described mechanism. Ortiz et al. (2016) noted that model fall productivity is likely biased high, particularly due to strong small phytoplankton blooms. These fall blooms are further enhanced during the cold period, due to relatively deeper mixed layer depths and enhanced vertical transport of sub-surface nitrate. The colder period is also defined by a comparatively weaker spring diatom bloom, leaving more nitrate in the summer surface layer for small phytoplankton (Figure 13). Model surface DIC does not contain a significant bias though, suggesting a reasonable seasonal biological drawdown of carbon. Surface NO3- is the only model variable that displays much greater variability compared to the observations (Figure 3), though model NO3- is also biased somewhat high compared to the observations, which is not what we would expect if model productivity is biased high. Ammonium inhibition of nitrate uptake could also lead to the higher nitrate concentrations, though this process is modeled to prevent the complete limitation of nitrate uptake at high ammonium concentrations (Gibson and Spitz, 2011).

Over the 10-year model hindcast, statistically significant increasing trends in surface ocean DIC and pCO2, along with a decreasing trend in ΩArag are apparent throughout most of the Bering Sea shelf. These are the long-term changes expected due to OA. That these trends emerge within a 10-year timeframe is surprising considering the highly variable coastal setting and substantial temperature variability. The transition from warm to cold water conditions acts to decrease surface ocean pCO2 values, due to increased gas solubility and primary production. Yet, surface ocean pCO2 values largely increase at a rate comparable to that expected due to increasing atmospheric CO2 (Figures 11, 12). For pCO2, the temperature driven signal and the atmospheric CO2 driven signal act in opposition. Conversely, both colder water temperatures and increasing atmospheric CO2 act to decrease ΩArag values. For ΩArag, the transition from a warm to cold regime acts to accelerate the decreasing trend. Comparing Figures 9C, 10C illustrates that the region of greatest temperature decrease in the southeastern middle shelf corresponds to the region of greatest ΩArag decrease. However, warm conditions returned to the Bering Sea starting in 2014 (Duffy-Anderson et al., 2017), thus this mechanism likely did not continue.

The Community Earth System Model (CESM) projects that, under RCP 8.5, annual average surface ΩArag will decrease by an additional 0.70 units by the end of the century and reach the saturation threshold (i.e., ΩArag < 1) by 2062 (Mathis et al., 2015a). If extrapolated forward in time, our model results suggest a faster timeline, with annual average ΩArag values < 1 by 2040. However, this estimate is highly uncertain due to the simplistic extrapolation method and likely biased by significant variability encapsulated within our model timeframe. Furthermore, observed and modeled spatial variability suggests that this timeframe may vary by up to several decades between different shelf regions. Future work will utilize ESM output to dynamically downscale OA projections using the Bering10K model (Hermann et al., 2016; Wallhead et al., 2017). This method will reduce the uncertainty of these critical threshold years and provide information on their spatial variance throughout the Bering Sea.


We use a regional model of the Bering Sea to illustrate how coastal biogeochemical processes and climate variability impact aragonite saturation state over a 10-year timeframe. Freshwater runoff from the Yukon River reduces surface ΩArag values to < 1 from December-May along the inner shelf region, particularly in Norton Sound. Conversely, substantial primary productivity increases ΩArag values to >2 along the outer shelf and shelf break regions, though ΩArag still decreases to ∼1.5 during winter months. A shift from relatively warm, low sea ice conditions (2003–2005) to relatively cold, high sea ice conditions (2010–2012) produces a significant increase in fall productivity and ocean carbon uptake throughout the outer shelf region. Surface DIC increases at a rate of 2.5–4.0 mmol/m3/year and surface ΩArag decreases at a rate of 0.025 – 0.04 units/year. These trends are statistically significant throughout most of the middle and inner shelf domains, though they are likely enhanced by the transition from the warm to cold period. Nonetheless, these trends are relatively greater than those observed for the open ocean, illustrating how high-latitude coastal shelf regions serve as hotspots of global change. Future work will utilize statistical downscaling to project these changes through 2100 in order to determine critical threshold years for persistent annual ΩArag undersaturation and variability in this projected change between the inner and outer shelf regions.

Author Contributions

DP and JM conceived the project. DP and DN were responsible for the model simulations, analysis, and producing the manuscript figures. JC provided the observational data. AH provided the model forcing fields. AH, SS, and GG provided guidance and support for model development. DP wrote the first draft of the manuscript. All authors provided feedback and contributed to the final manuscript.


Part of this research was performed while DP held a National Research Council Research Associateship award at the Pacific Marine Environmental Laboratory. Funding was provided by the NOAA Arctic Research Program and through the Joint Institute for the Study of the Atmosphere and Ocean (JISAO) NOAA Cooperative Agreement NA15OAR4320063 (2015–2020). Funding for DN was provided by the NOAA Hollings Scholarship Program. This is PMEL contribution number 4802 and JISAO contribution number 2018-0158. Model output data supporting the conclusions of this manuscript will be made available by the authors, without undue reservation, to any qualified researcher.

Conflict of Interest Statement

The authors declare that the research was conducted in the absence of any commercial or financial relationships that could be construed as a potential conflict of interest.


We would like to acknowledge high-performance computing support from Yellowstone (ark:/85065/d7wd3xhc) provided by NCAR’s Computational and Information Systems Laboratory, sponsored by the National Science Foundation. Use of the Yellowstone system was provided by a postdoctoral allocation for DP.


Banas, N. S., Zhang, J., Campbell, R. G., Sambrotto, R. N., Lomas, M. W., Sherr, E., et al. (2016). Spring plankton dynamics in the Eastern Bering Sea, 1971-2050: mechanisms of interannual variability diagnosed with a numerical model. J. Geophys. Res. Oceans. 121, 1476–1501. doi: 10.1002/2015JC011449

CrossRef Full Text | Google Scholar

Barton, A., Hales, B., Waldbusser, G. G., Langdon, C., and Feely, R. A. (2012). The Pacific oyster, Crassostrea gigas, shows negative correlation to naturally elevated carbon dioxide levels: implications for near-term ocean acidification effects. Limnol. Oceanogr. 57, 698–710. doi: 10.4319/lo.2012.57.3.0698

CrossRef Full Text | Google Scholar

Bates, N. R., Mathis, J. T., and Jeffries, M. A. (2011). Air-sea CO2 fluxes on the bering sea shelf. Biogeoscience 8, 1237–1253. doi: 10.5194/bg-8-1237-2011

CrossRef Full Text | Google Scholar

Bednarsek, N., Tarling, G. A., Bakker, D. C. E., Fielding, S., Jones, E. M., Venables, H. J., et al. (2012). Extensive dissolution of live pteropods in the Southern Ocean. Nat. Geosci. 5, 881–885. doi: 10.1038/NGEO1635

CrossRef Full Text | Google Scholar

Brewer, P. G., and Goldman, J. C. (1976). Alkalinity changes generated by phytoplankton growth. Limnol. Oceanogr. 21, 108–117. doi: 10.4319/lo.1976.21.1.0108

CrossRef Full Text | Google Scholar

Brown, Z. W., and Arrigo, K. R. (2013). Sea ice impacts on spring bloom dynamics and net primary production in the Eastern Bering Sea. J. Geophys. Res. Oceans 118, 43–62. doi: 10.1029/2012JC008034

CrossRef Full Text | Google Scholar

Budgell, W. P. (2005). Numerical simulation of ice-ocean variability in the Barents Sea region: towards dynamical downscaling. Ocean Dyn. 55, 370–387. doi: 10.1007/s10236-005-0008-3

CrossRef Full Text | Google Scholar

Carstensen, J., Chierici, M., Gustaffson, B. G., and Gustafsson, E. (2018). Long-term and seasonal trends in estuarine and coastal carbonate systems. Global Biogeochm. Cycles 32, 497–513. doi: 10.1002/2017GB005781

CrossRef Full Text | Google Scholar

Chavez, F. P., Pennington, J. T., Michisaki, R. P., Blum, M., Chavez, G. M., Friederich, J., et al. (2017). Climate variability and change: response of a coastal ocean ecosystem. Oceanography 30, 128–145. doi: 10.5670/oceanog.2017.429

CrossRef Full Text | Google Scholar

Cheng, W., Curchitser, E., Stock, C., Hermann, A., Cokelet, E., Mordy, C., et al. (2016). What processes contribute to the spring and fall bloom co-variability on the Eastern Bering Sea shelf? Deep-Sea Res. II 134, 128–140. doi: 10.1016/j.dsr2.2015.07.009

CrossRef Full Text | Google Scholar

Cross, J. N., Mathis, J. T., and Bates, N. R. (2012). Hydrographic controls on net community production and total organic carbon distributions in the eastern Bering Sea. Deep-Sea Res. II 6, 98–109. doi: 10.1016/j.dsr2.2012.02.003

CrossRef Full Text | Google Scholar

Cross, J. N., Mathis, J. T., Bates, N. R., and Byrne, R. H. (2013). Conservative and non-conservative variations of total alkalinity on the southeastern Bering Sea shelf. Mar. Chem. 154, 110–112. doi: 10.1016/j.marchem.2013.05.012

CrossRef Full Text | Google Scholar

Cross, J. N., Mathis, J. T., Lomas, M. W., Moran, S. B., Baumann, M. S., Shull, D. H., et al. (2014). Integrated assessment of the carbon budget in the southeastern Bering Sea. Deep-Sea Res. II 109, 112–124. doi: 10.1016/j.dsr2.2014.03.003

CrossRef Full Text | Google Scholar

Dai, A., Qian, T., Trenberth, K. E., and Milliman, J. D. (2009). Changes in continental freshwater discharge from 1948-2004. J. Climate 22, 2773–2791. doi: 10.1175/2008JCLI2592.1

CrossRef Full Text | Google Scholar

Danielson, S. L., Curchitser, E. N., Hedstrom, K. S., Weingartner, T. J., and Stabeno, P. J. (2011). On ocean and sea ice modes of variability in the Bering Sea. J. Geophys. Res. 116:C12034. doi: 10.1029/2011JC007389

PubMed Abstract | CrossRef Full Text | Google Scholar

Doney, S. C., Fabry, V. J., Feely, R. A., and Kleypas, J. A. (2009). Ocean acidification: the other CO2 problem. Annu. Rev. Mar. Sci. 1, 169–192. doi: 10.1146/annurev.marine.010908.163834

CrossRef Full Text | Google Scholar

Duarte, C. M., Hendriks, I. E., Moore, T. S., Olsen, Y. S., Steckbauer, A., Ramajo, L., et al. (2013). Is ocean acidification an open-ocean syndrome? Understanding anthropogenic impacts on seawater pH. Estuaries Coasts 36, 221–236. doi: 10.1007/s12237-013-9594-3

CrossRef Full Text | Google Scholar

Duffy-Anderson, J. T., Stabeno, P. J., Siddon, E. C., Andrews, A. G., Cooper, D. W., Eisner, L. B., et al. (2017). Return of warm conditions in the southeastern bering sea: phytoplankton – Fish. PLoS One 12:e0178955. doi: 10.1371/journal.pone.0178955

PubMed Abstract | CrossRef Full Text | Google Scholar

Eisner, L. B., Napp, J. M., Mier, K. L., Pinchuk, A. I., and Andrews, A. G. III. (2014). Climate-mediated changes in zooplankton community structure for the eastern Bering Sea. Deep Sea Res. Part II Top. Stud. Oceanogr. 109, 157–171. doi: 10.1016/j.dsr2.2014.03.004

CrossRef Full Text | Google Scholar

Evans, W., Mathis, J. T., and Cross, J. N. (2014). Calcium carbonate corrosivity in an Alaskan inland sea. Biogeoscience 11, 365–379. doi: 10.5194/bg-11-365-2014

CrossRef Full Text | Google Scholar

Fabry, V. J., McClintock, J. B., Mathis, J. T., and Grebmeier, J. M. (2009). Ocean acidification at high latitudes: the bellwether. Oceanography 22, 160–171.

Google Scholar

Feely, R. A., Sabine, C. L., Hernandez-Ayon, J. M., Ianson, D., and Hales, B. (2008). Evidence for upwelling of corrosive “acidified” water onto the continental shelf. Science 320, 1490–1492. doi: 10.1126/science.1155676

PubMed Abstract | CrossRef Full Text | Google Scholar

Fennel, K., Wilkin, J., Previdi, M., and Najjar, R. (2008). Denitrification effects on air-sea CO2 flux in the coastal ocean: simulations for the northwest North Atlantic. Geophys. Res. Lett. 35:L24608. doi: 10.1029/2008GL036147

CrossRef Full Text | Google Scholar

Gibson, G. A., and Spitz, Y. H. (2011). Impacts of biological parameterization, initial conditions, and environmental forcing on parameter sensitivity and uncertainty in a marine ecosystem model for the Bering Sea. J. Mar. Sys. 88, 214–231. doi: 10.1016/j.jmarsys.2011.04.008

CrossRef Full Text | Google Scholar

Goldman, J. C., and Brewer, P. G. (1980). Effect of nitrogen source and growth rate on phytoplankton-mediated changes in alkalinity. Limnol. Oceanogr. 25, 352–357.

Google Scholar

Haidvogel, D. B., Arango, H., Budgell, W. P., Cornuelle, B. D., Curchitser, E., Di Lorenzo, E., et al. (2008). Ocean forecasting in terrain-following coordinates: formulation and skill assessment of the Regional Ocean Modeling System. J. Comput. Phys. 227:3624. doi: 10.1016/

CrossRef Full Text | Google Scholar

Heintz, R. A., Siddon, E. C., Farley, E. V. Jr., and Napp, J. M. (2013). Correlation between recruitment and fall condition of age-0 pollock (Theragra chalcogramma) from the eastern Bering Sea under varying climate conditions. Deep-Sea Res. II 94, 150–156. doi: 10.1016/j.dsr2.2013.04.006

CrossRef Full Text | Google Scholar

Hermann, A. J., Gibson, G. A., Bond, N. A., Curchitser, E. N., Hedstrom, K., Cheng, W., et al. (2013). A multivariate analysis of observed and modeled biophysical variability on the Bering Sea shelf: Multidecadal hindcasts (1970-2009) and forecasts (2010-2040). Deep-Sea Res. II 94, 121–139. doi: 10.1016/j.dsr2.2013.04.007

CrossRef Full Text | Google Scholar

Hermann, A. J., Gibson, G. A., Bond, N. A., Curchitser, E. N., Hedstrom, K., Cheng, W., et al. (2016). Projected future biophysical states of the Bering Sea. Deep-Sea Res. II 134, 30–47. doi: 10.1016/j.dsr2.2015.11.001

CrossRef Full Text | Google Scholar

Hunt, G. L. Jr., Coyle, K. O., Eisner, L. B., Farley, E. V., Heintz, R. A., et al. (2011). Climate impacts on eastern Bering Sea foodwebs: a synthesis of new data and an assessment of the Oscillating Control Hypothesis. ICES J. Mar. Sci. 68, 1230–1243. doi: 10.1093/icesjms/fsr036

CrossRef Full Text | Google Scholar

Hurst, T. P., Fernandex, E. R., and Mathis, J. T. (2013). Effects of ocean acidification on hatch size and larval growth of walleye Pollock (Theragra chalcogramma). ICES J. Mar. Sci. 70, 812–822. doi: 10.1093/icesjms/fst053

CrossRef Full Text | Google Scholar

Hurst, T. P., Fernandez, E. R., Mathis, J. T., Miller, J. A., Stinson, C. M., and Ahgeak, E. F. (2012). Resiliency of juvenile walleye Pollock to projected levels of ocean acidification. Aquat. Biol. 17, 247–259. doi: 10.3354/ab00483

CrossRef Full Text | Google Scholar

Jolliff, J. K., Kindle, J. C., Shulman, I., Penta, B., Friedrichs, M. A. M., Helber, R., et al. (2009). Summary diagrams for coupled hydrodynamic-ecosystem model skill assessment. J. Mar. Syst. 76, 64–82. doi: 10.1016/j.jmarsys.2008.05.014

CrossRef Full Text | Google Scholar

Khatiwala, S., Primeau, F., and Hall, T. (2009). Reconstruction of the history of anthropogenic CO2 concentrations in the ocean. Nature 462, 346–349. doi: 10.1038/nature08526

PubMed Abstract | CrossRef Full Text | Google Scholar

Kroeker, K. J., Kordas, R. L., Crim, R., Hendriks, I. E., Ramajo, L., Singh, G. S., et al. (2013). Impacts of ocean acidification on marine organisms: quantifying sensitivities and interaction with warming. Glob. Change Biol. 19, 1884–1896. doi: 10.1111/gcb.12179

PubMed Abstract | CrossRef Full Text | Google Scholar

Laruelle, G. G., Cai, W.-J., Hu, Z., Gruber, N., Mackenzie, F. T., and Regnier, P. (2018). Continental shelves as a variable but increasing global sink of atmospheric carbon dioxide. Nat. Comm. 9:454. doi: 10.1038/s41467-017-02738-z

PubMed Abstract | CrossRef Full Text | Google Scholar

Lewis, E. R., and Wallace, D. W. R. (1998). Program Developed for CO2 System Calculations. Rep. BNL-61827. Oak Ridger, TN: U.S. Dep. of Energy, Oak Ridge Natl. Lab., Carbon Dioxide Inf. Anal. Cent.

Google Scholar

Liu, C. L., Zhai, L., Zeeman, S. I., Eisner, L. B., Gann, J. C., Mordy, C. W., et al. (2016). Seasonal and geographic variations in modeled primary production and phytoplankton losses from the mixed layer between warm and cold years on the eastern Bering Sea shelf. Deep-Sea Res. II 134, 141–156. doi: 10.1016/j.dsr2.2016.07.008

CrossRef Full Text | Google Scholar

Lomas, M. W., Moran, S. B., Casey, J. R., Bell, D. W., Tiahlo, M., Whitefield, J., et al. (2012). Spatial and seasonal variability of primary production on the Eastern Bering Sea shelf. Deep-Sea Res. II 65–70, 126–140. doi: 10.1016/j.dsr2.2012.02.010

CrossRef Full Text | Google Scholar

Long, W. C., Swiney, K. M., and Foy, R. J. (2016). Effects of high pCO2 on Tanner crab reproduction and early life history, Part II: carryover effects on larvae from oogenesis and embryogenesis are stronger than direct effects. ICES J. Mar. Sci. 73, 836–848. doi: 10.1093/icesjms/fsv251

CrossRef Full Text | Google Scholar

Long, W. C., Swiney, K. M., Harris, C., Page, H. N., and Foy, R. J. (2013). Effects of ocean acidification on juvenile red king crab (Paralithodes camtschaticus) and tanner crab (Chionoecetes bairdi) Growth, condition, calcification, and survival. PLoS One 8:e60959. doi: 10.1371/journal.pone.0060959

PubMed Abstract | CrossRef Full Text | Google Scholar

Mathis, J. T., Cooley, S. R., Lucey, N., Colt, S., Ekstrom, J., Hurst, T., et al. (2015a). Ocean acidification risk assessment for Alaska’s fishery sector. Prog. Oceanogr. 136, 71–91. doi: 10.1016/j.pocean.2014.07.001

CrossRef Full Text | Google Scholar

Mathis, J. T., Cross, J. N., Evans, W., and Doney, S. C. (2015b). Ocean acidification in the surface waters of the Pacific-Arctic boundary regions. Oceanography 28, 122–135. doi: 10.5670/oceanog.2015.36

CrossRef Full Text | Google Scholar

Mathis, J. T., Cross, J. N., and Bates, N. R. (2011). Coupling primary production and terrestrial runoff to ocean acidification and carbonate mineral suppression in the eastern Bering Sea. J. Geophys. Res. 116:C02030. doi: 10.1029/2010JC006453

CrossRef Full Text | Google Scholar

Mathis, J. T., Cross, J. N., Bates, N. R., Moran, S. B., Lomas, M. W., Mordy, C. W., et al. (2010). Seasonal distribution of dissolved inorganic carbon and net community production on the Bering Sea shelf. Biogeoscience 7, 1769–1787. doi: 10.5194/bg-7-1769-2010

CrossRef Full Text | Google Scholar

Mathis, J. T., Cross, J. N., Monacci, N., Feely, R. A., and Stabeno, P. (2014). Evidence of prolonged aragonite undersaturations in the bottom waters of the southern Bering Sea shelf from autonomous sensors. Deep-Sea Res. II 109, 125–133. doi: 10.1016/j.dsr2.2013.07.019

CrossRef Full Text | Google Scholar

Moreau, S., Vancoppenolle, M., Delille, B., Tison, J.-L., Zhou, J., Kotovitch, M., et al. (2014). Drivers of inorganic carbon dynamics in first-year sea ice: a model study. J. Geophys. Res. Oceans 120, 2121–2128. doi: 10.1002/2014JC010388

CrossRef Full Text | Google Scholar

Mueter, F. J., Bond, N. A., Ianelli, J. N., and Hollowed, A. B. (2011). Expected declines in recruitment of walleye Pollock (Theragra chalcogramma) in the eastern Bering Sea under future climate change. ICES J. Mar. Sci. 68, 1284–1296. doi: 10.1093/icesjms/fsr022

CrossRef Full Text | Google Scholar

Norton Sound Economic Development Corporation [NSEDC] (2016). Annual Report 2016. Anchorage, AK: Norton Sound Economic Development Corporation.

Orr, J. C., Fabry, V. J., Aumont, O., Bopp, L., Doney, S. C., Feely, R. A., et al. (2005). Anthropogenic ocean acidification over the twenty-first century and its impact on calcifying organisms. Nature 437, 681–686. doi: 10.1038/nature04095

PubMed Abstract | CrossRef Full Text | Google Scholar

Ortiz, I., Aydin, K., Hermann, A. J., Gibson, G. A., Punt, A. E., Wiese, F. K., et al. (2016). Climate to fish: synthesizing field work, data and models in a 39-year retrospective analysis of seasonal processes on the eastern Bering Sea shelf and slope. Deep-Sea Res. II 134, 390–412. doi: 10.1016/j.dsr2.2016.07.009

CrossRef Full Text | Google Scholar

Pan-Arctic River Transport of Nutrients, Organic Matter, and Suspended Sediments Project [PARTNERS]. (2010). Arctic River Biogeochemistry Data Set. Available at

Peters, W., Jacobson, A. R., Sweeney, C., Andrews, A. E., Conway, T. J., Masarie, K., et al. (2007). An atmospheric perspective on North American carbon dioxide exchange: CarbonTracker. Proc. Natl. Acad. Sci. U.S.A. 104, 18925–18930. doi: 10.1073/pnas.0708986104

PubMed Abstract | CrossRef Full Text | Google Scholar

Punt, A. E., Foy, R. J., Dalton, M. C., Long, W. C., and Swiney, K. M. (2016). Effects of long-term exposure to ocean acidification conditions on future southern Tanner crab (Chionoecetes bairdi) fisheries management. ICES J. Mar. Sci. 73, 849–864. doi: 10.1093/icesjms/fsv205

CrossRef Full Text | Google Scholar

Rodgers, K. B., Lin, J., and Frölicher, T. L. (2015). Emergence of multiple ocean ecosystem drivers in a large ensemble suite with an Earth system model. Biogeoscience 12, 3301–3320. doi: 10.5194/bg-12-3301-2015

CrossRef Full Text | Google Scholar

Rysgaard, S., Glud, R. N., Lennert, K., Cooper, M., Halden, N., Leakey, R. J. G., et al. (2012). Ikaite crystals in melting sea ice – implications for pCO2 and pH levels in arctic surface waters. Cryosphere 6, 901–908. doi: 10.5194/tc-6-901-2012

CrossRef Full Text | Google Scholar

Saha, S., Moorthi, S., Wu, X., Wang, J., Nadiga, S., Tripp, P., et al. (2010). The NCEP climate forecast system reanalysis. Bull. Am. Meteorol. Soc. 91, 1015–1057. doi: 10.1175/2010BAMS3001.1

CrossRef Full Text | Google Scholar

Saha, S., Moorthi, S., Wu, X., Wang, J., Nadiga, S., Tripp, P., et al. (2014). The NCEP climate forecast system version 2. J. Climate 27, 2185–2208. doi: 10.1175/JCLI-D-12-00823.1

PubMed Abstract | CrossRef Full Text | Google Scholar

Sarmiento, J. L., and Gruber, N. (2006). Ocean Biogeochemical Dynamics. Princeton, NJ: Princeton University Press, 526.

Google Scholar

Semiletov, I., Pipko, I., Gustafsson,Ö., Anderson, L. G., Sergienko, V., Pugach, S., et al. (2016). Acidification of East siberian arctic shelf waters through addition of freshwater and terrestrial carbon. Nat. Geosci. 9, 361–365. doi: 10.1038/NEGO2695

CrossRef Full Text | Google Scholar

Shchepetkin, A. F., and McWilliams, J. C. (2005). The regional oceanic modeling system (ROMS): a split-explicit, free-surface, topography-following-coordinate oceanic model. Ocean Modell. 9, 347–404. doi: 10.1016/j.ocemod.2004.08.002

CrossRef Full Text | Google Scholar

Siedlecki, S. A., Pilcher, D. J., Hermann, A., Coyle, K., and Mathis, J. T. (2017). The importance of freshwater to spatial variability of aragonite saturation state in the Gulf of Alaska. J. Geophys. Res. Oceans 122, 8482–8502. doi: 10.1002/2017JC012791

CrossRef Full Text | Google Scholar

Sigler, M. F., Stabeno, P. J., Eisner, L. B., Napp, J. M., and Mueter, F. J. (2014). Spring and fall phytoplankton blooms in a productive subarctic ecosystem, the eastern Bering Sea, during 1995-2011. Deep-Sea Res. II 109, 71–83. doi: 10.1016/j.dsr2.2013.12.007

CrossRef Full Text | Google Scholar

Stabeno, P. J., Kachel, N. B., Moore, S. E., Napp, J. M., Sigler, M., Yamaguchi, A., et al. (2012). Comparison of warm and cold years on the southeastern Bering Sea shelf and some implications for the ecosystem. Deep-Sea Res. II 6, 31–45. doi: 10.1016/j.dsr2.2012.02.020

CrossRef Full Text | Google Scholar

Steele, M., Morley, R., and Ermold, W. (2001). PHC: a global ocean hydrography with a high-quality arctic ocean. J. Climate 14, 2079–2087. doi: 10.1175/1520-04422001014<2079:PAGOHW<2.0CO;2

CrossRef Full Text | Google Scholar

Stow, C. A., Jolliff, J., McGillicuddy, D. J., Doney, S. C., Allen, J. I., Friedrichs, M. A. M., et al. (2009). Skill assessment for coupled biological/physical models of marine systems. J. Mar. Syst. 76, 4–15. doi: 10.1016/j.jmarsys.2008.03.011

PubMed Abstract | CrossRef Full Text | Google Scholar

Striegl, R. G., Dornblaser, M. M., Aiken, G. R., Wickland, K. P., and Raymond, P. A. (2007). Carbon export and cycling by the Yukon, Tanana, and Porcupine rivers, Alaska, 2001-2005. Water Resour. Res. 43:W02411. doi: 10.1029/2006WR005201

CrossRef Full Text | Google Scholar

Sutton, A. J., Wanninkhof, R., Sabine, C. L., Feely, R. A., Cronin, M. F., and Weller, R. A. (2017). Variability and trends in surface seawater pCO2 and CO2 flux in the Pacific Ocean. Geophys. Res. Lett. 44, 5627–5636. doi: 10.1002/2017GL073814

CrossRef Full Text | Google Scholar

Waldbusser, G. G., Hales, B., Langdon, C. J., Haley, B. A., Schrader, P., Brunner, E. L., et al. (2015). Saturation-state sensitivity of marine bivalve larvae to ocean acidification. Nat. Climate Change 5, 273–280. doi: 10.1038/NCLIMATE2479

CrossRef Full Text | Google Scholar

Wallhead, P. J., Bellerby, R. G. J., Silyakova, A., Slagstad, D., and Polukhin, A. A. (2017). Bottom water acidification and warming on the western eurasian arctic shelves: dynamical downscaling projections. J. Geophys. Res. Oceans 122, 8126–8144. doi: 10.1002/2017JC013232

CrossRef Full Text | Google Scholar

Wanninkhof, R. (1992). Relationship between wind speed and gas exchange over the ocean. J. Geophys. Res. Oceans 97, 7373–7382. doi: 10.1029/92JC00188/

CrossRef Full Text | Google Scholar

Weiss, R. (1974). Carbon dioxide in water and seawater: the solubility of a non-ideal gas. Mar. Chem. 2, 203–215. doi: 10.1016/0304-4203(74)90015-2

CrossRef Full Text | Google Scholar

Wiese, F. K., Van Pelt, T. I., and Wiseman, W. J. (2012). Bering sea linkages. Deep-Sea Res. II 6, 2–5. doi: 10.1016/j.dsr2.2012.03.001

CrossRef Full Text | Google Scholar

Xue, Z., He, R., Fennel, K., Cai, W.-J., Lohrenz, S., Huang, W.-J., et al. (2016). Modeling pCO2 variability in the Gulf of Mexico. Biogeoscience 13, 4359–4377. doi: 10.5194/bg-13-4359-2016

PubMed Abstract | CrossRef Full Text | Google Scholar

Zeebe, R., and Wolf-Gladrow, D. (2001). CO2 in Seawater: Equilibrium, Kinetics, Isotopes, Elsevier Oceanograph Series, 65. Amsterdam: Elsevier.

Keywords: ocean acidification, aragonite saturation, Bering Sea, climate variability, coastal carbon cycling, coastal biogeochemistry, freshwater inputs

Citation: Pilcher DJ, Naiman DM, Cross JN, Hermann AJ, Siedlecki SA, Gibson GA and Mathis JT (2019) Modeled Effect of Coastal Biogeochemical Processes, Climate Variability, and Ocean Acidification on Aragonite Saturation State in the Bering Sea. Front. Mar. Sci. 5:508. doi: 10.3389/fmars.2018.00508

Received: 12 July 2018; Accepted: 19 December 2018;
Published: 09 January 2019.

Edited by:

Andrea J. Fassbender, Monterey Bay Aquarium Research Institute (MBARI), United States

Reviewed by:

Raffaele Bernardello, Barcelona Supercomputing Center, Spain
Diane Lavoie, Department of Fisheries and Oceans (Canada), Canada

Copyright © 2019 Pilcher, Naiman, Cross, Hermann, Siedlecki, Gibson and Mathis. This is an open-access article distributed under the terms of the Creative Commons Attribution License (CC BY). The use, distribution or reproduction in other forums is permitted, provided the original author(s) and the copyright owner(s) are credited and that the original publication in this journal is cited, in accordance with accepted academic practice. No use, distribution or reproduction is permitted which does not comply with these terms.

*Correspondence: Darren J. Pilcher,