2.1. MPAS-Seaice zbgc

The MPAS-Seaice model evolves a set of physical ice variables (including enthalpy, salinity, ice concentration, ice and snow volume, and melt pond volume and area) which characterize the sea ice state and enable assessment of coupled interactions with other Earth system components. Version 1.0 is adapted from CICE (the Los Alamos sea ice model; Hunke and other (Reference Hunke, Lipscomb, Turner, Jeffery and Elliott2013) and Hunke and Dukowicz (Reference Hunke and Dukowicz2002)) for non-regular polygons (Golaz and others, Reference Golaz2019) and implemented on a B-grid with state scalars (sea ice concentration, temperature, salinity and biogeochemical tracers) on cell centers and velocity components defined on cell vertices. MPAS-Seaice uses the Icepack library, version 1.0.2 (Hunke and others, Reference Hunke2018) for its column physics and BGC.

The sea ice physics package and dynamical core is comprised of several submodules representing the following key sea ice processes: (1) mushy layer thermodynamics (with temperature and salinity evolution from Turner and others (Reference Turner, Hunke and Bitz2013)), which defines sea ice and snow growth and melt due to conductive, radiative, salinity and turbulent fluxes; (2) shortwave radiative transfer using the delta-Eddington multiple-scattering scheme of Briegleb and Light (Reference Briegleb and Light2007); (3) sea ice momentum which evolves the ice velocity field assuming an elastic-viscous-plastic rheology (Hunke and Dukowicz, Reference Hunke and Dukowicz1997); (4) incremental remapping transport which describes advection of ice concentration, volume, physical and biogeochemical tracers, and other state variables (Lipscomb and Hunke, Reference Lipscomb and Hunke2004); and (5) a ridging parameterization which transfers sea ice concentration between each of the five sub-grid thickness categories (Lipscomb and others, Reference Lipscomb, Hunke, Maslowski and Jakacki2007).

The zbgc (vertically resolved, biogeochemical) submodule uses information about the sea ice state, particularly the vertical distribution of temperature, salinity and internal shortwave radiation, and changes to the ice and snow thicknesses, to determine vertical boundary flux velocities and internal changes in brine tracer concentrations from sea ice physics processes. Biogeochemical exchanges are modeled as coupled non-linear reaction terms (Appendix: Biogeochemical Reaction Terms).

The zbgc submodule procedure is as follows. At each time step and in each of five thickness categories within a given horizontal grid cell, solve mushy layer thermodynamics to update vertical profiles of ice/snow enthalpy and salinity and compute thermodynamic changes in the ice/snow thicknesses. Next, compute microstructural characteristics (porosity, permeability, gravity drainage diffusivity, Darcy flow velocities) from the ice temperature and salinity profiles and ice boundary growth and melt rates after (Jeffery and others, Reference Jeffery, Hunke and Elliott2011). Use this information to solve for the change in brine level (Appendix: Brine Evolution) relative to the sea ice surface. Then determine the mobile and attached fractions for each biogeochemical tracer (section: Mobile and attached tracers) and solve the vertical transport equation for the mobile tracers (section: zbgc tracer equation) and melt losses of attached tracers. Finally, compute the coupled set of non-linear biogeochemical reaction terms for each brine tracer concentration (Appendix: Biogeochemical Reaction Terms). The remaining sea ice redistribution and horizontal transport processes are then performed while conserving the total brine volume of each biogeochemical tracer. The five sets of biogeochemical variables, one for each subgrid thickness category, are retained at each time step and saved as output data. All the results presented in this work are from the net horizontal grid cell tracer field computed as a category sum weighted by sea ice area concentration.

2.2. The sea ice algal ecosystem

