Self-Shading and Meltwater Spreading Control the Transition From Light to Iron Limitation in an Antarctic Coastal Polynya

Dotson Ice Shelf (DIS) in West Antarctica is undergoing rapid basal melting driven by intrusions of warm, saline Circumpolar Deep Water (CDW) onto the continental shelf. Meltwater from DIS is thought to influence biology in the adjacent Amundsen Sea Polynya (ASP), which exhibits the highest net primary productivity (NPP) per unit area of any coastal polynya in the Southern Ocean. However, the relative importance of iron and light in colimiting the spring phytoplankton bloom in the ASP remains poorly understood. In this modeling study we first investigate the mechanisms by which ice shelves impact NPP, then map spatio-temporal patterns in iron-light colimitation, and finally examine the environmental drivers of iron and light supply. We find that ice shelf melting leads to greater upper ocean iron concentrations, both directly due to release of iron from sediments entrained at the glacier bed, and indirectly via a buoyancy-driven overturning circulation which pulls iron from CDW to the surface. Both of these mechanisms increase NPP compared to experiments where ice shelf melt is suppressed. We then show that the phytoplankton self-shading feedback delays the bloom and reduces peak NPP by 80% compared to experiments where light penetration is independent of chlorophyll. Compared to light limitation, iron limitation due to phytoplankton uptake is more important a) later in the season, b)

Satellite ocean color measurements indicate that, per unit area, the Amundsen Sea Polynya (ASP) is the most productive coastal polynya in the Southern Ocean. Annual net primary productivity (NPP) has been estimated at 105 g C m −2 yr −1 , reaching a peak NPP of up to 2.5 g C m −2 day −1 at the height of the spring bloom (Arrigo & Van Dijken, 2003;Arrigo et al., 2012). The ASP is also noted for its close proximity to Dotson Ice Shelf, where satellite elevation measurements show rapid basal melting (Gourmelen et al., 2017). This concurrence of high melt rate with high productivity makes the ASP a suitable location to model the factors that dictate the magnitude and timing of phytoplankton blooms around Antarctica.
A phytoplankton bloom is light-limited when the flux of downwelling irradiance I reaching cells within the bloom is lower than that needed to maximize the rate of photosynthesis. At low irradiances the response of photosynthetic rate to increases in irradiance is approximately linear, with proportionality constant α. At higher irradiances this response becomes progressively less sensitive, as the photosynthetic reaction centers become saturated with photons. Attenuation by water molecules and other optically active constituents of the mixed layer means that the light available for phytoplankton varies significantly within a single phytoplankton community. There may be sufficient light for photosynthesis near the surface; however, the greater the depth over which phytoplankton are mixed, the less light is available on average for cells throughout the bloom. Hence there is an upper bound on the mixed layer depth, beyond which it cannot sustain a phytoplankton bloom. This critical depth D cr is the solution to the equation (1) derived by Sverdrup (1953). Here it is assumed that surface irradiance I 0 is low enough that the photosynthetic response is linear over the entire water column, and that light attenuation can be represented by a single constant k 0 . It is also assumed that wind-driven mixing removes any vertical gradients in biomass, and that the loss rate of phytoplankton λ is constant. In this context, light limitation arises due to a deep mixed layer (>D cr ) and/or turbid waters (high k 0 ) and/or weak surface irradiance (low I 0 ).
The surface irradiance at a given location is a product of solar angle, overlying cloud cover and, in the high latitudes, sea ice cover. Both the albedo and attenuation coefficients of sea ice are far larger than those of seawater, so the presence of sea ice greatly reduces the irradiance actually penetrating into the water column. Sea ice also has a strong influence on the seasonal cycle of mixed layer depth; an effect demonstrated in the Southern Ocean using a combination of seal, float and ship data (Pellichero et al., 2017). Brine rejection from newly formed sea ice destabilizes the water column, whilst sea ice melting provides a stabilizing layer of freshwater at the ocean surface. Where present, ice shelves are an additional source of freshwater to the water column. The mixed layer depth is then a function of these salinity forcings as well as wind-driven mixing and heat fluxes.
Dotson Ice Shelf (DIS) is the westernmost of a series of ice shelves in the Amundsen Sea. Sea ice is advected from east to west in front of these ice shelves by a strong coastal current, but Thwaites Fast-Ice Tongue (St-Laurent et al., 2017) prevents much of the sea ice formed upstream from entering the ASP. Thus in summer the ASP is consistently free of sea ice . Furthermore, mixed layers are generally shallow in summer, rarely exceeding 70 m Park et al., 2017), suggesting that it is instead water column turbidity which leads to light limitation in the ASP (Schofield et al., 2015).
When chlorophyll contributes significantly to overall light attenuation in the mixed layer via self-shading, k 0 can no longer be regarded as a constant. Instead k 0 will increase with increasing chlorophyll, leading to a self-shading feedback by which biomass near the surface limits the light available for deeper phytoplankton to photosynthesize (Vernet et al., 2008).
A phytoplankton bloom is iron-limited when the demand from photosynthesizing cells for bio-available iron exceeds the supply. Unlike light, the supply of iron at any given time is in part a product of how much has already been taken up by phytoplankton cells earlier in the bloom development. Numerical models can be used to construct iron budgets that show how physical and biogeochemical processes control iron supply over the course of a growing season. In general, iron availability in the upper ocean is expected to be at a maximum during winter when there is strong vertical mixing from deeper waters. Availability of iron then progressively declines as the bloom develops in summer. This gives rise to the carrying capacity hypothesis TWELVES ET AL.
10.1029/2020JC016636 2 of 28 (Hopkinson et al., 2013), whereby the inventory of iron just before bloom onset places a limit on the NPP which a given location can ultimately support. St-Laurent et al. (2017) applied this reasoning to the ASP within a purely physical model. A summertime drawdown of iron was designed to mimic biological uptake, and then different contributions to the wintertime iron inventory were assessed. Possible sources of iron to the ASP include sea ice, sediments, Circumpolar Deep water (CDW), icebergs, atmospheric dust, and glacial meltwater. Shipboard measurements in the ASP show a gradient of increasing iron concentration with increasing proximity to DIS, suggesting a source within the ice shelf cavity Gerringa et al., 2012). Accordingly St-Laurent et al. (2017) identified both sediments and glacial iron as important contributors to the wintertime iron inventory, with the contribution from the former mediated by the meltwater pump effect. This mechanism entrains iron-rich deep water into the buoyancy-driven overturning beneath the ice shelf, resulting in a redistribution of iron from depth to the upper ocean. A meltwater pump effect has already been identified as a key driver of high NPP around Greenland, where it serves to bring limiting macronutrients to the surface (Cape et al., 2019).
Further modeling of the ASP St-Laurent et al., 2019) has resolved the principal biogeochemical components of the dissolved iron budget: uptake, scavenging and remineralization. The westwards flowing coastal current in the Amundsen Sea was shown to be an important driver of seasonal iron cycles. Oliver et al. (2019) undertook extensive optimization of biogeochemical parameters, making use of datasets from the ASPIRE research cruise (Yager et al., 2012) (see values listed in Table 1). The resulting one dimensional model was able to provide a good fit to data for multiple stations, and reproduced many important features of the 2010/2011 bloom in the ASP. Importantly, a transition from light limitation to iron limitation was observed over the course of the season.
In the ASP, there is strong evidence for colimitation of productivity, whereby the relative importance of iron versus light varies with: • Time: As in Oliver et al. (2019), depletion of iron stocks over the growing season may lead to a temporal succession from light to iron limitation. This is difficult to verify in the ASP due to the narrow time window for which it is accessible to research cruises. However Arrigo et al. (2017) do observe this succession from light to iron limitation near the West Antarctic Peninsula. • Depth: Light availability decays exponentially with depth, whilst the iron concentration will often show a near-surface minimum. Hence iron limitation may dominate near the surface even whilst light limitation dominates deeper in the mixed layer. • Location: Large horizontal gradients in sea ice cover, mixed layer depth and external iron supply may lead to differences in limitation between different regions of the polynya. Alderkamp et al. (2015) measured greater iron limitation in the central ASP, despite its higher productivity compared to coastal waters.
In addition, the demands of phytoplankton for iron and light are codependent. Phytoplankton cells demand iron in part to maintain their photosynthetic apparatus (Strzepek & Harrison, 2004). When ambient iron concentrations are below those needed for efficient photosynthesis, and cells are thus in a state of iron stress, light limitation can become more severe (Geider et al., 1997). An important marker of iron stress in phytoplankton is the ratio of variable fluorescence to maximum fluorescence (Fv/Fm), a method used in the ASP by Alderkamp et al. (2015) and Park et al. (2017). More generally, the various light and nutrient requirements of phytoplankton communities are often intricately coupled rather than being independent of each other (Saito et al., 2008).
Within the consensus that the phytoplankton bloom in the ASP is co-limited by both iron and light, it remains important to study whether the overall productivity of the polynya-and therefore its capacity to sequester carbon (Yager et al., 2016)-is in fact more sensitive to one variable than the other. Park et al. (2017) conduct an intercomparison of two polynyas: the ASP and nearby Pine Island Polynya (PIP). Results show that compared to the PIP, the ASP exhibits both more severe iron stress and higher productivity. This apparent paradox can be resolved if the Amundsen Sea is primarily light-limited, since reanalysis data show lower cloud cover above the ASP, leading to surface irradiances up to 15% greater than in the PIP (Park et al., 2017). Conversely Arrigo et al. (2015) complete a continental scale survey of coastal polynyas and identify basal melting as the key driver of NPP, by supplying iron to phytoplankton communities where it is strongly limiting.
TWELVES ET AL.
10.1029/2020JC016636 3 of 28 The rapid basal melting of DIS is driven by intrusions of warm, saline CDW onto the continental shelf via Dotson Trough (Jacobs et al., 2012). These intrusions are in turn coupled to the wind and sea ice conditions at the shelf break, with variability at the shelf break thus leading to variation in basal melt rate (T. W. Kim et al., 2017). However, for the melting which takes place beneath the ice shelf to have an impact on biological processes in the polynya, meltwater must first undergo horizontal spreading away from the shelf (Schofield et al., 2015). Modeling is required to elucidate the link between melt rate and meltwater pathways, with recent results from Kimura et al. (2017) indicating that greater melting does not necessarily lead to higher meltwater concentrations away from the continent. Since melt rates are expected to increase under future climate change scenarios, it is important to understand the sensitivity of coastal polynya ecosystems to the changes in iron and light availability which will follow.
In this study, we use the Biology Light Iron Nutrients Gases (BLING) model developed by Galbraith et al. (2010) with no a-priori assumption of iron limitation, but with boundary and initial conditions sourced from the Biogeochemical Southern Ocean State Estimate (B-SOSE) (Verdy & Mazloff, 2017). We modify the BLING model to include a parametrization of self-shading based on that of Manizza et al. (2005), which can be turned on or off at runtime. In addition, we conduct experiments wherein we artificially suppress the TWELVES ET AL. , which allows us to make comparative statements about iron and light limitation on a regional scale. We include an ice shelf of fixed size and shape within the domain in order to produce both a meltwater pump effect and an external source of iron from the glacier itself. The size and shape of the polynya evolves according to a thermodynamic and windblown sea ice model.
We address four main research questions: 1. How do the meltwater pump and glacial iron supply impact productivity in the ASP? 2. What effect does self-shading have on the distribution and timing of the phytoplankton bloom? 3. How do iron and light limitation constrain NPP over the course of the summer growing season? 4. How does productivity in the ASP respond to changes in CDW intrusion and overlying cloud cover?
We structure this paper as follows. We present our methods in Section 2, our results in Section 3 and a discussion of these results in Section 4. In Section 5, we link our study to wider implications for biogeochemical cycles in the Southern Ocean; furthermore, we explain how our methods might inform future modeling studies. Finally in Section 6 we present our conclusions.

