Darkening of the Svalbard Fjords Waters Observed With Satellite Ocean Color Imagery in 1997–2019

Introduction

In the mid-1970s a significant change in the Arctic Ocean system began. Since then, systematic satellite observations have proved that the Arctic Ocean ice cover is shrinking, and in the following decades the trend has continued with a rate of −13.3% dec^–1 (Serreze and Stroeve, 2015). The main driving force for the decline of the Arctic Sea ice area and its thickness was an increase in global mean temperature, which in the Arctic is two times stronger than elsewhere (Weller, 2005; Cohen et al., 2014). Sea ice forms a boundary layer in polar and sub-polar marine basins, that reduces the heat, mass, and momentum exchange between the ocean and atmosphere (McPhee, 2017), alters the surface albedo and reduces radiative energy transfer to the underlying water, which affects the heating rate and primary production (Arrigo, 2017). Sea ice, especially in the Arctic, regulates global climatic processes-, because its albedo is drastically higher than the open water surface or ice-free land areas. Therefore, enhanced melting and ice cover retreat results in higher doses of radiation reaching sea waters which when absorbed, additionally enhances water temperatures accelerating the heating rate that is the so-called sea ice albedo feedback (Curry et al., 1995; Kashiwase et al., 2017). The consequences of global warming have impacted the whole Arctic Ocean and its terrestrial watershed leading to a cascade of changes in the physical environment, biogeochemical process, and its ecosystems and their functions (Smedsrud et al., 2016). The sea-ice decline in the Arctic Ocean has been accompanied by an increase of its productivity (Ardyna and Arrigo, 2020) that has been observed in eutrophic shelf area, oligotrophic Central Arctic Basin and exchange gateways in Bering and Fram Straits (Kohler et al., 2007) and Canadian Archipelago (Hill et al., 2018). Terrestrial systems in the vast watersheds surrounding the Arctic Ocean are subject to changes too (Parmentier et al., 2017). Increased precipitation and accelerating permafrost thawing affected coastal erosion and water runoff from land (Kipp et al., 2018). Terrestrial fluxes of heat, dissolved, and particulate material through melt-water discharge and riverine inflow contributed to an increase of productivity in continental shelf margins. Through the exchange gateways, particularly, the Fram Strait, the Arctic Ocean received steadily increased advection of heat and salt from mid-latitudes with inflowing Atlantic Water (Walczowski, 2014; Polyakov et al., 2017). Together with the inflow of warm and salty Atlantic Water a significant northward advance of temperate phytoplankton and changes of the planktonic organism size distribution toward smaller organisms, i.e., pico- and nanoplankton have been observed (Oziel et al., 2017; Neukermans et al., 2018).

Changes in the light availability in the ocean are not often discussed as one of the consequences of global climatic changes. The light regimes in the surface waters of the Arctic Ocean are controlled mainly by the sea ice albedo feedback, terrestrial fluxes, and marine primary production. As those three processes are strongly affected by global warming, the light availability in the Arctic Ocean has been significantly altered across the whole ocean basin. The amount of light at various wavelengths reaching the ocean surface through the atmosphere is regulated by the Sun’s activity, the Earth’s orbit geometry, and transmission through the atmosphere, which depends on cloud cover, the concentration of water vapor and other gasses, and concentration of scattering and absorbing aerosols (Dera, 1992). In the water column, the transmission of solar energy into deeper layers depends mainly on the spectral optical properties of water and optically significant water constituents. In the simplified model, the optical properties of seawater are a superposition of spectral scattering and absorption coefficients of water molecules themselves, suspended organic and inorganic particles, colored dissolved organic matter (CDOM), dissolved inorganic salts, and dissolved inorganic gasses (Woźniak and Dera, 2007). Spectral absorption properties of each of them are additive and combined constitute inherent optical properties of the water (IOPs). In the Nordic Seas, the IOPs variability is driven by sea ice and glacier ice-melt processes (Noël et al., 2020; Schuler et al., 2020), phytoplankton dynamics in spring and summer (Winter et al., 2014; Pavlov et al., 2015; Kowalczuk et al., 2019), and mixing of multiple water masses that contain different CDOM concentrations (Granskog et al., 2015; Pavlov et al., 2015; Makarewicz et al., 2018). Melting ice provides input of fresh water, but also inorganic suspended matter (Bélanger et al., 2013) which sub-glacial melt is particularly rich in (Holinde and Zielinski, 2016). Plankton communities are a primary source of organic particles and are tightly related to particulate carbon stocks in Fram Strait (Engel et al., 2019). Phytoplankton also contributes to the autochthonous production of dissolved organic matter in the West Spitsbergen Current (Makarewicz et al., 2018). Therefore, this is a region of particular environmental significance, but the remoteness, ice cover, and harsh meteorological conditions drastically limit the amount of in situ data on inherent (IOPs) and apparent optical properties (AOPs). Most of the papers have been published in recent years (see citations above) and only few in the first decade of 21st century (Stramska et al., 2003, 2006). Long-term changes and evolution of IOPs and AOPs is largely unknown in the Arctic Ocean, especially in the European sector. The scarcity of the long-term records could be partially addressed using available ocean color data time-series for this region.

