Relative Impact of Sea Ice and Temperature Changes on Arctic Marine Production

We use a modern Earth system model to approximate the relative importance of ice versus temperature on Arctic marine ecosystem dynamics. We show that while the model adequately simulates ice volume, water temperature, air‐sea CO2 flux, and annual primary production in the Arctic, itunderestimates upper water column nitrate across the region. This nitrate bias is likely responsible for the apparent underestimation of ice algae production. Despite this shortcoming, the model appears to be a useful tool for exploring the impacts of environmental change on phytoplankton production and carbon dynamics over the Arctic Ocean. Our experiments indicate that under a warmer climate scenario, the percentage of ocean warming that could be apportioned to a reduction in ice area ranged from 11% to 100%, while decreasing ice area could account for 22–100% of the increase in annual ocean primary production. The change to CO2 air‐sea flux in response to ice and temperature changes averaged an Arctic‐wide 5.5 Tg C yr−1 (3.5%) increase, into the ocean. This increased carbon sink may be short‐lived, as ice cover continues to decrease and the ocean warms. The change in carbon fixation from phytoplankton in response to increased temperatures and reduced ice was generally more than a magnitude larger than the changes to CO2 flux, highlighting the importance of fully considering changes to the marine ecosystem when assessing Arctic carbon cycle dynamics. Our work demonstrates the importance of ice dynamics in controlling ocean warming and production and thus the need for well‐behaved ice and BGC models within Earth system models if we hope to accurately predict Arctic changes.


Introduction
Through their impacts on water temperature and sea ice cover, changes to the global climate system are likely to have large implications for Arctic marine ecosystems and consequently for Arctic marine production and the carbon cycle. Temperature changes can impact phytoplankton growth rates, while sea ice, in addition to providing a habitat for algae, suppresses water column mixing and the air/sea exchange of gas or particles and attenuates solar radiation. Recent estimates of global surface air temperature (SAT) increases are around 0.1°C per decade for the 1998-2012 period (Huang, Zhang et al., 2017), while Arctic SAT has increased at a rate over six times greater (0.755°C per decade) during the same period. Coincident with higher temperatures, sea ice in the Arctic Ocean has undergone an unprecedented reduction in area and thickness (Cavalieri & Parkinson, 2012;Serreze & Meier, 2019;Stroeve et al., 2014;Stroeve & Notz, 2018). Correspondingly, open water in the Arctic has been increasing at a rate of ©2020. The Authors. This is an open access article under the terms of the Creative Commons Attribution License, which permits use, distribution and reproduction in any medium, provided the original work is properly cited.

10.1029/2019JG005343
Key Points: • Varying by region, reduced sea ice could account for 22-100% of the increase in annual ocean primary production under a warmer climate scenario • Changes to carbon fixation by primary producers may be an order of magnitude greater than changes to air-sea CO 2 flux under future warming • The HiLAT model replicates observed patterns of ice, surface temperature, and phytoplankton production in the Arctic but has a low nitrate bias Supporting Information: • Supporting Information S1 • Table S1 • Table S2 • Table S3 • Table S4 • Figure S1 • Figure S2 • Figure S3 • Figure S4 • Figure S5 • Figure S6 • Figure S7 Correspondence to: G. Gibson, gagibson@alaska.edu 0.07 × 10 6 km 2 yr −1 since the late 1990s, with the greatest increases in the Barents, Kara, and Siberian sectors, particularly over the continental shelf ; models suggest that this increase will continue in the future (Barnhart et al., 2015). While often presented as a linear trend, there is evidence that the increase in open water is actually occurring in a stepwise fashion (Goldstein et al., 2018).
While some of the Arctic Seas have seasonal patches among the most productive (>200 g C m −2 yr −1 ) in the world Rysgaard et al., 2001;Sakshaug, 2004;Springer et al., 1996;Walsh, 1989) at 1-17 g C m −2 yr −1 , annual primary production in the central Arctic Ocean is generally quite low (Fernández-Méndez et al., 2015;Gosselin et al., 1997). Reductions in sea ice area will expose an increasing fraction of the sea surface to solar radiation, thus increasing the available area and length of growing season for phytoplankton . The increase in open water can also result in an increase in wind induced mixing, which has the further potential to increase the supply of nutrients to the surface and to increase productivity of the region (ACIA, 2005;Anderson & Kaltin, 2001;Bates et al., 2006). A reduction in snow and ice thickness will allow increased light penetration and thus a more suitable habitat for the growth of ice algae. However, a substantial reduction in ice cover will inevitability reduce the viable substrate and thus the contribution of ice algae to total production (Gosselin et al., 1997;Watanabe et al., 2019).
Since the late 1990s, net primary production (carbon fixation by photosynthesis) in the Arctic has increased by 30% (Arrigo & van Dijken, 2015) and shifts in species composition (e.g., Grebmeier, 2012;Weydmann et al., 2014) in response to changes in timing of sea ice formation and retreat have already been observed. Arctic air temperature is projected to increase by 0.3-4.8°C by the end of the 21st century, depending on which Representative Concentration Pathway is considered (IPCC, 2013), which may lead to completely ice-free summers in the Arctic by the mid-21st century. As Arctic sea ice continues to contract, it is thought likely that primary production will increase, leading to a significant increase in the amount of carbon sequestered by the Arctic Ocean (ACIA, 2005;Bates et al., 2006;Brown & Arrigo, 2013). However, the magnitude and direction of these changes will be seasonally and regionally specific, as nutrient supply and wind speed will also play a role (Manizza et al., 2019, Tremblay et al., 2015. Vancoppenolle et al., 2013).
Understanding and predicting change to the Arctic marine ecosystem requires an integrated biogeochemical approach, connecting the relatively small Arctic Ocean to adjacent basins and adequately resolving vertical nutrient supply processes at regional and local scales. Earth system models promise to be a useful tool in this regard, but one of the largest sources of model uncertainty is located in the sea ice zone, where a subtle balance between light and nutrient limitations determines the change in primary production (Vancoppenolle et al., 2013). Here we take a regional approach to understanding the mechanistic relationship between ocean temperatures, sea ice, and biogeochemical processes controlling primary production in the Arctic. By analyzing a series of alternate model scenarios with a modern Earth system model (E3SMv0-HiLAT), we ascertain the proportion of changes to the Arctic environment and annual marine primary production that can be independently attributed to changes in sea ice and air temperature. While the relationship between the complex physical environment and ocean productivity is inherently nonlinear, an improved understanding of the relative roles of ice and temperature will help with our understanding of potential changes to productivity in the Arctic.