Physical Model
We use MITgcm to model ocean physics, with a domain forming an idealized version of Dotson Ice Shelf (DIS) and the ASP (Figure 1). We employ a Cartesian grid of 1 km horizontal resolution, extending 150 km in the zonal and meridional directions. We represent DIS with a static ice shelf bordered to the west and east by blocks of land representing Martin Peninsula and Bear Peninsula respectively. The vertical resolution is variable, ranging from 1 m at the surface to 64 m at the bottom of Dotson Trough. Vertical mixing using the K-profile parametrization (KPP) (Large et al., 1994) acts on temperature, salinity and biogeochemical variables. The trough extends from the north west corner of the domain to the southern boundary under the ice shelf. It has a base depth of 950 m, with sides sloping up to the 400 m sill outside the trough. There is an open northern boundary, whilst the zonal boundaries are kept periodic. This means that the model outputs represent the result of a series of adjacent ice shelves, analogous to Pine Island Glacier, Crosson Ice Shelf, Thwaites Glacier and DIS, which run from east to west along the coastline of the ASP.
We represent the ice shelf in an idealized fashion, specifically as a wedge with thickness increasing from 100 m at the cavity front to a maximum of 800 m at the southern boundary. The ice is grounded where it meets the sloping sides of the trough, but overlying the deepest part of the trough (where the bathymetry is at 950 m) the grounding line lies beyond the southern edge of the domain. Ice shelf thermodynamics is based on the three equation formulation of D. M. Holland and Jenkins (1999) for heat and salt fluxes across the ice-ocean boundary. This treatment fixes the ice-ocean interface at the local freezing temperature T B (°C), which in turn depends on the local salinity S B (psu) and pressure p B (Pa) via Both of these fluxes are influenced by the interface temperature. The MITgcm ice shelf package that we use generates meltwater by way of a freshening tendency in the grid cells adjacent to the shelf (Losch, 2008;Dansereau et al., 2014;P. R. Holland et al., 2008). Following Schodlok et al. (2012), we keep friction velocities constant, so that ice shelf melt is only dictated by temperature and not by ocean velocities at the interface. The freshening tendency at the top of each water column underneath the ice shelf acts as a source for meltwater tracers, two of which are used in the following investigations. These are sourced identically at the ice shelf, but only one of them (the "global" tracer) is allowed to pass through the periodic ocean boundaries, reentering the domain from the east. The other, "local" tracer, is relaxed to zero at periodic boundaries and thus allows us to isolate meltwater from DIS without the influence of upstream sources in the Amundsen Sea.
We use the sea ice package of Losch et al. (2010), with dynamics driven by wind fields and by circulation in the top level of the ocean model. A mutual drag acts between the ice floe and ocean surface. We impose a constant westwards sea ice velocity at the northern boundary to prevent artificial buildup of ice. The sea ice thermodynamics is based on the assumptions of a zero layer model (Semtner Jr, 1976), namely that the ice has uniform conductivity and zero heat capacity. We do not include precipitation in our model forcings, so the layer of snow with different conductivity and albedo as modeled in Zhang et al. (1998) is absent.
Work by Stammerjohn et al. (2015) suggests that polynya formation in the Amundsen Sea depends on the presence of the Thwaites Fast-Ice Tongue, and we found the same behavior in our preliminary investigations (not shown). Without any obstruction to westwards moving sea ice, the polynya failed to open in spring. Hence in all subsequent model runs we add a 2 m thick portion of ice shelf (outlined in red in Figure 1) at the western boundary, serving as an obstacle to sea ice reentering the domain from the east. This ice tongue is a negligible source of meltwater, but it succeeds in maintaining a small region of open water adjacent to the ice shelf from which the polynya can develop in spring.
At the northern boundary, the temperature and salinity are relaxed toward prescribed values on a one week timescale ( Figure 2). These values are informed by profile outputs from Kimura et al. (2017) on a transect ∼100 km from the front of DIS, but are kept constant over the course of the year. Therefore, we neglect the seasonal fluctuations in onshore CDW transport thought to influence cycles in basal melt rate . Surface forcing of the model is via monthly fields for 2 m air temperature and humidity, 10 m winds, longwave and shortwave downwelling radiation. Our monthly fields for temperature, humidity and radiation are based on the climatological forcings used in Petty et al. (2013), which in turn derive from NCEP-CFSv1 reanalysis data (Saha et al., 2006). In contrast, we treat winds as constant over time. This is in order to simplify our analysis of couplings between thermocline depth, melt rate and sea ice cover. Constant winds imply that there is no seasonal variation in upwelling and downwelling across the domain (apart from that variability associated with overlying sea ice cover). The winds are south-easterly, with magnitude zonally constant but decreasing with meridional distance from the ice shelf (T. W. Kim et al., 2017).  of the nine tracers included in BLING impact phytoplankton growth directly: nitrate, phosphate, and iron. Two tracers (alkalinity and dissolved inorganic carbon) represent carbon chemistry, whilst a further tracer represents carbon in the form of phytoplankton biomass. The remaining three tracers are for dissolved oxygen, dissolved organic nitrogen and dissolved organic phosphorous.

Biogeochemical Model
The most critical parameters describing our setup of BLING are listed in Table 1, with particular focus on those parts of the model where we have modified the existing BLING routines.

