Decline in Atlantic Primary Production Accelerated by Greenland Ice Sheet Melt

Projections of climate impacts on marine net primary production (NPP) are reliant on Earth System Models (ESMs) that do not contain dynamic ice sheets. We assess the impact of potential Greenland ice sheet meltwater on projections of 21st century NPP using idealized ESM simulations. Under an extreme melt scenario, corresponding to 21st century sea level rise close to 2 m, Greenland meltwater amplified the decline in global NPP from a decrease of 3.2 PgC/yr to a decrease of 4.5 PgC/yr, relative to present. This additional reduction in NPP predominately occurs in the North Atlantic subtropical and subpolar gyres, as well as Atlantic eastern boundary upwelling systems. Accelerated NPP declines are the result of both surface freshening and reductions in upwelling‐favorable winds enhancing phytoplankton nutrient limitation. Our findings indicate that including a dynamic Greenland ice sheet in ESMs could have large impacts on projections of future ocean circulation and biogeochemistry.


Introduction
Net primary production (NPP) in the oceans is the base of the marine food web and the primary source of energy to higher trophic levels (Chassot et al., 2010;Pauly & Christensen, 1995). NPP and the subsequent export of organic material to the deep ocean also influence atmospheric CO 2 concentrations and therefore represent an important climate feedback (Matear & Hirst, 1999). Relative changes in NPP and phytoplankton biomass are typically amplified in zooplankton (Chust et al., 2014;Stock et al., 2014), especially in the oligotrophic oceans . Accurately projecting future NPP is therefore critical to characterizing climate impacts on the marine food web and ultimately fisheries (Barange et al., 2014;Blanchard et al., 2012;Cheung et al., 2010), as well as constraining climate-carbon cycle feedbacks (Friedlingstein et al., 2006). ESM projections of how ocean primary production will respond to future climate change have a high degree of associated uncertainty (Frölicher et al., 2016;Taucher & Oschlies, 2011). Nonetheless, across large ensembles of ESMs, there is general agreement that global ocean NPP will decline under climate change (Bopp et al., 2013;Cabré et al., 2014;Steinacher et al., 2010). Under high emissions scenarios, 21st century global NPP is projected to decline by up to 20% (Bopp et al., 2013) with declines exacerbated on multicentennial timescales (Moore et al., 2018). Such global declines are typically driven by surface ocean warming in low and midlatitudes increasingly stratifying the upper ocean and reducing nutrient concentrations in the euphotic zone of nutrient-limited waters (Bopp et al., 2013;Laufkötter et al., 2015;Steinacher et al., 2010). These declines are generally somewhat offset by shoaling of the mixed layer depth and sea ice loss reducing light limitation in high latitudes (Bopp et al., 2013;Steinacher et al., 2010).
Across ESM ensembles similar processes appear to determine the sensitivities of NPP to climate change and El Niño/Southern Oscillation climate variability. Consequently satellite observations of the sensitivity of NPP to El Niño/Southern Oscillation have been used to constrain projections of low latitude NPP decline to 3 ± 1% per degree warming in equatorial sea surface temperature (SST) (Kwiatkowski et al., 2017). Critically however, projections of future NPP, including such emergent constraints, are reliant on ESMs that do not contain coupled ice sheet models for Greenland and Antarctica.
Ice sheet meltwater has the potential to influence ocean biogeochemistry by affecting mixed layer depth and ocean circulation (Vizcaíno et al., 2008) as well as through the provision of associated nutrients (Bergeron & Tremblay, 2014;Death et al., 2014). The current absence of dynamic ice sheets in ESMs is primarily a consequence of a mismatch between the high spatial resolution required by ice sheet models and that which is computationally permissible in atmospheric and ocean components of current ESMs (Lipscomb et al., 2013). It is further complicated by difficulties in ice sheet initialization under past climatic forcing due to scarcity of observations, in addition to the community focus on century-scale climate projections, which has resulted in the inclusion of dynamic ice sheets not being prioritized. Here we use simulations of an ESM to explore the impact of idealized flash Greenland melting events on projections of 21st century primary production.