In general, a typical Svalbard fjord would consist of three water layers: Surface Waters at the top (local), Intermediate Waters at intermediate depths (advected), and Winter Water at the bottom of the water column (Cottier et al., 2010). The near-surface, relatively fresh layer is formed mostly by sea ice and glacial-melt, and freshwater runoff discharged in the spring and early summer months. A gradual formation of the density gradient reduces vertical mixing and particle sedimentation. Particulate suspended and dissolved matter trapped in the surface layer accelerate water heating in the months when insolation is the most intense. A warmer surface layer may provide optimal temperature conditions for various organisms over a season that would support primary production. Although photosynthetically available radiation, “PAR,” intensity values reach a maximum in June, the highest air temperatures are lagged by about a month, and water temperatures are usually the highest in late July and early August (Calleja et al., 2017; Payne and Roesler, 2019). The vertical stratification stabilizes in summer, which is accompanied by shallowing of the mixed layer and accumulation of the Chla near the sea surface, which was clearly demonstrated in the results of Calleja et al. (2017). Calleja et al. (2017) compared vertical profiles of T, S, Chla measured at locations in the inner, middle and outer part of the Kongsfjorden in May, August, and October. They found that in the central part of the fjord in August, the thermo- and pycnocline were only several meters deep. Variability of the depth of the top mixed layer is a key-aspect regarding remote sensing.

Spectral IOPs and AOPs, and the relationships between them served as theoretical foundations of the ocean color remote sensing (OCRS). The first so-called “blue-green ratio” algorithms were developed for retrieval of chlorophyll a concentration (Morel and Prieur, 1977; Gordon et al., 1980; Gordon and Morel, 1983) during the “proof of concept” of the Coastal Zone Color Scanner (CZCS) mission launched in 1979. Since then, OCRS has evolved and has been used successfully for studying spatial and temporal dynamics of phytoplankton abundance, IOPs and AOPs, and oceanic primary production across a range of spatial and temporal scales (Blondeau-Patissier et al., 2014). After the outstanding scientific and operational success of the CZCS mission, a number of the next generation of ocean color sensors have been launched into the Earth’s orbit: Sea-viewing Wide Field-of-View Sensor (SeaWiFS), placed in orbit in 1997, was followed by the Moderate-Resolution Imaging Spectroradiometer (MODIS) deployed on the AQUA satellite, launched in 2002 and operated by United States NASA (Emberton et al., 2016). In 2002, the European Space Agency (ESA) also launched an Earth observation satellite. The Medium-Resolution Imaging Spectrometer (MERIS), in operation until 2012, was highly successful. The successor of MERIS is the Ocean and Land Color Instrument (OLCI). In 2010, the Climate Change Initiative (CCI) with the special Ocean Color section (OC-CCI) was launched to create data series of climate-quality by combining all available satellite measurements (Sathyendranath et al., 2012, 2019). A set of common spectral bands derived from the ocean color sensors with overlapping operation periods provided after merging a time series of coherent satellite observations spanned more than two decades (Mélin et al., 2016). We have used available data sets to analyze temporal changes in chlorophyll a and selected IOP and AOP values in West Spitsbergen fjords and adjacent oceanic shelf waters.

Fjords have been defined as aquatic critical zones (ACZs) (Bianchi et al., 2020) as these environments have been recognized as important carbon burial hot-spots sequestering the largest amount of organic carbon per unit of area (Smith et al., 2015). In those specific environments, especially in boreal and polar zones, integrated and accumulated effects of interactions between atmosphere, ocean, cryosphere, and surrounding landmass could be observed, which make those geomorphological units very sensitive to climatic changes (Bianchi et al., 2020). Effects of environmental changes could be observed visually as a progressive decline of water transparency in fjords, as the increase in terrestrial flux and autochthonous production of optically significant water constituents decreases water transparency and limits the amount of light in the water column, which is commonly called ‘darkening.’ This term was introduced about a decade ago when this effect was described for the first time along the Norwegian coast and has come to be used too, especially beyond the marine optics community (Aksnes et al., 2009). It was associated with enhanced riverine input of the CDOM and suspended particles due to increasing precipitation over terrestrial watersheds. The impact of climatic changes on water transparency in specific fjords in Norway and Greenland has been the subject of recent studies by Mascarenhas et al. (2017) and Mascarenhas and Zielinski (2018, 2019). To this day, very few papers have been published on seawater optical properties in Svalbard fjords and adjacent waters. Sagan and Darecki (2017) reported differences in IOPs distribution between the Hornsund and the Kongsfjorden and Pavlova et al. (2019) reviewed impacts of climatic changes on the light field in Kongsfjorden based on the existing data record. The local bio-optical and remote sensing algorithms for retrieval of chlorophyll a and suspended matter concentrations were developed for Kongsfjorden by Son and Kim (2018) and in adjacent coastal waters by Konik et al. (2020). The only long-term remote sensing study was published on the changes of IOPs in the Kongsfjorden (Payne and Roesler, 2019) in the context of modeling light availability for phytoplankton distribution in this fjord. No other studies of long-term changes in seawater optical properties have been conducted in the other fjords of the Svalbard Archipelago or the Nordic Seas. Despite the lack of direct observations, there are multiple lines of indirect evidence of ongoing long-term darkening of fjords. Among those is an already observed increased occurrence of the jellyfish Periphylla periphylla in Svalbard waters (Geoffroy et al., 2018), and shifts of kelp forest to shallower depths in Kongsfjorden likely associated with light limitation (Bartsch et al., 2016). Besides, the following manifestations of climate change are expected to contribute to the further darkening of coastal waters on Svalbard: accelerated melting of glaciers (Kohler et al., 2007), the anticipated increase in precipitations (Førland et al., 2011) in combination with accelerated thawing of permafrost (Isaksen et al., 2007) and greening of terrestrial vegetation (Vickers et al., 2016) will likely further elevate inputs of suspended particles and CDOM into fjords, thus making surface waters less transparent.

In this work, we aimed to fill the knowledge gap on the long-term changes of the water transparency in West Spitsbergen fjords and adjacent oceanic shelf, using ocean color data archive and developed locally optimized algorithms for retrieval of IOPs, AOPs, and chlorophyll a concentration. The long term assessment in fjord waters in Spitsbergen is important due to the high value of ecosystem services provided by the coastal areas, and presumable adverse effects of darkening and importance of this effect for the environmental management in Norwegian coastal waters (Tiller et al., 2017) and from the ecological and scientific standpoint for the global community.

Data and Methods