Phytoplankton Growth and Transport
The light saturated (per capita) photosynthesis rate Pc max (units of day −1 ) is adjusted by a dimensionless light limitation term to give the per capita growth rate μ: Here again I (units of W m −2 ) is the photosynthetically available radiation (PAR). The light saturation parameter I k (also W m −2 ) is sensitive to local iron concentration as where α now has units of (W m −2 ) −1 day −1 . I k is also sensitive to the light level to which phytoplankton are acclimated (Galbraith et al., 2010;Geider et al., 1997), henceforth termed the irradiance memory. This term I mem (W m −2 ) is defined in full in Section 2.3. The light saturated photosynthesis rate is itself a function of temperature T and the degree of nutrient limitation n lim : where μ max = 1.47 days −1 is the maximum growth rate under light and nutrient replete conditions at 0°C, and κ = 0.063°C −1 is a constant (Eppley, 1972). The nutrient limitation term is calculated according to Liebig's law of the minimum applied to nitrate (NO 3 ), phosphate (PO 4 ) and dissolved iron (Fe) concentrations: The saturation parameters for nitrate (K N = 2.05 mmol N m −3 ), phosphate (K P = 10.25 μmol P m −3 ) and iron (K Fe = 0.16 μmol Fe m −3 ) are constants.
In phytoplankton acclimated to very low light I mem in Equation 6 tends to zero. If the incident irradiance is then increased from zero, I will at first be much smaller than I k . Applying these two conditions and substituting Equation 6 into Equation 5, we see that α is therefore equivalent to the initial slope of the photosynthesis-irradiance curve. This can be broken down further as a product of a chlorophyll specific response α chl (in units of g C (mg Chl) −1 ) day −1 ) and the chlorophyll:carbon ratio θ dark (in units of mg Chl (g C) −1 ) at low light: When the iron concentration is low relative to the half saturation constant K Fe , and phytoplankton cells are undergoing iron stress, chlorophyll production is reduced and the cell photosystem is inhibited. Hence the values of θ dark and α chl in the model vary between minimum (min) and maximum (max) values as set out in Table 1, modulated by the local iron concentration: The impact of light saturation on the chlorophyll:carbon ratio is described in Geider et al. (1997) and Galbraith et al. (2010).
Increases in biomass due to phytoplankton growth are balanced by depletion due to grazing, which is based on a fixed per capita grazing rate λ = 0.19 days −1 . Changes in biomass of the different phytoplankton classes (large, small, and diazotroph) are treated separately at each time step. As described in Verdy and Mazloff (2017) and Dunne et al. (2005), the model is designed so that a bloom of large phytoplankton undergoes less severe grazing pressure relative to its developing size than an equivalent bloom of small phytoplankton. In the context of ASP modeling, the large phytoplankton class represents various diatom species and the small phytoplankton class represents the haptophyte Phaeocystis antarctica.
This study makes use of the advected phytoplankton tracer parametrization in BLING. A single biomass tracer is transported around by the physical model, whilst the fractions of large, small, and diazotroph phytoplankton are retained within BLING. Then at each time step the updated biomass tracer is re-partitioned according to the fractions calculated at the preceding time step. Our preliminary investigations (not shown) exhibited unrealistic accumulation of diazotrophs in the cold waters of the ASP. Therefore, in this study, the diazotroph fraction was explicitly removed from the biogeochemical model, leaving only the large and small fractions; these are treated identically to Verdy and Mazloff (2017).
Finally, NPP is calculated as the product of biomass and growth rate: where B lg and B sm represent large and small phytoplankton biomass respectively in units of g C m −3 ; hence NPP is measured in units of g C m −3 day −1 .

Iron Budget
In this section, we show the closure of the upper ocean iron budget in BLING. The rate of biological iron uptake (Fe upt ; units of μmol Fe m −3 day −1 ) varies with phytoplankton growth rate, phytoplankton biomass, and the iron:nitrate uptake ratio σ (units of μmol Fe (mmol N) −1 ): where the carbon:nitrogen ratio η = 0.081 g C (mmol N) −1 is constant. In contrast, σ is itself a function of the iron concentration: Here σ min = 0.014 μmol Fe (mmol N) −1 and σ max = 0.17 μmol Fe (mmol N) −1 are the minimum and maximum iron:nitrate ratios respectively, with K σ = 0.82 μmol Fe m −3 denoting the half saturation constant for iron uptake. A proportion of this iron uptake is then exported downwards through the water column in detrital form. This proportion is not constant; it is instead a function ϕ sink of temperature, biomass and the ratio of large to small phytoplankton fractions, as derived in Dunne et al. (2005) and adapted for BLING in Galbraith et al. (2010) and Verdy and Mazloff (2017).
The flux of sinking particulate iron is supplemented further by scavenging, whereby free dissolved iron becomes adsorbed onto particles and sinks in colloidal form. The equations used to parameterize scavenging in BLING are described in full in Galbraith et al. (2010). As this particulate iron sinks through the cell, a portion is remineralized into dissolved iron. Any particulate iron which is not remineralized contributes to a net imbalance between the particulate iron flux where Δr (measured in m) is the depth range covered by the cell. The remineralization length scale z remin (measured in m −1 ) is a function of oxygen (O 2 ) saturation: where γ POM = 0.12 day −1 is the decay timescale for particulate organic matter, K O2 = 20 mmol O 2 m −3 is the half saturation constant for oxygen and r min = 0.15 is a constant. The sinking speed w sink (units of m day −1 ) is parameterized identically to Galbraith et al. (2010). In acting as a sink of particulate iron, remineralization is also a source Fe remin (μmol Fe m −3 day −1 ) of dissolved iron: A further source of dissolved iron is recycling, whereby the microbial loop returns a quantity Fe rec (μmol Fe m −3 day −1 ) of the iron taken up by phytoplankton back into the dissolved pool. In BLING, this process is assumed to happen instantaneously, so that at each time step the quantity of iron recycled is simply equal to the difference between uptake and particulate export: Since iron is conserved within each cell, the sum of exchanges between the dissolved and particulate iron pools must be zero, unless there is an external flux from the physical model. Therefore, the conservation equation, holds throughout the domain. Here ΔDFe (units of μmol Fe m −3 day −1 ) is the overall change in dissolved iron concentration, whilst ΔDFe phys (also μmol Fe m −3 day −1 ) is the change in dissolved iron concentration due to advection, diffusion and vertical mixing calculated within the physical model.

Substituting Equations 17-19 into Equation 15
and re-arranging, we arrive at a description of the dissolved iron budget in the upper ocean: In this study, we do not include an air-sea flux of iron from dust. At the sea floor there are additional terms corresponding to exchange of iron with sediments, but these do not feature in the upper ocean budget. In the remainder of this paper, we will limit our analysis of the iron budget to biogeochemical processes only, thereby treating ΔDFe phys as a residual.

Irradiance and Irradiance Memory
Our implementation of phytoplankton self-shading in MITgcm makes use of the bio-optical model of Manizza et al. (2005) as in previous earth system modeling (G. Kim et al., 2015;Manizza et al., 2008). PAR is split into two bands of approximately equal power following Foujols et al. (2000). Attenuation coefficients k red (m −1 ) and k bg (m −1 ) for the red and blue-green components respectively are calculated as: are the same as in Manizza et al. (2005) and are based on the study of Morel (1988). This treatment does not resolve individual scattering or absorption processes, but assumes that the number of photons actually used in photosynthesis is small compared to those otherwise attenuated by chlorophyll. Thermodynamically, this implies that self-shading should contribute significantly to the heat budget in ocean cells close to the ocean surface, but we do not implement this feedback in the current study. This approach allows us to test different versions of the biogeochemical model whilst maintaining a fixed set of physical properties for the water column.
As phytoplankton undergo mixing within the euphotic zone, they sample a range of different light levels. Therefore there is generally a mismatch between the instantaneous irradiance available at a given depth and the irradiance memory (Galbraith et al., 2010;Geider et al., 1997). Furthermore, the latter is a property associated with individual phytoplanton cells, making it difficult to capture in a Eulerian model.
We proceed by considering two timescales associated with mixing in the euphotic zone, using insights from Kida and Ito (2017) and from our own Lagrangian one-dimensional modeling (not shown). The first is the mixing timescale, τ mix , which we define as the average time taken for a phytoplankton cell to move over the full range of depths in the mixed layer. The second timescale is the acclimation timescale, τ acc , which represents the speed at which photosynthetic apparatus adapt to new light conditions. Based on Morris (1980), we make the assumption within the mixed layer, and the assumption, below the mixed layer. This describes a mixed layer where between time t − τ acc and time t each phytoplankton cell has been transported over the full depth of the mixed layer multiple times; below the mixed layer each cell remains at approximately the same depth between time t − τ acc and time t. Therefore, the irradiance memory becomes a bulk property within the mixed layer, and I mem evolves according to the equation where I mld is the irradiance averaged over the depth of the mixed layer at time t.

