Regional-scale phytoplankton dynamics and their association with glacier meltwater runoff in Svalbard

. Arctic ampliﬁcation of global warming has accelerated mass loss of Arctic land ice over the past decades and lead (cid:58)(cid:58) led (cid:58) to increased freshwater discharge into glacier fjords and adjacent seas. Glacier freshwater discharge is typically associated with high sediment loads (cid:58)(cid:58)(cid:58) load (cid:58) which limits the euphotic depth, but may also (cid:58)(cid:58) aid (cid:58)(cid:58) to (cid:58) provide surface waters with essential nutrients, thus having counter-acting (cid:58)(cid:58)(cid:58)(cid:58)(cid:58)(cid:58)(cid:58)(cid:58)(cid:58)(cid:58)(cid:58) counteracting (cid:58) effects on marine productivity. In-situ observations from a few measured fjords across the Arctic indicate that glacier fjords dominated by marine-terminating glaciers are typically more productive 5 than those with only land-terminating glaciers. Here we combine chlorophyll a from satellite ocean colour, an indicator of phytoplankton biomass, with glacier meltwater runoff from climatic mass-balance modelling to establish a statistical model of summertime-phytoplankton (cid:58)(cid:58)(cid:58)(cid:58)(cid:58)(cid:58)(cid:58)(cid:58)(cid:58)(cid:58)