Study Area

We have conducted our analyses in the coastal waters of West Spitsbergen, the biggest island of the Svalbard Archipelago. We have focused on the fjords directly connected with the Greenland Sea: the Hornsund Fjord, the Van Mijenfjord with the Bellsund, the Isfjorden, and the Kongsfjorden (Figure 1). Water temperature and salinity within these fjords are controlled by the interaction of several water masses, typically dominated by the Atlantic Waters (AW) with modification of terrestrial and meltwater input as well as the contribution of cold and fresh Arctic Water (ArW) carried by the Sørkapp Current (SC). The AW carried by the West Spitsbergen Current (WSC) is the major heat source in the area that accelerates the ice melting and has a significant effect on the biogeochemical cycles and nutrient loads (Hegseth et al., 2019). The composition of water masses in the Hornsund fjord is strongly regulated by the SC that flows around the southern tip of Spitsbergen, bringing cold and fresh ArW from the Barents Sea. It propagates northward along the West Spitsbergen Shelf (WSS) parallel to the coast and the WSC, which results in reduced advection of AW into the fjord (Promińska et al., 2017). The other fjords on the west coast of Spitsbergen are more impacted by the WSC and AW. In particular, the seabed topography of the Isfjorden and the Kongsfjorden allows for easier penetration of AW into these fjords (Cottier et al., 2010; Skogseth et al., 2020; Bloshkina et al., 2021). In contrast, the Van Mijenfjord is a semi-enclosed water basin, separated from the Bellsund that connects it with the Greenland Sea by two elongated islands Akseløya and Mariaholmen positioned across the entrance of the fjord. It is an example of a fjord controlled by semi-diurnal tidal forcing with a strong vertical stratification. Relatively fresh and warm water, formed by runoff from glaciers, resides near the surface, whereas salt and cold water that flow into the fjord, remain in deeper layer (Støylen and Fer, 2014).

Satellite Data

As the main indicators describing quantitative changes in water transparency, we selected the variables routinely retrieved by OCRS and directly linked with the notion of ‘darkening.’ Of the three analyzed variables, chlorophyll a concentration, Chla, is a proxy for phytoplankton biomass. The total absorption coefficient at 443 nm, a_tot(443), overlaps with the main maximum of phytoplankton pigments absorption spectra and contains the contribution of other absorbing constituents in natural water: CDOM and suspended particles. The solar downward irradiance diffuse attenuation coefficient at 490 nm, K_d(490), describes the transmission of solar energy through the water column as the rate at which the solar downward irradiance is attenuated with increasing depth. It is one of the standard AOPs of the water column retrieved from satellite data. However, we based our analyses on the local algorithms and we used only remote sensing reflectance products at several wavelengths.

The satellite data used for the analysis were obtained from the EU Copernicus Marine Environment Monitoring Service (CMEMS¹). This was a long time series of daily mosaics merged from multiple optical satellite sensors (MODIS, VIIRS, SeaWiFS, MERIS, and OLCI). The merged product had been developed in the framework of the ESA Climate Change Initiative (ESA-CCI/C3S project) (Sathyendranath et al., 2012, 2019). We exploited the remote sensing reflectance, R_rs [sr^–1], product at six spectral bands: 412, 443, 490, 510, 555, and 670 nm from the Arctic region in a form of Level-3 OPTICS product of 1-km ground resolution (CMEMS id: OCEANCOLOUR_ARC_OPTICS_L3_REP_OBSERVATIONS_0 09_068) that covered years 1997 – 2019. These data set had been provided to CMEMS by Plymouth Marine Laboratory (McClain et al., 1995; CMEMS, 2020).

Remote Sensing Algorithms for Determination of Chlorophyll a Concentration and Total Absorption Coefficient

Available spectral remote sensing reflectance values were used to derive Chla, K_d(490), and a_tot(443) values at selected locations, and assess the darkening effects. We have used the empirical band ratio algorithms developed and validated in the Greenland Sea and the West Spitsbergen Shelf based on the in situ radiometric measurements conducted in the study area in summer seasons between 2013 and 2015 (Konik et al., 2020). Details of field work, methodologies of measurements and analytical methods that were used for determination of values of optically significant constituents and radiometric quantities needed for algorithm development and validation have been given in detail in studies by Makarewicz et al. (2018), Kowalczuk et al. (2019), and Konik et al. (2020). In brief, chlorophyll a concentration values were determined in collected water samples using the spectrophotometric method by Lorenzen (1967) and the formula by Strickland and Parsons (1972). The total absorption coefficient a_tot(λ) is the superposition of the spectral values of light absorption coefficient by CDOM a_CDOM(λ) and the light absorption coefficient by suspended particles a_p(λ). We omitted the influence of the light absorption coefficient after subtracting the pure water contribution to the light absorption coefficient, because it is almost constant in the natural waters and in typical for natural waters temperature ranges (Pope and Fry, 1997). The absorption coefficients by suspended particles were determined using Perkin-Elmer Lambda 650 double-beam spectrophotometer equipped with an integrating sphere. Details on sample handling, measurements protocol, and data processing have been included in Kowalczuk et al. (2017, 2019). The CDOM sample handling, measurements protocol, and data processing were described by Makarewicz et al. (2018). The OCRS algorithms, mathematical equations for determination of Chla and a_tot(443) in the Greenland Sea, validation details, and uncertainty assessments have been presented in detail by Konik et al. (2020).

Remote Sensing Algorithm for the Calculation of Solar Downward Irradiance Diffuse Attenuation Coefficient at 490 nm K_d(490)

The often-used spectral band, centered at 490 nm is located in the spectral range of maximum transmission (least attenuation) of solar energy in clear oceanic waters. At 490 nm, impact of phytoplankton is the least pronounced and it is possible to detect there more prominent signals by the other particles. In turbid coastal waters that are typical for the Arctic fjords, the spectral region of maximum transmission of solar energy can often be shifted toward the longer wavelengths, even to 560–580 nm. For this reason, the region-specific algorithm for retrieval of K_d(490) should be developed and applied in the Svalbard fjords.