Description of Experiments
To obtain our biogeochemistry boundary conditions we average the B-SOSE outputs over the period 2008-2012 on a monthly basis, for the region covered by 72.9°-74.8°S, 110°-115°W. We then transpose these fields onto our higher resolution vertical grid. We apply a mask to the 10 cells nearest the Northern boundary, and here relax BLING tracers toward their respective boundary conditions. The exception is iron, which we do not prescribe on the boundary above 600 m in order to prevent artificial relief from limitation in the iron-depleted waters of the domain interior.
We first spin up the physical model for 18 months to reach quasi-equilibrium with external forcings before BLING is enabled. The eight core BLING tracers (i.e. dissolved inorganic carbon, alkalinity, nitrate, phosphate, oxygen, iron, dissolved organic nitrogen and dissolved organic phosphorous) are initialized based on 2008-2012 B-SOSE outputs for the month of June. Starting from midwinter allows us to initialize the biomass tracer close to zero. With biogeochemistry enabled we spin up for a further year, taking our results from the year beginning 30 months after model initialization. We verified that the model was spun-up using time series in ocean heat flux, ice shelf melt rate and iron budget (not shown), which were approximately equal in the period with months 31-42 as in the period with months 19-30. In Table 2, we set out the 14 different experiments which we conduct.
The first question we address is whether the presence of the ice shelf has an impact on biology in the polynya. We answer this by evaluating first the role of the meltwater pump on the physical model, then by adding an iron tracer into the glacial meltwater. Since the meltwater pump effect arises from the positive buoyancy of freshwater underneath the ice shelf, we isolate its impact by carrying out an experiment (no_melt) where the freshening tendency from melting is suppressed. We achieve this suppression of the melt rate in the model code by explicitly setting the fluxes latent T Q and brine S Q to zero. The melt_pump experiment meanwhile has these fluxes calculated as per Equations 3 and 4. Thus the difference in outputs between melt_pump and no_melt shows a) the impact of the meltwater pump on circulation within the domain and b) the impact this change in circulation has on biogeochemistry in the polynya.
In melt_pump the iron concentration in meltwater is identically zero, but for the next experiment-gmw_ iron-we follow St-Laurent et al. (2019) in specifying a 20 μmol m −3 concentration for iron in glacial meltwater. This value is similar to that inferred from measurements in the Amundsen Sea Miles et al., 2016). In all other regards the setup of gmw_iron is identical to melt_pump. Hence the difference in outputs between gmw_iron and melt_pump shows the impact of iron originating from glacial meltwater on biogeochemistry in the polynya. Furthermore, for the gmw_iron experiment, we employ an additional iron tracer, which is relaxed to zero at the zonal boundaries by the same method as the local meltwater tracer. This allows us to quantify the importance of iron from DIS versus ice shelves further upstream.
Next we look to investigate the significance of iron limitation due to phytoplankton uptake and of light limitation due to phytoplankton self-shading. The flat_iron, fixed_zeu, and max_yield experiments all have strictly identical physics to gmw_iron, but implement different versions of the biogeochemical code. For the flat_iron experiment, we set σ from Equation 13 to zero, so that phytoplankton uptake no longer serves as a sink of dissolved iron in the model. Since recycling is calculated as a fraction of uptake, this is also fixed at zero. As a consequence, the upper ocean iron budget in flat_iron represents a balance between TWELVES ET AL.  remineralization, scavenging and physical transport processes only. We interpret the different NPP observed in gmw_iron as compared to flat_iron to be the result of increased iron limitation arising from phytoplakton uptake.
We investigate the impact of using the Manizza et al. (2005) self-shading scheme by comparison against an experiment (fixed_zeu) which employs the constant light extinction profiles used in Verdy and Mazloff (2017). Specifically, the chlorophyll-dependent k red and k bg coefficients used in gmw_iron are replaced in fixed_zeu with a single, chlorophyll independent attenuation constant k 0 . A spatially and temporally constant light extinction profile implies a constant euphotic depth (the depth at which PAR is 1% of its surface value). For the fixed_zeu experiment we use k 0 = 0.04 m −1 , giving a euphotic depth of 111 m. Thus the difference in irradiance profiles and euphotic depths from gmw_iron as compared to fixed_zeu shows the effect of the modeled chlorophyll concentrations on the polynya light environment. We interpret the different NPP observed in gmw_iron as compared to fixed_zeu to be the result of increased light limitation arising from phytoplankton self-shading.
We conclude this part of the study by investigating a hypothetical phytoplankton bloom which is neither strongly iron-limited nor strongly light-limited. For this max_yield experiment we both set σ = 0 μmol Fe (mmol N) −1 as in flat_iron and use k 0 = 0.04 m −1 as in fixed_zeu. The result is expected to be a much more productive bloom than in gmw_iron. In summary, these three experiments (flat_iron, fixed_zeu, and max_ yield) allow us to switch iron and light sinks on and off, helping us to infer the relative importance of iron and light limitation at different locations and timings within the bloom.
Finally, we look to investigate how the ocean, ice shelf and phytoplankton bloom respond to changes in environmental conditions. We do this by varying external forcing at the boundaries of the model domain ( Figure 2). In total, we perform eight sensitivity experiments in addition to gmw_iron, which serves as our base case. The relevant boundary conditions, plus their values in gmw_iron (and all other experiments named thus far) are as follows: 1. At the northern (open) boundary we impose a thermocline between 400 m and 600 m depth. Above the thermocline there is a layer of (cold and fresh) winter water with temperature increasing from −1.8°C at the top of the thermocline to −1.6°C at the surface. Salinity decreases from 34.1 psu at the top of the thermocline to 33.9 psu at the surface. Below the thermocline there is a layer of modified CDW, which is warm (0.6°C), highly saline (34.5 psu) and homogenous. Water properties are interpolated linearly over the depth of the thermocline between the two layers. 2. At the ocean surface we impose a seasonal cycle of downwelling shortwave irradiance. Based on the NCEP-CFSv1 re-analysis (Petty et al., 2013;Saha et al., 2006), the peak irradiance, applied in January, is 320 Wm −2 . Irradiance falls to a minimum of 10 Wm −2 in polar winter.
Observations from the Amundsen Sea (Dutrieux et al., 2014;Randall-Goodwin et al., 2015;Sherrell et al., 2015) indicate the magnitude of interannual variability in CDW intrusions onto the continental shelf.
Based on these studies we perform experiments (named with the prefix warm−) where the thermocline (and the halocline) at the northern boundary is situated between 250 m and 450 m depth, to approximate years with larger than average intrusion of CDW. Similarly we perform experiments (prefixed cold−) where the thermocline (and the halocline) is situated between 550 m and 750 m to approximate years with smaller than average intrusion of CDW onto the continental shelf. The prefix base− refers to experiments conducted with the "base case" CDW condition that is, a thermocline (and a halocline) between 400 m and 600 m. We expect to see a positive correlation between the quantity of heat transported onshore by CDW and the rate of basal melting from DIS.
Recent work by Park et al. (2017) attributed large differences in NPP between the ASP and PIP to differences in surface irradiance on the order of 15%. To examine whether changes in surface irradiance of this magnitude can have large impacts on productivity in the ASP, we perform experiments where the amplitude of the seasonal cycle in downwelling shortwave radiation was either increased or decreased by 15%. In the low irradiance experiments (suffixed -low) the surface irradiance peaks in January at 272 W m −2 ; in the high irradiance cases (suffixed -high) it peaks at 368 W m −2 . The suffix -med then refers to experiments with the original irradiance forcing. With these three different regimes for surface irradiance, in combination with the three different regimes for the northern thermocline, we have nine different iterations of the model (including the base case gmw_iron). The annual productivity which results from each of these different setups will indicate whether known variability in surface or deep ocean conditions has a greater impact on biology in the ASP. Further we will ascertain whether the hypothesized relationships between NPP and cloud cover (Park et al., 2017) and between NPP and basal melt rate  are reproduced in our model.