The zbgc model configuration used in this analysis has eight components (see Table 1 for symbols used in this text): three algal groups measured in units of nitrogen (diatoms, small plankton and Phaeocystis sp.), three macro-nutrients (nitrate + nitrite, silicate and ammonium), one micro-nutrient (dissolved iron) and dissolved organic nitrogen. We assume that phosphorus is generally not a limiting macro-nutrient in the sea ice (Deal and others, Reference Deal2011. This is consistent with the findings of Saenz and Arrigo (Reference Saenz and Arrigo2012) for Southern Ocean sea ice. However, we acknowledge observational evidence which suggests potential phosphorus limitation of sea ice algae in the Baltic Sea (Piiparinen and others, Reference Piiparinen, Kuosa and Rintala2010). Conversions to carbon, silica and iron are based on fixed prescribed ratios. All components are two-way coupled with MPAS-O BGC. A schematic of zbgc and biogeochemical interactions with the ocean are depicted in Figure 1. Our choice of ecosystem components is based on Deal and others (Reference Deal2011), which includes only the dominant biogeochemical controls on total sea ice algal chl-a: sea ice diatoms, nitrate, ammonium and silicate. In this study, we expand to eight fields to allow for potentially interesting interactions with the polar marine PP system. For example, though ice diatoms tend to dominate sea ice chl-a observations, we also include the prymnesiophyte, Phaeocystis sp. and small flagellates which can be present in abundance (Lizotte, Reference Lizotte2015), contribute significant fractions to the polar oceanic system and are important producers of dimethyl sulfide. In addition, though sea ice algae are generally not observed to be iron-limited (Lannuzel and others, Reference Lannuzel2016), we include dissolved iron in our component fields because sea ice melt is thought to be an important source of the micro-nutrient for Southern Ocean phytoplankton (Lannuzel and others, Reference Lannuzel2013). Finally, dissolved organic nitrogen is included in sea ice to track production and exchange with the ocean. This feature is essentially diagnostic with respect to sea ice production, though in future versions we plan to add sympagic remineralization.

Fig. 1. Schematic of coupled components in the E3SMv1.1 sea ice–ocean eco-dynamics. Boxes represent biogeochemical fields. Thin arrows are fluxes between components. Thick arrows between ice algal/phytoplankton groups indicate fluxes involving the entire group while the thin arrow represents a flux between silicate and diatoms. The thick open arrow indicates respective fluxes: ice diatoms exchange with ocean diatoms, small flagellates with small phytoplankton, and Phaeocystis sp. with its ocean counterpart. Only shown are ocean biogeochemical components involved in direct fluxes with the sea ice. MPAS-O BGC also includes inorganic carbonate chemistry, phosphate, zooplankton, sinking detrital pools, diazotrophs and coccolithophores assigned as implicit members of the small plankton group.

Table 1. Sea ice biogeochemical tracers

2.3. Mobile and attached tracers

One basic assumption of our ice biogeochemical model is that BGC tracers are only present in the sea ice brine. Consequently, the vertical domain for the zbgc model consists of only brine saturated sea ice. We do not, however, require that all sea ice have non-zero brine volume or that the brine surface level within the ice crystal matrix equal the ice surface. For Arctic sea ice under light snow conditions, the ice surface typically rises above the ocean surface, draining the upper ice of brine. The resulting equilibrium brine level then lies below the upper $\sim 10\percnt$ of the ice. By contrast, in many regions of the Southern Ocean, relatively heavy snow depresses the ice below the ocean surface, driving brine intrusions into the snow layer. In this case, the brine level lies above the ice surface, and so does the vertical extent of the zbgc model. Motion of the brine level (h _b) is fundamentally a result of gravitational forces and depends on the permeability of the ice, changes in the ice bottom and surface boundaries, and surface sources of water: snow and ice meltwater, sea water flooding and liquid precipitation. A full description of the brine level dynamics equation is given in the Appendix: Brine Evolution.

The extent to which bio-chemical tracers are advected by the brine depends upon whether the tracer is dissolved or particulate and whether it has some affinity for the ice crystal/fluid boundary. Ice algal diatoms are well known for their ability to attach to ice crystals by means of frustules. In addition, sea ice algae produce ice-active proteins in exopolysaccharids that potentially increase ice affinity (Boetius and others, Reference Boetius, Anesio, Deming, Mikucki and Rapp2015; Juhl and others, Reference Juhl, Krembs and Meiners2011). Indeed, ice diatoms not only resist brine motion through attachment, but actively migrate within the ice domain. Early estimates suggest diatoms move at speeds of 1.5 cm/d (Welch and Bergmann, Reference Welch and Bergmann1989), while more recent lab studies found speeds of 14.4 cm/d (Aumack and others, Reference Aumack, Juhl and Krembs2014; Krembs and others, Reference Krembs, Gradinger and Spindler2000).

In contrast, the dissolved salts, nitrate and silicate, are treated as purely mobile (Fripiat and others, Reference Fripiat2014) and their transport strongly governed by brine motion. This may be an oversimplification, however, as there is some evidence of enhanced nitrate levels in sea ice in the Norwegian Young sea ICE cruise (N-ICE2015; (Duarte and others, Reference Duarte2017; Fripiat and others, Reference Fripiat, Sigman, Massé and Tison2015)). Still, there is more definitive evidence that ammonium and organic matter produced by phytoplankton have an affinity for the ice matrix (Fripiat and others, Reference Fripiat2014; Müller and others, Reference Müller2013; Verdugo and others, Reference Verdugo2004). This is likely due to the prevalence of colloidal matter including nano- and microgels which adhere to the ice crystals. Similarly, Lannuzel and others (Reference Lannuzel2016) report that dissolved iron is actively extracted from ocean water during sea ice growth and incorporated into the sea ice matrix through adsorption onto particles such as organic ligands.

In general, we assume that the brine concentration of bio-chemical tracers may be present in both a purely mobile phase, c _m, which advects with the brine, and an attached (or stationary) phase, c _s, which adheres to the ice crystals or adsorbing particles. Algal diatom locomotion is treated separately and discussed below. The mobility of small flagellates is assumed to be ineffective against ice melt losses. Transfer between phases is governed by prescribed adsorption/desorption timescales and the assumption that ice growth enables adsorption while ice melt facilitates desorption.

Our treatment is based on the linear adsorption model of Harmand and others (Reference Harmand, Rodier, Sardin and Dodds1996). Sea ice biogeochemical tracers adhere to ice crystals with a time constant of τ_ret, and release with a time constant τ_rel, i.e.

(1)$$\eqalign{{\partial c_{\rm m}\over \partial t} & = -{c_{\rm m}\over \tau_{{\rm ret}}} + {c_{\rm s}\over \tau_{{\rm rel}}}\cr {\partial c_{\rm s}\over \partial t} & = {c_{\rm m}\over \tau_{{\rm ret}}} - {c_{\rm s}\over \tau_{{\rm rel}}}}$$

where c _m is the mobile phase tracer concentration and c _s is the attached phase concentration. The time constants are modeled as crude functions of the ice growth and melt rates defined as changes in the ice thickness (dh/dt). All tracers except algal diatoms follow the simple rule: when dh/dt ≥ 0, then τ _rel → ∞ and τ _ret is finite. For dh/dt < 0, then τ _ret → ∞ and τ _rel is finite. In this way, ice growth promotes adsorption while ice melt enables transitions to the mobile phase. Currently, there is a lack of observational and experimental data to motivate a more sophisticated parameterization. Even with this approach, there are potentially 2N additional unknown parameters representing the adsorption/desorption time constants for an N component biochemical ecosystem. To alleviate the proliferation of unconstrained model parameters, we assume all multi-phase tracers adsorb/desorb in the same manner and specify only two time constants.

The adsorption of tracers is also bound by a saturation concentration which is determined from geometric properties of the ice crystals and adsorbing particles (assumed to be spheres). Johnson and others (Reference Johnson, Blue, Logan and Arnold1995) define the maximum number of adsorbate that can attach to a given collector during porous transport, S _d:

(2)$$S_{\rm d} = {S_{\rm c}\over P_{\rm d}} = {4 \pi r_{\rm c}^2\over \pi r_{\rm d}^2}$$

where S _c is the surface area of the collector crystal and P _d is the projected area of the adsorbate (e.g. diatom). For sea ice diatom attachment to ice crystals, we use a diatom radius of r _d = 4.7 μm corresponding to a cell volume of 435 μm³, which is appropriate for Nitzschia frigida (Olenina and others, Reference Olenina2006), and an ice crystal radius of r _c = 3000 μm (Weeks, Reference Weeks2010). This is then multiplied by the ice crystal number density times the fraction of diatom surface area available for attachment (Johnson and others, Reference Johnson, Blue, Logan and Arnold1995). We calculate the ice crystal number density for each vertical layer as $\lpar 1-\phi \rpar \delta h\pi r_{\rm c}^2/V_{\rm c}$, where V _c is the ice crystal volume assumed to be spherical, δh is the ice thickness in a vertical grid layer, and $\delta h \pi r_{\rm c}^2$ is the volume of an ice layer column of radius r _c. When the stationary concentration reaches the saturation concentration, adsorption stops and the additional material remains in the mobile phase.

Diatoms are distinct from other tracers in that they can move independently. Although we do not model motility as an advective mechanism, we assume that diatoms actively maintain their relative position within the ice, i.e. bottom (interior, upper) algae remain in the bottom (interior, upper) ice, unless melt rates exceed a threshold, v _algae. See Table 2 for parameter values.

Table 2. MPAS-Seaice zbgc reaction parameters

^a Parameters for 1) diatoms, 2) small flagellates and 3) Phaeocystis sp., respectively