We developed a locally optimized algorithm for retrieving the K_d(490) values based on four field campaigns conducted in the Svalbard fjords, mainly in the Hornsund fjord and the Kongsfjorden in 2009, 2010, 2016, and 2017 (Figure 2A). More than one hundred profiles of spectral downward irradiance E_d(λ) and upward radiance L_u(λ) have been made with the use of the high-performance free-fall aquatic profiler C-OPS (Compact Optical Profiling Systems, Biospherical Instruments, Inc.). C-OPS radiometer equipped with 19 wavebands (from 305 to 765 nm) was mounted in a free-fall backplane, allowing deployment far from the ship’s shadow. The instrument was equipped with a reference irradiance sensor mounted on the deck to measure incident irradiance and to control the stability of ambient light during deployment. The remote sensing reflectance, R_rs(λ), was calculated as a ratio of the upwelling radiance just above the water surface, L_w(λ), to downwelling irradiance measured above the water, E_s(λ). The water-leaving radiance L_w(λ) was obtained from the upwelling radiance estimated just below the water surface, L_u(z = 0^–,λ), and propagated through the water-air interface according to standard procedures described in the literature (Mueller et al., 2003). The K_d(490) for the subsurface layer was calculated as the local slope of ln[E_d(z,490)] measured within a depth interval spanning a few meters within the surface layer, depending on the extent to which the surface layer was homogeneous:

The optimal regression between K_d(490) and R_rs band ratio (R² = 0.85) has been found for R_rs(490)/R_rs(670) ratio, as presented in Figure 2B and Equation 1. The algorithm has been validated based on in situ radiometric measurements. The calculated mean normalized bias MNB (systematic error) and the normalized root-mean-square RMS error (random error) were 9 and 52%, respectively. We used this algorithm in our study to assess the long-term changes of the K_d(490) in the Svalbard fjords.

Statistical Analyzes

To assess the long-term variability of selected IOPs and AOPs in the fjords, we selected several locations labeled from L1 to L13 (Figure 1). The choice of locations allowed for capturing: open waters, the mouths of the fjords and the central parts of the fjords. There was at least one location placed within the fjord and one location on the West Spitsbergen Shelf near the mouth of a fjord. In order to minimize the noise influence we used median values of remotely sensed Chla, K_d(490), and a_tot(443) that were extracted from each processed map at 12 selected locations from a window of 3 pixels × 3 pixels (9 km² area) (Figure 1). We used values where at least 5 out of 9 pixels were valid. Only results from locations with characteristic patterns or trends are presented in the paper. Results from all locations are included in the Supplementary Figures 1–6. The L4 location was excluded. Since this position had no valid data in relatively long periods, we feared that this could potentially bias our trend analysis. The full picture of the amount of valid data may be found in the Supplementary Figure 7.

Regarding the log-normal distribution of Chla, a_tot(443) and the K_d(490) (Campbell, 1995), all statistical analyses were performed using median metrics for all parameters. Statistical analyses were performed using R-CRAN software [libraries: ggplot2 (Wickham et al., 2020), ggformula (Kaplan and Pruim, 2020)]. Boxplots depict median values with lower and upper hinges corresponding to the first and third quartiles, and whiskers extend to 1.5 ^∗ Inter-Quartile Range (IQR). The trend significance was assessed based on the lm fitting function (R Core Team, 2020). The regime shifts in the summer levels of Chla and K_d(490) values were identified using the changepoint library (Killick and Eckley, 2014; Killick et al., 2016). It identifies points within a data set where the statistical properties of a time-series change, and it is used in a variety of fields like finances or genetics, but it has been gaining recognition also in the climatology and oceanography (Giamalaki et al., 2018; Nordli et al., 2020). A time series is partitioned into smaller segments where data distribution remains homogeneous, and the beginnings and ends of these segments are characterized by an abrupt change of the mean value, which are sometimes called “tipping points” (Jackson et al., 2005). The algorithm used here is a combination of an optimal partitioning and pruning algorithm denoted “PELT the Pruned Exact Linear Time” (Killick et al., 2012; Jackson et al., 2005).

Results

Monthly Variability

The dark polar night and almost permanent cloud cover limit the number of available and usable satellite images to the period between March and September. Throughout the whole data series spanning 1997–2019, the broadest range of Chla concentration values were recorded mostly in April when phytoplankton blooms start in this region. It was visible particularly in the open waters at locations L1, L5, and L9 (Figure 3 and Supplementary Figure 1). In the Hornsund fjord at locations L12 and L13, there was also a substantial variability in March (Supplementary Figure 1). The spring bloom peaked in May when the local maxima of Chla were observed at most locations. It was followed by a slight Chla decrease in June and another increase with a stationary phytoplankton growth phase until the end of summer, sometimes sloping down in the last summer month (September). It is worth noting that at the locations outside the fjords, maximum Chla values were observed in May (L1, L5, and L11). On the contrary, in the central parts of the fjords, Chla reached maximum concentrations later in the year, usually around July – August (L2, L3, L6, L7, L8, L12, and L13) (Figure 3 and Supplementary Figure 1). The exception is L10, where both global and local maxima were noted only in July and August.