Impact of Ice Shelf on Physics and Biogeochemistry
In order to compare our melt_pump results with observations, we examine a transect along the front of the ice shelf in January (Miles et al., 2016;Randall-Goodwin et al., 2015). At the ice shelf front, the hydrography resembles the northern boundary condition, with a warm and saline CDW layer below approximately 600 m (Figures 3a and 3b). Near the surface the effect of summertime warming by shortwave radiation is visible, with temperatures as much as 0.87°C above the subsurface minimum of −1.75°C. This surface warming is weakly visible in the shelf front transect by Miles et al. (2016) (their Figure 5a) and more strongly in their transects away from the ice shelf (their Figures 3a and 4a). At the western edge of the cavity we observe a strong warming anomaly (Figure 3c) relative to the depth averaged temperature, as in Randall-Goodwin et al. (2015) (their Figure 7a). This is due to ice shelf melting, which also drives the outflow velocities in Figure 3d. The maximum modeled outflow velocity is 0.14 m s −1 , whilst the maximum inflow velocity (at the eastern edge of the cavity) is 0.38 m s −1 . Compared to Randall-Goodwin et al. (2015) (their Figure 7b) the modeled outflow is both weaker and more diffuse, whilst the modeled inflow is stronger.
Buoyant meltwater exiting the cavity rises through the water column and is pulled to the west by the strong boundary current. Figure 4a shows the column depth of the local meltwater tracer, which is concentrated near the coast and is almost entirely confined to the western half of the domain. Meanwhile, Figure 4b shows the additional meltwater which enters the domain from the east via the zonal periodic boundary, calculated by subtracting the local meltwater tracer from the global meltwater tracer. This latter figure shows that our model can produce large quantities of meltwater in the eastern half of the domain, but only due to recirculation via the periodic zonal boundaries. This in turn suggests that ice shelves upstream of DIS may dominate the meltwater supply to the eastern ASP. However the lack of a strong sink for the global meltwater tracer (in contrast to the local meltwater tracer, which has a strong sink at the zonal boundary) means that, unlike Figure 4a, the distribution in Figure 4b has not reached a steady state. Henceforth we therefore use only the local meltwater tracer to quantitatively compare against data from the ASP.  slight underestimate. However, the spatial patterns remain similar, with meltwater focused in a small region extending from the western edge of the shelf. Randall-Goodwin et al. (2015) also show the isopycnal at which each meltwater maximum occurs (their Figure 12b), whilst we plot the depth of the meltwater maximum in Figure 4d. The accordance between these measures is less close, principally because in our model the accumulation of meltwater at depth is generally greater than at the near-surface peak. The small region which does show a near-surface peak in Figure 4d corresponds to the most intense region of surface meltwater. The near-seafloor peak is the result of meltwater reaching neutral buoyancy at depth and spreading along the contours of the bathymetry.

Randall
In order to assess the impact of ice shelf melting on the physical state of the ASP, we compare the melt_pump and no_melt experiments. Figure 5a shows that sea surface temperature (SST) peaks in February in both experiments, with only a small (0.1°C) decrease in maximum SST due to melting. The April peak in melt rate two months later suggests that the downwelling of warmed surface waters is responsible for melting at the front of the ice shelf, a mechanism previously observed at the Ross Ice Shelf (Stewart et al., 2019). The cycles in sea ice coverage and spatially averaged mixed layer depth with and without melt are shown in Figure 5b. The presence of a melt-driven circulation reduces maximum wintertime sea ice coverage from 88% to 79% of the ocean domain. The impact of ice shelf melt on horizontally averaged mixed layer depth is most visible at the start and end of summer, where it is likely due in part to the aforementioned differences in sea ice distribution between no_melt and melt_pump.
Unlike the meltwater tracer, the iron tracer is periodically depleted on a seasonal basis due to biogeochemical processes, and thus does not accumulate over successive years. Therefore we continue to include iron sourced from upstream ice shelves in integrating the global iron tracer to produce the time series in Figure 6. The importance of iron supplied from upstream ice shelves in our results ( Figure 16) reflects previous modeling studies ( St-Laurent et al., 2019). The near-surface iron pool follows a strong seasonal cycle in the melt_pump case, with a peak in November. The wintertime iron inventory is depleted by 74% over the course of the bloom, reaching a minimum in April. With gmw_iron there is an additional (glacial) source of iron, resulting in a 49% increase in wintertime iron inventory. Similar to melt_pump, this is then depleted by 72% to an April minimum.
The annually integrated NPP reaches 41 g C m −2 yr −1 in the melt_pump case, peaking at 0.39 g C m −2 day −1 in early December ( Figure 6). This compares to annual and peak NPP values of 1 g C m −2 yr −1 and 0.01 gC m −2 day −1 in the no_melt case, with the peak occurring in February. Thus the meltwater pump brings forward the spring bloom by around 2 months and causes a 40-fold increase in productivity. This is attributable to a similarly large increase in the pool of dissolved iron available within the top 100 m. Time series of NPP for melt_pump and gmw_iron show similar patterns with the initial December peak followed by a secondary peak in January. The additional iron available in gmw_iron causes peak NPP to increase by 28% to 0.51 g C m −2 day −1 and annual NPP by 34% to 55 g C m −2 yr −1 . That the regions of highest NPP map onto the regions of highest iron concentrations is detailed in Section 3.4, and can be seen in Figure 16. Results from the gmw_iron experiment shows high iron concentrations and high NPP closer to the coast, within a relatively narrow coastal current.

Impact of Self-Shading
The inclusion of the Manizza et al. (2005) formulation for self-shading in the model leads to reduced summertime light penetration. However there is an overall deepening in the annual-and spatially-averaged euphotic depth from 111 m to 118 m (Figure 7). This is due to lower light attenuation in winter, months when the water column is free of chlorophyll. The spatially averaged euphotic depth calculated in gmw_iron shoals from a maximum of 135 m in winter to a minimum of 95 m in summer. Therefore, in fixed_zeu, the euphotic depth is underestimated by as much as 18% in winter and overestimated by up to 17% in summer. There is also significant variation in summertime euphotic depth across the domain when self-shading is enabled, with a range of over 100 m between different locations. There are no differences in mixed layer depth between gmw_iron and fixed_zeu due to the identical physics of these experiments, and mixed layer is consistently shallower than the euphotic depth.
The reduced summertime light availability due to self-shading causes a reduction in December phytoplankton growth rate throughout most of the water column in gmw_iron as compared to fixed_zeu (Figure 8a). However above approximately 40 m self-shading leads to an increase in phytoplankton growth rate. Similarly self-shading leads to a shoaling of the December deep chlorophyll maximum (DCM) (Figure 8b) from 70 m in fixed_zeu to 30 m in gmw_iron. In both experiments the large phytoplankton fraction dominates biomass at and above the DCM with the smaller fraction dominating below (Figure 8c). December NPP follows a similar profile to growth rate, with gmw_iron more productive near the surface (Figure 8d). However a profile of annual integrated NPP (Figure 8e) shows that on a yearly basis self-shading leads to a reduction in NPP throughout the entire water column. The anomaly in NPP between gmw_iron and fixed_zeu is represented in a Hovmoller plot (Figure 8f), which shows the emergence of a positive surface anomaly following a strong negative anomaly earlier in the season. Hence the higher growth rates in fixed_zeu early in the season lead to greater depletion of iron near the surface, but this depletion in turn leads to lower surface growth rates late in the season. Figure 9 demonstrates the impact that self-shading has on the vertical distribution of biogeochemical tracers and on their seasonal cycles. The reduced uptake of nutrients in the upper ocean in gmw_iron leads to elevated nitrate and reduced oxygen close to the surface, with contrasting patterns at depth. The relative scarcity of iron means that iron concentrations are more sensitive to the degree of uptake in the system, so that the anomalies in iron concentration between fixed_zeu and gmw_iron can be greater than an order of magnitude near the surface.

Iron-Light Colimitation
In order to understand the seasonal cycling of iron, and how this impacts iron-light colimitation, we consider the four key processes which BLING calculates in the upper ocean: uptake, remineralization, scavenging and recycling (Equation 20). Following this approach, the budgets in Figure 10 do not include the advection and diffusion processes that are calculated in the physical core of the MITgcm or the vertical mixing calculated using KPP. In the melt_pump and gmw_iron cases, the primary influences on the biological tendency of iron are the removal of dissolved iron by phytoplankton uptake and the addition of dissolved iron via remineralization of organic material. Both uptake and remineralization peak during December. Contributions from scavenging and recycling are smaller, with the latter showing a peak around February, after the peak of the bloom. Uptake is 51% higher in gmw_iron as compared to melt_pump, TWELVES ET AL.   in line with the increased supply of iron from the ice shelf. For the flat_ iron experiment iron uptake and recycling were suppressed as outlined in methods, but there is still a small biological tendency due to scavenging and remineralization processes.
In Figures 11a, we examine the relative importance of iron and light in dictating NPP using time series from gmw_iron, fixed_zeu, flat_iron, and max_yield. The latter shows annual NPP an order of magnitude greater than that for gmw_iron, illustrating the large combined impact of iron uptake and phytoplankton self-shading in reducing the magnitude of the spring bloom. For all four experiments NPP remains close to zero until mid-October due to low light levels characteristic of Antarctic winter. At this point the time series diverge, with both max_yield and fixed_zeu beginning to show an increase in productivity. The spring bloom does not commence until approximately one month later for flat_iron and gmw_ iron, both of these being experiments where self-shading is enabled. As the bloom progresses, depletion of the initial iron pool due to phytoplankton uptake becomes cumulatively more severe. Hence NPP in the flat_ iron case first exceeds that in gmw_iron, then eventually surpasses that for fixed_zeu. As the bloom declines, the flat_iron experiment remains more productive than gmw_iron due to ongoing iron depletion in the latter.
We identify the time at which flat_iron becomes more productive than fixed_zeu as a transition from light to iron limitation. Figures 11b shows how the timing of this transition varies spatially, with flat_iron first surpassing fixed_zeu in the center of the domain. Meanwhile close to the peninsulas fixed_zeu remains more productive than flat_iron for most of the year.