2.4. zbgc tracer equation

Purely mobile biogeochemical tracers are tracers which move with the brine and thus, in the absence of biochemical reactions, are transported like salinity. The flux-conserved quantity for purely mobile tracers is the column integrated bulk tracer concentration in the ice column, i.e. C = h _bϕ[c], where h _b is the brine level, ϕ is the porosity and [c] is the in situ brine tracer concentration.

The continuity equation is (Jeffery and others, Reference Jeffery, Hunke and Elliott2011 Eq. (A3)):

(3)$${\partial C\over \partial t} = h_{\rm b} {\partial \over \partial Z}\left\{\left(w + D{\partial\over \partial Z}\right)[ c ] \right\}+ h_{\rm b}\phi R$$

where w is the Darcy velocity (A3), D = D _MLD + ϕD _m is a combination of the effective molecular diffusivity (D _m) scaled by the porosity and the mixing length diffusivity (D _MLD of Jeffery and others (Reference Jeffery, Hunke and Elliott2011)), and R([c]₁, [c]₂, …[c]_N) are biogeochemical reaction terms described in detail in Appendix: Biogeochemical Reaction Terms. In general, the biogeochemical reaction terms are non-linear functions of the N ecosystem tracer brine concentrations.

The vertical coordinate, Z ∈ [Z _bot(t), Z _top(t)], defines a domain bounded by the ice bottom and brine level surface, and so (3) represents a moving boundary problem. The biogeochemical module solves this equation on a fixed domain z ∈ [0, 1] through the coordinate transformation z = (Z _top − Z)/(Z _top − Z _bot). The resulting equation (Jeffery and others, Reference Jeffery, Hunke and Elliott2011 eqn (C16)):

(4)$${\partial C\over \partial t} = {\partial \over \partial z}\left\{\left(v\phi + w_{\rm d} + {D\over h_{\rm b}}{\partial\over \partial z}\right)[ c ] \right\}+ h_{\rm b}\phi R\lpar [ c ] \rpar $$

includes an additional advection term v:

(5)$$v = \lpar z-1\rpar {d z_{{\rm top}}\over {\rm d}t}- z{{\rm d} z_{{\rm bot}}\over {\rm d}t}$$

which explicitly accounts for changes in the sea ice thickness from melt and ice growth.

Solutions to (4) are found using a flux-corrected and positive definite finite element Gelerkin discretization (Brenner and Scott, Reference Brenner and Scott2005; Jeffery and others, Reference Jeffery, Elliott, Hunke and Lipscomb2016).

2.5. Tuning and model parameters

MPAS-Seaice zbgc is part of ICEPACK, the CICE/Consortium column package (Hunke and others, Reference Hunke2018). Parameter values were first tuned in 1-D sea ice only configurations of the Arctic (N-ICE2015 forcing data set; Duarte and others, Reference Duarte2017), and initial tuning of the tracer transport physics was conducted using the Southern Ocean ISPOL forcing data set (Jeffery and Hunke, Reference Jeffery and Hunke2014). Parameters were adjusted using the fully coupled simulations in the HiLAT-E3SM.v0 configuration (Hecht and others, Reference Hecht2019) and tested over 10-year periods during ~200 years of ocean–sea ice spin-up of E3SMv1.1. In all 3-D model configurations, the primary tuning criterion was agreement between observations and simulated ice-integrated chl-a using field measurements from six locations: three in the Arctic (N-ICE2015 (Duarte and others, Reference Duarte2017); Point Barrow (Lee and others, Reference Lee, Whitledge and Kang2008); and BEST (Gradinger and others, Reference Gradinger, Iken and Bluhm2012)) and three in the Southern Ocean (ISPOL (Lannuzel and others, Reference Lannuzel2013); Marguerite Bay (Massom and others, Reference Massom2006); Prydz Bay (Laybourn-Parry, Reference Laybourn-Parry2012)).

For simplicity and to reduce the number of tunable parameters, this version of the model assumes that new sources of nitrogen to the sea ice (primarily nitrate, though there is a small ammonium flux from the ocean) dominate sea ice PP. In other words, in situ remineralized nitrogen is restricted to ammonium sources from the grazing of algae and their natural mortality, and not remineralization of the dissolved organic nitrogen pool, which is difficult to validate. As a further simplification, nitrification rates are set to zero.

As previously discussed, nitrate and silicate are treated as purely mobile tracers. All other tracers attach with an hour timescale and release on a 90-day timescale. Long release timescales allow for high retention of biomass throughout the spring, resisting melt-water flushing of biogeochemical material. Short adsorption timescales allow for maximal seeding of the initial population during sea ice formation and minimize loss of newly produced algal biomass from brine exchange with the ocean. Given that the time-step for MPAS-Seaice is 30 min, 1 h represents a near lower bound for the attachment timescale. Loss of ice mass, however, removes both attached and mobile tracers. The exceptions are the diatoms, the primary algal group, which actively migrate in the sea ice to retain their relative position with respect to the ice bottom. We use a moderate rate of 0.86 cm/d, which is half the rate assumed by Lavoie and others (Reference Lavoie, Denman and Michel2005) and well within the observed migration rates (Welch and Bergmann, Reference Welch and Bergmann1989; Aumack and others, Reference Aumack, Juhl and Krembs2014; Krembs and others, Reference Krembs, Gradinger and Spindler2000).