The ESM and Biogeochemical Model Configuration
Numerical simulations were conducted with the Institute Pierre Simon Laplace low resolution Earth System Model (IPSL-CM5A-LR; Dufresne et al., 2013) which contributed to the Coupled Model Intercomparison Project Phase 5. The ocean component has an ORCA tripolar grid with 2°horizontal resolution and 31 vertical levels (ORCA2; Madec, 2015) and is coupled to the Pelagic Interactions Scheme for Carbon and Ecosystem Studies (PISCES) ocean biogeochemical model . PISCES describes the biogeochemical cycles of carbon and the main ocean nutrients (N, P, Fe, and Si). The model contains 24 prognostic tracers including 2 phytoplankton types (nanophytoplankton and diatoms) and 2 zooplankton size classes (microzooplankton and mesozooplankton).
In the PISCES biogeochemical model, the response of NPP to internal climate variability and long-term climate change is predominately a consequence of changes in phytoplankton growth rates (μ) (Kwiatkowski et al., 2017). In the version of the PISCES biogeochemical model that was coupled within IPSL-CM5A-LR, μ is the product of a maximum phytoplankton growth rate (μ max ), temperature limitation (T f ), nutrient limitation (N lim ), and light limitation (L lim ): where T f is based on Eppley (1972), N lim is defined as the most limiting nutrient limitation term for nanophytoplankton and diatoms, respectively: and L lim is parameterized, based on Geider et al. (1997), as where i is the phytoplankton group, α is the initial slope of the photosynthesis-irradiance curve, θ Chl is the chlorophyll to carbon ratio, I PAR is photosynthetically active radiation, and f(mxl) represents additional limitation that occurs when the mixed layer depth exceeds the depth of the euphotic zone.

Greenland Ice Melt Simulations
The IPSL-CM5A-LR model was run under historical  and Representative Concentration Pathway 8.5 (RCP8.5; 2006-2100) radiative forcing, with RCP8.5 representing a business-as-usual 21st century emissions scenario (Riahi et al., 2011). Four additional freshwater hosing simulations, analogous to large and rapid Greenland ice melt, were conducted under identical 2006-2100 climate forcing. Freshwater flow of 0.05 Sv (GrIS50), 0.075 Sv (GrIS75), 0.1 Sv (GrIS100), and 0.2 Sv (GrIS200; 1 Sv = 10 6 m 3 /s) was added continuously from 2006 to 2100, corresponding to additional 21st century sea level rise of approximately 43, 65, 86, and 173 cm, respectively. For context, the largest of these simulations represents the loss of approximately a quarter of the Greenland ice shelf. Model freshwater inputs were released in surface waters of regions of North Atlantic deepwater formation (45 to 65°N, 45°W to 5°E), coincident with observations of Greenland meltwater (Bamber et al., 2012;Kjeldsen et al., 2015). The idealized meltwater was assumed to lack nutrients and associated sediments. Similar simulations of freshwater hosing, conducted with the IPSL-CM5A-LR model, have previously shown that Greenland ice melt can reduce West African monsoon rainfall, increasing the vulnerability of the Sahelian agroecosystem (Defrance et al., 2017).

NPP Response to Greenland Ice Melt
The IPSL-CM5A-LR model projects global NPP declines of 3.2 PgC/yr (9%) over the 21st century under RCP8.5, with declines of 0.7 PgC/yr (17%) in the North Atlantic and 0.2 PgC/yr (4%) in the South Atlantic ( Figure 1). Under the most extreme Greenland ice sheet melt scenario (GrIS200), projected global NPP declines increase by 1.3 PgC/yr over the 21st century, relative to RCP8.5. The overwhelming majority of these additional NPP reductions is confined to the Atlantic basin, with declines enhanced by 0.8 PgC/yr in the North Atlantic and 0.5 PgC/yr in the South Atlantic. The magnitude of additional NPP decline increases with the intensity of ice sheet melt scenario (Figure 1). At the end of the 21st century, there are additional NPP declines of 0.3, 0.4, and 0.9 PgC/yr for the GrIS50, GrIS75, and GrIS100 scenarios, relative to RCP8.5, respectively. Although projected NPP is influenced by climatic variability, for the majority of GrIS simulations, the suppression of NPP increases throughout the duration of simulations. However, under the most extreme scenario (GrIS200), suppression of NPP in the North Atlantic plateaus around 2060 ( Figure 1).
Projected global declines in NPP are 2.0 PgC/yr (7%) for nanophytoplankton and 1.2 PgC/yr (19%) for diatoms under RCP8.5. The greater relative decline in diatom NPP is attributed to the higher nutrient half-saturation constants of diatoms and consequently their greater sensitivity to nutrient limitation (Bopp et al., 2005). Under the GrIS200 simulation, Greenland meltwater enhances the NPP declines of nanophytoplankton and diatoms by 1.1 and 0.2 PgC/yr, respectively.