Sensitivity to Thermocline Depth and Surface Irradiance
The depth of the thermocline and the strength of surface irradiance both impact ice shelf melt rates.     cords with our earlier highlighting of the coupling between summertime heating of the surface layers and frontal ice shelf melt; quantitatively we find that increases in this frontal ice shelf melt rate are proportional to the square of the respective increases in SST.
The maps in Figure 13 show spatio-temporal patterns in the timing of spring sea-ice retreat for each of the nine sensitivity experiments. For the base− and warm− cases, the sea ice cover retreats from both the North and South, with a band through the center of the ocean domain the last to retain ice cover greater than 15%. This pattern is different in the cold− experiments subject to a deepened thermocline, with a more heterogenous distribution of retreat date. Across all CDW scenarios, the result of increased surface irradiance is earlier sea-ice melt and thus earlier opening of the polynya at the start of summer.
Both peak and annual NPP show a sub-linear response to surface irradiance, with a stronger sensitivity to thermocline depth ( Figure 14, Table 2). In the (lowered melt rate) cold− cases, there is a single peak in NPP, as opposed to the double peak observed in base− and warm− experiments, with the cold− peak up to 20% higher than the base− peak. This elevated NPP from cold− boundary conditions is in spite of a reduction in horizontally averaged iron concentration across the near-surface (<100 m) waters of the domain. In particular, the annual productivity of the polynya increases from 55 g C m −2 yr −1 in base_med to 66 g C m −2 yr −1 in the cold_med experiment. In the warm− cases, where the thermocline is shallow compared to base− experiments, NPP is again increased by around 20%, though this time in conjunction with an increase in iron content. For all three ocean conditions, increased irradiance causes small increases in NPP and small reductions in upper ocean iron concentration (due to increased phytoplankton uptake).
We seek to explain the NPP time series by examining the impact of variations in melt rate on physical variables relevant to phytoplankton growth ( Figure 15). Most significantly, meltwater concentrations increase in the central portion of the ASP, whilst decreasing close to the coast, in cold_med as compared to base_med. Mixed layer depth responds most strongly near the coastline, with a slight shallowing visible in both the cold_med and warm_med cases, as compared to the base_med experiment. Reduced melt rate in cold_med produces a weak sea surface warming close to the shelf, probably due to the reduced inputs of freshwater; meanwhile in the warm_med case the strong warming is likely due in part to changes in upwelling.
We investigate further the differences in spatially averaged upper ocean iron and NPP in Figure 16 with maps of their respective distributions. In base_med and warm_med (Figures 16b and 16c), iron levels in October (before the bloom) are highest close to the coast and in front of the peninsulas. In contrast, the cold_med experiment (Figure 16a) shows significant levels of iron in front of the ice shelf in the center of the domain. In Figures 16d-16f, the distribution is plotted for the local iron tracer only, so as to remove the impact of upstream meltwater via the zonal periodic boundaries. In these figures the different pathway of iron from the ice shelf cavity in the cold_med case can be seen clearly. These same plots also show a gradual strengthening of the modeled gyre as the thickness of the CDW layer is decreased. The maps of annually integrated NPP in Figures 16g-16i show that a shallowing of the thermocline (and thickening of the CDW layer) in warm_med strengthens productivity in the coastal current, but does not radically change the distribution of the bloom. However the deepening of the thermocline (and thinning of the CDW layer) in cold_med leads to a qualitative change in NPP distribution as compared to base_med, with the bloom now focused in the center of the domain rather than in front of the peninsulas.

Impact of the Ice-Shelf on Productivity
Our results indicate that the high primary productivity observed in the ASP is critically dependent on the basal melting of adjacent ice shelves, with annual NPP increasing from 1 g C m −2 yr −1 to 41 g C m −2 yr −1 based on the meltwater pump alone. This is in agreement with previous investigations of the meltwater pump, in the Amundsen Sea et al., 2020). These studies found large decreases in NPP when iron concentrations at the shelf front were set to zero. Significantly, we arrive at this conclusion despite using a novel approach to isolate the effect of the meltwater pump, with melting of the ice shelf prevented in the physical model by suppressing heat fluxes. This in turn suggests that the importance of the meltwater pump in driving high NPP is not an artifact of individual model setups, but instead follows naturally from the iron and freshwater concentrations measured near DIS (Randall-Goodwin et al., 2015).
The set of experiments presented here does not allow us to distinguish between iron originating off the continental shelf in CDW and iron sourced from sediments on the continental shelf. There remain large uncertainties in specifying the end-member concentration of dissolved iron in glacial meltwater, compounded by the uncertainty surrounding the partial bio-availability of the particulate iron pool . In our experiments, we follow St-Laurent et al. (2019) in using a 20 μmol m −3 concentration for iron in meltwater, and observe a 49% increase in maximum wintertime iron inventory as a result. This in turn drives a 51% increase in iron uptake, indicating a near-linear response as expected for an iron-limited system where Fe ≪ K Fe (cf. Equation 8). The accompanying increase in NPP is smaller-only 34%-due to the variable iron-to-nitrate uptake ratio, as discussed in Section 4.3.  use of zonal periodic boundaries permits us to account for the influence of upstream ice shelves, which is significant-as expected from St-Laurent et al. (2019). A comparison of the total and local meltwater tracers demonstrates that a large proportion of the meltwater present in the domain is sourced from the eastern periodic boundary, corresponding to an upstream source. This is reflected in the relative values of local and global iron tracers. The presence of a deep meltwater maximum in Dotson Trough (Figure 4) raises the possibility that some of the iron which is melted out from the ice shelf may become re-entrained in the inflow of modified CDW into the ice shelf cavity. Thus a portion of the iron pulled up from depth by the meltwater pump may itself be glacial in origin.
Melting of the ice shelf is shown to be an important driver not only of iron distribution but also of upper ocean circulation, sea ice cover, and vertical mixing. Our sensitivity experiments, conducted with no change in wind forcing, show that melt rate dictates the strength of the coastal current which counters the wind driven gyre (see Section 4.4). This, in combination with surface irradiance, largely sets the timing of the seaice retreat in spring. Finally, the mixed layer depth shows sensitivity to melt rate, though this sensitivity is mostly confined to the region where the coastal current is strong.