The three algal groups are distinguished in the following ways. Ice diatoms are the only group with a demand for silica. We use a silicate uptake half-saturation constant of 4 mmol/m³ (Lavoie and others, Reference Lavoie, Denman and Michel2005). The algal growth parameters follow Jin and others (Reference Jin2006) for diatoms: maximum specific algal growth rates of 1.44 1/d, light limitation factor of 0.8 m²/d and light inhibition factor of 0.002 m²/d. For estimates of prymnesiophyte growth parameters, Arrigo and others (Reference Arrigo, Mock and Lizotte2010) report maximum specific growth rates of 0.35–0.71 1/d for P. antarctica, which grows well in Southern Ocean sea ice and open ocean. We use the intermediate value of 0.63 1/d and a light inhibition factor of 0.04 m²/W. There are insufficient observations to constrain growth parameters for the small flagellate group. However, we assume a distinct niche from the other algal groups by specifying the lowest maximum specific growth rate (0.41 1/d) but relatively weak light inhibition (0.001 m²/W). We also assume that small plankton and Phaeocystis sp. experience higher grazing pressure as observed in pelagic phytoplankton communities (Yang and others, Reference Yang, Ha and Kang2015) in addition to mortality, while ice diatoms losses occur through mortality alone. This assumption produces ice algal assemblages dominated by diatoms in both poles (Lizotte, Reference Lizotte2015; Arrigo and others, Reference Arrigo, Worthen, Lizotte, Dixon and Dieckmann1997; Arrigo, Reference Arrigo2003).

The most important tuning parameters are predation rate (f_graze), modeled as a fraction of production; mortality rate (m_pre); the mixing length parameter (l _sk) which governs the strength of gravity drainage fluxes for biogeochemical tracers; and adsorption/retention timescales. Table 2 lists biogeochemical reaction parameters and chosen values. During the tuning process, we found that modeled sea ice algal chl-a tended to be biased low, particularly in the Arctic. As a result, our chosen parameter values skew toward (1) maximizing algal access to ocean nutrients (i.e. high gravity drainage fluxes; l _sk = 20 m), (2) enhancing chl-a retention in the sea ice (short adsorption and long retention timescales) and (3) reducing predation and mortality of algal groups. We understand now that biases in the climate state (discussed further in Simulation results and Discussion) contribute significantly to these underestimations.

3.1. HIST-forcing and CNST-forcing

For the majority of our analysis, we use the CNST-forcing simulation because it represents a relatively stable 157-year time-series of natural system variability. The HIST-forcing simulation, though more representative of the atmospheric forcing conditions on Earth in the last 40 years of the simulation, has a too strong aerosol-cloud effect in response to industrial aerosol emissions which produces enhanced initial cooling of the Arctic (also observed in the E3SMv1.0-DECK simulations of Golaz and others (Reference Golaz2019)). The result is sharp increases in Arctic sea ice volume and extent relative to CNST-forcing beginning in the early 1970s. Figure 2 shows five-year running means of integrated sea ice area and volume for both simulations. Decreasing trends in Arctic ice area and volume in HIST-forcing begin in the mid to late 1970s with values returning to pre-industrial levels over the final 20 years of simulation. Simulated Arctic sea ice area exceeds mean observations (Cavalieri and others, Reference Cavalieri, Parkinson, Gloersen and Zwally1996) by ~1 million km² for CNST-forcing, and between 0.5 and 3 million km² for HIST-forcing. The additional sea ice is primarily an artificial extension in the Labrador Sea likely due to insufficient resolution in ocean circulation around the southern tip of Greenland (Burrows and others, Reference Burrows2020). Southern Ocean sea ice area is more consistent with SSMI observations (Cavalieri and others, Reference Cavalieri, Parkinson, Gloersen and Zwally1996) particularly for HIST-forcing over the final 30 years. Unlike the Arctic, industrial aerosol concentrations do not play a significant role in Southern Ocean ice extent or volume. Extent and volume trends in HIST-forcing are weaker but notably decreasing from the early 1980s. Over the final 30 years, Southern Ocean CNST-forcing extent exceeds HIST-forcing by ~1.2 million km².

Fig. 2. Five-year running means of hemispheric integrated (a,b) sea ice area (million km²), (c,d) sea ice volume (1000 km³) and (e,f) net sea ice algal primary production (PP, TgC/a) for CNST-forcing (blue) and HIST-forcing (red) simulations. Northern hemisphere averages are shown in the left column while Southern hemisphere are on the right. Only grid cells with at least 15% sea ice concentration are included in the averaging. Also shown is the 1979–present mean sea ice area (black dashed line, (Cavalieri and others, Reference Cavalieri, Parkinson, Gloersen and Zwally1996)) and PIOMAS estimate of Arctic sea ice volume (black solid line, (Schweiger and others, Reference Schweiger2011)).

Mean annual sea ice PP in the Arctic is on the low side of previous estimates even with too extensive mean sea ice. This suggests that both HIST-forcing and CNST-forcing likely underestimate per area PP rates. HIST-forcing mean values fluctuate between ~5.8 and 9 TgC/a while CNST-forcing has generally lower values ranging between 5 and 7 TgC/a. Previous field estimates vary by an order of magnitude with 6 TgC/a (Subba-Rau and Platt, Reference Subba-Rau and Platt1984) and from 9 to 73 TgC/a (Legendre and others, Reference Legendre1992). Modeling studies estimate 10.1 TgC/a (Deal and others, Reference Deal2011) to 36 TgC/a (Dupont, Reference Dupont2012). Increases in primary production in HIST-forcing relative to CNST-forcing occur during the years of high Arctic sea ice area in the 1970s and 1980s. However, we also see slightly higher net sea ice PP per ice area in the last 30 years of HIST-forcing (570–650 mgC/m²/a) relative to CNST-forcing (500–610 mgC/m²/a). This suggests that the regions in HIST-forcing of more extensive sea ice are also regions of greater sea ice primary production density.

Five-year running mean values of Southern Ocean sea ice primary production range between 11 and 20 TgC/a. These values are again low compared to previous estimates from field observations 63–70 TgC/a (Legendre and others, Reference Legendre1992) though more consistent with the range in model estimates: 23.7 (Saenz and Arrigo, Reference Saenz and Arrigo2014) to 35.7 TgC/a (Arrigo and others, Reference Arrigo, Worthen, Lizotte, Dixon and Dieckmann1997). CNST-forcing and HIST-forcing estimates are similar, and differences generally reflect variability in sea ice area. Average PP per ice area is similar in both simulations, ranging from 1400 to 1700 mgC/m²/a.