Differences between median values of K_d(490) for a given month were around two orders of magnitude and were the smallest among analyzed parameters. In most cases, the strongest variability was observed in April except L10 in the Van Mijenfjord and L12 and L13 at the Hornsund (Figure 3 and Supplementary Figure 3). At the first location, K_d(490) values increased gradually until July and started sloping downward in August. The latter two locations in the Hornsund fjord are the most interesting due to the fact that the highest variability was observed there in March when there were also observed the highest median values of K_d(490) during analyzed months (March–September). The seasonal changes at other locations in the Isfjorden and the Kongsfjorden depend on whether they are outside or inside the fjords. At the West Spitsbergen Shelf locations (L1, L5, and L11) (Figure 3 and Supplementary Figure 2) or the ones open to enhanced water exchange with the shelf, like L9, the highest medians of K_d(490) were observed in May, whereas at L2, L3, and L10 they were noted in July – August. In the Isfjorden and the Hornsund fjord the July – August local maximum is less visible, but in the inner parts of these two fjords, the higher K_d(490) values were noted later in the year (L6 – L8 and L12 – L13) (Figure 3 and Supplementary Figure 3).

The a_tot(443) changes partly follow the patterns of Chla with the high variability of values in April, local maxima in May, and a rather unimodal variability with the maximum values in August at L10. However, the local maxima in May could be clearly identified mainly at the open water locations (L1, L5, L9, and L11) or with higher exchange of water with the Greenland Sea (L3). In the inner parts of the Isfjorden, the Kongsfjorden or the Hornsund fjord no maximum in May was observed (L6 – L8 and L12 – L13) (Figure 3 and Supplementary Figure 2). The a_tot(443) gradually increased there until September (L6, L7, and L12) or slightly sloped downward in September (L2, L8, and L13) (Figure 3 and Supplementary Figure 2).

Due to the presence of the two local maxima in most analyzed cases – the one around May and the other in July – August, we decided to divide trend analyses and conduct them for the two seasons, separately. The low number of data points in March and April and very high variability of values of considered parameters would significantly disturb the multiyear trend analysis. Therefore, we decided to limit long-term trend analysis only to the period from the beginning of May until the end of September. We estimated separately multiyear trends in spring that covered May – June and in the summer that covered the months of July – September.

Long-Term Trends

Results of trend analysis for all selected locations for all considered optical parameters, Chla, K_d(490), and a_tot(443) are presented in Table 1. We have only included the trend estimates that were statistically significant at p-values < 0.05. The trend analysis revealed a clear seasonal differentiation (Figure 4). It should be noted that for the majority of selected locations, the trends estimates for Chla and a_tot(443) were statistically insignificant. Negative trends for Chla, K_d(490), and a_tot(443) in spring were identified at the mouth of the Kongsfjorden (L3) and in the Van Mijen Fjord (L10) (Figures 1, 4). The negative K_d(490) trends were also calculated in the Isfjorden (L6 – L8), at the northernmost location on the shelf (L1) and near the mouth of the Hornsund fjord (L12). The only positive trend in spring was observed for a_tot(443) in the central Isfjorden at L8 (Table 1 and Figure 4). The opposite pattern was observed in summer. Clearly, the positive trends in the Isfjorden, the Bellsund area and Hornsund fjord were identified. In total, an increase was confirmed at 10 and 9 out of 12 selected locations in the case of Chla and a_tot(443), respectively (Figure 4). At the locations on the shelf (L1, L5, and L11), significant positive trends for Chla and a_tot(443) were identified, and a positive trend of K_d(490) was confirmed at L1. In the case of the K_d(490), we identified positive and statistically significant trends at 4 out of 12 locations. Positive summer trends had weaker statistical significance compared to those calculated for Chla and a_tot(443), and were found at the northern location L1 on the shelf, in the inner part of Isfjorden (L7), in the Bellsund (L9), and the outer part of the Hornsund fjord (L12) with p-values < 0.01 for all of them (Table 1). In the inner part of the Hornsund fjord (L13) and on the shelf (L5), no statistically significant trends for K_d(490) were identified. However, in the recent years at L13, we have observed evidently higher variability of K_d(490) values what might indicate beginning of some regime shifts on the long term. Examples of trends calculated for Chla and a_tot(443) are characterized by the highest value of determination coefficient, R², that clearly present positive trend of linearly increasing parameter values within the inner part of the Hornsund (L13) and on the shelf adjacent to the entrance to the Isfjorden (L5) were presented on (Figure 5).

To expand information gained with trend analysis, we computed temporal distribution of anomalies of selected optical parameters. The anomalies were calculated as the difference between the median value of each month and the median value of all available data in given spring or summer season in the period between 1997 and 2019. At all locations, the amplitudes of K_d(490) anomalies were the smallest, which is consistent with the smallest variability of K_d(490) in the time series dataset (Figure 3). Positive anomalies were stronger than negative. The exceptionally high values were noted in the Isfjorden in 2013 and the Hornsund fjord in 2014 (Figure 6 and Supplementary Figure 3). Since 2010, we have observed at most of the analyzed locations the predominance of positive anomalies with higher amplitudes at the majority of locations, especially in the Kongsfjorden (L2), in the outer part of Isfjorden (L6 and L7), and in the Hornsund fjord (L12 and L13). A similar pattern was identified at the West Spitsbergen Shelf locations like L1, L5, and L9, except L11, where strong year-to-year changes from positive to negative anomalies occurred (Figure 6).

The trend analysis applied to the whole data series in years 1997 – 2019 gives us an overall view on the long-term changes in the Svalbard area. In order to get a more detailed view on the processes and changes throughout the analyzed time period, we plotted year-to-year anomalies and applied the changepoint detection PELT algorithm. Positive anomalies were predominant for both Chla and a_tot(443), especially in the Kongsfjorden (L2), the Isfjorden (L6 – L8), the Hornsund fjord (L12 and L13) and even in the Van Mijen fjord (L10). However, within the latter fjord, positive Chla anomalies were present along with negative a_tot(443) anomalies, all of higher amplitudes. This unusual situation could be explained by greater contributions of other constituents than phytoplankton pigments to the total absorption coefficient. On the shelf (L1, L5, and L11) also positive anomalies persisted after 2010, but at the northernmost shelf locations (L1 and L5), negative anomaly in 2019 was spotted (Figure 6 and Supplementary Figures 4,6).