Light and Vertical Mixing
The vertical structure of the phytoplankton bloom is primarily a product of light attenuation and vertical mixing in the upper ocean. By using the Manizza et al. (2005) formulation, we demonstrate that phytoplankton self-shading has a large impact on light penetration and consequently on the euphotic depth TWELVES ET AL.   Figure 7). Furthermore, there is large spatial and temporal variability in the self-shaded euphotic depth, a phenomenon which cannot be captured using a single attenuation constant throughout the model. We find that self-shading leads to a shoaling of the December mean DCM by 40 m (Figure 8b) and a reduction in annual NPP from 115 g C m −2 yr −1 to 55g C m −2 yr −1 , compared to the model run with constant euphotic depth. Self-shading also leads to a delay in bloom onset, which in turn results in enhanced late-season surface NPP for gmw_iron compared to the fized_zeu experiment (discussed in section 4.3). Moreover, these changes in productivity impact the distribution of nitrate and oxygen within the model. Changes in nitrate export due to self-shading were identified in Manizza et al. (2005) and Kim et al. (2015) as impacting on productivity in the macronutrient-limited lower latitudes. Meanwhile the reduced oxygen concentrations above 200 m and increased concentrations below illustrate how changes in light environment may impact the wider polynya ecosystem.
The distribution of biomass shows that the dominant phytoplankton type varies with depth ( Figure 8c). This is a result of the differential treatment of loss rate between small (P. antarctica) and large (diatom) phytoplankton classes in BLING; furthermore, it is likely that we underestimate the differences in community composition with depth. Biomass is mixed through the water column as a single tracer, so that variability in community composition tends to be smoothed out by the physical model.
Vertical mixing is also important in determining how well adapted phytoplankton are to their light environment. As outlined in Section 2.3 we make the assumption that the irradiance memory is a function of the average light throughout the mixed layer. This implies that near-surface phytoplankton are acclimated to efficiently utilize lower light levels than those which they are exposed to. Thus their already high photosynthesis rates are enhanced (Schofield et al., 2015). Conversely, phytoplankton at the bottom of the mixed layer are exposed to light levels lower than those to which they have acclimated, so here the already low rate of photosynthesis is suppressed further. The contrary assumption would remove this enhancement/ suppression of photosynthesis at the surface/base of the mixed layer, and thus likely give a flatter profile of productivity with depth. However, data from Argo floats in the Southern Ocean indicate that it may not be possible to define a single time-scale for mixing of bio-optical properties (Carranza et al., 2018). Instead they are mixed episodically following storms, becoming restratified in the days that follow. As a result, the validity of our assumption can vary in time and space. Meanwhile, shipboard measurements from the Amundsen Sea show that diatoms are better able to make use of variable light conditions than P. antarctica . This indicates that the degree of photoacclimation may be species-dependent as well as time-dependent over the course of the bloom.
TWELVES ET AL. In this study, we use KPP to conduct vertical mixing. This follows previous modeling in the ASP  but diverges from B-SOSE, which uses the GGL_90 mixing scheme (Gaspar et al., 1990). We find mixed layer depths of less than 20 m throughout the domain except adjacent to the ice shelf front where there is a strong meltwater outflow from below 100 m (as discussed above in Section 4.1). These depths were calculated using the second derivative of the density profile with depth, and were used as diagnostics as well as for homogenizing the irradiance memory as described above. Our mixed layer depths are shallow compared to Alderkamp et al. (2015) and Park et al. (2017), who find mixed layers of up to 70 m and 50 m, respectively, in the ASP. This may in part be due to the fact that these observations use a density threshold to calculate mixed layer depth rather than the second derivative of density. If instead the mismatch is due to an underestimate of wind-driven mixing in our model, then this may result either from our choice of mixing scheme or from our choice of wind forcing. Conversely, our model produces euphotic depths in excess of those measured by Park et al. (2017), even when we include the impact of self-shading. A more accurate light field for the ASP might be obtained by tuning coefficients in the bio-optical model to better represent local conditions, or by allowing attenuation to vary with other water constituents such as colored dissolved organic matter or detritus (Dutkiewicz et al., 2015). The combination of shallow mixed layers and deeper euphotic zones would be expected to reduce light limitation, yet our model does not overestimate productivity. Instead we find annual NPP in the range 50-150 g C m −2 yr −1 depending on location within the domain, which is in line with Yager et al. (2012) (their Figure 3) as well as most stations modeled in Oliver et al. (2019). However, depending on our boundary conditions, the pattern of NPP often differs substantially (see Section 4.4).
While attenuation due to phytoplankton can potentially impact the upper ocean heat budget (Manizza et al., 2005), this feedback is not yet implemented in BLING within MITgcm. Instead, the total shortwave irradiance used for heating and the PAR relevant for biology are treated independently in the ocean interior. This means that we allow the mixed layer depth to impact the euphotic depth via changes in chlorophyll distribution, but the euphotic depth cannot in turn feedback onto the mixed layer depth. Implementing the full biophysical feedback is expected to be particularly important in regions where the euphotic depth exceeds the mixed layer depth, as is observed throughout our model and by Alderkamp et al. (2012) in parts of the ASP. Our suite of sensitivity experiments already demonstrates that changes in shortwave heating can impact not only the sea ice cover but also the frontal melting of the ice shelf, further motivating the inclusion of this feedback in future ice-ocean modeling.

Iron-Light Colimitation
The balance between iron and light limitation in the ASP varies with location and with time. In all experiments, the phytoplankton bloom cannot initiate until the sea ice cover has largely retreated. This accords with observations, with low humidity and strong katabatic winds preventing the formation of the melt ponds which permit under-ice blooms in the Arctic (Horvat et al., 2017). However, a comparison of experiments with identical sea ice cycles in Figure 11 shows that the sea ice cover is not the only control on the timing of bloom initiation. The inclusion of self-shading delays the onset of the bloom in gmw_iron by around a month compared to fixed_zeu. This implies that after sea-ice begins to retreat, the early growth of the bloom is slowed by the negative feedback between biomass and light penetration. The peak NPP reached in gmw_iron of 0.51 g C m −2 day −1 is 80% lower than the 2.47 g C m −2 day −1 reached in fixed_zeu. Up to this point, flat_iron does not diverge significantly from gmw_iron. This demonstrates that light limitation is the primary control in the early part of the bloom leading to the initial peak.
The timing of the initial gmw_iron peak is around a month earlier than the peak generally observed in the ASP. It is followed by a second peak at the start of January, by which point flat_iron has surpassed fixed_zeu, indicating a transition from a predominantly light-limited to a predominantly iron-limited system. Nonetheless, for the entirety of December and January both fixed_zeu and flat_iron are more productive than gmw_iron, so that the bloom can be regarded as co-limited by iron and light. The severe iron limitation that occurs in gmw_iron as compared to flat_iron is a cumulative result of uptake by phytoplankton. Moreover, by December the bloom in fixed_zeu has been substantially more productive than gmw_iron and has thus exacerbated the depletion of the iron pool. Nearest the surface where light limitation is least severe, this al-lows for higher growth rates and greater productivity in the self-shaded experiment gmw_iron despite being exposed to less intense light than in fixed_zeu (Figure 8d).
The seasonal succession from light to iron limitation reflects previous modeling results in the ASP , as well as observations from the West Antarctic Peninsula (Arrigo et al., 2017) and the Weddell Sea (von Berg et al., 2020). We use the transition from light to iron limitation described above on a spatially averaged basis to examine differences across the domain. Alderkamp et al. (2015) show that there are clear differences in the characteristics of phytoplankton across the ASP, in particular with respect to distance from ice shelves. Growth rate limitation due to iron deficiency is found to be more severe in the central ASP compared to stations in the midst of the meltwater outflow. Our results in Figure 11b reproduce some of this spatial variability, with an earlier transition to iron limitation found in the center of the domain away from meltwater inputs. In front of the peninsulas where concentrations of the meltwater tracers are highest, light limitation remains the dominant control long after the bloom has peaked-in some locations for the entirety of the growth season.
An important distinction between our model setup and that used in Oliver et al. (2019) and St-Laurent et al. (2019) is in our use of an iron-to-nitrate uptake ratio which increases with iron concentration. As a result, the 51% higher annual iron uptake in gmw_iron compared to melt_pump (Figure 10d) is only associated with a 34% increase in productivity. More broadly, our iron budgets are dominated by the balance between uptake and remineralization, with smaller contributions from recycling and scavenging. When we remove the phytoplankton iron sink as in flat_iron, this is equivalent to the addition of a continuous flux compensating for phytoplankton uptake. The result is a large increase in the carrying capacity of the ASP, from 54 g C m −2 yr −1 to 116 g C m −2 yr −1 . This result is consistent with the conclusion of Alderkamp et al. (2015) that a continuous resupply of iron is necessary to increase the carrying capacity of the bloom.

Importance of Melt Rate and Sea Ice Cover
In our sensitivity experiments, the perturbation of downwelling short-wave irradiance is an idealized representation of the variation in cloud cover over the ASP and PIP as reported in Park et al. (2017), whilst the raising and lowering of the thermocline depth at the northern boundary is an attempt to represent interannual variability in the amount of CDW observed on the Amundsen Sea shelf (Dutrieux et al., 2014;Jenkins et al., 2016). We observe that both sets of forcing impact on productivity, but there are important qualitative differences in the response. The response of peak and annual NPP to perturbations in irradiance is sub-linear and symmetric with respect to the sign of the perturbation, i.e. NPP increases or decreases with the irradiance change with roughly similar magnitudes, regardless of thermocline depth. It is worth noting that this response is not solely due to a change in PAR, but represents a convolution of signals from PAR, modified sea ice retreat and advance, surface stratification due to heat input, and modified late-season melt input due to downwelling of solar-heated water.
Raising and lowering of the thermocline yields changes in peak and annual NPP which are highly nonlinear and greater in magnitude than changes driven by surface irradiance. This is in part due to the fact that the melt response to thermocline change is not symmetric, due both to the complexity of sub-ice shelf circulation and to Ekman downwelling at the ice-shelf front (discussed below; see also Kim et al. [2017]). However, this is not sufficient to explain the NPP response, as wintertime iron inventory does increase monotonically with the shoaling of the thermocline.
The response to elevated iron levels when the thermocline is raised to 450 m is a modest increase in peak and annual NPP attributable at least in part to a strengthening of observed patterns of iron concentration and sea-ice opening (Figures 13 and 16). A lowering of the thermocline to 750 m yields lower iron levels on a spatially integrated basis, and yet here too productivity is enhanced. This response is likely due to a largescale change in the mode of circulation-a transition from a strong, westward-flowing, coastally-intensified current to a circulating gyre, transporting iron to the center of the polynya (Figure 16). This change in the strength of the coastal current can be attributed at least in part to the reduced melt rate. In all experiments, the ice-shelf melt leads to a meridional density gradient at the ice-shelf front due to the rising, buoyant, melt-laden waters, and this gradient induces a boundary current. This boundary current is weakened when melt rate is lowered, allowing transport of melt water, and therefore iron, away from the coast. A similar coupling between coastal current and meltwater distribution is seen in the studies of Nakayama et al. (2014) and Kimura et al. (2017) (their Figure 13). We conclude that the weakening of the coastal current in years with reduced CDW intrusion allows iron to spread into the region of the domain where it is most limiting (Figure 17), thus driving higher productivity.
It is important to note that in the experiments shown the ferricline-i.e. the point above which zero concentration of iron was imposed at the northern boundary-is held constant, even when the thermocline was modified. This approach allows us to isolate the effect of the melt pump on productivity. Thus the warm− experiment results likely give a lower bound on the effect of increased ocean heat flux. However, an additional cold_med experiment was run with a lowered ferricline (results not shown); it showed qualitatively the same pattern of productivity as the results described above.
Another factor that we do not consider in our sensitivity experiments is the effect of the imposed wind forcing. In all our experiments there is a southeasterly wind stress forcing, zonally constant but decreasing northward. Such a wind stress pattern leads to Ekman upwelling over most of the domain (T. W. Kim et al., 2017). As such upwelling can lead to an input of cyclonic vorticity (Hughes, 2005); this may explain the gyre-like circulation which develops in the absence of a strongly buoyancy-forced boundary current. The strength of this upwelling (along with details of the background stratification) may then quantitatively affect the transition from a boundary current regime to a gyre circulation regime, and thus impact the nonlinear response of NPP to thermocline depth.
It is not currently known whether the mechanism identified above-lowered melting due to a deeper thermocline leading to increased transport of iron to more iron-limited regions-is important for productivity in coastal polynyas in the Amundsen Sea. Melt variability is not as wellknown as that of bloom strength or sea ice cover due to the need for in situ measurements (Jenkins et al., 2018). Moreover, the results of Arrigo et al. (2012) suggest that productivity is strongly determined by the mean open water area in the polynya, which is strongly modulated by wind variability, as well as other factors such as the presence of grounded icebergs . In this study we use a set of independent sensitivity experiments to study the consequences of different regimes of CDW intrusion on the continental shelf. In the future a multi-year model with interannual variability may be necessary to disentangle these signals from those arising due to wind forcing at the surface (Jacobs et al., 2012).