Sea ice PP observations are relatively scarce and not sufficient for ice algal model validation on a global scale. Rather, sea ice core measurements of algal chl-a offer the most extensive comparison. Figures 3 and 4 show seasonal cycles of column integrated sea ice algal chl-a from HIST-forcing and CNST-forcing simulations at several locations in the Arctic and Southern Ocean, respectively, contrasted with core measurements. The shaded regions in the Arctic figure indicate the range in simulated monthly mean over the 157-year CNST-forcing time-series and the final 20 years of HIST-forcing. Individual cores are indicated as well as observational monthly means and standard deviations. For the Arctic, core observations are from (a) the Bering Sea (Gradinger and others, Reference Gradinger, Iken and Bluhm2012), (b) Baffin Bay (Michel and others, Reference Michel, Nielsen, Nozais and Gosselin2002), (c) the Norwegian Sea (Duarte and others, Reference Duarte2017), (d) Point Barrow (Lee and others, Reference Lee, Whitledge and Kang2008 and C. Deal, pers. comm.), (e) the Lincoln Sea (Lange and others, Reference Lange2015), (f) the Canadian Arctic (Melnikov and others, Reference Melnikov, Kolosova, Welch and Zhitina2002) and (g,h) the Beaufort Sea (Horner and Schrader, Reference Horner and Schrader1982; Alou-Font and others, Reference Alou-Font, Mundy, Roy, Gosselin and Agusti2013, respectively). Simulated chl-a is most consistent with observations in the Bering Sea and Point Barrow, which is perhaps not surprising given that both these locations were used in model tuning. Modeled chl-a is within the range of observations in the Baffin Bay and Beaufort Seas, though with some discrepancies in the seasonal cycle. In the Norwegian Sea, model results suggest values as high as 10 mg/m², nearly an order of magnitude larger than observations. However, the reported measurements were taken from a newly formed lead and perhaps do not generally represent sea ice in this region (Duarte and others, Reference Duarte2017). In the Lincoln Sea and Canadian Arctic, the model overpredicts chl-a concentrations, though by only 1–3 mg/m². This is because observations in these locations indicate very low chl-a (< 3 mg/m²). In fact, the peak in observations is more consistent with the modeled minima. Notably, for the Canadian Arctic, the change in chl-a over the seasonal cycle is similar to observations. This may indicate that modeled losses of algae in the spring and summer are too weak. We note here that low simulated primary production is potentially incongruous with reasonably well simulated chl-a. We discuss the implications of this result in detail in the Discussion section.

Fig. 3. Monthly simulated sea ice algal chlorophyll-a concentrations (mg/m²) in CNST-forcing (blue; averaged over 157 years) and HIST-forcing (red; averaged over the last 20 years) plotted against year day for eight locations (a–h) in the Arctic. Shaded regions indicate the range in monthly means. Also shown are core observations: both individual measurements (open circles) and monthly means (solid diamonds with ±1 std). Data references are in text.

Fig. 4. Monthly simulated sea ice algal chlorophyll-a concentrations (mg/m²) in CNST-forcing (blue line; averaged over 157 years) and HIST-forcing (red line; averaged over the last 20 years) plotted against day of year for all grid points in six regions (a–f) (as defined in Meiners and others (Reference Meiners2012)). Shaded regions represent the range in the monthly mean. Also shown are time-series averages of the CNST-forcing monthly mean chlorophyll-a computed at grid cells corresponding to observed core locations for each region (light blue squares with 1 std error bars). Core observations are depicted as black symbols; both individual core measurements (open circles) and monthly means (solid diamonds with ± 1 std) are shown. Data references are in text.

The Southern Ocean observations (Fig. 4) are from the ASPECT data set (Meiners, Reference Meiners2013) and represent a compilation of field campaigns. The red (HIST-forcing) and blue (CNST-forcing) curves and shaded regions are averages over all grid cells in six regions defined as sectors south of 50°S and bounded by the following latitudinal lines (Meiners and others, Reference Meiners2012): (a) the Western Weddell (300 to 315°E), (b) Eastern Weddell (− 45 to 20°E), (c) Indian Sector (20 to90°E), (d) Pacific Sector (90 to 160°E), (e) Ross Sea (160 to 230°E) and (f) Bellingshausen/Amundsen Seas (230 to 300°E). The shaded regions are the range in monthly means over the time-series. Also shown in light blue symbols with 1 std error bars is the monthly mean CNST-forcing climatology, using only grid cells corresponding to core measurement locations for each region. In general, the two simulations are quite similar. In addition, simulated mean climatology generated from grid cells at core locations is consistent with that using all grid points except perhaps in the Indian Sector, where temporal coverage of observations is sparse. Agreement between model and observations is better in the Southern Ocean than the Arctic, though admittedly averaging over a wider spatial region smooths over heterogeneity in the data. Notably, simulated climatology from grid cells at core locations (light blue squares) is not more consistent with observations than the regional average except for perhaps during a couple months of the year in the Pacific Sector and Western/Eastern Weddell. Overall, modeled monthly mean Southern Ocean chl-a is somewhat overestimated except perhaps in the Bellingshausen/Amundsen sector.

With respect to the modeled pre-industrial and historical time periods, differences in sea ice chl-a at observed locations and sectors are smaller than differences between simulations and observations. Several locations show possible increases in peak biomass over the last 20 years of HIST-forcing compared to pre-industrial: the Norwegian and Lincoln Seas, Baffin Bay, the Beaufort and perhaps the Pacific Sector, Ross Sea and Bellingshausen/Amundsen sectors. However, all increases are well within the range of CNST-forcing monthly means (Figs. 3 and 4). For this reason, we focus the remainder of our analysis on the CNST-forcing simulation with the understanding that the primary discrepancies between model and observations are likely not due to the transient forcing conditions of the industrial historical period.

3.2. Ocean and sea ice state

Based on the results from the previous section, the E3SMv1.1 model simulates low values of net annual sea ice PP in both polar regions. In order to understand the processes behind these values, it is critical to first assess the quality of the modeled polar climate state as it pertains to the sea ice algal ecosystem. Figures 5 and 6 show contours of the mean surface ocean nutrients, nitrate and silicate, respectively, in the Arctic and Southern Oceans. The upper row depicts seasonal means in Arctic winter/austral summer (JFM, January–March) and summer/austral winter (JAS, July–September). Figure 5a is notable for the low values in surface nitrate during winter and early spring, particularly in the Chukchi and Canadian Arctic waters. The second row are differences between World Ocean Atlas observations (Garcia and others, Reference Garcia2009) and modeled values. The large regions of red shading in winter and summer in both figures and in both poles indicate systemic underestimations in modeled surface ocean macro-nutrients.