The change-point method introduces some generalization to the anomalies view by determining some patterns and regimes observed throughout the time of the analyzes (Figure 7). Moreover, identified segments follow the anomalies groups and show broader patterns. In the central Isfjorden (L7) the first shift of the mean Chla values was in 1998, at the beginning of our data series, when the mean Chla increased from about 0.3 to 0.8 mg m^–3. It was followed by a long stationary period. Since 2017, there was another abrupt change to mean values reaching 1.4 mg m^–3, which is supported by 4 years characterized by clearly positive Chla anomalies (Figure 6). The negative Chla anomalies in the Bellsund (L9) contributed to the first segment of relatively low mean Chla levels around 0.4 mg m^–3 until 2001 when an evident increase to the level of 1 mg m^–3 was detected and has persisted since then. Larger variability was observed in the Hornsund fjord (L12) and at L1, the northernmost shelf location. At L1, the Chla values were about 0.3 mg m^–3, similar to concentrations at L7. In 2002, the mean concentrations increased to about 0.6 mg m^–3, but in 2006 Chla concentrations dropped to about 0.4 mg m^–3. However, it remained higher than at the beginning of the time series. In 2012, another shift took place and the mean Chla values around 0.7 mg m^–3 have persisted since then, despite some negative anomalies in 2015 and 2019 (Figures 6, 7). In the Hornsund fjord (L12), we observed a stepwise increase in Chla values from around 0.5 mg m^–3 to over 1 mg m^–3 reached in 2012 and persisting since then. However, a particular drop was revealed in the years 2008–2012, notably pronounced in the 2008 – 2009, which are also years of certain negative Chla anomalies (Figure 6).

The mean values of K_d(490) varied less in time than Chla and mostly in the second half of the analyzed period of time (2012–2019). The one exception was a peak observed in 2006 at L1, which was extreme. The K_d(490) doubled in comparison to the mean value in the years 1997–2012 (Figure 7B). Since 2012 the mean K_d(490) has been about 20% higher than at the beginning of the time series with the mean values around 1.2 m^–1, but within that last period, numerous strong positive anomalies were observed (Supplementary Figure 6). In the Hornsund fjord, despite clearly negative anomaly in 2009, situation seemed to be balanced until 2013, followed by a number of positive anomalies which are reflected in the results of change-point analysis since 2013. The last 2 years 2018–2019, especially a decrease of K_d(490) values in 2019 was detected as a downward shift, but depending on the next years it might turn out to be only short-term like the one in 2015 (Figures 6, 7). In the Isfjorden (L7), the average K_d(490) was 0.15 m^–1 with evident reduction of K_d(490) values in years 2008–2009 and 2014–2015, and a peak in 2013. Since 2015, the average K_d(490) has remained at the higher level of 0.2 m^–1. An upward shift of K_d(490) average by 0.1 m^–1 was first observed in the Bellsund (L9) in 2010. Since then, the average has remained above 0.2 m^–1, except for the years 2014 and 2015.

Discussion

Monthly variability in most cases showed two local maxima of the analyzed water characteristics during each year. The first maximum was typically observed in May, and it was particularly pronounced at the shelf locations and in the outer parts of the fjords. The other peak, was observed in late July – early August. This peak was particularly distinct in the central parts of the fjords. The increase of Chla and a_tot(443) in the Svalbard fjords during analyzed years (1997 – 2019), in particular, is associated with the spring phytoplankton bloom that is additionally stimulated by the advection of the warm AW inflow. The only positive trend in the spring (May – June) months was identified in the Isfjorden, where there is the strongest influence of the AW due to favoring seabed topography (Cottier et al., 2010). In previous publications, regarding phytoplankton blooms timing in the open waters of Fram Strait, it was acknowledged that blooms started in this region in May, reached maxima in June, and lasted through July (Cherkasheva et al., 2014). Satellite observations in this paper revealed that phytoplankton blooms within the Svalbard fjords start, in general, a month earlier – in April and peak in May. Although we do not have in situ measurements collected in April – June, the timing of the identified blooms in May is consistent with the findings of other studies conducted in the Svalbard fjords, e.g., the Kongsfjorden (Calleja et al., 2017; Ericson et al., 2019; Hegseth et al., 2019). The late summer maximum is probably of a more complex nature that will be further discussed. Anyway, this clear division combined with the awareness of the seasonal differences in the water column structure was a basis for analyzing available time series in the two separate seasons – spring (May – June) and summer (July – September).

A shallow mixed layer is typically favorable for accumulation of particulate matter within it. Payne and Roesler (2019) reported a systematic increase of the backscattering coefficient b_b(440) and absorption coefficient by colored detrital material (particulate and dissolved) a_dg(440) values from May to August in the Kongsfjorden, which was attributed to an increase in stratification. Higher quantities of optically active sea-water constituents become trapped in a confined mixed layer, i.e., within light penetration depth and satellite sensors detection range. It was estimated back in the 1970s that ocean color satellite sensors measure light that penetrates water column to the depth above which 90% of the diffusely reflected irradiance originates z₉₀ (except specular reflectance) (Gordon and McCluney, 1975). According to Jerlov (1976) in coastal waters light penetration depth z₉₀ usually does not exceed 10 m and in homogeneous water bodies such as the top mixed layer of oceanic waters, e.g., near the mouth of the fjord, the z₉₀ may range up to around 20 m (Gordon and McCluney, 1975; Szeligowska et al., 2020). This means that within the fjords analyzed here, only changes of the IOPs in the maximum top 20 m would be visible in the satellite images. Therefore, the summer season is a time when the probability of obtaining the most complete picture of the primary productivity is the highest. Thus, in order to observe and identify any long-term changes we focused on the summertime when strong stratification ensures accumulation of the optically significant sea water constituents within the z₉₀.