Method
Version 4 (CICE4) with the newly released CICE5 model, which includes new physics (i.e., mushy layer dynamics) critical for simulating biogeochemical processes in sea ice. Additionally, explicit and vertically resolved sea ice biogeochemistry was incorporated. A more extensive description of the E3SMv0-HiLAT has been previously documented (Hecht et al., 2019) and is not repeated here.
The data atmosphere to drive the ocean-ice-ecosystem model is Version 2 of the Coordinated Ocean-ice Reference Experiments (CORE2, Griffies et al., 2009). This global reanalysis data set is an interannual hindcast spanning the 62-yr period from 1948 to 2009 (Large & Yeager, 2009). The complete atmospheric data set comprises air temperature, air density, U and V wind speed, and specific humidity at 10 m, as well as surface upwelling shortwave flux, surface downwelling shortwave flux, surface downwelling longwave flux, precipitation, and sea level pressure. This atmospheric forcing data set has a 6-hr resolution in most of its components. The data-runoff model used the Dai and Trenberth Global River Flow and Continental Discharge Dataset (Dai et al., 2009;Dai & Trenberth, 2002). This river runoff data set has a monthly mean resolution from 1948 through 2007 and climatology for subsequent years. The runoff climatology is based on monthly means for a span of 5 yr from October 1999 through September 2004. Both the ocean and sea ice models were run on a global grid with a Greenland pole and a nominal 1°resolution (gx1v6).
The Biogeochemical Elemental Cycling (BEC) model implemented within the HiLAT model has been described in detail previously (Moore et al., 2002(Moore et al., , 2004Wang et al., 2015). In brief, this model is fully coupled with the ocean component and includes four key phytoplankton functional groups (diatoms, diazotrophs, phaeocystis, and small phytoplankton), a single zooplankton group, and two detrital pools. The nonsinking detrital pool largely represents dissolved organic matter (DOM), but it can also be considered to include colloidal and small nonsinking particulate organic matter. The large particulate detrital pool represents sinking particulate organic matter. BEC nutrient variables include carbon, nitrogen, phosphorous, iron, and silica. The single zooplankton compartment "grazes" on each of the phytoplankton groups, as well as the large particulate detritus, and has been parameterized to encompass the actions of both the microzooplankton and larger zooplankton.
The sea ice biogeochemistry (Deal et al., 2011;Jin et al., 2006), simulated within the CICE model, is based on the 2-D PhEcoM model. In this version, biogeochemical processes occur throughout the eight-layer ice interior, while vertical transport between layers and at the ice-water interface is modeled using the scheme of Jeffery et al. (2011). This submodel simulates ice algae, nitrate, ammonium, silicate, dissolved iron, and DOM, all of which exchange with the ocean and BEC models at a 30-min coupling frequency.
The net change in phytoplankton carbon biomass due to photosynthesis (g C m −3 day −1 ) for each primary producer in the model is given by where PCphoto is the C-specific rate of photosynthesis (day −1 ), lamda (g C (g N) −l ) represents the cost of biosynthesis, Vnc (g N (g C) −1 day −l ) is the total nitrogen uptake by phytoplankton, Rref (day −1 ) is the respiration rate, and phytoC is the phytoplankton carbon biomass (g C m −2 ). A full explanation of each of model variables and equations is provided in Moore et al. (2002) and is not repeated here.
Oceanographic initial conditions (temperature and salinity) were from Version 2 of the Polar Science Center's hydrographic climatology, a gridded climatology based on Steele et al.'s (2001) original version merging World Ocean Atlas data with the regional Arctic Ocean Atlas. As described in Wang et al. (2015), initial distributions of nutrients, inorganic carbon, and alkalinity in the sea were based on the World Ocean Atlas and Global Data Analysis Project databases. The model was initialized with no ice or ice associated variables-rather, these were allowed to spin-up at run time. The baseline HiLAT model was run for two full cycles of the 62-yr CORE2 forcing data set to ensure adequate model spin-up. Model perturbation experiments were initialized on 1 January 1960 (Model Year 137), following the twelfth year of the baseline model's third CORE2 cycle, and then run through to the end of 2009, the last year of forcing available, for a total of 50 yrs. Perturbation experiments were designed to explore the impact of sea ice on the ice and ocean ecosystem in the Arctic.
Perturbation experiments. To force a sea ice reduction, such that the resultant impact on the ice and ocean ecosystems in the Arctic could be determined, two approaches were taken: (1) increase the baseline CORE2 air temperature and (2) artificially reduce sea ice by modifying sensitive parameters. For the first experiment, the baseline air temperature in the CORE2 data set was modified using projections by the Intergovernmental Panel on Climate Change (IPCC) Working Group I atlas of global and regional climate projections (IPCC, 2013). The atlas provides a predicted change in global mean SAT relative to the early 21st century (1986)(1987)(1988)(1989)(1990)(1991)(1992)(1993)(1994)(1995)(1996)(1997)(1998)(1999)(2000)(2001)(2002)(2003)(2004)(2005) reference period. For its fifth assessment report (AR5), the IPCC adopted alternative greenhouse gas concentration trajectories, or Representative Concentration Pathways (RCP  Figure 1) illustrate that the temperature anomalies in the Arctic can be much higher, approaching 8°C in the fall and winter.
To discern the relative impact of sea ice changes on the ice-ocean ecosystem, without the compounding effects of atmospheric temperature changes, a second set of experiments was conducted using the baseline CORE2 forcing data set but modifying three model parameters that have been found to influence the sea ice area and volume predicted by the CICE model (Urrego-Blanco et al., 2016). Because of the complex, nonlinear nature of the interactions between the ice-ocean ecosystem, it is difficult to predict the exact response of the ice cover to parameter changes and thus to reduce the ice by the exact amount as seen in the RCP4.5 atmospheric warming experiment. By performing both a moderate and a more extreme ice reduction experiment, we were able to develop an appreciation of the relative role that ice plays in controlling primary production in the Arctic and the strength of the feedback between the temperature and the ice. The parameters modified were thermal conductivity of snow (ksno), maximum melting snow grain size (rsnw_mlt), and the sigma coefficient for snow grain (R_snw). To maximize impact of each parameter on ice volume, we adjusted baseline parameters to values at the far end of the known range of variation (Urrego-Blanco et al., 2016). Reducing ksno increases the insulation effect of snow, resulting in smaller heat loss and thus reduced ice growth. For this parameter we used a value of 0.03 W m −1 k −1 -an order of magnitude less than the baseline value of 0.3 W m −1 k −1 (range = 0.03-0.65). Maximum melting snow grain size (rsnw_mlt) can have a large influence on sea ice volume and snow albedo. Increasing rsnw_mlt results in a decreased albedo and increased absorption of solar radiation. We used rsnw_mlt = 3,000 μm in place of the default value of 800 μm (range = 250-3,000 μm). Decreasing R_snw has the effect of increasing grain size in the delta-Eddington radiation scheme, decreasing albedo and increasing absorption of solar radiation, resulting in increased sea ice melt rate and decreased sea ice volume. In our experiments we use a value of −2 for this parameter, in place of the default 0.48 (range = −2 to 2). The ksno is known to be a highly sensitive parameter, causing large variation in sea ice volume relative to the other parameters, particularly in fall and winter. We completed one model run in which only rsnw_mlt and R_snw were modified (Scenario S1), and a separate model run in which only ksno was modified (Scenario S2). This enabled the comparison of the ecosystem response to a relatively gentle and a relatively dramatic ice reduction.
Analysis. The environmental model output variables considered included spring sea ice area and thickness, sea surface temperature (SST), summer mixed layer depth (MLD), summer mixed layer nitrate, and average photosynthetically active radiation (PAR) in the upper layer of the ocean (midpoint at 5 m below the surface) during the growing season, all of which had a monthly temporal resolution. Snow depth over our analysis regions could also potentially vary, as snow cover is intrinsically dependent on ice area (as a substrate), upward heat flux through the ice, and on the atmospheric temperature that can melt it. Snow cover can also impact ice thickness by restricting heating and cooling, and impact the penetration of PAR to the surface ocean. Thus, to present a fuller picture of the change in the environment, we also include regionally averaged spring snow depth in our analysis. As a measure of water column stratification, we use potential energy anomaly (ϕ; Simpson et al., 1977;Simpson & Bowers, 1981) to compute the Simpson stability index, determined from the density profile, according to where h is the depth of the water column, ϕ is the amount of energy per unit volume required to completely mix the water column, ρ is the density, g is gravity, and dz is layer depth. The ϕ is 0 in fully mixed conditions and increases with stratification. For each model run, we computed the average summer stability over each region. Each of the environmental analysis variables selected was thought to play an important role in determining annual production.
Annual phytoplankton production (Prod o ) was computed by summing daily average, depth integrated productivity (prod o ) for each calendar year period, essentially an approximation to the annual integration. The mean and standard deviation of annual production were then computed for the final 10 yrs. of the model simulation. Correspondingly, annual production by the ice algae (Prod I ) was computed from monthly average values.
For each region (Figure 2), variables were averaged over all associated grid points to produce time series of regional monthly averages. For this analysis, "spring" was considered to be March, April, and May; "summer" was June, July, and August; and the "growing season" was March through September. The change in regionally averaged environmental variables and annual primary production between the baseline model run and each model scenario was determined by computing differences in model projections for the last 10 yrs. (2000)(2001)(2002)(2003)(2004)(2005)(2006)(2007)(2008)(2009) of the 50-yr model run. This analysis period ensured that the model had the maximum time possible to adjust to our perturbations in forcing (40 yrs.) while allowing sufficient years to explore interannual variability. For the same 10-yr period, the proportion of daily phytoplankton productivity and annual primary production contributed by each phytoplankton group, the day of peak phytoplankton biomass, and the bloom start date (BSD) were also computed. Adapting the cumulative biomass approach outlined in Brody et al. (2013) for each year examined, BSD was assumed to be the day at which regionally averaged cumulative phytoplankton biomass exceeded 15% of the total annual biomass. Mean and standard deviation were then calculated for the 10-yr period.

Results
Model performance. The HiLAT models ability to simulate ice volume, ocean temperate, ocean nitrate concentration, annual phytoplankton and ice algae primary production, chl-a biomass, and CO 2 is presented in the supporting information section. Model performance with respect to ice, temperature, and nitrate, environmental variables that can impact ocean primary production, is summarized here in brief. The ice volume simulated by the baseline HiLAT model varied annually (17.3 ± 6.7 · 10 3 km 3 ) and compares favorably to sea ice volume calculated using the Pan-Arctic Ice Ocean Modeling and Assimilation System ( Figure S1; Schweiger et al., 2011;Zhang & Rothrock, 2003) that incorporates reanalysis data. Model estimates of both summer and winter upper water column temperature generally had a cool bias relative to observations ( Figure S2). Excluding the Barents and the Bering Sea, the average summer model bias was −0.7°C. The Eurasian Basin had a small positive (warm) simulated bias (0.4°C), while the Barents Sea had a relatively large (−3.7°C) cool bias. The winter temperature bias of model estimates over most analysis regions was generally less than −0.5°C. The Barents Sea also had a low bias (−3.2°C) relative to observations in this season, as did the Bering Sea (−1.7°C). Despite the differences, model estimates and observations overlapped within one standard deviation in both summer and winter.
The broad spatial pattern of upper water column nitrate throughout the Arctic prior to the onset of spring phytoplankton growth ( Figure S3a) matches reasonably well with the patterns seen in consolidated observations (Codispoti et al., 2013). The summer regional average NO 3 concentration for all analysis regions as simulated by the baseline HiLAT model ( Figure S3b) has a low bias, averaging~2 mmol N m −3 across regions. However, with the exception of the Eurasian basin, where the bias exceeded 3 mmol N m −3 the modeled and observed nitrate concentrations did overlap within one standard deviation. With a bias exceeding 5 mmol N m −3 , the HiLAT simulated average NO 3 concentration for the prebloom period in the Beaufort, Barents, and Eurasian Basins was notably less than observation-based estimates. The number of observational data points for these regions is relatively small (Table S2) and the spatial coverage was quite heterogeneous. The bias in the Nordic seas, where observational data were more plentiful and the coverage more even, was closer to 2 mmol N m −3 and modeled and observed nitrate concentrations overlapped within one standard deviation. In the Bering Sea and the Canadian Basin model estimates were even closer to observations, with a bias of only 1.0 and 0.03 mmol N m −3 respectively. The ocean-atmosphere CO 2 flux and the annual phytoplankton production in most regions were in good agreement with observations, but the model appears to underestimate ice algae production ( Figure S4-S6).
Sea ice. The extent of sea ice over the Arctic Ocean for each of the model scenarios is visualized in Figure 3. Summer sea ice area in the baseline model averaged 5,720,400 km 2 (Figure 3a, iv). Increasing rsnw_mlt and decreasing R_snw (Scenario S1) reduced summer sea ice area by 22% to 4,440,900 km 2 (Figure 3a, i). Reducing ksno (Scenario S2) caused a larger reduction (39%) in ice area, down to 2,222,300 km 2 (Figure 3a, ii). An intermediate (30%) reduction in summer sea ice area resulted from an increase in air temperature (RCP4.5 experiment, Figure 3a, iii). Despite the differences in ice area, the overall spatial pattern of ice cover was similar in each model experiment; the reduction in overall ice area was due to the reduction in fraction of the model grid cell covered. Simulated median ice cover fraction over the entire Arctic was 0.6, 0.4, 0.2, and 0.3 in the baseline, S1, S2, and RCP4.5 experiments, respectively. In each scenario, the greatest ice cover was seen off the northern coast of Greenland and the Canadian Arctic Archipelago.
Upper water column temperature. All three ice reduction experiments resulted in an increase in annual average water temperature in the upper 25 m, with respect to the baseline run ( Figure 3b), though the strength and the extent of the warming was notably different between model scenarios. With only a small reduction in ice cover and no modification to atmospheric temperature forcing (S1), temperature changes across the Arctic Ocean were positive, but generally less than 0.25°C. A more dramatic reduction in ice cover (S2 model scenario, Figure 3b, ii) resulted in a stronger (~0.5-1°C) and more extensive warming of the Arctic Ocean that extended southward into the northern Bering and western Nordic Seas region. The increased atmospheric temperatures in the RCP4.5 scenario caused extensive ocean heating throughout the Arctic, North Pacific, and North Atlantic (Figure 3b, iii); however, the warming of the central Arctic Basins was actually more extensive (by~0.25°C) in the S2 experiment (Figure 3b, iv).
Primary production. Changes to phytoplankton production in response to ice and temperature forcing were reached within the first few years of the 50-yr model simulation ( Figure S1). With few exceptions, annual phytoplankton production, relative to the baseline model run, increased across the Arctic in each of the three ice reduction scenarios ( Figure 3c and Table 1). With a moderate reduction in the ice cover (S1), there was only a small increase in annual production (~1-10 g C m −2 ) over the Arctic shelves and basins. The more dramatic reduction in ice, in the S2 model scenario, resulted in production anomalies reaching~5-15 g C m −2 across the Arctic, with larger anomalies in the Gulf of Anadyr (~58 g C m −2 ) and the western Greenland Sea, and the northern edge of the Barents Sea (~30 g C m −2 ). The regions of high positive change across the Arctic generally aligned with the retreating ice edge. Coincident with the increased production across much of the Arctic, a small reduction in production (~2.5-5 g C m −2 ) occurred in the Labrador Sea, a region that did not support notable summer sea ice. In the atmospheric warming scenario (RCP4.5), these negative production anomalies were more extensive, occurring south of Greenland and the Labrador Sea, and stronger, indicating a larger reduction in production. Despite these regions of negative production anomalies, positive production anomalies extended into the southern Bering Sea and the Gulf of Alaska under the warming scenario. As discussed in the regional analysis section, the relative importance of ice and ocean temperature for primary production varied by region. Only the S2 ice reduction scenario resulted in regional differences in annual phytoplankton production that could be deemed significant (i.e., simulated increase was greater than two standard deviations). The significant increases were in the East Siberian Sea, the Laptev Sea, and the Canadian and Eurasian Basins.
Neither diazotrophs nor phaeocystis was significant contributors (<1%) to total annual phytoplankton production in any of the regions examined and so are not discussed further here. The relative contribution from diatoms and small phytoplankton to daily productivity varied by region and time of year ( Figure 4 and Table 1). In both the Bering Sea and the Nordic Seas, phytoplankton productivity increased in the spring (late April), with bloom days of 115 and 119 for the two regions, respectively, and was dominated by diatoms. At bloom peak (Day 133 and Day 127 in the Bering and the Nordic, respectively), diatoms contributed 99% of the daily production. Small phytoplankton took over as the dominant producer in mid-July (Day 198) in the Bering Sea and in early September (Day 255) in the Nordic Seas, though diatoms contributed more than 65% Figure 3. Arctic sea ice cover, ocean temperature, annual primary production, and annual CO 2 flux in each model scenario. (a) Summer fraction of model grid cell covered in sea ice in each sea ice reduction experiment S1 (i), S2 (ii), and RCP4.5 (iii), and the baseline model run (iv). (b) Annual average temperature anomaly in the upper 25 m relative to the baseline run (Temp EXP -Temp BASE ) for S1 (i), S2 (ii), and RCP4.5 (iii). Red colors indicate warmer waters than the baseline run, while blue colors indicate cooler waters than the baseline run. Panel (b, iv) shows the temperature difference between the S2 and the RCP4.5 experiment (Temp RCP4.5 − Temp S2 ). (c) Difference (Prod EXP − Prod BASE ) in annual average phytoplankton production relative to the baseline model run for experiments S1 (i), S2 (ii), and RCP4.5 (iii). Red colors indicate higher production than the baseline run, while blue colors indicate lower production than the baseline run. of total annual phytoplankton production. In the Barents Sea, although the day of peak productivity slightly lagged (17 days) behind the Nordic Seas region, phytoplankton productivity exhibited similar trends with respect to timing and diatom contribution. In the other, more interior Arctic seas and basins, the bloom did not start until the end of May/beginning on June. Except for the East Siberian Sea, diatoms still dominated during the initial bloom phase, but by the productivity peak in mid-July to early August (Day 193-218), small phytoplankton contributions had increased, reaching a maximum in mid to late September (Day 253-268). In the Eurasian Basin, diatoms contributed 63% to total annual phytoplankton production, though in other interior seas their contribution was less than 52%. Throughout the seasonal cycle in the East Siberian Sea, diatom contribution to the total phytoplankton productivity was never more than 53%. Contribution to total annual phytoplankton production in the region was only 19%, with the remainder coming from small phytoplankton.
In general, the changes to ice and temperature in our three model scenarios did not have a notable impact on the relative contribution of the phytoplankton groups to annual production. However, the diatom contribution was significantly less in the Arctic Basins under the more dramatic ice reduction scenario (S2) compared to the baseline run (Table 1). Although the three ice reduction scenarios generally resulted in an earlier regionally averaged bloom day, none of the timing shifts was significant. The day of peak productivity was Note. Regional means and one standard deviation are for the 1999-2009 period. The maximum percentage of productivity (prod o ) contributed by the diatoms (D) and small phytoplankton (PS), the day on which the maximum contribution occurs and the annual contribution of diatoms to annual production, and the bloom day and day of peak productivity in each model run are also shown. The shading indicates when the values from the model scenarios do not overlap with values from the baseline simulation.
also generally shifted earlier under each scenario. The most notable changes were in the Laptev Sea, the Kara Sea, and the Canadian Basin, where simulated average peak productivity day was 20-30 days earlier than in the baseline model run, although all predictions overlapped within two standard deviations.
Ocean-atmosphere carbon flux. Changes to ocean-atmosphere CO 2 flux resulting from each model experiment are shown in Figure 3d and Table 2. A small reduction in ice (S1) resulted in only a small increase (<0.5 mmol CO 2 m −2 day −1 ) in the flux of CO 2 into the ocean, primarily over the Barents Sea in the vicinity of the receding ice edge. A more dramatic reduction in ice (S2 and RCP4.5) resulted in an increase in the CO 2 flux into the ocean (~1-3 mmol CO 2 m −2 day −1 ) in the western Nordic Seas between Iceland and Svalbard, over much of the Barents Sea, and the northern Bering Sea-regions that were ice covered for at least part of the year in the baseline experiment. This increase in CO 2 flux into the ocean was more pronounced in the atmospheric warming scenario (RCP4.5) and wider spread in the central Arctic Basin. In the eastern Nordic Seas and the Southern Bering Sea, regions that were generally not ice covered in the baseline model run, the CO 2 flux into the ocean was reduced slightly in the S2 scenario and more noticeably in the RCP4.5 scenario (−5.1 and −13.9 Tg C yr −1 for the two regions, respectively). Despite this change, the overall direction of CO 2 flux in the Nordic Seas region remained positive into the ocean. However, under the RCP4.5 atmospheric warming scenario, the Bering Sea region exhibited a notable reversal in the direction of CO 2 flux from a 7.6 Tg C yr −1 influx to a 6.5 Tg C yr −1 efflux, indicating a regional shift from a sink to a source of CO 2 . More broadly, the atmospheric warming scenario (RCP4.5) resulted in an Arctic Ocean-wide net increase in CO 2 flux into the ocean of 5.5 Tg C yr −1 (or +3%). This overall increase was essentially equal

Journal of Geophysical Research: Biogeosciences
to the net increase over the Arctic Basins (5.7 Tg C yr −1 ), though this was a much larger percentage change (+35%) for the basins.
Regional analysis. The following regional analysis relates the changes in annual net primary production in each model simulation to environmental changes. All three ice reduction experiments resulted in an increase in annual phytoplankton production relative to the baseline run, while the response of annual ice algal production was more variable (Figure 5 and Tables 1 and S4). Over all three scenarios, the correlations between regional changes in Prod o and environmental variables were highest with annual SST (r = 0.7) and spring ice area (r = −0.6). The correlation to ice depth, snow depth, and PAR was approximately −0.5, while the correlation to MLD, water column stability and summer mixed was only −0.05, −0.24, and −0.12, respectively.

Inflow shelves: The Chukchi and Barents Seas.
In the Chukchi Sea, the S1 scenario resulted in only a small reduction in ice area (4,624 km 2 ) and ice thickness (3 cm). Associated with these changes, the regionally averaged SST increased by 0.2°C, snow depth decreased by 1 cm, and PAR increased by 3 W m −2 . The increase in phytoplankton production in response to these changes was only 1.9 g C m −2 . At~7 g C m −2 (or~19%), the increase in phytoplankton production was not notably different between the atmospheric warming experiment (RCP4.5) and the more dramatic ice reduction experiment (S2). Reducing spring ice cover over the region by 29,234 km 2 and the ice thickness by 102 cm (S2) led to an average surface warming of 0.6°C, a reduction in snow depth by 3 cm, and an increase in PAR of 10 W m −2 . Atmospheric warming (RCP4.5) reduced ice area by only 10,805 km 2 , but resulted in a slightly larger (0.9°C) increase in ocean temperature averaged across the region, a reduction in ice thickness of only 17 cm, a reduction in snow depth of 15 cm, and an increase in PAR of 7 W m −2 across the region. In the Barents Sea, the largest ice area reduction (332,980 km 2 ) and production increase (14 g C m −2 /23%) were associated with the atmospheric warming scenario, resulting in an ocean temperature increase of 1.4°C. Ice depth was reduced by only 5 cm, while snow depth was reduced by 15 cm and the corresponding PAR increase was 7 W m −2 (Table S4). Comparison with the S2 scenario suggests that a similar PAR increase (6 W m −2 ) and production increase (12 g C m −2 /19%) in the Barents resulted from a much smaller ice loss (152,720 km 2 ) and larger reduction in ice thickness (47 cm). These changes drove a 0.5°C water temperature increase and an 11-cm reduction in snow depth. Ice algae production in the Chukchi Sea increased by 34% (0.3 g C m −2 ) in the atmospheric warming scenario but decreased by twice that amount (69%/0.6 g C m −2 ), with further ice reduction and loss of habitat (S2). In the Barents Sea, a small reduction in ice thickness and area (S1) led to a notable increase in ice algae production (21%), while further ice area reduction was associated with a decrease in ice algae production. In this region, a larger reduction in ice algae production was seen in S2 relative to RCP4.5, despite smaller ice loss in the former scenario-likely because ice thickness was insufficient to support ice algae Note. Positive numbers indicate flux of CO 2 from the atmosphere into the ocean. The change in the flux, relative to the Baseline model, that was seen in the S2 and RCP4.5 scenarios is also presented; in this instance, a positive number indicates an increase in the influx relative to the baseline scenario, whereas a negative number indicates a reduced influx to the ocean. "Arctic" refers to the combined ocean region above 66.56°, and "Basins" refers to the Canadian and Eurasian Basins combined. For comparison, the corresponding regional change in net annual primary production (Prod o ), computed from the annual average primary production in Table 1, is also shown.
growth over parts of the region. In both inflow regions, water column structure does not appear to be a primary driver of the increase in production, as MLD and stability changes were not consistent with simulated production changes, and nitrate concentration in the mixed layer was slightly reduced in all experiments.
Interior shelves: Kara, Laptev, ESS, Beaufort. Each of our model scenarios (S1, S2, and RCP4.5) resulted in an increase in phytoplankton production, ranging from 0.8 to 11 g C m −2 (3-47%) across the four Arctic interior shelves (Tables 1 and S4). The increased production appears to be primarily due to the increase in light availability that comes with a reduction in ice thickness and snow depth. The reduction in ice area relative to the baseline experiment varied by region and scenario but only ranged from 1% to 6% over each of the interior Figure 5. Regional average magnitude and direction of change in annual phytoplankton (Prod o ) and ice algae production (Prod I ), and accompanying environmental variables, in response to three ice reduction Scenarios S1 (white), S2 (gray), and RCP4.5 (black). Environmental variables examined include spring ice area (Ice Area), spring ice thickness (Ice h ), snow depth (Snow h ), annual average sea surface temperature (Temp), summer mixed layer depth (MLD), summer water column stability (Stability), photosynthetically available radiation during the growing season (PAR), and summer mixed layer nitrate concentration (NO 3 ). The regional average differences were calculated for the last 10 yr of the model run (2000)(2001)(2002)(2003)(2004)(2005)(2006)(2007)(2008)(2009), over the regions as defined in Figure 2. The computed differences for all regions are presented in supporting documents Table S4. The correlation coefficients (r Prodo ) between changes in Prod o and changes in each of the environmental variables are also shown.

Journal of Geophysical Research: Biogeosciences
shelves. In each region, the largest reductions in ice thickness (74-121 cm) were seen in the S2 scenario, while the largest reduction in snow depth (7-14 cm) was in the RCP4.5 atmospheric warming scenario. At 0.6°C and 0.8°C, the increase in surface temperature in S2 and RCP4.5 across the interior shelves was quite similar. Reduced ice presence was more influential than atmospheric warming in raising the water column temperature in the Laptev ( Figure 5) and ESS, while in the Kara and Beaufort Seas, water column temperature increases were more responsive to atmospheric warming than a reduction in ice cover. On each interior shelf, larger PAR increases were seen in the S2 scenario (10-13 W m −2 ) than in the warming scenario RCP4.5 (7-8 W m −2 ). Summer mixed layer nitrate concentration decreased in each region and scenario, while the MLD and stability response was mixed-suggesting none of these environmental variables were primary factors in controlling the production increase.
In contrast to phytoplankton production, annual ice algal production generally decreased across the interior shelves with each ice reduction scenario. Depending on the scenario, the decrease ranged from 15% to 89% across the ESS, Laptev, and Beaufort. In the Kara Sea, while an 85% reduction in ice algae production resulted from the more dramatic ice reduction scenario (S2), there was no notable reduction in production with a more gentle ice reduction (S1), and there was an increase (19%) in production in the atmospheric warming scenario (RCP4.5).
Central Arctic Basins. The phytoplankton production response to environmental changes in the two Arctic Basins (Canadian/Eurasian) was quite similar ( Figure 5 and Tables 1 and S4), with both regions showing a 4 g C m −2 increase with moderate ice reduction (S1), a 10/13 g C m −2 increase with a more dramatic ice reduction (S2), and a 6/7 g C m −2 increase in response to an increased atmospheric warming (RCP4.5).
The larger production response in the S2 scenario relative to the RCP4.5 scenario appears to be in response to the more dramatic reduction in ice thickness (120/108 cm), and to the resulting 15 W m −2 increase in PAR and surface temperature e (~0.7°C) in the S2 scenario. Snow depth decrease by 16/18 cm in RCP4.5, a larger decrease than seen in S2 (9/14 cm), but the PAR increase was half as much (~7 W m −2 ). A small reduction in ice area and thickness (S1) could account for 0.2°C (50%) of the 0.4°C increase in ocean temperature seen with atmospheric warming (RCP4.5). Across scenarios, with the increase in surface ocean warming, stability of the water column increased and MLD decreased, which would have the effect of retaining phytoplankton in the well-lit upper water column. Production increases were also seen, despite reduced summer mixed layer nitrate levels, suggesting this region was not primarily nutrient limited.
Marginal, subarctic seas: Bering Sea (inflow) and Nordic Seas (outflow). The RCP4.5 atmospheric warming experiment resulted in a 9.4 g C m −2 (8%) increase in annual average Bering Sea production. This increase coincided with a 1°C warming at the sea surface and a small (~5 cm) deepening and an increase in the stability of the mixed layer. The average mixed layer summer nitrate concentration also increased slightly (0.7 mmol N m −3 ). The S1 and S2 scenarios, including the modification of ice parameters and no temperature forcing, resulted in a < 0.1°C surface temperature increase, <0.1 mmol N m −3 nitrate increase, and <2% increase in production. In this Bering Sea region, there was no notable change in regionally averaged ice thickness in S1, only a 1-cm reduction in thickness in RCP4.5, and a 22-cm reduction in S2. Snow depth decreased by 10 cm in RCP4.5 but by only 3 cm in S2 and <1 cm in S1. In all three scenarios, the increase in PAR was <2 W m −2 . Ice algae production increased (0.1 g C m −2 ) in both the S1 scenario, which had a relatively small reduction (30,884 km 2 /8%) in ice area, and in the RCP4.5 scenario (0.18 g C m −2 ), which had the largest reduction in ice area (158,470 km 2 /38%). The S2 scenario, which had an intermediate ice

Journal of Geophysical Research: Biogeosciences
area reduction (27%) gave rise to a large reduction in ice algae production (0.42 g C m −2 ). In the Nordic Seas region, the largest increase in phytoplankton production (9.4 g C m −2 ) was seen in the RCP4.5 scenario. The reduction in ice area (120,040 km 2 ) and ice depth (11 cm) was less than in the S2 scenario (146,020 km 2 and 71 cm), but the average surface temperature of the Nordic Seas increased by 0.93°C-nearly double the temperature increase (0.48°C) in the S2 ice reduction scenario.
In the RCP4.5 warming experiment, the average surface temperature of the Nordic Seas increased by 0.93°C, and the production increased by 9.4 g C m −2 . The temperature increase in the S2 ice reduction scenario was less than half this (0.48°C), but the production increase was greater (12 g C m −2 ). In concert with a relatively small increase in temperature (0.1°C) and PAR (0.94 W m −2 ), the production increase in S1 was only 2.5 g C m −1 . The MLD shoaled and stability decreased in the S2 scenario, while MLD deepened and stability increased in the RCP4.5 scenario. The summer mixed layers nutrient increase was small and equivalent in both experiments. The reduction in ice thickness (5, 71, and 11 cm) and snow Figure 6. The predicted change in production (indicated with an asterisk) that would be expected from the ice loss in the warming scenario (RCP4.5), if we assume a linear relationship (dashed line) between ice area loss and production increases simulated by S1 (white circle) and S2 (gray circle). Actual production change simulated in RCP4.5 (black circle) is also shown.

10.1029/2019JG005343
Journal of Geophysical Research: Biogeosciences depth (2, 13, and 22 cm) in the S1, S2, and RCP4.5 scenarios only caused relatively small increases in PAR (0.9, 3.8, and 2.89 W m −2 , respectively). The change in ice algae production was small (<0.1 g C m −2 ); however, in the case of the S2 scenario, this small change was equivalent to a 63% reduction in total production in the region.
Summarizing the impact of temperature and ice area on production. Previous research has shown that net primary productivity in the Arctic increases proportionally with the amount of open water, with rates of net primary production per unit area remaining virtually unchanged (Arrigo & Van Dijken, 2011). As such, here we assume a linear relationship between ice area and annual phytoplankton production changes to understand, at least to the first order, the relative impact of ice presence and atmospheric warming on production throughout the Arctic Ocean ( Figure 6 and Table 2).
Assuming a simple linear regression, changes to regionally averaged ice area and associated ocean production changes seen in the S1 and S2 scenarios were used to predict the change in production that we would expect in the RCP4.5 warming scenario due to ice loss alone-that is, By comparing ΔProd P to the production that was actually simulated (ΔProd S ) in the RCP4.5 warming scenario ( Figure 6) we can calculate what percentage of these production changes can be apportioned to a reduction in ice area (Table 3), and thus what fraction is due to the additional warming.
Note. ΔProd P and ΔSST P respresent Prod O and SST predicted from the S1 ans S2 scenarios (i.e. Eq 2).
In the Bering, Chukchi, East Siberian, Laptev, and Beaufort Seas, the production increase in the RCP4.5 scenario was greater than predicted from the simulated ice area reduction. Changes in ice area could account for only 21.5% of the production changes in the Bering Sea, but for as much as 85.5% in the Laptev Sea. In all other regions, the production predicted based on loss of ice area was greater than the production that was actually simulated in the warming scenario, suggesting ice plays a more important role in controlling production in these regions than the air temperature. Generally, the regionally averaged production increase in the atmospheric warming scenario was less than 5 g C m −2 yr −1 more than could be expected from ice reduction alone. In the Bering Sea, the additional heating of the surface ocean because of atmospheric warming resulted in an additional 7 g C m −2 yr −1 beyond the production increase that we could associate with the loss of ice area.
If we similarly assume a linear relationship between ice area and SST, we can begin to quantify the importance of increased atmospheric warming relative to diminishing ice cover, and the resultant increased heat exchange between the ocean and atmosphere, in causing surface ocean heating in each Arctic region. Figure 7 illustrates the change in SST that we would predict based on the reduction in ice area and surface heating seen in Scenarios S1 and S2, along with the surface ocean temperature increase that was actually simulated in the atmospheric warming scenario (RCP4.5). In the Canadian and Eurasian Basins, the SST increase that we can predict from ice area reduction more than accounted for the temperature increase observed in the RCP4.5 atmospheric warming scenario. In all other regions, the temperature increase that was simulated in the RCP4.5 scenario exceeded any temperature increase that could be predicted solely from a reduction in ice area. In the Laptev, Kara, Barents, Beaufort, and East Siberian Seas, ice reduction could account for an average of 43% of the simulated SST increase in the RCP4.5 scenario. The additional 0.37°C increase can thus be attributed to the atmospheric warming. In the more marginal Chukchi and Nordic Seas, the ice area reduction could account for~37% of the SST increase seen in RCP4.5, with the additional increase in SST with atmospheric warming averaging 0.56°C. A reduction in ice area in the Bering Sea could only account for 11% of the surface warming seen in the RCP4.5 scenario, thus the atmospheric warming was responsible for an additional 0.89°C increase in SST.

Summary and Discussion
Cold and covered in ice for much of the year, the Arctic Ocean is a complex, often light-limited, environment, historically considered a region of low biological production and not a

Journal of Geophysical Research: Biogeosciences
the Arctic to increase proportionally with the amount of open water and with rates of net primary production per unit area remaining virtually unchanged. Open water can result from ice melt due to increased atmospheric temperatures and also from changes in atmospheric and ocean circulation. In a positive feedback loop, increased open water will allow additional heating of the surface ocean and increased ice melt. Given the speed of change in the Arctic (Huang et al., 2017;Serreze & Meier, 2019), understanding the relative controls of atmospheric temperature and sea ice on production in this region is of great interest. By analyzing a series of alternate ice reduction scenarios with a modern Earth system model (E3SMv0-HiLAT), we have determined to the first order the percentage of annual marine primary production that can be independently attributed to changes in sea ice and air temperature. All experiments were conducted with the HiLAT model in a partially coupled G-configuration, meaning that atmospheric conditions do not evolve in response to flux changes between the atmosphere and the ocean. Without two-way coupling, some potentially important dynamics, such as changes to wind forcing and air-sea/ice gas flux interactions, could have been missed. Regardless, validation of the baseline HiLAT model found it Figure 7. The SST that would be expected (indicated with an asterisk) from the ice are loss in the warming scenario (RCP4.5), if we assume a linear relationship (dashed line) between ice area loss and changes to SST simulated by S1 (white circle) and S2 (gray circle). The actual change to SST simulated in RCP4.5 (black circle) is also shown.

10.1029/2019JG005343
Journal of Geophysical Research: Biogeosciences capable of simulating observed patterns in upper water column temperature, ice volume, and primary production in the Arctic, making it a suitable tool for our investigation. Our 10-yr analysis period spanned both positive and negative phases of the Arctic Oscillation index (AMAP, 2012). Thus, by considering mean annual production response to ice and temperature perturbations, our results have provided a fairly good summary of how a region can be broadly expected to respond to reduced ice and increased atmospheric temperature. Our results illustrate that the relative control of atmospheric temperature and ice on annual Arctic production are not homogeneous but vary by region.
In the Arctic Basins and the Nordic and Barents Seas, a reduction in ice area could account for all of the SST increase that was simulated in an atmospheric warming scenario. In these regions, the model actually simu-lated~0.37°C more surface warming when ice was reduced, than when additional atmospheric heating was applied, presumably because of the additional ocean surface area that was available for interactions (i.e., heat exchange and wind mixing) with the atmosphere. Conversely, reduction of ice area over the Pacific Arctic Region (i.e., the Chukchi, East Siberian, and Beaufort Sea) and adjacent Laptev Sea could account for less than 50% of the simulated surface ocean warming. In the Bering Sea, a region with high year-to-year variability in spring sea ice extent (Stabeno et al., 2012(Stabeno et al., , 2017 ice area reduction only accounts for~11% of the SST increase. Correspondingly, the relative control of ice and atmospheric temperature on phytoplankton production is very different for the Pacific Arctic Region, compared to the central Arctic Basins, the Nordic Sea, and the adjacent Kara and Barents Seas. In the Pacific Arctic Region, by imposing a reduction in ice area (increasing open water) without increasing atmospheric temperature, we could explain only 22-86% of the simulated increase in production observed, given an equivalent ice reduction under an atmospheric warming scenario. This suggests that the remaining increase in regionally averaged annual net production was due to an increase in productivity, caused by temperature-dependent factors, such as the growth rate or surface heat-driven reductions in ice and snow depth or porosity that allowed additional light penetration under the ice. Throughout the remainder of the Arctic, we found that a reduction in ice area, which would correspondingly increase the open water available for production, was more than sufficient to explain

10.1029/2019JG005343
Journal of Geophysical Research: Biogeosciences the regional increase in production observed under our atmospheric warning scenario. An increase in primary productivity during the spring bloom was recently observed in both the Barents and Kara Seas (Renaut et al., 2018). While this appears to contradict our findings for annual changes in production, it merely highlights the importance of time frame when considering the response of Arctic production to a changing environment.
Annual net phytoplankon production increased throughout the Arctic in each of our model scenarios, with the largest regionally averaged annual production increases (exceeding >20 Tg C) in the Nordic, Barents, Bering, and Arctic Basins under our atmospheric warming scenario. Slagstad et al. (2015) found that the greatest increase in annual primary production would occur along the Eurasian perimeter of the Arctic Ocean (up to 40 g C m −2 ) and in the Barents Sea and Kara Sea (40-80 g C m −2 ). We also found some of the largest increases in annual primary production to occur in the Barents, in the vicinity of the retreating ice edge. However, these patches of higher production were relatively narrow, and the region wide average production increase was much smaller at~6-14 g C m −2 . While not a primary focus of our study, this suggest that subregional ice and temperature dynamics will also be important. In a study focused on the controls of ice and temperature in the Barents Sea, Dong et al. (2020) found that the spring phytoplankton bloom occurred earlier and had a higher magnitude in warm versus cold years in the northern part of the Barents Sea but found no such association in the southern Barents Sea. Subregional heterogeneity has also been documented for the Bering Sea, where Brown and Arrigo (2013) found that annual net primary production during warm years of early sea ice retreat was enhanced by 40-50%, compared to years with late sea ice retreat in the southeastern Bering Sea. Across all three model scenarios, regional SST and ice area changes were most highly correlated with changes to annual phytoplankton production. On these scales, there appears to be little correlation between the annual production and the summer water column structure (i.e. MLD and stability) or summer mixed layer nitrate concentration. This is not to say that these environmental drivers are not important to overall Arctic production; rather, any controlling influence they have must be happening on subregional and subseasonal time scales, or during a time period not considered.
The HiLAT model has four categories of phytoplankton, though phaeocystis and diazotrophs were not of significant concentrations in our study area. Diatoms dominated annual primary production in the Eurasian Basin and the Bering, Nordic, and Barents Seas in our baseline simulation, but contributed~50% or less in the other Arctic regions. Morán et al. (2010) found that temperature changes alone could account for 73% of the variance in the relative contribution of small cells to total phytoplankton biomass in the eastern and western temperate North Atlantic Ocean and predicted a gradual shift toward smaller primary producers in a warmer ocean across a diverse range of environmental conditions. Under our more dramatic ice reduction scenario, the HiLAT model suggests the relative contribution from small phytoplankton and diatoms to the regionally averaged annual production in the Arctic will remain relatively stable. The exception is the prediction of a significant decline in the contribution of diatoms to annual production in the Arctic Basins. Such a shift toward smaller phytoplankton cells has actually already been observed in the Canadian Basin (Li et al., 2009). With reduced ice and increased atmospheric warming, our experiments did not reveal a notable shift in bloom day but did find a notable shift (generally earlier) in the peak production day. This shift was significant, and by as much as 3 to 4 weeks, in the Laptev, Kara, and Canadian Basins. Kahru et al. (2011) has previously found a significant trend toward earlier Chl maximum day (by up to 50 days) in about 11% of the Arctic, in areas where summer ice concentration has decreased. Our results suggest that this trend toward earlier peak production will spread across the Arctic, as spring and summer ice continues to recede. Given ongoing sea ice decline and warming ocean temperatures, it is possible that phytoplankton species will adapt accordingly, such that future Arctic assemblages comprise cells better able to take advantage of the warmer temperature and increased light availability. While the HiLAT model can simulate shifts between its four phytoplankton categories, as with most global Earth system models, it is not currently well equipped to capture such a species shift, or how it would impact annual production.
The baseline HiLAT model estimate for net annual Arctic Ocean CO 2 uptake captured magnitude and spatial patterns of the annual average air-sea flux previously determined for this region (MacGilchrist et al., 2014;Manizza et al., 2019;Schuster et al. 2013;Yasunaka et al., 2018). Cai et al. (2010) challenged the common assumption that a more ice-free Arctic Ocean would directly translate into a larger CO 2 sink simply because of the increased area for air-sea gas exchange. More recently, a reduction to the Arctic CO 2 sink, despite decreasing ice area, has been noted (Manizza et al., 2019). Unlike phytoplankton production, which exhibited a mean annual increase in each Arctic analysis region in response to ice area reduction atmospheric warming, we found the response from the CO 2 flux to be more variable. The balance between ice area and surface water temperature in controlling the overall direction of CO 2 flux in the Arctic appears to be quite delicate. Under our warming scenario, the HiLAT model simulated an additional net influx of 5.7 Tg C yr −1 , relative to the baseline run, in the central Arctic Basins. However, in the more dramatic ice reduction experiment, which had larger ice loss and greater surface warming, the change to the net CO 2 influx was much smaller (1.3 Tg C yr −1 ). This suggests that while the Arctic Basins may initially become a larger sink for CO 2 , as ice area continues to recede, this change will be a relatively short-lived. The interior shelf seas-that is, the Chukchi, Kara, Laptev, ESS, and Beaufort Seas-each exhibited only a very minimal change in annual net CO 2 flux in both the warming and ice reduction scenarios. None of these regions are currently considered significant CO 2 sources or sinks, and it appears this will likely remain true in the foreseeable future. With a predicted CO 2 outflux close to the present influx, our atmospheric warming simulation indicates that the Bering Sea will become a CO 2 source as the climate warms. This flux reversal appears to be driven directly by the temperature, rather than ice changes, as no such flux reversal was seen when we forced ice reduction without increasing atmospheric temperature. The change in carbon fixation due to increased net primary production with an atmospheric warming scenario was generally an order of magnitude greater than the change in air-sea CO 2 flux.
While the HiLAT model could replicate observed spatial patterns of nitrate in the upper water column throughout the Arctic and Subarctic, our analysis indicates that the model has a persistent low bias of nitrate over this domain. It is important to consider the sparsity of nutrient data available for validation, but the bias appears to be notably large (~5 mmol N m −3 ) over the Arctic shelf seas during the prebloom period, when surface waters are usually replenished with nutrients. Sufficient nitrate supply is a fundamental requirement for most marine primary producers, and adequately simulating nutrient fluxes into the Arctic will be important for fully understanding how production responds to ice reduction and temperature increases in a future warmer Arctic (Arrigo & Van Dijken, 2011). It is likely that the low nitrate bias is also the cause of the model's lower than observed ice algae primary production estimates. Fixing such a nitrate shortfall would potentially enable an additional~14 g C m −2 of production. It is possible that the nitrate bias caused us to underestimate the importance of the summer water column structure and nitrate concentration. Although, given the relative strength of the correlations between production and SST and ice area, it seems likely that these other environmental factors would be of secondary importance.
The HiLAT model simulates temporal patterns in ice primary production, similar to other studies (Deal et al., 2011;Gradinger, 1999) and captures observed ice algae chl-a concentrations, though production is underestimated by an order of magnitude when compared with past estimates. This underestimation is consistent with the model's low surface nitrate bias. Unlike oceanic phytoplankton, which can thrive at the subsurface nutrient maximum, sea ice algae are dependent on the upper ocean for essential nutrient fluxes, and so they exhibit greater sensitivity to ocean surface biases. At times, the primary productivity of ice algae can significantly exceed water-column productivity (Gosselin et al., 1997;Gradinger, 2009), and ice algae production could constitute an important part of the overall production in the Arctic-that is, 16-22% in the Barents Sea (Hegseth, 1998). Thus, fully comprehending HiLAT model shortcomings with respect to mixed layer nutrient concentrations will be important for understanding not only the overall response of production to climate change but its partitioning between the ice and ocean. HiLAT model estimate for annual Arctic (all analysis regions except the Bering Sea) ice-algae primary production (2.4 Tg C) were in line with estimates from the UAF-G model but significantly less than estimates from the other four model configurations (Watanabe et al., 2019). The HiLAT model and the UAF-G model are both based on the biogeochemical module (Jin et al., 2012;Moore et al., 2004), which has been incorporated into CESM.
The authors suspect the nitrate bias stems from underrepresentation of nitrate transported into the Arctic through the shallow Bering Strait. The Pacific inflow to the Arctic comprises three water masses with very different nutrient properties (Codispoti et al., 2013). Anadyr Water can have nitrate concentrations exceeding 25 mmol N m −3 , Bering Shelf Water has moderate prebloom nitrate concentrations (~15 mmol N m −3 ), and Alaska Coastal Water has low nutrient concentrations (prebloom nitrate <10 mmol m −3 ), such that the average prebloom nitrate concentration is~20 mmol N m −3 (Codispoti et al., 2009(Codispoti et al., , 2013. While the HiLAT model adequately simulates the nitrate concentration in the Bering Sea on a regional scale, and at 0.73 ± 0.31 Sv also reproduces the observed mean transport through the Bering Strait (Roach et al., 1995;Woodgate, 2018), the model presently underestimates the concentration of nutrients in the Bering Strait ( Figure 8). The nitrate concentration over the Northern Bering shelf is primarily dependent on in situ regeneration on the shelf, rather than from the seasonal entrainment of Bering slope waters (Granger et al., 2013). Thus, it seems likely that it is this regenerative process within the model that needs development and tuning. More fully assessing the nature and extent of the nitrate bias within the Arctic region should be a focus of model improvement in subsequent model iterations.

Concluding Remarks
The controlling influence of ice and temperature on Arctic ecosystem dynamics will likely vary through the year; however, our study has provided an improved understanding of the importance of these physical drivers in controlling production on an annual scale. We have shown that the relative control of ice and temperature over annual primary production and CO 2 flux is heterogeneous across the Arctic and that under a warmer reduced ice scenario, the additional carbon fixation by phytoplankton may be an order of magnitude larger than any change in the CO 2 flux. A fuller appreciation of the Arctic carbon dynamics will require a fully coupled atmosphere-ocean model with improved representation of nitrate dynamics.

Data Availability Statement
The E3SMv0-HiLAT source code will be made publicly available at github.com/lanl/E3SMv0-HiLAT website, once the review process is completed. The supporting processed model output for this manuscript is available at https://doi.org/10.5281/zenodo.3901591.