Fig. 5. Upper row are mean surface ocean nitrate concentrations (mmol/m³) in CNST-forcing. Columns (a,c) are January–March averages while (b,d) are averages over July–September. The second row shows differences between observations (Garcia and others, Reference Garcia2009) and model results.

Fig. 6. Upper row are mean surface ocean silicate concentrations (mmol/m³) in CNST-forcing. Columns (a,c) are in January–March averages while (b,d) are averages over July–September. The second row are differences between observations (Garcia and others, Reference Garcia2009) and model results.

Biases in mixed layer depth (MLD) suggest that insufficient mixing may be partly to blame. Figure 7a indicates underestimations of up to 80 m in MLD in winter North Atlantic waters and, to a lesser degree, in the North Pacific. Both regions are important for transporting nutrient-rich waters into the Arctic, particularly through the Barents Sea Opening and the Bering Strait (Torres-Valdés and others, Reference Torres-Valdés2013). Insufficient winter mixing may contribute to a reduced resupply of upper ocean nutrients to Arctic waters.

Fig. 7. Upper row is ocean mixed layer depths (m) in CNST-forcing for the Arctic (a,b) and Southern Ocean (c,d). Averages are (a) January–March, (b) annual, (c) June–August and (d) annual. The second row are differences between observations (Holte and others, Reference Holte, Talley, Gilson and Roemmich2017) and model output.

Figures 8 and 9 contrast the winter and summer sea ice concentration and thickness in CNST-forcing and observations (Comiso, Reference Comiso2012; Yi and Zwally, Reference Yi and Zwally2009, respectively). Figure 8a shows clearly the additional sea ice tongue of the winter simulation extending into the Labrador Sea. This feature is responsible for the overestimation in total Arctic sea ice area discussed in the previous section. Figures 9a and b are perhaps more revealing. Underestimations in simulated ice thickness along the Alaskan shelf, Canadian Archipelago and in the Barents sea are compensated by excessive ice in the central Arctic and Siberian shelf. Together these features indicate a displaced Beaufort High and insufficient transport of the Bering Sea inflow to the Canadian Arctic.

Fig. 8. Upper row is sea ice concentration for the Arctic (a,b) and the Southern Ocean (c,d) in CNST-forcing simulations. Averages are (a,c) January–March and (b,d) July–September. The second row plots are differences between observations (Comiso and others, Reference Comiso, Cavalieri, Parkinson and Gloersen1997) and model output.

Fig. 9. Upper row is sea ice thickness (m) in CNST-forcing for the Arctic (a,b) and Southern Ocean (c,d). Averages are (a,d) February–March and (b,d) October–November. The second row figures are differences between observations (Yi and Zwally, Reference Yi and Zwally2009) and model output.

There are also biases in Southern Ocean sea ice concentration (Figs 8c and d) and thickness (Figs 9c and d) that impact sea ice algal distributions. In general, modeled winter sea ice is thicker close to the Antarctic continent and less extensive, except in the Bellingshausen/Amundsen and Ross sectors. Summer sea ice, particularly along the western Penninsula and Bellingshausen/Amundsen sector, is far less extensive and much thinner. The exception is in the Weddell and Ross seas where simulated summer sea ice is generally thicker along the coast.

3.3. Net sea ice primary production

Sea ice algal production is controlled by light and nutrient availability and mediated by temperature and salinity. After the polar night, variable snow cover may reduce internal shortwave penetration, limiting production particularly in the bottom ice layers. Restricted access to upper ocean nutrients also places limitations on ice algal production. This might be a result of ice physical processes or simply low nitrate and/or silicate concentrations in the upper ocean.

In order to assess the impact of the upper ocean nutrient biases on simulated ice algal production, it is useful to define upper bounds for sea ice algal growth rates based on observational variables of the polar climate state. We then contrast with bounds defined equivalently but using CNST-forcing observables. Observations–estimates of a maximum bound for annual mean algal growth can be derived using monthly climatology of (1) gridded observations of downwelling shortwave radiation (CORE 2; Large and Yeager, Reference Large and Yeager2009), (2) upper ocean nitrate and silicate concentrations (World Ocean Atlas; Garcia and others, Reference Garcia2009) and (3) sea ice extent (SSMI; Comiso and others, Reference Comiso, Cavalieri, Parkinson and Gloersen1997). We assume that the maximum nitrate/silicate concentration available for ice algal growth is given by the surface ocean concentration. Similarly, the maximum photosynthetically available radiation (PAR) is the mean monthly incident PAR at the ice surface. We then compute monthly mean ice algal growth rates for gridded Arctic and Southern Ocean locations within the 15% sea ice extent contour using the limitation functions defined in Appendix: Algal Growth and Nutrient Uptake.

The resulting maximal bound for annual average algal growth rate in the Arctic is plotted in Figure 10 for (a) model based (CNST-forcing) and (b) observation based. It is apparent that from the perspective of sea ice algae, our mean modeled polar environment allows for at most 35% of the maximum specific algal growth rate (μ _max) in a very narrow region of the Sea of Okhotsk and generally far less than 20%. Notably, the Beaufort, Canadian Archipelago and central Arctic see mean values less than 5% of μ _max. In contrast, observations suggest that mean annual growth rates as high as 60% of μ _max are possible in the Sea of Okhotsk, and typical values in the Bering, Chukchi and Canadian Archipelago range from 20 to 50% of μ _max.

Fig. 10. Upper row is the bound for mean algal growth rate (in units of maximum growth rate, μ _max) which depends on ocean surface nutrients, sea ice extent and incoming shortwave radiation for (a) CNST-forcing and (b) observations (references in text). The second row indicates the frequency (fraction of months per year) that the growth bound is determined by surface ocean nitrate concentrations. Light blue line indicates the mean 15% sea ice concentration.

Also shown in Figure 10 (second row) are contour plots of the frequency (fraction of months per year) with which the maximal bound is determined by ocean nitrate concentrations. Given that light necessarily limits growth during the polar night, values of 0.5, which indicate that ocean surface nitrate limits the growth bound for 50% of the year, are extreme. The contours of nitrate limitation frequency based on model output indicate extreme nitrate limitation of the maximal bound throughout the algal growing season. This is reflected in depressed maximal growth rates. The observationally based metric also indicates significant nitrate limitation though not as strong. However, this does not necessarily imply that the computed growth rates from observations are severely depressed. Indeed, nitrate may limit the maximum growth rate in the observational metric, but concentrations are still sufficient to allow moderate (20% of μ _max) to high (55% of μ _max) algal growth rates.