The most interesting result was the consistent increase of all analyzed parameters in July – August in the long-term at numerous locations (Figure 4). On the one hand, a high correlation between K_d(490), Chla, and a_tot(443) implies that the phytoplankton (with Chla as a proxy) controls the IOPs variability to a large extent in the area of study. Our observations seem to coincide with the positive trends in long-term AW inflow and glacier-melt runoff during the analyzed period (1997 – 2019) which are the two possible reasons, potentially explaining the positive trends in IOPs. The increase of the AW impact on the water masses within the Svalbard fjords, phytoplankton bloom dynamics, and acceleration of the melting processes were reported in some of the recent publications (Pavlov et al., 2013; Promińska et al., 2017; Hegseth et al., 2019; Bloshkina et al., 2021). The presence of the warmer AW may cause faster ice-melting from the bottom and higher heat content at the starting point which means that the temperature optima of the planktonic organisms may be reached or exceeded earlier in the year. It may result in earlier blooms onset or phenology shifts, or both. On the other hand, nutrients are the limiting factor for primary productivity, and nutrient levels remain low during summer, starting from June (Calleja et al., 2017; Ericson et al., 2019). However, particular attention should be given to glacier meltwaters considered an important local nutrient source (Stachnik et al., 2019). Microbial activity at the glacier surfaces and within the cryoconite holes where even nitrogen-fixing cyanobacteria are found, can provide a substantial source of nutrients into fjords (Telling et al., 2011; Anesio et al., 2017). Despite the fact that composition of the glacier and fjord bacterial communities are entirely distinct, a significant relationship between nutrients (NO₃^– and SiO₄^2–) and glacier melt-water bacterial communities within the fjords was already confirmed (Thomas et al., 2020). This suggests that acceleration of the glacier retreat may bring excess loads of nutrients and photosynthetic microbes that stimulate primary productivity, and may be responsible for the higher values of Chla and a_tot(443) observed here in July – August. Previous research also revealed that, in the Hornsund fjord, the contribution of the total precipitation over land (excluding winter snowfall), total precipitation over the fjord area, and melting of the snow cover over unglaciated areas are roughly less than 36% of the total freshwater input. In contrast, glacier melt-water runoff and frontal ablation of tidewater glaciers are predominant in the area (Błaszczyk et al., 2013). A coupled energy balance–subsurface model also confirmed this tendency for the whole of Svalbard and predicted a further increase of the glacier-melt runoff in the coming years (van Pelt et al., 2019). This means that the nutrient pool may rise for some time in the Svalbard area, it supports our observations of positive Chla and a_tot(443) trends and forecasts further possible increases at least as long as the glaciers are still present.

In spring, we mainly observed negative trends in the Svalbard area during 1997 – 2019 (Figure 3). However, Chla and a_tot(443) seemed to decrease only in the Van Mijen Fjord (L10) and around the mouth of the Kongsfjorden (L3). The negative trends of K_d(490) are visible also within the fjords or the entrance of the fjords. Due to higher year-to-year variability and more possible sources of satellite data contamination, conclusions should be drawn more carefully for this season. Firstly, we should consider the role of adjacency effects. Their impact may be twofold, and in both cases, they are closely related to the albedo. Typically, at locations where albedo is high, particularly due to the presence of snow or ice, photons are scattered into the atmosphere and contaminate registered radiance. It often results in overestimating water-leaving radiance in the neighboring areas (Bélanger et al., 2005). Analogical effects are observed near the clouds, because of their high albedo. It was proved that during the melting season cloud cover above the Arctic intensifies (He et al., 2019). Even though the snow and ice cover have decreased in the Kongsfjorden, a significant increase in humidity and the cloud cover in that area were confirmed with the highest variability in May – July (Maturilli et al., 2019). The impact of land topography may also play a role here. The narrow, elongated shape and steep walls of the fjords enhance scattering and multiply the adjacency effects (Rozwadowska and Górecka, 2017). The other aspect is the high year-to-year variability that in spring may stem from the fact that dynamics of the three factors controlling the environment in this area which are the AW inflow, the snow and ice cover, and glacier melt-water discharge, has been changing recently (Hegseth et al., 2019). In the Kongsfjorden, ice cover has been declining in the last 10 years (Pavlova et al., 2019). The estimated first snow-free day at Svalbard as a whole is around 29th of June, so both the sea ice and land snowmelt may intensify around April – June (Vickers et al., 2020). Except for winter, spring was the season when the fastest increase of air temperatures was identified based on the observations at Ny-Ålesund BSRN site at a rate of 1K ± 1.4K per decade (Maturilli et al., 2019). Processes in the water column cannot be analyzed separately from the atmospheric conditions and land influence. The air temperatures govern the pace of the ice and snowmelt in the whole area, and therefore govern the amount of the meltwater runoff from the land into the fjords. In spring, fresh meltwater forms a thin near-surface layer due to the low salinity. This is a relatively transparent layer (0–10 m) which is optically clear in comparison to the underlying water masses (Granskog et al., 2015), because meltwater effectively dilutes the CDOM (Kowalczuk et al., 2017, 2019). Although the meltwater could contain higher loads of the suspended particulate material (Bélanger et al., 2013), this does not contribute significantly to the decreased transparency, especially at the beginning of the melting time. Approximately half of the ice-free area in Svalbard falls beyond the vegetation season (Vickers et al., 2016). The ice and snow water themselves contains a depleted level of nutrients (Ardyna et al., 2011) which does not allow for intense phytoplankton growth. When it predominates the water column above the first optical depth, it leads to underestimating values of remotely sensed Chla, a_tot(443) or K_d(490), because subsurface chlorophyll a maxima usually form below ocean color sensor detection depth (Ardyna et al., 2013; Kowalczuk et al., 2019). This may be the main reason why we observe negative trends in spring. Smaller amounts of the sea ice combined with intense melting of it may result in earlier and more pronounced near-surface diluted layer observed from the satellites.