Introduction
The Arctic cryosphere is experiencing rapid transitions due to Arctic amplification of global warming. Climate change is reflected in changing oceanic and atmospheric circulation patterns, permafrost degradation, decline in sea ice thickness and extent, as well as shrinking glaciers (AMAP, 2017, IPCC, 2019. Over the past few decades, glaciers and ice caps in the 25 Arctic have retreated and lost mass at accelerating rates (e.g., Hugonnet et al., 2021), including glaciers in Svalbard (Schuler et al., 2020). A long-term trend of increased mass loss is also observed for the Greenland ice sheet, despite of a temporary slowdown of mass loss in -2017(IMBIE Team, 2019. Ice mass loss in form of glacial meltwater runoff or frontal ablation, i.e. iceberg calving and submarine melt, constitutes a significant source of freshwater being discharged into glacial fjords and adjacent seas (Bamber et al., 2018). This glacier freshwater discharge has implications for the physical oceanographic 30 conditions (Straneo and Cenedese, 2015;Carroll et al., 2017) and the biogeochemistry of water masses (Wadham et al., 2013;Hopwood et al., 2016), which affects the biological productivity in the fjords and the ocean (Etherington et al., 2007;Juul-Pedersen et al., 2015;Arendt et al., 2016;Meire et al., 2016Meire et al., , 2017Calleja et al., 2017;Kanna et al., 2018;Hegseth et al., 2019;Hopwood et al., 2018Hopwood et al., , 2020. Arctic marine ecosystems display strong seasonal cycles in productivity and functioning due to pronounced seasonality of 35 environmental variables such as solar radiation, sea ice concentration, sea-surface temperature and salinity, as well as terrestial freshwater input (Sakshaug, 2004). Marine primary production, i.e. the generation of phytoplankton biomass, ultimately depends on the availability of light and the supply of essential, 'limiting' nutrients (Sakshaug, 2004;Hopwood et al., 2020).
Seasonal changes in any of these factors lead to periods of high or low primary production (Sakshaug, 2004;Arrigo and van Dijken, 2015). A characteristic 'phytoplankton spring bloom' follows the rapid increase in incoming solar radiation after the 40 polar night, combined with high initial nutrient levels and the development of a weak stratification (e.g., Sakshaug, 2004;Juul-Pedersen et al., 2015;Meire et al., 2016). The persistence of sea ice, with or without snow cover, may delay the penetration of light into the water column and thus the phytoplankton spring bloom (Sakshaug, 2004;Rysgaard and Nielsen, 2006). Stratification relies on a positive gradient in potential water density with depth, which is controlled by salinity and temperature. Stratification during spring bloom is due to freshwater input, mainly from melting of sea ice, as well as solar heating.
Subglacial discharge drives buoyant upwelling of plumes near the calving front of tidewater glaciers, which lead to entrainment of large volumes of ambient seawater from all depth levels, thereby supplying nutrient-depleted surface layers with nutrients from nutrient-rich deep water layers (Meire et al., 2017;Kanna et al., 2018;Hopwood et al., 2018Hopwood et al., , 2020. A study by Hopwood et al. (2018) suggests that this 'nutrient pump' may provide the euphotic zone with two orders of magnitudes more nutrients than what is directly supplied by the glacial meltwater. Glacier runoff may also enhance the general estuarine circulation within 75 fjords and embayments, which is considered to have positive effects on biological productivity (Rysgaard et al., 2003;Juul-Pedersen et al., 2015;Meire et al., 2016). Down-fjord katabatic winds facilitate export of brackish/low-density surface water out of the fjord and lead to a compensating return flow of nutrient-rich saline water at depth (Svendsen et al., 2002;Cottier et al., 2010;Straneo and Cenedese, 2015;Spall et al., 2017;Sundfjord et al., 2017). In either case, positive effects of glacier runoff on primary productivity are expected to occur only where suspended particles have settled deeper into the water column 80 and light conditions in surface waters become more favorable (Etherington et al., 2007;Lydersen et al., 2014;Halbach et al., 2019).
In Godthåbsfjord, a sub-Arctic tidewater glacier fjord in SW Greenland, Juul-Pedersen et al. (2015) observed a secondary 85 peak in primary production, or 'summer bloom' that coincided with substantial runoff from the Greenland Ice Sheet. This summer bloom may be of similar magnitude, or even exceed the spring bloom. Similar findings are available from Glacier Bay, Alaska (Etherington et al., 2007). On Svalbard, glacier runoff is known to affect the distribution and species composition of phytoplankton (Piquet et al., 2014;van de Poll et al., 2018), but it is a matter of debate whether or not glacier runoff facilitates higher productivity during summer (Halbach et al., 2019). 90 capturing a snapshot during summertime. This highlights the need for innovative long-term monitoring programs of proglacial marine-ecosystems (Straneo et al., 2019). In addition, efforts should be taken to up-scale local in-situ observations in space and time. This can be achieved by the application of modelling approaches and/or satellite remote sensing.
This study aims to investigate the overall effects of glacier runoff on phytoplankton dynamics and marine primary productivity in Svalbard, focusing on a regional, rather than local scale. We utilize a 10-year timeseries of glacier runoff from 100 high-resolution climatic mass balance simulations of all glaciers in Svalbard for the timeperiod 2003-2013 (Aas et al., 2016) and chlorophyll a concentrations from satellite ocean-colour, an indicator of phytoplankton biomass (Moses et al., 2009;Matrai et al., 2013;Kahru et al., 2014;Lee et al., 2015). Chlorophyll a products and other physical ocean variables, including sea surface temperature (SST) and sea ice fraction (SIF), are available through the Copernicus Marine Environment Monitoring Service (CMEMS). We use a statistical model to identify significant associations of chlorophyll a with runoff, while accounting 105 for the potentially confounding effects of physical ocean and sea ice variables that may covary with runoff. We focus on the summer melt period, from mid-June to September, anticipating that this period follows the termination of the spring bloom.
Specifically, we investigate whether there are significant associations between runoff and chlorophyll a in coastal waters around Svalbard, and if there are spatial variations in association strength, e.g. with respect to regional characteristics or distance to coast.

Research region
The Svalbard archipelago in the Eurasian Arctic is bordered by the Barents Sea to the east, the Greenland Sea to the west and the Arctic Ocean to the north (Fig. 1). The climate in Svalbard is relatively warm, given its high Arctic location. This is due to the West Spitsbergen Current (WSC), an extension of the North Atlantic Current, which transports warm Atlantic water up north along the West Spitsbergen Shelf (Svendsen et al., 2002;Walczowski and Piechura, 2011, Fig. 1a). The eastern 115 side of Svalbard is dominated by the East Spitsbergen Current (ESC), which transports cold Arctic water clockwise around the southern tip of Spitsbergen (Loeng and H., 1991;Svendsen et al., 2002). It continues northwards on the West Spitsbergen Shelf, forming a coastal current, which is subsequently freshened by the export of brackish surface water from the fjords (Svendsen et al., 2002;Nilsen et al., 2016, Fig. 1a).
From 1971 to 2017, Svalbard has experienced strong atmospheric warming by 3-5 • C (Hanssen-Bauer et al., 2019), evident 120 in all seasons, but most pronounced during winter and spring (Nordli et al., 2014). Strong atmospheric warming is attributed to a general decline in sea ice and an increase in sea-surface temperatures (Isaksen et al., 2016). Climate projections under medium to high emission scenarios indicate that air temperatures may rise by 7-10 • C until 2071-2100, as compared to 1971-2000, which may lead to a five-fold increase in glacier mass loss (Hanssen-Bauer et al., 2019). 57% (34000 km 2 ) of the total land area on Svalbard is covered by glaciers and ice caps, and 68% of the glacierized area 125 drains into the ocean through tidewater glaciers with a total calving-front length of ∼740 km (Nuth et al., 2013). The degree of glacier coverage and the size of individual glaciers reflects the general climatic gradient across Svalbard. Glaciers in the southern and western parts, characterized by relatively warm atmospheric and oceanic conditions, are generally smaller than glaciers in the north-eastern parts of Svalbard, where colder climatic conditions prevail. Consequently, the total glacier coverage is lower in the southern and western parts, with a minimum in the dry central parts of Spitsbergen (Nuth et al., 2013). Overall, 130 glacier in Svalbard have been loosing mass since the 1960's, with a pronounced increase in mass loss since the 2000's (Schuler et al., 2020). A compilation of available mass balance assessments for the period 2000-2019 reveals a total mass balance of -8 ± 6 Gt a −1 , of which -7 ± 4 Gt a −1 are attributed to the climatic mass balance and -2 ± 7 Gt a −1 to the poorly constrained frontal ablation, i.e. iceberg calving and submarine melt (Schuler et al., 2020). The climatic mass balance simulation by Aas et al. (2016), from which we extract glacier runoff, is included in this reconciled mass balance estimate. For the period 2003-135 2013, Aas et al. (2016) found a mean annual mass balance of about -8.7 Gt, which is well within the error margins of the consensus estimate by Schuler et al. (2020).
Fjords in Svalbard are affected by terrestrial freshwater discharge, on one hand, and the exchange of water masses with the adjacent shelf, on the other hand (Svendsen et al., 2002;Cottier et al., 2005;Nilsen et al., 2016;Sundfjord et al., 2017). Glacier ablation constitutes the major component of the terrestrial freshwater discharge into Svalbard fjords (Pramanik et al., 2018;van 140 Pelt et al., 2019). During the summer melt season, glacier runoff enters the fjord in the form of surface runoff and subglacial discharge, in addition to iceberg calving and submarine melt. This freshwater mixes with ambient fjord water to form a layer of brackish surface waters, its thickness typically decreasing from the head towards the mouth of the fjord (Svendsen et al., 2002).
The exchange of water masses between the fjords and shelf depends on stratification and wind-stress, as well as the presence or absence of a topographic barrier, e.g. in form of a shallow sill at the fjord mouth (Cottier et al., 2010). The dominating wind 145 field in Svalbard fjords is down-fjord, due to katabatic winds and orographic steering of the large-scale wind-field (Svendsen et al., 2002;Cottier et al., 2005). This drives brackish surface water out of the fjord and a compensating inflow of Atlantic water from the shelf, thereby stimulating estuarine circulation and vertical mixing of water masses (Svendsen et al., 2002;Cottier et al., 2010;Sundfjord et al., 2017). Changes in atmospheric circulation patterns since the early 2000's have caused repeated overflow of the WSC onto the West Spitsbergen Shelf and inflow of warm saline Atlantic water masses into some 150 of the major fjords, with implications for water mass composition and heat content, significantly reducing sea ice production during wintertime (Cottier et al., 2007;Nilsen et al., 2016). Svalbard fjords can be regarded as broad fjords, i.e.fjord circulation is influenced by rotational dynamics or 'Coriolis' effects (Svendsen et al., 2002;Cottier et al., 2010). For our regional-scale assessment of glacier-runoff effects on phytoplankton dynamics and marine primary production, we consider 14 primary drainage basins or hydrological regions of Svalbard (Fig. 1a), following the most recent Svalbard glacier 155 inventory (Nuth et al., 2013;König et al., 2014). The identification system follows Hagen et al. (1993), where the first digit represents one out of five major areas: (1) Spitsbergen, (2) Nordaustlandet, (3) Barentsøya, (4) Edgeøya, (5) Kvitøya, the 5 https://doi.org/10.5194/bg-2021-181 Preprint. Discussion started: 20 July 2021 c Author(s) 2021. CC BY 4.0 License. latter of which is not included in this study. The second and third digits indicate the primary and secondary drainage basins, respectively. For each hydrological region, we distinguish between different marine zones, defined by their distances from the coast, namely 0 to 10 km, 10 to 20 km and 20 to 50 km. The innermost zone contains most of the fjords, which typically have 160 a width of less than 20 km. The outer regions beyond 10-km distance from the coast extend into the open ocean. Along the western and northern side of Spitsbergen, the 50 km offshore-distance contourline corresponds approximately with the shelf edge. In addition to the primary hydrological regions, we consider one subregion near the research hub of Ny Ålesund in NE Spitsbergen (15). The Kongsfjorden-Krossfjorden system consists of two secondary drainage basins, Kongsfjorden (155), and Krossfjorden (156) and serves as a key site for interdisciplinary studies on glacier-ocean interactions, focusing on physical 165 oceanographic conditions in response to glacier runoff (Svendsen et al., 2002;Cottier et al., 2005;Sundfjord et al., 2017;Torsvik et al., 2019) and their implications for the marine ecosystem (Lydersen et al., 2014;Piquet et al., 2014;Calleja et al., 2017;Halbach et al., 2019;Hegseth et al., 2019).

Climatic glacier mass balance and meltwater runoff
170 We extract regional glacier meltwater runoff from a 10-year simulation of the climatic mass balance of all glaciers in Svalbard, later referred to as glacier runoff or simply runoff. The coupled atmosphere-glacier model was run over the time period September 2003 to September 2013 (Aas et al., 2016). The glacier model computes the climatic mass balance (CMB), i.e. the mass fluxes at the surface of the glacier, mainly due to deposition of snow during the accumulation season (typically October to May) and surface melt followed by runoff during the ablation season (typically June to September). The CMB model is implemented into the Weather Research and Forecasting model (WRF), which provides precipitation and other meteorological variables to the CMB model, required to compute the climatic mass balance, considering the surface energy balance. WRF is a mesoscale atmospheric model (Skamarock and Klemp, 2008). In Svalbard it has been applied to study boundary layer processes (Kilpelainen et al., 2011(Kilpelainen et al., , 2012 and atmosphere-land interactions over both tundra (Aas et al., 2015) and glaciers (Claremar et al., 2012;Aas et al., 2016). Coupled model simulations were run over all of Svalbard at 3-km horizontal resolution 180 using sea-surface temperature and sea ice concentration from the Operational Sea Surface Temperature and Sea Ice Analysis (OSTIA) and ERA-Interim climate reanalysis data as boundary conditions. Results were validated against field observations of meteorological conditions and in-situ measurements of snow accumulation and surface-mass balance across the archipelago (Aas et al., 2016).
For grid cells covered by glaciers, the land-surface scheme of WRF was replaced by a modified version of the CMB model of 185 (Mölg et al., 2008(Mölg et al., , 2009, specifically adjusted for Arctic conditions (Aas et al., 2016). The model simulates the development of multi-year snowpacks and their transition into firn and ice. The CMB model employs meteorological variables generated by WRF, near-surface temperature, humidity, pressure, wind speed and incoming radiation to solve the surface energy balance and determine the energy available for melt. Solid precipitation along with surface and subsurface melt then yield the column-specific mass balance over 17 layers down to 20 m depth. Variables are computed at a 20 seconds temporal resolution and are then aggregated into daily values.
Daily glacier runoff is determined as the difference between a production and a retention term of liquid water at or near the glacier surface. Production of liquid water is given as the sum of surface melt, internal melt and rain (liquid precipitation).
Meltwater retention is the sum of internal refreezing within the snow and firn, superimposed ice formation, i.e. water refreezing on top of impermeable ice, and liquid water storage or more precisely, the change in liquid water content. Meltwater production 195 is highest at lower glacier elevation, but not restricted to the ablation area. At higher elevation within the accumulation area, locally produced meltwater may be stored in the snow and firn column, thus reducing or preventing runoff. Runoff from each region is first computed in absolute terms (Gt; Fig. 1b), and then normalized by the associated area of the sea (km 2 ), up to a defined distance from the coast (10, 20 or 50 km). This yield specific runoff received by the sea in terms of mm water equivalents (RUNOFF, in mm w.e.), i.e. the same units as used for expressing precipitation amounts or specific glacier mass 200 balance. Note that our CMB model does not include a scheme for transport and routing of meltwater. The exact location of meltwater input to the fjords and ocean is therefore unknown. However, this does not compromise our regional-scale analysis, where all glacier runoff generated within a primary hydrological region, drains into the same associated fjord system or adjacent sea.
Mean specific climatic net mass balance of Svalbard glaciers for the period 2003-2013 was negative, -257 mm w.e. yr −1 , 205 which corresponds to a mean annual mass loss of about 8.7 Gt (Aas et al., 2016). Interannual variability in climatic mass balance is large, and dominated by a high variability in summer ablation. This is closely reflected in the annual cumulative runoff curves for the various hydrological regions (Fig. 1b). Regional glacier runoff is a function of the total regional glacier area and regional specific ablation. On average, Svalbard-wide specific glacier ablation and thus total annual glacier runoff amounted to 919 mm w.e. and 31.2 Gt, respectively, with a minimum in summer 2008 (673 mm w.e.; 22.9 Gt) and a maximum 210 in summer 2013 (1508 mm w.e.; 51.3 Gt).

Ocean data
Chlorophyll a concentration (CHL, in mg m −3 ) in near-surface waters was quantified using satellite data from the European Space Agency (ESA) Ocean Colour Climate Change Initiative (CCI). We used Arctic reprocessed version L3 data obtained from the Copernicus Marine Environment Monitoring Service (CMEMS), providing 8-day means of merged, bias-corrected 215 remote sensing reflectance at 1-km resolution from 1998 to 2014 (http://marine.copernicus.eu). This product merges reflectance data from SeaWiFS, MODIS-Aqua and MERIS sensors by realigning the spectra to that of the SeaWiFS sensor. Chlorophyll a is estimated from the OC5ci algorithm, which is a combination of two ocean colour algorithm for chlorophyll retrieval. The

Statistical analysis
All data (CHL, RUNOFF, SST, MLD, SIF) were first aggregated into regional time-series with the same 8-day temporal resolution as CHL. For each of the 14 hydrological regions (plus one sub-region), we constructed three time-series of different spatial scale and near-shore influence: 0-10 km, 10-20 km and 20-50 km distance from land. Main emphasis is on 0-10 km from land, as this covers the major fjord systems where we expect largest potential RUNOFF effects.

235
To test if associations between RUNOFF and CHL were statistically significant we restricted the data to late summer (June 13 to October 15, i.e. annual 8-day periods 21 to 36). This period includes the main glacier summer-melt period (mid-June to September) and is expected to start after termination of the phytoplankton spring bloom. For each region and spatial scale we considered the following generic model:

240
Here log(CHL r,t ) is the natural logarithm of CHL in region r (and a given distance interval from land) at time t, α r is the intercept, β r is the auto-regressive effect of CHL in the previous time step, c r is a row vector with coefficients for environmental effects, e r,t is a column vector with the environmental covariate values, ε r,t is a normal and independently distributed error term with variance σ r 2 /n r,t and n r,t is the number of CHL observations that were averaged to calculate CHL r,t . By weighting the error variance with sample size, region-time combinations with few CHL observations, e.g., due to cloud cover, have less 245 influence on results than region-time combinations with many observations.
To determine which environmental variables to include for each region, we used a two-step approach. We first found the best model without RUNOFF, using data for all years 1998-2014 (whereas RUNOFF was only available from September 2003 to September 2013). Variables were selected by step-wise adding terms if leading to lower value of the Akaike Information Criterion corrected for small sample size, AIC C (Hurvich and Tsai, 1989). The AIC C helps to find the best trade-off between 250 the goodness-of-fit of a model and the simplicity of the model; a model with lower AIC C is preferred over a model with higher AIC C . Terms only marginally significant (P > 0.05) were removed from the model. Nine candidate variables were considered at this step: (1) SST r,t , (2) SST r,t−1 , (3) SST r,t −SST r,t−1 , (4) log(MLD r,t ) , (5) log(MLD r,t−1 ) , (6) log(MLD r,t )− log(MLD r,t−1 ), (7) SIF r,t , (8) SIF r,t−1 and (9) SIF r,t − SIF r,t−1 . The difference variables for SST and MLD were included as possible indicators of mixing of deeper nutrient-rich water masses into the surface layer. We then added RUNOFF and 255 RUNOFF t−1 to the model selected in the first step, but only if leading to lower AIC C (for the reduced period with RUNOFF data) and only if the association was significant at P < 0.05.
To assess if key model assumptions were met, we checked if residuals were independent and approximately normal distributed. Specifically, Pearson residuals (i.e., residuals standardized to unit standard deviation) from the final model for each region were explored for independence by plotting the autocorrelation function and the partial autocorrelation function and for 260 approximate normality by plotting quantile-quantile normal plots. The final model for each region was uncorrelated in time and approximately normal distributed with a possible exception of region 22 in the analysis for 0-10 km from coast, which showed indications of unequal variance. We also checked if results were strongly influenced by a few outlying observations.
Outliers were identified as residuals more than 3.3 × standard deviations away from zero, which is expected to occur by chance for 1 out of 1000 normal distributed cases, i.e. for about 2-3 of the >2000 observations analysed. If outliers were identified, 265 we refitted the model with the outliers removed and report significant changes in results, but kept the outliers in the presented model. 13 residuals distributed among 10 regions were identified as outliers in the analysis for 0-10 km from coast, and similar number for other distances from coast. Removing these outliers had little influence on parameter estimates for RUNOFF effects (all the coefficients remained statistically significant at P < 0.05). All statistical analyses were performed using the R programming environment (R Core Team 2016).

Results
We first present regional associations of CHL with glacier runoff (Sec. 4.1), before moving on to associations with physicalocean and sea ice variables (Sec. 4.2). Interpretation of these results will be discussed in the following section (Sec. 5).
Our statistical model identifies the environmental variables that best explain the observed regional summertime CHL ( A positive association also exist for the subregion of Kongsfjorden/Krossfjorden (155), whereas no significant association exists for NW Spitsbergen (15) as a whole. Positive associations are mainly restricted to within 10-km distance from the coast, indicating that the RUNOFF effect on CHL is mainly limited to within the fjords. Fjords in Svalbard have a maximum width of typically less than 20 km and are thus entirely covered by this range. Beyond 10-km distance from the coast, as well as  specific runoff rates that exceed those in the other regions by 46% and 69%, for mean specific runoff rates over ten subse-305 quent summers from 2004 to 2013 and specific mean-annual-peak runoff rates, respectively. Our statistical model suggests that an increase in specific runoff of 10 mm w.e. 8-days −1 raises summertime chlorophyll a concentrations in these regions by 5.2% to 20.0%, or 9.3% on average, with a standard deviation of 4.6% (Tab. 1). During the annual peak discharge we estimate that runoff increases chlorophyll a by 13.1% to 50.2% or 28.4 ± 13.5% on average, compared to situations with no runoff.

310
There are both negative and positive associations with CHL and any of the physical ocean and sea ice variables, although only for a limited number of regions. Concerning sea ice variables, the current sea ice fraction (SIF) has little association with CHL ( Fig. 2d). However, there is a delayed positive association of CHL with SIF in northern Svalbard, mainly within 10 km from the coast (regions 15, 16, 23; Fig. 2e), but also 10-20 km (16) and 20-50 km (21), while CHL is negatively associated with a change in SIF at 0-10 km and 10-20 km (regions 12,15,17,21,24,31,41;Fig. 2f).

315
Moving on to sea-surface temperature (SST), current SST has a few positive associations at 20-50 km distance from the shore (regions 12, 14 and 17) and negative association north of Nordaustlandet at 0-10 and 10-20 km distance from the coast (24, 25; Fig. 2g). There is a positive delayed association of CHL and SST along the entire west coast of Spitsbergen at 0-10 and/or 10-20 km distance from the coast (12,13,14,15; Fig. 2h), as well as in Hinlopen straight off Northeast Spitsbergen (17).
There is a negative instantaneous association of CHL with SST north of Nordaustlandet (25). The association of CHL with a 320 change in SST is negative all around Edgeøya (31) and Barentsøya (41), as well as western Nordaustlandet (23) and weakly positive in the outer region of NE Spitsbergen (17), at 20-50 km distance from the coast (Fig. 2i).
Mixed-layer depth shows some positive association with CHL at the outer regions along the west coast of Spitsbergen (13, 14, 15) and Hinlopen (17; Fig. 2j). The delayed association between CHL and MLD is negative in two northern regions (16, 21) within 10-km from the coast and positive at 10-20 and 20-50 km for Isfjorden (14) and E Spitsbergen (11), respectively 325 (Fig. 2k). The change in MLD has a few both positive and negative associations (Fig. 2l).

Discussion
We first discuss the observed associations of summertime CHL with any of the environmental variables and provide physical and biological explanations. We start with the associations of summertime CHL with RUNOFF (Sec. 5.1), before moving on to ocean and sea ice variables which point at the effect of persistent sea ice coverage, and the influence of the West

330
Spitsbergen Current (Sec. 5.2). We then describe the seasonal evolution of chlorophyll a in relation to environmental variables (Sec. 5.3). Finally, we discuss challenges related to the use of remotely sensed chlorophyll a as a proxy of phytoplankton biomass (Sec. 5.4).

Glacier runoff effects on marine primary production
Our study suggests that the overall effect of glacier runoff on marine primary production is positive for 7 out of 14 hydrolog-335 ical regions, including the major fjord-systems in Svalbard. We find that regions that display significant positive associations between CHL and RUNOFF have a 26% higher mean summertime chlorophyll a, and a 19% higher mean annual maximum chlorophyll a, than regions without such association (Tab. 1). Regions which display a positive association between CHL and RUNOFF are characterized by a highly variable fraction of tidewater-glacier drained area, ranging from 40.2% for Isfjorden to 87.4% for Southern Spitsbergen, with a regional mean of 62.3 ± 21.0%. This is slightly less, than the corresponding mean 340 value of 66.4 ± 21.0% in the other regions.
Field observations across the Arctic show that glacial fjords dominated by tidewater glaciers have generally higher productivity than those dominated by land-terminating glaciers (Hopwood et al., 2020). Marine terminating glaciers are generally thought to enhance marine primary production, while land-terminating glaciers are thought to obstruct it (Meire et al., 2017).
Consequently, one might expect that regions with a high fraction of tidewater glaciers yield significant positive associations 345 between CHL and RUNOFF, whereas regions with a low fraction of tidewater glaciers, yield weaker positive, or potentially negative associations. However, we do not find a clear relationship between the fraction of tidewater glaciers and the sign or strength of associations between CHL and RUNOFF (Tab. 1). This indicates that a fraction of tidewater glaciers above ∼40% is sufficient to provide upwelling of subglacial discharge plumes capable of stimulating regional-scale marine primary production.
Alternatively, other mechanisms by which glacier runoff stimulates marine primary productivity may play a role.

350
While our method allows us to assess the overall effect of glacier runoff on regional-scale phytoplankton dynamics, it does not reveal the very mechanism(s), by which the effect is achieved. We suggest that the positive association between CHL and RUNOFF could be explained by the several processes, which may act independently or in combination, dependent on regional characteristics: (1) buoyant upwelling of subglacial discharge plumes at the calving front of tidewater glaciers (a few tidewater glaciers may be sufficient to fuel primary production in the entire fjord system); (2) glacier runoff may enhance the general 355 estuarine circulation; and (3) glacier runoff provides a direct source of limiting nutrients. The first two points are considered indirect effects and the third a direct effect of glacier runoff on marine primary production.
Considering the first mechanism, buoyant upwelling of subglacial discharge plumes is associated with the entrainment of large volumes of ambient sea-water from deep to intermediate depth. This process is considered to deliver significant quantities of nutrients to surface waters (Svendsen et al., 2002;Meire et al., 2017;Kanna et al., 2018;Hopwood et al., 2018). These 360 nutrients are first expected to enhance primary productions some distance away from the glacier front, where light conditions become more favorable as progressively more suspended particles have settled deeper into the water column (Etherington et al., 2007;Halbach et al., 2019;Hopwood et al., 2020). Glacier erosion rates, the amount and size of suspended particles and thus glacier runoff effects on light regime is controlled by the glacier bedrock lithology as well as subglacial drainagesystem configuration and total discharge (Halbach et al., 2019). Glaciers in Svalbard are grounded at shallow depth compared to 365 glaciers in Greenland. Entrainment factors are therefore expected to be significantly smaller for Svalbard than for Greenland, as they scale with the depth at which subglacial discharge enters the water column (Hopwood et al., 2020). Nevertheless, Halbach et al. (2019) found nutrient upwelling in Kongsfjorden to be significant source of nutrients to the euphotic zone, as comparably small discharge volumes were sufficient for the plume to reach the surface (Slater et al., 2017) and plumes were present for a long timeperiod during summer (How et al., 2017). In addition, upwelling of ammonium released from the shallow seafloor of 370 Kongsfjorden was found to be a significant source of bioavailable nitrate (Halbach et al., 2019).
The second mechanism concerns the estuarine circulation, driven by down-fjord katabatic winds, which facilitate the export of relatively fresh or 'brackish' surface waters out of the fjord (e.g., Svendsen et al., 2002). This outflow of surface waters will induce a compensating return flow of warm and saline water masses from the shelf area at intermediate depth (Svendsen et al., 2002;Cottier et al., 2010). Sundfjord et al. (2017) used a high resolution ocean-circulation model, forced with glacial freshwater 375 discharge to simulate water exchange in Kongsfjorden, Svalbard. Simulations revealed that glacial freshwater discharge drives a strong outflow in the upper surface layer and a significant compensating inflow of Atlantic water in the upper 15-20 m, which was enhanced in times of peak discharge. The volume flux was strongly influenced by the local wind field. Vertical mixing by wind stress and tidal forcing provides a mechanism of bringing nutrients from intermediate water into the euphotic zone where they become available for phytoplankton, fueling primary production. Svalbard fjords are considered broad fjords, where 380 rotational 'Coriolis' effects play a role (Svendsen et al., 2002;Cottier et al., 2010). These rotational dynamics may contribute to vertical mixing of surface and intermediate depth waters, thereby enhancing the effect of the general estuarine circulation on nutrient availability in surface waters.
The third candidate mechanism concerns the direct fertilization of seas by nutrients contained in glacier runoff. In light of the reported low concentrations of nutrient in glacier meltwater, compared to ambient seawater (Halbach et al., 2019; Hopwood the role of subglacial discharge plumes, we cannot exclude that also the enhancement of the general estuarine circulation contributes to the observed positive effect of glacier runoff on marine primary productivity.

5.2
The role of ocean and sea ice variables on summertime CHL

390
The northern regions of Svalbard show a positive delayed association of CHL with SIF (Fig. 2e). This suggests high CHL in response to previously high SIF. The exact timing and breakup of sea ice is highly variable. It depends not only on the initial sea ice extent, thickness and stability, but also wind conditions and wave action, as well as sea ice conditions further offshore (Cottier et al., 2010). In northern Svalbard, oceanic pack ice can prevent sea ice from being exported out of the fjord, thus extending the sea ice season (Cottier et al., 2010). This is expected to lead to a significant delay of the phytoplankton spring 395 bloom. Presence of sea ice in the previous 8-day time period in the summer months in this region is thus an indication of hydrological spring conditions. This interpretation of a late spring bloom is supported by a negative association of CHL with changes in SIF, meaning that chlorophyll a is increasing when sea ice coverage is decreasing (Fig. 2f). The latter association is, however, not restricted to northern Svalbard, but significant also for other regions in Svalbard.

Advection of water masses of Atlantic origin 400
Similar as for the sea ice variables, we found delayed associations of CHL with SST and with changes in SST. A delayed positive association with SST is revealed along the entire west coast of Spitsbergen (Fig. 2h). This may indicate the influence of the WSC, flowing along the West Spitsbergen shelf and spilling onto the shelf. Note that the 50-km offshore-distance aligns approximately with the shelf edge along the western and northern side of Spitsbergen, indicating that variations in overflow of the West Spitsbergen current may affect the outer region (20-50 km). High SST points at the advection of warm Atlantic water, 405 which is also characterized by high salinity and nutrient content, thus capable of enhancing primary production and hence, CHL. The importance of warm saline Atlantic water for fjord and shelf water masses and the marine ecosystem was previously reported by Hegseth and Tverberg (2013) and Nilsen et al. (2016).
Around Edgeøya, a strong negative association of CHL with a change in SST coincides with the positive association of CHL with RUNOFF (Fig. 2i,c). Cooling SST may be associated with meltwater spreading out on the surface away from the coast, 410 meaning that the association of CHL with this variable and RUNOFF may reflect the same process. The negative association of CHL with change in SST might also be caused by increased stratification and nutrient limitation due to solar heating.
Vertical mixing is closely linked with the mixed-layer depth (MLD). The generally positive associations between MLD and CHL along the west coast is possibly caused by advection of Atlantic water onto the shelf, leading to increased vertical mixing as evident in a deepening of the MLD. Vertical mixing increases the supply of essential nutrients to surface water layers, 415 thereby increasing primary production as indicated by high CHL (Fig. 2j). multiply more slowly because they get access to less light (Sakshaug et al., 2009). Deepening of MLD can also have a dilution effect on near-surface phytoplankton biomass (e.g. Behrenfeld and Boss, 2014). These phenomena could explain the negative 420 associations between MLD and CHL in some northern regions.

Phytoplankton dynamics during the productive season
Our timeseries of chlorophyll a, glacier runoff, as well as physical ocean and sea ice variables allows us to put the summer bloom into a larger temporal context. We will discuss phytoplankton dynamics in Svalbard over the entire productive season, which lasts from about April to September, and compare our findings to those from other regions. Investigating fjords in SW

425
Greenland, Juul-Pedersen et al. (2015) were able to divide the productive season into three distinct phases: the spring bloom (April/May; phase 1), a transition period with low primary production (June; phase 2) and the summer bloom (July-August, phase 3).
To investigate whether these three phases can be identified in Svalbard, we average monthly means of all relevant variables over the period 2003-1013 (Fig. 3). The spring bloom typically occurs in May (Fig. 3a), coincident with increased solar 430 insulation, sea ice breakup (Fig. 3c) and initialization of a weak stratification, in line with phase 1 of Juul-Pedersen et al. (2015).
Stratification (shallow MLD; Fig. 3e) seems to be dominated by solar heating (increasing SST; Fig. 3d). Significant runoff starts in June, when stratification is already established (Fig. 3b,e), but CHL has declined from its spring-bloom value, indicative of nutrients depletion (Phase-2 in Juul-Pedersen et al. (2015)). Runoff during the later summer, i.e. July and August, coincides with a second period of high CHL (Phase-3; Fig. 3a,b), in some cases exceeding the monthly mean values during spring bloom.

435
Note that the spring bloom typically only lasts for a short time, i.e. one 8-day period, during which concentrations can be several times larger then what is reflected in the monthly mean. Peak values of CHL during summer may be lower, but more persistent, resulting in monthly mean values similar or larger than those during spring time. For regions that show a positive association between CHL and RUNOFF (e.g. regions 11,12,13,14,16), monthly mean CHL during summer (July-August) typically match or exceed that during spring bloom (May), with a minimum in June, in line with phytoplankton dynamics 440 described by Juul-Pedersen et al. (2015).
In NE Svalbard and Nordaustlandet (regions 17, 21-24), the 10-year monthly mean SIF is around 40-50% in June and 20% in July. Several regions in northern Svalbard showed a delayed association of CHL with SIF (regions 15, 16, 23; Fig. 2e) that indicates a delayed spring bloom. In this case, two separate production phases cannot be distinguished, at least at monthly temporal resolution. Instead, CHL during spring is low and steadily increases towards a maximum in July (e.g. regions 17, MLD is already established when significant glacier runoff starts in July. Generally lower CHL in June than May suggests that phytoplankton may be nutrient limited when glacial melting sets in. The peak meltwater discharge coincides with elevated

Challenges and uncertainties of satellite-based surface chlorophyll a products
Although remotely sensed chlorophyll a is a commonly used proxy of phytoplankton biomass, there are several limitations to this approach. Firstly, data sampling relies on sufficient daylight, clear skies and largely sea ice free conditions, as ocean colour sensors cannot detect ice-algae or phytoplankton cells beneath sea ice (Arrigo, 2014). For Svalbard, chlorophyll a 455 observations are typically limited to late March to early September. In the beginning and end of the acquisition period, spatial sampling is generally poor, due to the persistence of sea ice and limited day light (low sun angles). Spatial sampling is also poor under cloudy conditions, typical for Svalbard during summertime. The variable sampling intensity was accounted for in the statistical analysis, as 8-day time periods and regions with many satellite observations of CHL were given more weight in the analysis than periods and regions with few observations. Secondly, although the algorithm used to estimate CHL from 460 surface reflectance accounts for the possible presence of inorganic particles, bias from inorganic particles originating from glacial meltwater cannot be ruled out. Some fjords of Svalbard are heavily influenced by suspended sediments from terrestrial or subglacial runoff, which influences ocean colour significantly (e.g., Dowdeswell et al., 2015;How et al., 2017). Thirdly, sub-surface maxima of chlorophyll a, as may occur in summer situations, are easily missed by satellite sensors, because data retrieval is restricted to the upper layer of the water column down to the 1% photosynthetically available radiation (Lee et al., 465 2007). It should therefore be kept in mind that our results show what happens in near-surface layers, and not the entire water column. Furthermore, phytoplankton can rapidly respond to reduced light availability, for example due to suspended matter, by increasing the chlorophyll a concentrations in their cells (Finkel, 2001;Finkel et al., 2004). It is therefore uncertain whether possible increased chlorophyll a concentrations at high meltwater runoff also reflects increased phytoplankton biomass. Further verification of remotely sensed chlorophyll a as a proxy of phytoplankton biomass in complex Arctic waters is required to gain 470 more confidence in the results from our statistical analysis. This can only be achieved by in-situ observations, extensive in both space and time, including simultaneous measurements of phytoplankton biomass, glacier runoff and nutrient concentrations in different water masses.

Conclusions
We investigated the effect of glacier runoff on regional-scale phytoplankton dynamics in Svalbard by combining chlorophyll a 475 from satellite ocean colour with glacier mass-balance modelling. Statistical analysis of regional timeseries revealed significant positive association of CHL and RUNOFF for 7 out of 14 primary hydrological regions. Our results suggest that the overall effect of glacier runoff on marine primary production in these regions is positive, despite counter-acting effects of glacier runoff on the availability of light and essential nutrients, both of which are required for an increase in phytoplankton biomass. We find that regions that display significant positive associations between CHL and RUNOFF have a 26% higher mean summertime 480 chlorophyll a, and a 19% higher mean annual maximum chlorophyll a, compared to regions without such association. Our analysis suggests that an increase in specific runoff of 10 mm w.e. 8-days −1 raises summertime chlorophyll a concentrations by 5.2% to 20.0%, or 9.3% on average, with a standard deviation of 4.6%. During the annual peak discharge the effect is even larger, when 13.1% to 50.2% of the increase in chlorophyll a, or 28.4 ± 13.5% on average, is associated with glacier runoff.
Glacier runoff thus facilitates a secondary phytoplankton bloom in July to August, typically following a spring bloom in May 485 and a minimum in June, in line with in-situ observations from Greenland (e.g. Juul-Pedersen et al., 2015). In terms of monthly mean CHL, the magnitude of the summer bloom is similar or may even exceed that of the spring bloom.
Observations across the Arctic have shown that glacier fjords dominated by tidewater glaciers are typically more productive than those with only land-terminating glaciers (Meire et al., 2017;Hopwood et al., 2020). Here we find positive association of CHL with RUNOFF regardless of a highly variable glacier coverage, ranging from 35.5% to 81.2%, as well as glacier area 490 drained through tidewater glaciers, ranging from 40.2% to 87.4%. This indicates that upwelling-effects of nutrients from subglacial discharge plumes at a few tidewater glaciers may be sufficient to fuel regional-scale primary production. Alternatively, other mechanisms, such as enhanced estuarine circulation may play a role. Estuarine circulation relies on a relatively fresh surface layer which may originate from runoff of both land and marine-terminating glaciers and down-fjord wind, which are typical for Svalbard fjords (Svendsen et al., 2002;Cottier et al., 2005;Sundfjord et al., 2017).

495
The association of regional-scale CHL with RUNOFF is typically restricted to the major fjord systems, and within 10-km distance from the coast. For regions characterized by open coastal conditions, and beyond 10-km distance from the coast, the relationship between glacier runoff and marine primary production vanishes. Since our statistical model also accounted for physical-ocean and sea ice variables, we are able to identify other factors which influence chlorophyll a in Svalbard during summer. These factors include sea ice conditions, especially in northern Svalbard, pointing at the influence of persistent sea 500 ice and late sea ice breakup. Furthermore, associations of CHL with SST and MLD along the West Spitsbergen shelf indicate the role of the West Spitsbergen Current, i.e. the advection of warm saline and nutrient-rich water masses of Atlantic origin.
To our knowledge, this study is the first to link large-scale chlorophyll a from satellite ocean colour, an indicator of phytoplankton biomass, with glacier runoff from glacier mass balance modeling. Our method can be applied on a regional to pan-Arctic scale, thereby complementing valuable in-situ observations which are only available from a few sites and often of 505 short duration, thus not capturing inter-seasonal to inter-annual variability.  Svendsen et al. (2002) where the red arrow delineates the West Spitsbergen Current (WSC) and the blue arrow the Arctic Coastal Current (ACC), originating as East Spitsbergen Current (ESC). (b) Regional timeseries of annual cumulative glacier runoff extracted from climatic mass-balance simulations by Aas et al. (2016).

Predictor variable
Partial e ect on log(CHL t ) for each region Figure A1. Partial effects of environmental variables on chlorophyll a, CHL, within 10 km from the coast. Each row shows the model (Eq. 1) for one hydrological region (Tab. 1). Each panel shows the relationship between a predictor variable (x-axes) and CHL (y-axes), with lines showing estimated partial effects, and points showing partial residuals. Blank panels imply that the variable was not selected. Asterisks show statistical significance at levels 5% ( * ), 1% ( * * ) or 0.1% ( * * * ). Predictor variable Partial effect on log(CHL t ) for each region Figure A2. Partial effects of environmental variables on chlorophyll a, CHL, within 10 to 20 km from the coast. Each row shows the model (Eq. 1) for one hydrological region (Tab. 1). Each panel shows the relationship between a predictor variable (x-axes) and CHL (y-axes), with lines showing estimated partial effects, and points showing partial residuals. Blank panels imply that the variable was not selected. Asterisks show statistical significance at levels 5% ( * ), 1% ( * * ) or 0.1% ( * * * ). Predictor variable Partial effect on log(CHL t ) for each region Figure A3. Partial effects of environmental variables on chlorophyll a, CHL, within 20 to 50 km from the coast. Each row shows the model (Eq. 1) for one hydrological region (Tab. 1). Each panel shows the relationship between a predictor variable (x-axes) and CHL (y-axes), with lines showing estimated partial effects, and points showing partial residuals. Blank panels imply that the variable was not selected. Asterisks show statistical significance at levels 5% ( * ), 1% ( * * ) or 0.1% ( * * * ).