Geophysical Research Letters
Atlantic basin. However, the greatest enhancements of NPP suppression occur in the North Atlantic subtropical gyre and African Eastern boundary upwelling regions, where NPP under GrIS200 can be reduced by 30-70 gC·m -2 ·yr -1 relative to RCP8.5 (Figure 2).
In the first 15 years of the GrIS simulations, the enhanced suppression of NPP predominately occurs near Greenland and in the North Atlantic subpolar gyre, while there is a small region of NPP increase in the Canary Current System ( Figure S1). From 2020 onward, enhanced suppression of NPP increasingly occurs in the North Atlantic subtropical gyre and the Atlantic eastern boundary upwelling regions of both the Northern and Southern Hemispheres ( Figure S1). The Canary Current region that initially experiences enhanced NPP increasingly exhibits NPP suppression, while small regions of NPP suppression and enhancement are present near the Amazon estuary. This temporal progression and spatial pattern of enhanced NPP declines relative to RCP8.5 are consistent across the four GrIS simulations ( Figure S2).

Phytoplankton Growth Rate Limitation
Given that L lim in the PISCES model is not independent of other phytoplankton growth rate limitation terms (equation (5)), we confine our analysis here to the impact on T f and N lim . Further analysis of the response of L lim is provided in the supporting material. The influence of Greenland meltwaters on nanophytoplankton T f and N lim is shown in Figure 3. T f and N lim are depth integrated and weighted by the vertical concentration of mean historical (1986-2005) nanophytoplankton biomass. As such, the T f and N lim anomalies reflect changes in growth limitation terms that impact the vertically integrated NPP of each grid cell. Greenland ice melt is shown to have limited impact on the temperature limitation of nanophytoplankton growth rates.
In the North Atlantic subpolar gyre, T f increases slightly (<15%) representing temperature increases that act to enhance phytoplankton growth rates, whereas in the rest of the North Atlantic, T f declines slightly (<15%).
In the Southern Hemisphere, T f generally exhibits increases of a similar order of magnitude.
The impact of Greenland ice melt on nanophytoplankton nutrient limitation is considerably greater than its impact on temperature limitation. N lim decreases substantially (>30%) in the North Atlantic subtropical gyre and the Atlantic eastern boundary upwelling regions of the Northern and Southern Hemispheres. This represents an increase in nutrient limitation that acts to suppress nanophytoplankton growth rates in these regions (Figure 3b). In contrast with changes in temperature limitation, the regions of greater nutrient limitation are strongly coincident with the areas of enhanced NPP decline. As such, the dominant driver of enhanced NPP reductions in the GrIS simulations is greater nutrient limitation, with minor additional influence due to changes in temperature and light limitation ( Figure S3). The spatial pattern of anomalies in diatom growth limitation terms is broadly similar to that of nanophytoplankton ( Figure S4). The magnitude of anomalies differs however, due to the growth parameterization and biomass distribution of each plankton functional type.