We acknowledge that the maximal growth rate bound as defined above, though useful for contrasting modeled and observed polar states, is a crude measure of sea ice algal growth. Certainly, variability in snow and ice thickness, melt-pond area and in situ nutrient concentrations play significant roles in determining sympagic algal growth rates. Similarly, nitrate limitation of the maximal growth bound does not mean that nitrate is necessarily limiting to the actual algal growth rate. However in the Arctic at least, we do find a significant correlation between the mean net annual ice algal PP of the fully coupled simulation and the maximal growth bound. Conversely, no such pattern correlation is found between mean annual snow thickness and Arctic sea ice primary production. Figure 11a is the total annual sea ice PP averaged over 157 years of CNST-forcing. There are notable similarities with Figure 10a, and indeed the pattern correlation coefficient with the maximal growth rate bound is remarkably high, 0.73 with a t-value of 175. We now use this relationship to improve our estimate of net Arctic PP given the observation-based maximal growth bound. The regressed estimate is shown in Figure 11b. Peak values are similar to the modeled result, but there is far more production along the Canadian shelf and Archipelago, Central Arctic and Bering Sea. Total Arctic PP for the regressed estimate is 60.7 TgC/a, about an order of magnitude larger than the simulated result (Fig. 2) and consistent with the higher estimates of Legendre and others (Reference Legendre1992). The regressed value may be viewed as a model derived estimate of global Arctic sea ice PP that adjusts for ocean nutrient biases in the simulated state. However, it assumes that the pattern correlation between the maximal growth rate bound and PP continues to hold under a state of reduced nitrate limitation.

Fig. 11. Mean total annual Arctic primary production (gC/m²/a) in CNST-forcing (a) is significantly correlated with the maximal growth bound of Figure 10. Figure (b) is an estimate of Arctic primary production using the regression coefficient from (a) (r = 0.73) and the observations based on maximal growth bound. Total integrated Arctic ice algal PP from (b) is 60.7 TgC/a. Light blue line indicates the mean 15% sea ice concentration.

A similar analysis for the Southern Ocean is not nearly so useful. Figure 12 shows (a) model-based and (b) observation-based maximal ice algal growth rate bounds for the Southern Ocean. The model-based bound is generally lower than that of observations, ranging between 20 and 80% of μ _max, while observations suggest values exceed 50% μ _max within the 15% sea ice concentration contour. Certainly, the Southern Ocean modeled polar state appears to be more conducive to ice algal growth than in the Arctic. However, unlike the Arctic, there is no pattern correlation between the maximal growth bound as defined here and net Southern Ocean annual ice PP (r = 0.05). This is because ocean macro-nutrients do not limit the maximal growth bound in the Southern Ocean for either modeled or observed states for any significant fraction of the year (not shown). Rather, incident PAR is always limiting. We find that contrasting algal growth rates based on surface PAR between model and observations, without information about snow and sea ice thickness, is simply too crude a measure. However, the extreme light limitation of the growth bound suggests that snow thickness patterns may provide some clues to regional variability in Southern Ocean ice PP.

Fig. 12. Contours of the upper bound for algal growth rate (in units of maximum growth rate, μ _max) based on mean annual ocean surface nitrate, silicate and incident PAR from (a) CNST-forcing and (b) observations (references in text). Light blue line indicates the mean 15% sea ice concentration.

Figure 13a shows the CNST-forcing mean total annual sea ice PP in the Southern Ocean. High production zones generally follow a thin band along the coastline. Figure 13b shows the simulated mean annual snow thickness. There appear to be two regimes: (1) a negatively correlated coastal zone where decreases in snow thickness coincide with increased PP and (2) the pack ice, where snow thickness and net PP in a grid cell (weighted by sea ice area fraction) both correlate positively with sea ice concentration.

Fig. 13. (a) Mean total annual primary production in CNST-forcing for the Southern Ocean is uncorrelated with the maximal algal growth bound of Figure 12. However, there are regions, particularly in the coastal zones, where (b) mean snow thickness (m) is negatively correlated with PP. Light blue line indicates the mean 15% sea ice concentration.

The scatter plots of Figure 14 support these relationships. Figure 14a contains only Southern Ocean sea ice grid points in the deep ocean (ocean depths > 1200 m) which we refer to as the pack ice region. Net annual ice algal PP per sea ice area fraction (note, this is now unweighted by sea ice area fraction) is plotted against mean sea ice thickness and color coded to indicate the mean sea ice concentration of each grid point. For most of the pack ice, snow thickness is positively correlated with sea ice concentration and weakly positive or uncorrelated with primary production. The exception is perennial pack ice where mean annual sea ice concentrations exceed 70–80%. In contrast, for the near coastal sea ice (ocean depths < 1200 m) in Figure 14b, negative correlations between PP and snow thickness are clearly evident when mean sea ice concentrations exceed ~60%. For seasonal ice of the coastal regions, snow thickness again becomes positively correlated with sea ice concentration, and correlations with PP are much weaker. We see similar but notably weaker correlations between sea ice thickness and ice PP (not shown), though this is likely a result of the positive correlation between snow and sea ice thicknesses.

Fig. 14. Net annual ice algal primary production per sea ice concentration in a grid cell is plotted against mean sea ice snow thickness and color coded to indicate the mean sea ice concentration. Figure (a) includes only grid cells with ocean depths exceeding 1200 m and (b) includes near coastal grid points with ocean depths < 1200 m. Black lines of equivalent slope (− 12.8 gC/m²/a/m of snow) are shown for reference. Figure inserts are scatter plots of sea ice area concentration versus snow thickness for (a) the sea ice pack and (b) near coastal grid cells. Also shown are correlation coefficients significant at the 99% confidence level.

3.4. Temporal variability in sea ice primary production and chlorophyll-a

In this section, we compare modeled seasonal cycles of sea ice primary production and ice algal chl-a in CNST-forcing with observations. Figures 15 and 16 contrast monthly average PP. For the Arctic, we use the compiled results from Leu and others (Reference Leu2015) plotted as gray shaded bars to denote the time-span and range of observed values. Southern Ocean sea ice PP measurements (Arrigo and others, Reference Arrigo, Mock and Lizotte2010) are very sparse, spanning only three months of the year. Again we use gray shaded bars to indicate the span of observed values. For model results, we compare seasonal cycles of PP using output from the data locations only (blue lines and shading) and for the entire polar region with at least a 15% sea ice concentration (pink lines and shading). Shaded regions indicate the range in monthly mean simulated PP. Also shown is the mean and standard deviation for each month over the CNST-forcing time-series.