Implications and Recommendations
Our study confirms that seasonal changes in ocean optics due to self-shading have a large impact on the highly productive phytoplankton blooms found near Antarctic ice shelves. Furthermore we expect that extending the treatment of self-shading so that it impacts the attenuation of incoming irradiance within the physical model (as in Manizza et al. (2005)) will have a measurable impact on the heat budget within the top 50 m of the ocean. This heat budget has been identified as an important driver of high-latitude ice-ocean systems (Hahn-Woernle et al., 2020); our results indicate that upper ocean heating may affect not only the sea ice cover but also the frontal melting of ice shelves. However, as in Manizza et al. (2005) and Kim et al. (2015), we use a parametrization of self-shading informed by bio-optical data from lower latitudes, rather than data specific to the region of interest (Morel, 1988). Previously, efforts to quantify the TWELVES ET AL.
10.1029/2020JC016636 24 of 28 Figure 17. Schematic diagram to demonstrate spatial variability in iron and light limitation. In the coastal current, most of the upper ocean is primarily light-limited, most of the time. In the central polynya, iron limitation extends deeper and comes earlier in the summer. Therefore, the central polynya is more sensitive to an increase in iron concentrations. attenuation of light by chlorophyll in Antarctic polynyas have been more focused on constructing algorithms to extract NPP from satellite ocean color data (Arrigo et al., 2008). In the future, there is set to be a large increase in the availability of in-situ data from the Amundsen Sea, including from autonomous platforms (Meredith et al., 2016). Thus, there will be a significant opportunity to improve the characterization of bio-optical feedbacks in this important region of the Southern Ocean, leading to better estimates of NPP both from satellite data and from models.
In this study, we use a novel methodology to study iron limitation within a polynya. By comparing against an experiment where phytoplankton growth is not allowed to deplete iron, we demonstrate the large role that phytoplankton uptake of iron plays in the decline of the bloom. This in turn emphasizes the importance of the iron-to-nitrate uptake ratio early in the bloom development. The greater this ratio in the early (iron-replete) stages of the bloom, the more severe will be the iron limitation in the later stages of the bloom. Therefore, an accurate representation of iron uptake is needed to successfully reproduce the dynamics of phytoplankton blooms in polynyas such as the ASP. We also show the importance of advection from upstream in supplying iron to the ASP, although it remains difficult to isolate different sources (CDW, glacial, and benthic) within a single model experiment, especially when different sources may exhibit differing degrees of bio-availability. Future modeling could make use of multiple passive tracers for iron (as with our "local" iron tracer), or follow an adjoint modeling approach Song et al., 2016), to examine the sensitivity of the phytoplankton bloom to iron from different sources.
Highly productive Antarctic polynyas such as the ASP act as sinks for carbon, iron, and macronutrients. Nitrate, phosphate, and silicic acid from the ACC are supplied to continental shelf systems via onshelf incursions of CDW, where they may be taken up in phytoplankton blooms. If not used, they will remain in the water column and be transported northwards, eventually returning to lower latitudes along the thermocline (Gruber et al., 2019;Sarmiento et al., 2004) or becoming sequestered in Antarctic Bottom Water (Marinov et al., 2006). Our analysis of self-shading shows that changes in the euphotic depth lead to changes in nitrate concentrations deeper in the water column. A shallower euphotic depth (i.e. increased self-shading) leads to a smaller nitrate sink, which in turn means that more nitrate is returned to lower latitudes. We also demonstrate that the buoyancy of iron-rich meltwater leaving the ice shelves is a strong control on the magnitude and spatial distribution of the bloom. When less buoyant meltwater is able to spread into more iron-limited regions of the domain, the iron contained within it is more efficiently taken up and the polynya becomes a stronger sink for carbon and nutrients. Primary production in the lower latitudes, and the ecosystems that this production supports, are thus likely to be sensitive to the changing availability of iron and light in Southern Ocean polynyas (Henley et al., 2020;Moore et al., 2018). Recent modeling (Moorman et al., 2020) has shown that over longer timescales a strengthening of coastal currents around Antarctica may suppress ice shelf melting; our study suggests a separate mechanism by which these strong coastal currents could impact the global climate system, by suppressing the iron supply to phytoplankton. We recommend that self-shading and three-dimensional meltwater inputs should be included in future modeling of the Southern Ocean to improve estimates of macronutrient export and biological carbon uptake at regional and larger scales.

Conclusion
The ASP is one of several highly productive polynyas around West Antarctica, most of which are close to fast melting ice shelves. This study demonstrates that (1) ice shelves play a crucial role in driving high NPP, and (2) there is strong connectivity between the polynya and ice shelves further upstream, due to the presence of a strong coastal current. Using a novel methodology we find qualitatively the same behavior as in St-Laurent et al. (2017): the melting of the ice shelf increases iron availability in the polynya both due to a meltwater pump effect and due to release of iron entrained at the glacier bed. However, the strong westwards transport of iron in the coastal current  suggests that changes to ice shelves immediately to the east of DIS may have as large an impact on biology in the ASP as changes to the DIS itself.
Our results show that both iron and light availability impact productivity in the ASP, and that it is their combined effect which controls annual NPP. However, we find strong spatial and temporal variability in iron-light colimitation. As in Oliver et al. (2019), we find that the phytoplankton bloom is primarily light-limited in its early stages, and we are able to attribute a large part of this light limitation to the self-shading feedback which arises due to chlorophyll in the water column. The inclusion of the chlorophyll dependent light penetration formula from Manizza et al. (2005) reduces peak NPP by 80%. Later in the season there is a transition toward iron limitation in most of the domain, as the pool of dissolved iron is depleted by phytoplankton uptake. Furthermore, we find that the central portion of the domain is the first to become strongly iron-limited, with light limitation remaining the dominant control on NPP at progressively later dates with increasing proximity to the ice shelf.
In this study we find that ice shelf melt rate dictates meltwater dispersal, and that this in turn exerts the dominant control on the spatial distribution of the phytoplankton bloom. We conclude that the future viability of Antarctic polynyas as biological carbon sinks may be subject to a trade-off between increased iron leaving the ice shelf cavity and stronger coastal currents preventing this iron from reaching the phytoplankton communities where it is most limiting.

Data Availability Statement
MITgcm code can be accessed publicly at mitgcm.org; the modified code and the input folders used in our model setup are available at the open-access repository 10.5281/zenodo.4268847; data used to plot figures in this manuscript are also available at the open-access repository 10.5281/zenodo.4268839, as well as being provided as supporting information. B-SOSE outputs are available at http://sose.ucsd.edu/bsose_solu-tion_Iter105.html.