Stratification and Upwelling Intensity
The reduction in phytoplankton growth rates and therefore NPP, in response to enhanced nutrient limitation in the GrIS simulations, can be attributed to changes in both stratification and the intensity of eastern boundary upwelling. The impact of increased stratification, as a consequence of surface freshening, is evident in shoaling of MLD max . The MLD max in the GrIS200 simulation shoals by up to 20 m, relative to RCP8.5, in the North Atlantic subpolar and subtropical gyres by the end of the 21st century ( Figure 4a). This results in less nutrient-rich deep waters being mixed into the euphotic zone and, in areas of typically high nutrient limitation such as the subtropical gyre, explains the enhanced NPP declines.
The enhanced nutrient limitation and subsequent declines in NPP that occur in the Atlantic eastern boundary upwelling systems are the result of changes in both stratification and upwelling. In the Benguela Current region, there is strong coherence between declines in upwelling-favorable winds (τ upw ) and declines in the vertical supply of NO 3 in the GrIS200 simulation ( Figure 4). The reduction in NO 3 supply, which enhances nutrient limitation and drives the observed declines in NPP in this region, can be explained by reduced vertical advection (upwelling), alongside reductions in the vertical mixing of nutrients into the upper ocean. In the Canary Current region, τ upw generally increases in the GrIS200 simulation relative to RCP8.5 (Figure 4b), with vertical advection predominately increasing ( Figure 4d). However, south of the Canary Islands the vertical supply of NO 3 , which is driven by both advection and mixing, declines. This is the location of greater nutrient limitation and declining NPP. North of the Canary Islands, the vertical supply of NO 3 shows limited increases, and there is reduced nutrient limitation; however coincident increases in temperature limitation offset any potential impact on NPP.

Mechanisms of Enhanced NPP Decline Under Greenland Ice Melt
It is somewhat expected that the simulation of Greenland ice melt results in enhanced regional declines in NPP. Greater phytoplankton nutrient limitation, in response to surface warming that intensifies stratification, is typically the dominant driver of projected 21st century NPP declines, across ESMs (Bopp et al., 2013;Steinacher et al., 2010). Our simulations of Greenland ice melt freshen local surface and upper ocean waters, thereby increasing the gradient of vertical density profiles and further enhancing thermally driven increases in stratification. This shoals the MLD max and reduces the mixing of nutrient-rich deep waters into the upper ocean.
It is, however, slightly counterintuitive that such a large proportion (38%) of enhanced Atlantic basin NPP declines occur in the Southern Hemisphere and in particular in the eastern boundary upwelling regions.  anomalies, which reduce rainfall associated with the West African monsoon. Here we show that Greenland ice melt also acts to reduce τ upw in the Benguela Current region, decreasing the vertical supply of nutrientrich deep waters and in turn reducing NPP both in the coastal region and more generally across the South Atlantic. The reduction in τ upw in the Benguela Current region is due to the development of low-level positive temperature anomalies over the South Atlantic and negative temperature anomalies over the adjacent African continental landmass ( Figure S5). This acts to reduce land-ocean sea level pressure gradients and consequently decrease upwelling-favorable winds.
It should be noted that climate change has often been hypothesized to intensify eastern boundary upwelling, due to greater continental warming and intensification of land-ocean atmospheric pressure gradients (Bakun, 1990). Recent ESMs, however, generally project an equatorward reduction in τ upw and poleward intensification of τ upw associated with the poleward extension of high-pressure atmospheric cells (Rykaczewski et al., 2015). In the RCP8.5 simulation, the IPSL-CM5A-LR model exhibits a similar poleward displacement of τ upw in the Atlantic ( Figure S6). However, the additional simulation of Greenland ice melt modifies this, particularly in the Southern Hemisphere where the poleward displacement of τ upw is intensified (Figure 4).

Caveats and Limitations
The Greenland ice sheet meltwater simulations presented herein are, by design, highly idealized. The ocean circulation response in IPSL-CM5A-LR is likely to differ in higher resolution models that resolve Greenland boundary currents, which might otherwise divert meltwaters away from deep convection regions (Weijer et al., 2012). Interactions with Antarctic Ice Sheet loss, which has been projected to enhance Southern Ocean primary production (Death et al., 2014;Person et al., 2019), are not explored. Moreover, the magnitude of simulated freshwater discharge is notably high compared to both recent observations of the Greenland contribution to global sea level rise (0.77 mm/yr; Bamber et al., 2018) and process-based estimates of the contribution to sea level rise under RCP8.5 (6-33 cm in 2100; Vizcaino et al., 2015;Fürst et al., 2015;Aschwanden et al., 2019). The duration of simulations is also relatively short compared to the likely multicentennial timescales of future ice loss. Given limited resources, our approach maximizes signal-to-noise ratios and reveals ocean-climate dynamics that might not be apparent under lesser freshwater forcing, without computationally prohibitive large-ensemble or long-duration simulations. Nonetheless, although unlikely, even the most extreme meltwater simulation is arguably plausible given paleorecords of meltwater pulse events that resulted in >40 mm·yr -1 sea level rise (Deschamps et al., 2012).
In addition, given the uncertainty associated with ESM ensemble NPP projections (Bopp et al., 2013), the impact of Greenland ice melt on projected NPP is likely to exhibit some degree of ESM dependence. Indeed, diverse ocean biogeochemical models coupled to the same physical ocean model can project a wide variety of NPP responses to identical physical ocean forcing (Kwiatkowski et al., 2014). As such, our results should be validated with alternative ESM and ocean biogeochemical model configurations.