Fig. 15. Comparison of monthly average Arctic sea ice primary production in CNST-forcing and observations (Leu and others, Reference Leu2015). Gray shaded regions denote the range of observed values. Colored shaded regions are the range in monthly averaged simulated sea ice primary production for grid cells matching observed locations (blue) and for all grid cells (red). Symbols indicate the mean value over the time-series while error bars denote 1 std. The black square and horizontal bar indicate one reported observation spanning March.

Fig. 16. Comparison of monthly average Southern Ocean sea ice primary production in CNST-forcing (blue and red) and observations (gray; Arrigo, Reference Arrigo2003). Gray shaded regions denote the range of observed values. Colored shaded regions are the range in monthly averaged simulated sea ice primary production for grid cells matching observed locations (blue) and for all grid cells (red). Symbols indicate the mean value over the time-series while error bars denote 1 std.

The underestimations in Arctic PP are now clearly evident and are reflected in a reduced and shifted production peak. Observations indicate a June peak, while modeled peak production is in April. This is consistent with severe nutrient limitation in the modeled ocean. The Southern Ocean results, on the other hand, are fairly consistent with observations particularly in the late spring and summer. However, the winter comparison points to possible underestimations in primary production in some regions of the Southern Ocean.

In addition to sea ice column integrated production and chl-a biomass, our model was designed with what we hoped were the important physical processes needed to resolve the vertical ice algal communities: the upper/freeboard communities abundant in the Southern Ocean, interior ice communities, and the bottom or skeletal layer algae common in both the Arctic and Southern Ocean. Figures 17 and 18 show model comparisons with observations of mean sea ice chl-a in each of three vertical regions. For this comparison, we conservatively interpolated vertically resolved observational data and our eight-grid cell vertical output to a 120-layer common vertical grid of normalized depth before averaging into three vertical sections: upper, interior and bottom. Southern Ocean observations are from the ASPeCT data set (Meiners, Reference Meiners2013). Measurements of sectioned core data in the Arctic are less common, in part because PP is found, or assumed to be, primarily in the bottom 3–5 cm. However, measurements from all three vertical sections were obtained from the Bering Sea Ecosystem Study (Gradinger, Reference Gradinger2009), N-ICE2015 (Duarte and others, Reference Duarte2017), Point Barrow (C. Deal pers. comm.) and the Canadian Archipelago (Lange and others, Reference Lange2015).

Fig. 17. Mean monthly algal chl-a concentrations averaged in the upper (top), interior (middle) and bottom (bottom) thirds of the sea ice column. Measurements and modeled output were first interpolated conservatively to a common 120-layer vertical grid. Symbols and error bars indicate the observed monthly mean concentrations and ± 1 standard for each vertical section. Modeled means and standard deviations are from grid cells corresponding to observed locations (black lines and shading) and all grid cells with sea ice $\gt 15\percnt$ for a given polar region (blue lines and shading). References for observations are given in the text.

Fig. 18. Observations and model comparison of the mean fraction of chl-a in a given vertical section of sea ice. See Figure 17 for description of symbols and shading. References for observations are given in the text.

Modeled chl-a concentration averages for vertical sections compare well with observations in both the Arctic and Southern Ocean. This is especially true in the bottom and interior ice, but less so in the upper ice. In particular, Southern Ocean upper ice concentrations are underestimated, though by less than 1 mg chl-a/m³, during the winter months. The March peak in Arctic chl-a in observations is due to the Bering sea region which is disproportionately represented in the profile data set. Our model indicates an April peak in Arctic mean chl-a. However, this later peak is also evident in mean simulated chl-a computed from the grid locations of corresponding core observations. This suggests model biases in the Bering Sea, possibly in the timing of sea ice retreat.

Figure 18 is more revealing. Here the fraction of chl-a in a given vertical section is computed before averaging. Upper ice communities in the Southern Ocean comprise, on average, between 12 and 35% of the ice core during an annual cycle. Our model means are between 10 and 19%, underestimating upper ice communities in the Southern Ocean. Modeled chl-a, rather, favors the bottom ice layer, attributing 50–70% of total chl-a biomass. Observations indicate a lower range (40–60%) for 11 months of the year. Observations of Arctic chl-a profile fractions exhibit greater variability in seasonality than in the simulated climatology. This is in part expected given the paucity of data. Simulated mean profiles do not well represent observed profiles with very high bottom ice concentrations (e.g. February and April). The model tends to spread out this biomass among the layers. This is not surprising given that the model uses only eight grid levels and numerical diffusion is high. Simulated monthly mean bottom ice fractions range between 60 and 70% between February and March while observations suggest values > 90% in March to ~55% in June. Note that we have likely underestimated the contribution of observed ice cores with very high bottom chl-a, given that we did not use data if upper and interior ice measurements were not made or not reported.

Given the observed relationships between modeled ocean surface nitrate (via the maximal growth rate bound) and PP in the Arctic, and between snow thickness and PP in the Southern Ocean, it is perhaps not surprising that we also see correlations in the CNST-forcing time-series anomalies (annual value minus its time-series mean) of these observables. Figures 19a and b show contours of the Pearson correlation coefficient between anomalies of ice algal annual net PP and annual mean surface nitrate in the Arctic and Southern Ocean, respectively. Only significant correlations at the 95% confidence level are shown. Strong positive correlations are evident in much of the western and central Arctic. Not surprisingly, Southern Ocean PP anomalies are poorly correlated with surface nitrate anomalies except just beyond the mean 15% concentration ice edge. Here, positive anomalies in algal PP result in additional draw-down of nitrate and hence significant negative correlations with ocean surface concentrations.

Fig. 19. Contours of the Pearson correlation coefficient (r) between the anomalies of (a,b) ice algal annual net primary production (PP) and annual mean ocean surface nitrate (NO3) and the anomalies of (c,d) ice algal annual net primary production and annual mean snow thickness (hs). Only significant values at the 95% confidence level are shown. The red line is the mean sea ice 15% contour.

In both the Arctic and Southern Ocean, snow thickness anomalies are negatively correlated with anomalies in annual net algal primary production (Figs 19c and d). This is particularly true in regions of modeled higher production such as the Eastern Siberian Sea and Sea of Okhotsk and in the near coastal band of the Southern Ocean.