With the start of the vegetative season the water composition may change. Land vegetation is a potential source of CDOM, suspended particulate material or nutrients, especially in summer months when flora thrives. The previous analyses revealed an increasing trend of the vegetation growth in Svalbard during the last decades (Vickers et al., 2016). It is particularly worth noting that the previous studies conducted by Vickers et al. (2020) on the timing of the snow disappearance is consistent with the inter-annual anomalies of all three analyzed water properties. In the years when the first snow-free days were observed earlier, in most cases negative anomalies in spring and positive anomalies in summer were observed, whereas later snow cover retreat resulted in the opposite patterns. This may imply that earlier start of the melting period may cause intensified ice and snow melt resulting in formation of a deeper or more diluted top layer in spring which may be followed by enhanced vegetation growth and glacier melt that provide extra nutrients inputs. It is reflected in negative anomalies in spring, followed by positive anomalies in summer. The long-term increasing trends of the glacier melt and vegetation growth in summer by the other researchers support this explanation, as well as justify the increasing summer trends observed here in the long term (Vickers et al., 2016; van Pelt et al., 2019).

Apart from air temperatures the AW water is also a huge heat store that gives up its heat and stimulates the melting processes of the sea ice as well as near-shore snow and ice cover. An evident manifestation of the AW impact was the extreme inflow to the Kongsfjorden (L2 and L1 on the shelf in the fjord vicinity) in 2006 (Głuchowska et al., 2016). Positive summer anomalies of all analyzed water characteristics were observed in that year at the locations L2 and L1, with a particularly strong anomaly of K_d(490) on the shelf at L1 (Supplementary Figure 6) which was detected using change-point analysis as an exceptional year (Figure 7). In the inner Hornsund fjord the AW was recognized for the first time in summer of 2014 and has been observed there with varied frequency and intensity since then (Promińska et al., 2017). The year 2014 was suggested to be a breakthrough year for the Hornsund fjord (Głuchowska et al., 2016) and it was identified as such here based on the change-point analysis of K_d(490) where we may see a regime shift between 2013 and 2014 (Figure 7). The other years with some AW intrusions that turned both the Hornsund fjord and the Kongsfjorden into the ‘warmer’ mode were 2002, 2006, and 2013 when we observed positive anomalies (Supplementary Figures 4–6) (Promińska et al., 2017). On the other hand, strong influence of the SC current on the Hornsund environment (L12) was reflected in our research. In the years when the SC inflow into the Hornsund was indicated by Promińska et al. (2017) such as 2001, 2004, 2008, and 2011 we observed negative summer anomalies of Chla, a_tot(443) and K_d(490) values (Supplementary Figures 4–6). A general message of the comparison of anomalies and time patterns at various locations around Svalbard is the clear upward shift of the Chla starting from the open-water locations at the beginning of 2000s, and an even more pronounced increase of all three parameters Chla, a_tot(443) and K_d(490) around 2012 – 2013 that manifested in strong positive anomalies in summer months and significant changes of the regimes that probably contributed the most to the increasing trends identified on the long term (Figures 6, 7 and Supplementary Figures 4–6).

Conclusion

Based on the analysis of the time series of satellite images in years 1997 – 2019, evident positive trends were revealed in the Kongsfjorden, the Isfjorden, the Hornsund fjord and the Bellsund. Even though there were differences between the trend estimates and their significance, we confirmed positive summer trends (July – September) of a_tot(443), Chla and K_d(490) at 10, 9, and 4 locations out of 12 consequently. We hypothesize that these changes can be attributed to changes and interplay in other environmental factors, such as hydrography, glacier and river runoff, precipitations, and atmospheric warming. It also should be noted that the enhanced amplitudes of anomalies observed in recent years are in line with observations and models showing that growing variability in the system can precede its regime shift (Spanbauer et al., 2014; Bury et al., 2021). Moreover, between 1997 and 2019 at the majority of locations we observed a regime shift toward darker waters. Our work provides the first documented record of ‘darkening’ in Svalbard fjords and adjacent shelf waters, and call for further investigation of these potentially profound environmental and ecosystem changes in the area. This work sets the baseline for consideration of the coastal darkening in the environmental management plans on Svalbard. It should be also considered in ecological modeling aiming at prediction of environmental changes related to global warming.

Data Availability Statement

The original contributions presented in the study are included in the article/Supplementary Material, further inquiries can be directed to the corresponding author.

Author Contributions

MK, PK, AP, MD, and SS: conceptualization, formal analysis, and writing–review and editing. MK and MD: methodology and data curation. MK: software and visualization. MK, PK, AP, and MD: investigation. PK, MD, and SS: resources. MK, PK, and MD: writing–original draft preparation. PK: project administration. PK, AP, MD, and SS: funding acquisition. All authors contributed to the article and approved the submitted version.

Funding

This study was supported by the Svalbard Environmental Protection Fund (Svalbards miljøvernfond), project nr. 19/21 - Monitoring Darkening of Svalbard Fjords from space.

Conflict of Interest

The authors declare that the research was conducted in the absence of any commercial or financial relationships that could be construed as a potential conflict of interest.

Publisher’s Note

All claims expressed in this article are solely those of the authors and do not necessarily represent those of their affiliated organizations, or those of the publisher, the editors and the reviewers. Any product that may be evaluated in this article, or claim that may be made by its manufacturer, is not guaranteed or endorsed by the publisher.

Acknowledgments

We thank the crews of R/V Lance and R/V Oceania and colleagues for the help on board during collecting the sea-truth data.

Supplementary Material

The Supplementary Material for this article can be found online at: https://www.frontiersin.org/articles/10.3389/fmars.2021.699318/full#supplementary-material

Footnotes

1. ^ http://marine.copernicus.eu/

References