Glacial Nutrients and Subglacial Discharge
Perhaps the greatest limitation of our model results is that we assume freshwater fluxes do not contain nutrients. There is emerging evidence that Greenland meltwater is a potential source of bioavailable iron (Bhatia et al., 2013;Hawkings et al., 2014) and N-rich dissolved organic matter (Lawson et al., 2014) to the North Atlantic. Furthermore, models of subglacial Greenland discharge have shown that the depth of glacier grounding lines can result in the entrainment of limiting nutrients into discharge plumes (Cape et al., 2019;Hopwood et al., 2018). Our idealized surface hosing simulations do not capture interactions between grounding line depth, subsurface discharge, and plume nutrient entrainment, which are likely to influence projections of localized primary production. Although we currently lack the data to constrain estimates of these nutrient fluxes, they could potentially compensate for stratification-driven reductions in euphotic zone nutrients. Indeed, simulations of Antarctic Ice Sheet meltwaters rich in iron have shown enhanced NPP in the Southern Ocean (Death et al., 2014;Person et al., 2019). The importance of iron for NPP in the North Atlantic is considerably less than in the Southern Ocean however, where iron is the main limiting nutrient. As such, the extent of Southern Ocean iron fertilization seen in model simulations is unlikely to be replicated in the North Atlantic with comparable Greenland simulations. Moreover, a fertilization effect due to Greenland meltwaters is likely to be highest locally and has limited impact on projected NPP declines in the subtropical gyre. It should almost certainly not influence the NPP declines in upwelling regions, which are driven by ocean-atmosphere teleconnections. Nonetheless, further sensitivity studies should be conducted to assess the aggregate impacts on basin-scale NPP.

Phytoplankton Stoichiometry
A further caveat of this study is the use of an ocean biogeochemical model (PISCES) with fixed C:N:P phytoplankton stoichiometry. Such Monod formulations are computationally efficient and therefore used as standard in ESMs, despite observations that phytoplankton exhibits varying degrees of stoichiometric plasticity (Moreno & Martiny, 2018). NPP projections in ocean biogeochemical models that incorporate variable C:N:P phytoplankton stoichiometry can be less sensitive to enhanced nutrient limitation under climate change (Kwiatkowski, Aumont, Bopp, & Ciais, 2018), while reductions in zooplankton growth efficiency can enhance trophic amplification of biomass declines in such models . Future work should assess the extent to which variable phytoplankton stoichiometry may moderate the impact of Greenland ice melt on primary production.

Conclusions
Using a coupled ESM, this study explores the potential impact of rapid Greenland ice sheet melt on 21st century projections of NPP, under a business-as-usual emissions scenario. The simulation of surface ocean freshening, in Atlantic waters of Greenland meltwater discharge, results in additional suppression of projected 21st century NPP. This enhanced NPP suppression predominately occurs in the North Atlantic subtropical and subpolar gyres, as well as Atlantic eastern boundary upwelling systems, and is shown to be largely due to greater phytoplankton nutrient limitation. Projected increases in phytoplankton nutrient limitation are the result of both enhanced stratification in the North Atlantic and atmospheric teleconnections that reduce upwelling-favorable winds, particularly in the Benguela Current region. Although highly idealized, our simulations indicate that dynamic ice sheets could have dramatic impacts on projections of ocean circulation and biogeochemistry. As such, their incorporation in ESMs should be considered a future priority.