Correlograms showed some degree of periodicity in estimated chl-a in most tracks, with apparent spatial autocorrelation length scales of 20–100 m for each of the transect runs (various colors; **Figure 8**). While the highest correlation is seen at the lower length scales (e.g., <20 m), longer length scales are evident in some of the transects (e.g., peaks of 50 and 110 m in the burgundy trace). The fact that these peaks were not evident in all of the data is an indicator of the spatial heterogeneity of these ice algal communities.

## DISCUSSION

The main purpose of this study was to develop a methodology to assess a compact, multispectral optical sensor, configured for AUV deployments for estimating spatial variability of sea ice algal communities. While similar approaches using an ROV and under-ice sleds have been demonstrated (e.g., Lange et al., 2016, 2017; Meiners et al., 2017; Lund-Hansen et al., 2018), this is the first attempt to our knowledge to make such irradiance measurements using an AUV. It is well documented that sea ice algae are optically active and make a significant contribution to the spectral attenuation in ice (e.g., Mundy et al., 2007; Ehn and Mundy, 2013; Melbourne-Thomas et al., 2015) making this type of remote sensing possible. Our data showed that it was possible to develop NDI and EOF correlations though with low correlations (R 2 = 0.31 and 0.39, respectively, N = 60) between predicted and observed chl-a, with broad 95% prediction intervals. These relationships are less reliable than previous studies, which developed predictive relationships with R 2 -values of up to 0.89 (see summary in Melbourne-Thomas et al., 2015) possibly due to the presence of platelet ice increasing the complexity of the ice structure where predictive relationships are weaker (e.g., R 2 = 0.07; Wongpan et al., 2018). When using the averaged values for each site to generate the NDI for use with the AUV, correlation was increased to R 2 = 0.41 (N = 19). A study based on radiance measurements below the platelet-rich ice in McMurdo Sound showed an even lower R 2 of 0.07 for correlation of NDI (471, 416) and chl-a concentrations (Wongpan et al., 2018). The relatively poor performance of the Wongpan et al. (2018) sampling was attributed to a relatively small range in chl-a that corresponded to a large change in optics, and the difficulty in

sampling the fragile platelet ice, similar to our experience. Clearly this type of ice and location is relatively difficult for remote sensing. The difference in within-station variability between observed and required chl-a, and the high variance as observed by Wongpan et al. (2018) of measured chl-a may indicate the main reason for the inability of our EOF and NDI models to provide a robust prediction of biomass. Firstly, it is possible that biomass was playing a minor role in attenuation of light within the sea ice. This seems unlikely, as the spectral qualities of under-ice light showed distinctive shapes typical of sea ice algae, with rapid attenuation in the blue (440 nm) and red (675 nm), affected by chl-a and carotenoids (Fritsen et al., 2011).

A second possibility is that the fragile nature of the platelet ice that was present on site resulted in incomplete and variable collection of biomass. This would be consistent with the high variability in observed chl-a compared to observed under-ice spectra in station replicates. To address this possibility, we undertook the simple modeling exercise that compared the concentrations of chlorophyll observed with those required to reproduce the observed spectra. Estimation of missing parameters using optical models has previously been used by Mundy et al. (2007), who derived a ∗ λ from under ice spectra, ice properties and chl-a biomass. In our case, we have independent estimates of a ∗ λ , and use these to estimate required chl-a to achieve observed spectra. This approach is limited by the sensitivity of such models to scattering (Ehn and Mundy, 2013). The scattering environment for ice with high amounts of embedded

and extending platelets, as seen at Cape Evans, is poorly constrained and quantitative estimates of required chl-a may be difficult to achieve.

A third option thus relates to the footprint size of the irradiance sensor and the ice corer. While the core samples have a diameter of 75 mm, the TriOS spectroradiometer, would have an effective footprint diameter of 350–500 mm when positioned 100–150 mm below the ice (Nicolaus et al., 2010). The difference in sampling footprint size may be important in sea ice of the type described here, if a high degree of small-scale variability

in biomass is related to the heterogenous matrix of the sub-ice platelet layer. The low variance in the required chl-a compared to measured chl-a for replicates at each station supports the view that, by sampling over a larger area, optical measurements are integrating variance, whereas under the conditions encountered at Cape Evans cores are not. Therefore, we suggest that when conducting these types of field observations under heterogeneous ice properties (i.e., a high probability of platelet ice) that multiple replicate cores are extracted from a region equivalent to the bottom-ice footprint size of the spectroradiometer.

Noise in the relationship between chl-a and NDI could be generated if there were substantial differences in the absorption spectra of ice algae across the sampling sites. However, while we did find significant differences between two groups of sites, the selected wavelengths did coincide with the selected NDI bands and differences were minor. These differences were systematic in that the two groups corresponded to sites with low (<5 µM photons m−2 s −1 ) and very low (<1 µM photons m−2 s −1 ) irradiances, and the shifts in carotenoid complements between these groups are consistent with a shift from energy harvesting (fuco) in low-light environments, to photoprotection (DDX) in higher light environments. Similar light-dependent carotenoid shifts have previously been reported in sea ice communities (Alou-Font et al., 2013). Robinson et al. (1997) also found that the diatoxanthin:chl-a ratio increased at higher irradiances in McMurdo Sound.

The shifts in absorption spectra and pigmentation of sea ice algal at Cape Evans occurred even though the highest light experienced by ice algae was low in absolute terms (e.g., average under-ice PAR was 2.61 µmol photons m−2 s −1 ). Models suggested that 1 to 23 µmol photons m−2 s −1 at a median incident irradiance of 720 µmol photons m−2 s −1 penetrated to the ice-algal layer, depending on snow cover. It is within, or slightly above, the range of previously reported critical minimum under-ice irradiance levels for algal growth of 2 to 9 µmol photons m−2 s −1 (e.g., Horner and Schrader, 1982; Gosselin et al., 1985, 1986; Lange et al., 2015). Active ice algae photosynthesis have recently been measured at 0.17 µmol photons m−2 s −1 (Hancke et al., 2018). In locations where ice and snow conditions allow a greater range of under-ice irradiance, more extensive divergence of pigmentation and absorption properties may be expected. We note that, while not significant in our study, this could potentially confound the development of NDIs based on carotenoid-specific wavebands, especially where there was a correlation between biomass, irradiance and carotenoids.

Our study region was dominated by a snow free ice cover with small snow patches (qualitatively observed to be <5 m) which did not alter the bottom-ice light levels below the snow patches due to the horizontal propagation/scattering of light from the surrounding snow-free regions. Because our study site was dominated by a snow free ice cover, there is a smaller decrease in the spatial coverage of sufficient bottom ice light levels for ice algal growth than expected. However, this finding has perhaps more significant implications for regions with thicker snow cover and less coverage of snow-free/thin-snow regions. In thicker snow regions the horizontal propagation of light may substantially increase the areal coverage of suitable light levels at the ice bottom around these sparse snow-free/thin-snow features. A future hypothesis to test is whether the scales of spatial variability in the algae compare to observed scales of snow cover.

If our failure to produce a robust NDI for use with the AUV data is due to the mismatch between footprint sizes of the ice corer and the spectroradiometer, our procedural difficulties relate to our calibration exercise rather than the performance of the AUV-mounted instrument. The AUV itself collected repeatable transects of multispectral optical data. NDIs can provide robust estimates of ice-algal biomass (Mundy et al., 2007; Fritsen et al., 2011; Melbourne-Thomas et al., 2015) and while our estimated chl-a concentration may have very wide confidence intervals, within the limitations of the large footprint of the AUV sensors, then patterns of variability of biomass seen at scales of several decimeters to meters appears to be reliable. Within this constraint, the observed pattern of spatial variation of biomass, with its patchiness on scales of 50–100 m, provides insight into the dynamics of sea ice algal communities and the design of sampling campaigns that need to be carried out to obtain area-representative samples. At Cape Evans, for example, the length scale of autocorrelation suggests that coring would need to be undertaken at ∼10 m intervals along 200 m transects to obtain adequate representation. The mechanisms underlying this broad scale pattern of biomass are as yet unclear. However, our data suggest that at present no single approach is adequate to fully understand multi-scale variability in sea ice algal biomass.

## CONCLUSION

In this study, we endeavored to determine whether AUV-mounted sensors could obtain systematic, transect-based information on sea ice algae from their optical signature. To achieve this, we attempted to develop a predictive model linking the under ice irradiance spectra to the biomass of bottom-ice algae and then applied this to AUV-derived data. Such approaches are now well accepted, but, while we were able to show a statistically significant relationship (R 2 = 0.41), this had a higher prediction error than other models produced using similar approaches in other ice-covered environments. We found some variability in the absorption characteristics of ice algae due to shifting pigment complements but conclude that the most likely challenge to development of a robust algorithm may be the small-scale variability of ice algal biomass. This is a feature of sea ice algae in general, but is particularly acute in our example as the under-ice environment at Cape Evans provided a physically complex structure, where variable amounts of platelet ice were frozen into first year sea ice and formed part of the ice-algal habitat.

Optical methods to estimate ice-algal biomass produced a smaller variance than coring, which we suggest reflects its ability to integrate biomass over a larger footprint than coring. Thus, our results suggest that remote assessment of biomass based on optical properties may, with suitable calibration for local algal attenuation, provide an estimate of the ice algal resource over spatial scales that is more robust

than coring approaches where small-scale variability of ice algal biomass is present. We used this approach, coupled to an AUV, to demonstrate what appear to be autocorrelated patterns of sea ice algal biomass on length scales of 10's to 100's m.

While sea ice imposes logistical challenges on the use of underwater technologies, the development of sensor arrays and the platforms that carry them offer a way forward to better understand sea ice algae communities. The development of such non-invasive techniques to measure environmental conditions under ice at scales commensurate with the inherent variability of these algal communities allows new understanding that would not be achievable with traditional methods and this alone represents a step change in our approach to studying these communities.

## AUTHOR CONTRIBUTIONS

All authors listed have made a substantial, direct and intellectual contribution to the work, and approved it for publication.

## REFERENCES

## FUNDING

This work was supported by the Air New Zealand through the New Zealand Antarctic Research Institute (Project NZARI-2014- 3). In addition, the authors would like to thank the ongoing support of the Danish Council for Independent Research (Project DFF – 1323-00335) and the Australian Research Council's Special Research Initiative for the Antarctic Gateway Partnership (Project ID SR140300001). Griffiths University and the Defense Science and Technology Group (Australia) also contributed through equipment support. BL was partly funded during this study by a visiting fellowship from the Natural Sciences and Engineering Research Council of Canada (NSERC).

## ACKNOWLEDGMENTS

We thank the Antarctica New Zealand for the provision of logistic support to ensure that this research was possible. Final thanks goes to Rowan Frost (Australian Maritime College Search) and the K081 Team for all of their work to make this project a success.

and transmittance of sea ice. Cold Reg. Sci. Technol. 62, 14–28. doi: 10.1016/j. coldregions.2010.03.001

**Conflict of Interest Statement:** 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.

The reviewer PD declared a past co-authorship with one of the authors BS to the handling Editor.

Copyright © 2019 Forrest, Lund-Hansen, Sorrell, Bowden-Floyd, Lucieer, Cossu, Lange and Hawes. This is an open-access article distributed under the terms of the Creative Commons Attribution License (CC BY). The use, distribution or reproduction in other forums is permitted, provided the original author(s) and the copyright owner(s) are credited and that the original publication in this journal is cited, in accordance with accepted academic practice. No use, distribution or reproduction is permitted which does not comply with these terms.

# Carbon Dynamics During the Formation of Sea Ice at Different Growth Rates

#### Daniela König1 \*, Lisa A. Miller 2 , Kyle G. Simpson2 and Svein Vagle2

1 Department of Environmental Systems Science, Institute of Biogeochemistry and Pollutant Dynamics, ETH Zurich, Zurich, Switzerland, 2 Institute of Ocean Sciences, Fisheries and Oceans Canada, Sidney, BC, Canada

Controlled laboratory experiments have shed new light on the potential importance of brine rejection during sea-ice formation for carbon dioxide sequestration in the ocean. We grew ice in an experimental seawater tank (1 m3 ) under abiotic conditions at three different air temperatures (−40◦C, −25◦C, −15◦C) to determine how different ice growth rates affect the allocation of carbon to ice, water, or air. Carbonate system parameters were determined by discrete sampling of ice cores and water, as well as continuous measurements by multiple sensors deployed mainly in the water phase. A budgetary approach revealed that of the initial total inorganic carbon (TIC) content of the water converted to ice, only 28–29% was located in the ice phase by the end of the experiments run at the warmest temperature, whereas for the coldest ambient temperature, 46–47% of the carbon remained in the ice. Exchange with air appeared to be negligible, with the majority of the TIC remaining in the under-ice water (53–72%). Along with a good correlation between salinity and TIC in the ice and water samples, these observations highlight the importance of brine drainage to TIC redistribution during ice formation. For experiments without mixing of the under-ice water, the sensor data further suggested stronger stratification, likely related to release of denser brine, and thus potentially larger carbon sequestration for ice grown at a colder temperature and faster growth rate.

## Keywords: sea ice, growth rates, CO2 system, brine drainage, tank experiment

## INTRODUCTION

While the polar regions are generally considered an important sink for atmospheric carbon dioxide (CO2), the influence of sea ice on air-sea carbon dynamics is unclear. Originally assumed to be an impermeable lid hindering air-sea gas exchange (e.g., Stephens and Keeling, 2000; Tison et al., 2002), sea ice has since been recognized to play an active role in the marine carbon cycle, as both biological and physico-chemical processes take place within and around the ice (Vancoppenolle et al., 2013). The sea-ice zone is believed to act as a sink for atmospheric CO2 on an annual basis (Rysgaard et al., 2011), although the direction and magnitude of the carbon flux differ between the different stages of the seasonal sea-ice growth and decay cycle and also appear to depend on the methods used (Miller et al., 2015).

During the period of ice formation, the conceptual model developed by Rysgaard et al. (2011) predicts CO2 release to the atmosphere. When ice is formed, the majority of ions and gases present in the original sea water are expelled from the crystal lattice, whereby liquid brine intrusions with high salinities and gas concentrations are formed between ice crystals (Petrich and Eicken, 2016).

#### Edited by:

Jeff Shovlowsky Bowman, University of California, San Diego, United States

#### Reviewed by:

Soren Rysgaard, Aarhus University, Denmark Sebastien Moreau, Norwegian Polar Institute, Norway Daiki Nomura, Hokkaido University, Japan

#### \*Correspondence: Daniela König

<EMAIL_ADDRESS>
#### Specialty section:

This article was submitted to Cryospheric Sciences, a section of the journal Frontiers in Earth Science

Received: 27 July 2018 Accepted: 03 December 2018 Published: 18 December 2018

#### Citation:

König D, Miller LA, Simpson KG and Vagle S (2018) Carbon Dynamics During the Formation of Sea Ice at Different Growth Rates. Front. Earth Sci. 6:234. doi: 10.3389/feart.2018.00234

**73**

A small part of the brine supersaturated in CO2 is then predicted to be expelled to the sea-ice surface and to release CO2 to the atmosphere (Rysgaard et al., 2011). Small, positive CO2 fluxes out of young, newly formed sea ice have indeed been detected in the field (Else et al., 2011; Geilfus et al., 2013; Barber et al., 2014; Nomura et al., 2018) and in artificial sea ice (Nomura et al., 2006; Geilfus et al., 2016; Kotovitch et al., 2016). However, while the CO2 flux out of the young ice itself is likely positive (upward), the overall CO2 flux during ice formation in the marginal sea ice zone may be negative, depending on the presence of open water, including leads and polynyas, which may enhance CO2 uptake from the air (e.g., Anderson et al., 2004; Loose et al., 2011).

Sea-ice formation also has the potential to export carbon to deeper water layers (Rysgaard et al., 2007, 2011). Carbon dioxideenriched brine produced during sea-ice formation could sink to levels below the surface mixed layer and lead to a net annual uptake of atmospheric CO2 of up to 33 Tg C yr−1 , according to Rysgaard et al. (2011). This uptake could be further enhanced by the precipitation of calcium carbonate (CaCO3) minerals, namely ikaite. During the formation of ikaite, CO2 is produced, and the ratio of total alkalinity (AT) to dissolved inorganic carbon (DIC) of the brine is reduced (DIC refers here to the dissolved fraction of TIC; we use TIC for the remainder of this manuscript since our samples were unfiltered and may have contained solid carbonates). Assuming that the CO2 enriched brine sinks to deeper water levels while the ikaite crystals stay in the ice, this would lead to an enrichment of surface water alkalinity during ice melt (Rysgaard et al., 2007). This, in turn, could result in a higher uptake of atmospheric CO2, on the order of 83 Tg C yr−1 , globally (Rysgaard et al., 2011). This effect might be more pronounced at lower ice temperatures, since ikaite precipitation is more likely (in an open system) at the higher brine salinities encountered at lower temperatures (Papadimitriou et al., 2013). Differences in ice temperatures have thus been suggested to explain different AT to TIC ratios observed at multiple field locations (Rysgaard et al., 2007, 2009). However, the importance of ikaite precipitation for atmospheric CO2 uptake is debatable, as is transport below the mixed layer. The latter almost certainly does not occur in all sea-ice environments (Loose et al., 2011), and according to a recent modeling study (Moreau et al., 2016), the fraction of DIC exported to depth could be no more than 2% (about 4 Tg C yr−1 ) of all DIC rejected during sea ice growth. This fraction was found to be even lower for a global warming scenario in a similar study conducted by Grimm et al. (2016) due to changes to the oceanic overturning circulation. However, Parmentier et al. (2013) also stressed the importance of open water in the sea-ice cover (i.e., leads and polynyas) for the formation of dense water needed for export to depth, which may not have been adequately represented in their model. Thus, as such open water areas may be more prevalent in a warming future (Parmentier et al., 2013), the export of carbon may also be enhanced.

The influence of the ice growth rate on TIC export from the forming ice has received little attention so far, but some insight might be gained from desalination studies, given that inorganic carbon is mostly present in dissolved form in seawater. Our understanding of sea-ice brine drainage processes is rapidly evolving, and mushy layer theory has been suggested as a useful model to explain the underlying mechanisms of sea-ice desalination (Feltham et al., 2006; Notz and Worster, 2009). This theory predicts gravity drainage to be the dominant desalination process during ice growth (Notz and Worster, 2009). Gravity drainage results from the temperature gradient across the ice, as high salinity brine formed at the lower temperatures in the upper layers of the ice overlies less saline, and thus less dense brine formed at the warmer temperatures closer to the sea ice-water interface. This unstable brine density profile eventually leads to convective overturning within the ice, when the density gradient and the ice thickness are large enough to provide the required potential energy to overcome the energy loss due to dissipation during convection (Wettlaufer et al., 1997). Thus, mushy layer theory predicts a delayed onset of desalination, which has been observed in laboratory ice growth studies (Wettlaufer et al., 1997; Notz and Worster, 2009).

The impact of different ice growth rates on gravity drainage is not straightforward. While the dependence of gravity drainage on the temperature gradient would suggest more effective desalination at colder ambient temperatures, desalination is hindered by lower ice permeability at these colder temperatures, which can limit the convection to the ice layers close to the icewater interface (Notz and Worster, 2009). Both laboratory and numerical studies found that for constant boundary conditions (i.e., constant cooling), an inverse relationship between salt flux and growth rates holds, i.e., increased bulk ice salinity at faster growth rates as a result of the decreased permeability (Cox and Weeks, 1975; Wettlaufer et al., 1997; Griewank and Notz, 2013). Similarly, results from field studies comparing salinity profiles of ice cores to their growth rate history suggest that faster ice growth leads to higher retention of salt in sea ice (Nakawo and Sinha, 1981; Eicken, 1992; Gough et al., 2012). Assuming the TIC redistribution during ice formation is dominated by brine drainage, a higher retention of TIC in ice at faster ice growth rates is thus predicted.

The correlation between TIC and salinity could be disturbed by multiple abiotic and biological processes. Apart from biological transformations of TIC into organic substances, the precipitation of ikaite or the nucleation of CO2 bubbles both lead to inorganic forms of carbon that could migrate within ice differently from dissolved salts (Vancoppenolle et al., 2013). Furthermore, both processes are dependent on the in-situ temperature and salinity conditions and are possibly interconnected: increased CO2 partial pressure (pCO2) at high salinities could lead to the formation of gas bubbles, which in turn could promote the precipitation of ikaite by lowering the CO2 concentration in the brine (Papadimitriou et al., 2014). Therefore, TIC and salinity decoupling could be more pronounced in colder, and thus more saline, brine. Nevertheless, experimental and model studies on newly forming sea ice have found that ikaite precipitation and CO2 outgassing to the atmosphere are small compared to the efflux of dissolved inorganic carbon with brine to the underlying water (Rysgaard et al., 2007; Sejr et al., 2011; Moreau et al., 2015; Kotovitch et al., 2016). This indicates that the transport in dissolved form is predominant, and the TIC signal should largely follow salinity. However, past laboratory experiments have been conducted at a single, relatively high temperature, and might not be representative of faster ice growth conditions.

In the context of the large uncertainties in the air-ice and ice-water carbon fluxes, transformation and transport processes within the ice, and the impact of different ice growth rates, further insight into carbon dynamics during sea-ice formation is needed. Therefore, we conducted laboratory experiments simulating ice formation at three different growth rates with discrete and continuous measurements of CO2 system parameters in ice and water. In addition, we estimated the air-ice CO2 flux using a carbon budget approach. In addition to discrete sample analyses, the installation of sensors in the ice and water, including thermistors, Dopbeam sonars, and a water salinity probe, allowed us to further investigate brine drainage processes.

## MATERIALS AND METHODS

## Ice Growth Experiments

The ice growth experiments were conducted in a walk-in freezer at the Institute of Ocean Sciences in Sidney, BC, which is equipped with three separate refrigeration condensers. The condensers can be set to keep the ambient temperature of the freezer stable anywhere between 0◦C and −45◦C, apart from regular defrost intervals, when the temperature can increase by up to 10◦C over a period of 1 h (in the case of the coldest experiments). Three ambient temperatures (−15◦C, −25◦C, and −40◦C) were chosen to cover a wide range of ice growth rates within a temperature range observable in nature, and to be sufficiently cold to reach the desired thickness. The ambient conditions in the cold lab were monitored by a weather station (König, 2017), mean air temperatures and their standard deviations are listed in **Table S2**.

The cold lab contained a water tank (ca. 105 × 95 × 99 cm; **Figure 1A**) holding natural, open ocean seawater, which was collected on June 16, 2016, at Ocean Station P20 in the Northeast Pacific (50◦N, 145◦W) at a depth of 5 m. To minimize biological activity, the experiments were run in the dark and the water was filtered through a 1µm filter and sterilized with a UV lamp for 24 h before the start of the experiment and whenever new seawater was added to the tank. Whereas in general, the seawater was re-used for the different experimental runs and additional seawater was only added on two occasions, tap water was added to the tank after each experiment, to counteract the increase in salinity due to evaporation and the removal of ice cores. This had a minor impact on the chemical composition of the seawater, observable as a slight increase in TIC relative to AT. The tank was insulated on the sides, to prevent ice growth on the walls, and open at the top. To enable water sampling from underneath the ice, a sampling system with adjustable heating was installed in the tank, connecting multiple Tygon R tubes from within the tank to a port outside the freezer (**Figure 1A**).

The first set of experiments was conducted with stirring to prevent stratification and ensure that all deployed sensors would measure seawater representative of the whole body of water. For a second set of experiments, the water was not stirred (except during water sampling to avoid sampling stratified layers) to allow brine plume detection, as well as to test the impact of

mixing on the results of the previous runs. Each experiment was stopped (i.e., water and ice were sampled) once the ice thickness had reached about 20 cm. The stirring system was made up of a pump connected to PVC pipes arranged in a rectangle at the bottom (**Figure 1A**), which generated an observable vortex through the entire depth of the tank. In addition, two pumps were installed vertically at a depth of ca. 30 cm in order to prevent stratification.

An overview of all successful experimental runs is shown in **Table 1**, including ambient temperatures and whether or not the tank was stirred. A more comprehensive overview, including all (attempted) runs is given by König (2017).

## Water and Ice Melt Samples

In order to calculate a carbon budget, seawater and sea-ice melts were analyzed for TIC content, AT, and salinity. In addition, the seawater was analyzed for pH.

For each run, water was sampled before the ice formed ("preice" samples), and once the ice had reached a thickness of 20 cm ("under-ice" samples). All pre-ice samples were taken once water temperatures were close to or below 0◦C, and while the tank was well-mixed. For some runs, two pre-ice samples were collected to determine the impact of the rapid cooling before each run on the carbon dynamics. For the stirred runs, we collected some preice samples while the water was cooling down relatively quickly, which may have contributed to poor precision between replicates.

TABLE 1 | Overview of experimental runs.

In the case of the unstirred runs, water was stirred while cooled relatively slowly to about −1 ◦C before pre-ice samples were collected. Only after sampling was the stirring system turned off, and the water left to stratify. For the under-ice samples of these runs, the stirring system was turned on at least 2 h before sampling, to ensure adequate mixing of the stratified water column. For all sampling periods, the sampling port was heated beforehand. In the case of the colder runs, the heating power had to be increased, and the heating period extended. In general, 0.5– 2 h of heating was required to open the port and allow sampling. The system was flushed before sampling to remove residual water from the sample port tubing. In the case of under-ice sampling, a significant number of bubbles formed in the sample tubing due to decreased solubility caused by decreased pressure and warming of the water during sampling. While we tried to avoid bubble entrainment during sampling, it could account for some of the differences between replicates sampled during these periods.

Samples were collected in the following order: salinity, TIC and AT (duplicates), and pH (quadruplicates). In the case of the last three runs, additional salinity samples were collected after TIC and AT, as salinity appeared to change slightly during the under-ice water sampling (possibly due to some residual stratification in the tank, despite the stirring). Salinity samples were collected into 200 ml borosilicate bottles rinsed three times with sample before filling. Samples for TIC and AT were tapped into 250 ml borosilicate glass bottles, following a short initial rinse and an overflow of at minimum one bottle volume. Subsequently, the headspace was adjusted, samples were poisoned with 100 µl of saturated mercuric chloride, and the bottles closed using a greased glass stopper and electrician's tape, in accordance with standard operating procedures (Dickson et al., 2007). All TIC and AT samples were stored in the dark at 4 ◦C before analysis. Samples taken for pH analyses were collected directly into 10-cm spectrophotometric cells, which were rinsed, flushed for approximately 15 s, and closed without a headspace using two Teflon R stoppers, as described by Dickson et al. (2007). The cells were immediately placed in a 25◦C incubator for 1 h before analysis.

Ice cores were collected after the under-ice water sampling and melted for subsequent chemical analyses. For all runs, we took four ice cores from approximately the same locations in the tank (König, 2017), except for Run 6, when only three cores were retrieved. The sampling spots were chosen so that different regions of the tank were sampled to best capture any horizontal variability. The ice was cored manually using a Kovacs Mark II corer (9 cm diameter). After a short visual inspection and determination of core length inside the cold lab (still at the pre-set temperature), the cores were placed in gas-impermeable ALTEF R bags which were then taken outside the lab and immediately sealed with C-clamps. The bags were closed outside the cold lab, because they became brittle and the clamps were difficult to close at low temperatures. The bags were evacuated immediately after sealing, using a hand-held pump to minimize excessive suction and potential loss of gases trapped inside the core, as suggested by Miller et al. (2015). The cores from the final three experiments were weighed after evacuation so as to determine the ice density, using an Ohaus GT4100 balance.

After evacuation, the cores were placed into a dark container and melted overnight at room temperature. All cores took approximately 24 h to melt completely. Once fully melted, the ice melts were mixed and allowed to warm up until all visible solid particles were dissolved, so as to assure homogeneity. The high reproducibility of TIC and AT in the ice melt sub-samples (**Table S2**) confirms that this procedure was indeed adequate. Before sampling the melts, accumulated air inside bags was removed using a syringe. The air volume differed between the melts (ca. 40–200 ml), possibly due to natural heterogeneity, inconsistent evacuation, or, in the case of the highest air volume, a leaky bag or valve leading to uptake of air during melting. However, TIC, AT, salinity, and the ratios between these parameters were no different in the ice melt sample with the largest air volume (200 ml) than those from the replicate cores without a large head-space. Ice melt samples were sub-sampled by connecting a short Tygon R tube to the valves of the bags. First, two TIC/AT samples were tapped, preserved, and sealed like the water samples, except for a slightly smaller overflow (0.5– 1 bottle volume) due to the limited ice melt volumes. The rest of the ice melt was used for 1-2 salinity samples, depending on the remaining volume.

#### Chemical Analyses

Total inorganic carbon was measured coulometrically on a SOMMA system (Johnson et al., 1993), run in combination with a UIC coulometer, using standard operating protocols (Dickson et al., 2007). To ensure the accuracy of the measurements, certified reference materials (CRM batch 152 and 160, provided by Andrew Dickson, Scripps Institute of Oceanography) were measured daily, and all values were corrected accordingly. The analytical precision, calculated from differences between duplicate samples, was about 1µmol/kg for water samples, and 0.5µmol/kg for ice melt samples.

Total alkalinity was determined by open cell, potentiometric titration, using a custom built system following the procedure outlined by Dickson et al. (2007). The accuracy of this system was ensured by calibrating against reference materials (CRM batch 152, 160, and 161) at the beginning and end of each day. The analytical precision (determined as for TIC) was about 4µmol/kg for water samples, and 3µmol/kg for ice melt samples. Both TIC and AT analyses were conducted within 4 weeks of collection.

Water sample pH was determined spectrophotometrically within 2 h of collection with an Agilent 8453 spectrometer, using m-cresol purple as an indicator dye, following Dickson et al. (2007). Each absorption measurement was conducted four times and averaged. All values were corrected for a potential pH perturbation by the dye addition. The temperature within the cell was recorded immediately after each spectrophotometric measurement using a hand-held temperature probe (Fluke 5611T). Since some of the under-ice water salinities were outside the range for which the dye pKa was defined by Clayton and Byrne (1993), pH was calculated (on the free pH scale) following Miller et al. (2011), using the parameterization by Millero et al. (2009). The pH values were then converted back to the total scale using the HSO− 4 acidity constant of Dickson (1990), assuming the sulfate concentrations of the samples were proportional to salinity.

Salinity of water and ice melt samples was determined from conductivity using a Guildline 8400B salinometer, after at least 24 h of acclimatization at the laboratory temperature, which was kept constant at 1◦ -2◦C below the salinometer bath temperature (24◦C). The instrument's accuracy was confirmed by the daily measurement of a certified standard (IAPSO Seawater Standard). The salinity analyses were conducted within maximum 6 weeks of collection, and the analytical precision was 0.0003 (on the practical salinity scale).

#### Carbon Budget

We used a budgetary approach to estimate how TIC in the initial, pre-ice seawater was distributed between the air, water, and ice by the end of each experiment. For this purpose, the absolute TIC content (i.e., moles of carbon) of the formed ice (absTICice) and drained brine (absTICdrain; calculated from differences in TIC concentration between pre-ice and under-ice water) was compared to the absolute TIC content (absTICsw, ini) of the volume of seawater which was converted into ice. Any differences between these three values were assumed to originate from exchange with air:

$$absTIC\_{sw, \, ini} = \, \, \, absTIC\_{ict} + \, absTIC\_{drain} + \, absTIC\_{air}$$

The net TIC exchanged with the air over the course of the experiment (absTICair) was thereby not determined directly and could be either positive or negative, as carbon could be released to or taken up from the air, respectively.

To calculate the absolute TIC content of water that eventually froze and the resulting ice, the TIC content measured for the corresponding water (TICsw, ini) and ice melt samples (TICice), both in µmol/kg, was multiplied by the final ice mass (mice), so as to consider only the pre-ice seawater from which the ice had formed. To account for TIC drained from the ice to the under-ice seawater, the difference between the TIC content in pre-ice and under-ice water samples (TICsw, fin − TICsw, ini) was multiplied by the total mass of the under-ice water, which was determined by subtracting the final ice mass from the initial water mass (msw.ini).

$$T\text{TC}\_{\text{sw},ini} \* m\_{\text{icc}} = T\text{TC}\_{\text{icc}} \* m\_{\text{icc}} + (T\text{TC}\_{\text{sw},fin} - T\text{TC}\_{\text{sw},ini}) \* 1$$

$$(m\_{\text{sw},ini} - m\_{\text{icc}}) + absT\text{TC}\_{\text{air}}$$

The initial water mass was calculated from seawater density, determined from salinity and in-situ temperature according to Millero and Poisson (1981), and the initial water volume based on water depth after sampling and tank dimensions, corrected for estimated sensor volumes (König, 2017). Due to rather large uncertainties regarding ice volume and density, the ice mass was derived from a salinity budget set up in a similar way as the above outlined carbon budget. For this purpose, we replaced the TIC concentrations with the salinities measured in water and ice melt samples in the budget, assuming these practical salinity values represent the mass concentration of salt in water. As we assumed the exchange of salt with air during ice growth to be negligible (thereby neglecting evaporation and sublimation), we set the air fraction of the budget to zero, and could thus determine the ice mass.

The overall uncertainties in the carbon budgets were propagated from analytical errors in TIC and salinity and from estimated errors in the tank volume, seawater density, and sensor volumes. Details of the error analysis are given by König (2017).

#### Sensors

The sensors installed in the seawater tank with relevance for this study are described in **Table 2** and shown in **Figure 1B**. The three temperature probes were installed at different depths to confirm that the tank was well-mixed, and the majority of the other sensors were installed in the middle of the tank, with an extra pH sensor at the bottom. Dopbeam sonars were deployed at the bottom of the tank. Additional details about sensor configurations during each experiment, including additional sensors installed, are given by König (2017).

#### Sensor CO2 System

The CONTROS HydroC pCO2 probe (Fietzek et al., 2014) we deployed was connected to an external pump (SBE 5M, Sea-Bird TABLE 2 | Overview of sensors deployed in the seawater tank. Additional sensors installed but not discussed in this paper are described in König (2017).

TABLE 3 | Average ice growth rates of all six runs, calculated from the runtime (duration) of each experiment and the average final ice thickness.

Scientific) for all experiments except Run 8, when the connection cable to the pump was broken, giving a longer sensor response time for that experiment. The sensor was programmed to zero every 12 h (24 h in the case of Run 8) (Fietzek et al., 2014), and pCO2 values calculated from pH and TIC of the pre-ice water samples were used to correct for a (generally small) offset in the pCO2 curve (König, 2017). For pH measurements, we deployed two SeaFET sensors (Martz et al., 2010). As for pCO2, the pH data were calibrated using the measurements of discrete water samples. Further details about signal correction and data processing are given by König (2017).

#### Other Parameters

Three RBR TR 1050 probes were installed in the tank at different depths to measure water temperatures. These data were used to check for stratification, correct the SeaFET pH data, and for the sensor CO2 system calculations. In addition, we obtained temperature data from an O2 optode (Aanderaa Data Instruments AS, model 4330) deployed in the tank to monitor potential contaminating biological activity.

Seawater salinity was measured by a conductivity sensor (Aanderaa Data Instruments AS, model 4319). This inductive measurement was likely influenced by the sensors and tank walls around it, and therefore, we applied an offset to the data using the discrete water samples (König, 2017).

A pressure sensor (RBR Solo Pressure) was installed in the tank to monitor the potential pressure increase in the water during ice formation and for the CO2 system calculations. This sensor was calibrated against ambient pressure, as described by König (2017).

#### Ice Thickness and Brine Release (Dopbeam Sonars)

Three 20 MHz single-beam Dopbeam sonars (Sontek) were deployed in the tank, facing upwards, to survey the growing ice. These sonars record the backscatter of the emitted sonar signal, as well as a velocity profile (Vagle et al., 2012) and were primarily installed to test whether they could provide information on the timing and extent of brine released from the growing ice layer, but we also used them to measure the thickness of the growing ice. This was possible thanks to the density differences between the high-salinity brine, the growing ice, and the underlying water layer, which lead to the reflection of the emitted sonar signals at these interfaces. Further information on the backscatter data collection is given by König (2017).

## CO2 System

The aqueous CO2 system in the tank was fully characterized using water sample data and the program of van Heuven et al. (2011), with temperature and pressure from sensors deployed in the tank, the borate-to-salinity ratio of Uppström (1974), and KSO4 dissociation constants of Dickson (1990). Total silicate and phosphate were assumed to be negligible for the calculations. We used the carbonic acid dissociation constants of Cai and Wang (1998), because they gave to the smallest discrepancies between calculated and measured values of the over-determined carbonate system (König, 2017).

We found that after the initial sampling point (before the experiments were started), calculated and measured alkalinity in the water diverged. We suspect that alkalinity was contaminated by an unknown process (possibly electrochemical) associated with the complex configuration of sensors (some of which were new) in the relatively small-volume tank during the initial trial runs. The relative contamination decreased when more seawater was added to the tank. As this AT increase occurred in a single step at the beginning of our experiments, the relative changes in AT during each experiment are still useful.

Carbonate system calculations were also conducted with sensor pH (SeaFET 245) and pCO2 data (HydroC), which allowed us to compute high-resolution TIC and AT time series. These calculations were only implemented for the stirred runs, since for the unstirred experiments the water was likely stratified, and therefore the sensors did not evaluate the same water.

## RESULTS AND DISCUSSION

## Average Growth Rates

Except at the coldest temperature, the average growth rates were lower for the stirred than the unstirred runs (**Table 3**), and the largest difference was observed between the slow growth experiments at −15◦C. These findings are not surprising, as the stratified water during unstirred runs allowed for quicker freezing due to decreased heat transfer from below. Ice growth during Runs 6 and 11 (coldest conditions, −40◦C), appeared to be faster with stirring, but we think this is the result of a delayed start of the main condenser units during the unstirred experiment (Run 11; **Table 1**), which led to higher temperatures (−30◦C) during the first hours of the experiment.

## Fate of TIC, AT, and Salinity During Ice Growth

#### Water and Ice Melt Samples

The relationship between salinity, AT, and TIC for all water and ice melt samples is illustrated in **Figures 2**, **3**. We included AT in this analysis, despite the potential contamination (section CO2 system), since this contamination appeared to be consistent for all discussed runs, and therefore the change in the AT signal in ice and water can still help identify potential ikaite precipitation in the ice.

The excellent fit of the linear regression to the data points in **Figure 2** shows that both TIC and AT are closely coupled with salinity. This indicates that the TIC (and AT) allocation in the airice-water system is strongly dominated by brine processes, and that gas exchange or precipitation processes within the ice are relatively minor. However, closer examination reveals additional differences between the different sample types and experiments (**Figure 3**). For this comparison, ratios between the variables (e.g., TIC to salinity) were calculated for each data point and subsequently normalized to the ratio of the pre-ice water sample from the same experiment, accounting for dilution (in the case of ice melts) and evaporation effects. For all experiments, changes in the parameter ratios were much more distinct for the ice melt samples than for the under-ice water. Even between ice cores from the same experiment, large discrepancies could be observed, particularly in the stirred runs. The TIC-to-salinity ratio in ice appeared to decrease slightly in most experiments (up to 5% decrease compared to pre-ice water), while AT seemed to generally increase in the ice, relative to salinity and TIC (increases up to 4 and 9%, respectively). The decoupling of the TIC and salinity signals implies a release of CO2 from the ice, either to

the air above or water below, but through a different mechanism from the salinity release. No significant increase in TIC relative to salinity was identified in the under-ice water, suggesting release to the atmosphere. However, the change in the under-ice water TIC concentration relative to salinity would have been very small, and might have not been detectable within our sampling and analytical uncertainties.

than the symbols. Analytical data are given in Table S2.

The detected increase in ice melt AT compared to both TIC and salinity implies precipitation and retention of CaCO3 minerals (e.g., ikaite) in some of the sampled ice. Since such minerals dissolve during the melting of the ice, they contribute twice as much AT as TIC to the resulting solution. The observed relative decrease in AT in some of the under-ice water samples might also be related to this process, as it leads to an increase in TIC (relative to AT) in the remaining brine, which can then influence the underlying water by drainage (Jones and Coote, 1981). When examining the melting ice cores, white particles, which could have been CaCO3, were visible for all cores from all experiments. No quantitative analysis of these particles was conducted, and we cannot confirm whether they were indeed CaCO3 and not some other salt precipitates or contamination. Even though ikaite may have precipitated during all our experiments, its impact on the carbon distribution was small, given the small relative increase in the AT to TIC ratio in the ice melt samples (mostly less than 5%, apart from one sample with a 10% increase).

Comparing the different experiments, it is difficult to find a clear pattern between the deviations in ice carbonate chemistry and the ice growth rates (**Table 3**). Nevertheless, it appears as if the AT to TIC ratio in ice is higher for experiments with slower ice growth (**Figure 3C**), indicating that CaCO3 precipitation may be enhanced at slower growth rates. In particular, the experiment with (by far) the slowest ice growth rate displayed an increase

in the AT to TIC ratio in the ice melts of up to 9.4%. This is contrary to expectations of increased ikaite precipitation at lower temperatures/higher brine salinities, but might be caused by the insufficient removal of dissolved CO2 from the ice during faster ice growth, which is required to achieve ikaite supersaturation (Papadimitriou et al., 2013). The apparent increase in CaCO3 precipitation at slower ice growth could therefore be related to the increased permeability at these higher temperatures, possibly leading to more efficient removal of CO2 by outgassing, or simply due to the longer duration of the run, which allowed more CO2 to escape. This is further supported by pCO2 values measured

Uncertainties were calculated from deviations between replicates.

in the ice during this run, which were lower than for other runs (König, 2017). Conversely, apparently smaller ikaite precipitation in ice samples in the other runs may simply be due to their shorter runtimes, which may not have been sufficient to allow substantial ikaite precipitation. Papadimitriou et al. (2014) reported that 1–3 weeks, depending on ice temperature, is needed for ikaite precipitation to reach steady state in sea ice.

Despite the generally high AT in our seawater due to contamination, the absolute AT to TIC ratios in the ice melt samples of our experiments (average of 1.09 ± 0.02, with a maximum of 1.16) are generally at the lower end of ratios measured in the field, where values of up to 2.2 have been reported (Rysgaard et al., 2007, 2009). However, our values are comparable to those measured in young or newly formed natural and artificial ice (Rysgaard et al., 2007, 2009; Fransson et al., 2015; Geilfus et al., 2016), which further highlights the importance of incubation time for ikaite precipitation. Some ice melts from the two colder, stirred runs (Run 5 and 6) show an inverse pattern from that of the warm run (Run 8), i.e., a relative increase in ice TIC combined with a decrease in AT. This might have been caused by either increased retention of CO2 in the ice (e.g., in the form of bubbles), or unidentified deficiencies in sampling and/or analyses. Transport of ikaite from the ice to the underlying water, as suggested by Geilfus et al. (2016), would also explain the relative decrease in ice AT compared to TIC. However, as no corresponding increase in underwater AT was detected, and since the permeability in ice was likely too low at our experimental temperatures, we believe this process is not responsible for the relative decrease in ice AT we detected.

#### Time Series Variations

While the water samples only provide point measurements at the beginning and end of the experiments, DIC and AT calculated from sensor pH and pCO2 allow us to examine changes in the water over the entire experimental period and to compare their development to those of salinity and ice thickness (**Figure 4**). Comparing TIC, AT, salinity, and ice thickness (**Figures 4A,B**), the parameter curves appear to be very similar. At the fastest ice growth rate (Run 6; **Figure 4**, left column), all four parameters increase nearly linearly over the whole experiment, while for the two experiments at slower ice growth rates, the increase in all three parameters was initially faster and then slowed over the course of the run. It is worth noting that for salinity, the increase was somewhat delayed for all three experiments compared to ice thickness. This was not the case for the TIC and AT, which appeared to increase even before the surface was ice covered. While an increase in TIC during this period is conceivable due to increased uptake of CO2 from the air, and has indeed been observed during a similar experiment by Kotovitch et al. (2016), a significant increase in AT is unrealistic. Furthermore, no such increase was observed for the pre-ice water samples, and it is therefore likely an artifact of the CO2 system calculations, coupled with the different response times of the sensors. We were only able to estimate the response time for the HydroC sensor and calculations with time-lag corrected sensor data did not provide useful insight (König, 2017).

FIGURE 4 | Temporal development of water TIC, AT (A), water salinity, and ice depth (B; fitted curves from the three Dopbeam sonars) during the three stirred experimental runs as well as the corresponding water sample measurements. TIC and AT were calculated from the HydroC pCO2 and the (corrected) external pH of SeaFET 245. Water AT increased before the experiments started, possibly due to some electrochemical process associated with the sensors and equipment in the tank at initial deployment; note that seawater was added to the tank between runs 6 and 8 (Table 1), reducing the relative impact of the alkalinity contaminant. (C) Ice growth rate (cm/d) derived from Dopbeam-measured ice depth and daily rates of change in salinity (dS/S per day) and TIC (dTIC/TIC per day). Note different vertical and horizontal scales. Black vertical lines indicate water sampling.

To better compare the temporal changes in TIC, salinity, and ice thickness, rates of change were calculated over time intervals adapted to the duration of the run (**Figure 4C**): 1.5 h for Run 5, 45 min for Run 6, 6 h for Run 8. For TIC and salinity, normalized change rates were calculated, e.g., dS = dS/S, using the data curves, which were averaged over the corresponding time interval for this purpose. As the calculated TIC curve was still very noisy, it was further averaged over four adjacent data points. For ice thickness, the change rate (i.e., growth rate) was derived from curve fits to the Dopbeam sonar data.

The ice growth rates decreased significantly over the course of each experiment, due to the insulating effect of the overlying ice (**Figure 4C**), by more than a factor of five for the slower two experiments and by half for the fastest experiment. Salinity and TIC appeared to increase more steadily over the course of the experiments. This implies a decoupling between ice growth and desalination rates, which can also be observed by comparing the timing of the peak desalination rate, which occurred later than the fastest ice growth rate in all three runs. The discrepancy between peak ice growth and desalination was largest for the fastest ice growth experiment, where the maximum increase in salinity was observed when the ice thickness was 5–6 cm. In the slower two experiments, salinity increased most at ice thicknesses of 0.5–2 cm (Run 8) and 1.5–3 cm (Run 5). These findings are consistent with mushy layer theory, which predicts a delayed onset of brine drainage (Wettlaufer et al., 1997; Notz and Worster, 2009). However, while the peak in the desalination rate appears to be delayed, salinity did begin increasing as soon as the first ice had formed, thereby contradicting the mushy layer predictions. This instantaneous increase could be related to the constant mixing of the tank during these three experiments, which could have disrupted the ice-water boundary layer, as described by Loose et al. (2011).

The rate of change in TIC is harder to interpret, due to the aforementioned difficulties with the sensor CO2 system calculations, which led to unrealistic increases in AT and TIC at the beginning of the experiments. Focusing only on the ice growth period (i.e., the time after the last pre-ice water sample was taken), it appears as if TIC was increasing more rapidly than salinity during initial ice growth, but more slowly after that. While this could be an artifact of CO2 system calculations, it might also be due to a dissimilar behavior of the CO2 gas and the dissolved salts, as was observed for conservative gas tracers in similar experiments by Loose et al. (2009), who suggested that the differences were related to the faster diffusion of gas away from the ice-water interface. In our case, this process may have been further supported by a slight, additional cooling of the under-ice water (i.e., from approximately −1 ◦ to −2 ◦C; König, 2017) and the concomitant increase in gas solubility.

## Brine Drainage at Different Growth Rates

The data collected during the experiments without stirring in the tank provide insight into brine release processes, and the temperature dependent differences we observed have potentially important consequences with regard to increasing air temperatures in polar regions. These findings are mainly based on the backscatter data obtained from the Dopbeam sonars, which indeed seemed able to detect brine release from the ice and were further supported by changes in other parameters.

As indicated by the rates of change in the salinity signal (**Figure 4C**), the maximum desalination from the ice occurred in the early stages of each (stirred) run, but after the maximum ice growth rate had been reached. Similarly, the backscatter signal from the Dopbeam sensors shows a relatively undisturbed, clearly defined ice-water interface up until an ice thickness of about 5– 10 cm, after which brine draining out of the ice starts to cloud the strong backscatter signal from the ice, indicating periods of relatively strong brine release (**Figure 5A**). This delayed onset of the maximum brine drainage is in agreement with predictions of the mushy-layer theory, as discussed above.

The blurring of the ice-water interface due to brine drainage is more pronounced for the two coldest experiments, either due to more intense outflow of brine or higher density differences between brine and underlying water, causing a more efficient reflection of the sonar pulses. The notion that brine released during the colder experiments is denser than for the warmer run is further supported by areas of intense backscatter in the water further away from the ice-water interface, as well as temperature measurements in the water column.

While for the warmest unstirred experiment (Run 12), the relatively weak backscatter signal attributed to brine drainage is restricted to the uppermost centimeters of the water column, we find that for the colder two runs, the sonar signal is also scattered back at deeper depths. Furthermore, the temperature signals of sensors installed at three different depths within the water column suggest that during colder runs, the tank was more stratified (**Figure 5B**). For the colder two runs, abrupt changes in the temperature curve of the bottom sensor coincide with the apparent brine drainage patterns we observed in the backscatter data, suggesting that during the intense brine release period at about 5–10 cm ice thickness, some plumes might have been sufficiently dense to reach the bottom of the tank. No such changes were observed during the warmest run, implying a lack of depth-penetrating brine plumes during this experiment. This finding has potentially far-reaching implications for carbon

FIGURE 6 | Carbon budgets for stirred runs (A) and unstirred runs (B). Bars show how much of the TIC initially in the water ultimately converted into ice ends up in the final under-ice water, ice, and air, respectively (Section Carbon Budget). Error bars were computed from uncertainties in water and ice TIC, initial water mass, and ice mass (propagated from salinity measurements).

transport to lower water layers with the so-called "brine pump" (section Implications for future sea ice carbon export, below).

However, as no salinity sensor was deployed at depth, the evidence for more effective vertical brine transport at faster ice growth rates (or colder air temperatures) is circumstantial. Nevertheless, given the potential impact on carbon sequestration of this brine transport, we suggest that further experiments should be conducted with additional probes (as well as in deeper tanks) to resolve the salinity profile across the water column.

## Carbon Budget

The results of the carbon budget calculations are shown in **Figure 6**. For both sets of runs, the majority of TIC (53–72%) was released to the under-ice water, whereas the rest (28–47%) remained in the ice, and exchange with air was negligible. This is in agreement with previous mass balance studies, which found

FIGURE 7 | Correlations between the TIC fraction in water (blue) or ice (red) and the average ice growth rate (A) or the ambient temperature (B) of each experiment. Values from stirred runs are marked with round, those from unstirred runs with triangular, symbols. Values measured by previous studies are shown in black, but have not been used for the regression. Black diamonds show values presented by Rysgaard et al. (2007), squares are the (average) results of Nomura et al. (2006). The air temperatures shown in panel B are average values measured during each run, as listed in Table S1.

that of all TIC released from the ice, only a small fraction escaped to the air (<0.1 to 1.1%), while the majority (89.9 to >99.9%) was transported to the underlying water (Nomura et al., 2006; Rysgaard et al., 2007; Sejr et al., 2011).The uncertainty was larger for the air than for water and ice fractions, because the air fraction was estimated from the results of the other fractions. In general, uncertainties in the carbon budget were dominated by the heterogeneity of the ice, which lead to uncertainties in ice TIC and salinity.

Both stirred and unstirred runs show very similar TIC redistribution patterns: in the fastest ice growth experiments, the ice phase contains nearly half of the initial TIC (47 and 46% for stirred and unstirred runs, respectively) but only less than a third in the slowest growth experiments (28 and 29%, respectively). This is in agreement with expectations that higher growth rates lead to less efficient desalination of the ice, and that TIC largely follows salinity. Carbon budget numbers reported in previous laboratory studies for similar growth rates appear to be in line with the values observed for our experiments (**Figure 7**).

Interestingly, the TIC redistribution appears to be less dependent on the growth rate than the ambient temperature at which the experiments were run. The slowest experiment (Run 8), for instance, shows a very similar distribution to Run 12, which was run at the same temperature, despite a more than three times lower growth rate. Thus, there is a better (linear) correlation between TIC fractions in water and ice with ambient temperature than with growth rate (**Figure 7**). Perhaps this reflects the importance of a temperature gradient within the ice in retaining brine and gases.

## Implications for Future Sea Ice Carbon Export

Due to large uncertainties in future sea-ice extent and volume (e.g., Meier, 2016; Stammerjohn and Maksym, 2016), meaningful estimates of the total future carbon export from sea ice are difficult. However, the implications of our findings for future potential sea-ice conditions are summarized in **Figure 8**. In the baseline, "present day" case the majority of TIC released from the sea ice remains in the surface layer, with only a small fraction exported to depth, as predicted in the modeling studies of Grimm et al. (2016) and Moreau et al. (2016).

In cases where warmer air temperatures result in a solid ice cover that is warmer and thinner the expected increase in TIC release (per unit of ice) is canceled out by the decrease in total ice production. Furthermore, export of TIC to the deep ocean is reduced due to the less sufficiently dense brine.

On the other hand, a thinner ice cover also leads to a more mobile, broken ice cover during the winter, with more leads and possibly polynyas. In this case, sea-ice and brine production may actually increase. In addition, that brine and the associated TIC may be exported deeper into the water column, thereby enhancing deepwater formation and carbon sequestration.

It is important to note that these simple future scenarios neglect crucial impacts such as changes to ocean circulation, or the related biological systems. However, they again stress the importance of open water within the sea ice cover, which may not only enhance gas exchange (Loose et al., 2009; Else et al., 2011), but also facilitate export of carbon out of the surface layer (Grimm et al., 2016; this study).

## CONCLUSION AND PERSPECTIVES

Sea ice growth experiments conducted in a cold laboratory showed that there is a connection between ice growth rates and the redistribution of TIC between the ice and the water underneath. This differential redistribution was evident from discrete water and ice samples, and also visible in continuous time series recorded by underwater sensors. A budgetary approach showed that at higher ice growth rates, relatively more of the initial TIC remains in the ice phase than at slower growth rates (46–47% compared to 28–29%), and the TIC redistribution appeared to depend more on the ambient temperature at which the ice formed than on the ice growth rate. For all experiments, most TIC was rejected into the under-ice water (53–72%, depending on the temperature), while exchange with the air was insignificant, relative to our experimental uncertainties.

The TIC redistribution appeared to be dominated by brine drainage processes, as indicated by a good correlation between TIC and salinity in the ice and water. Minor deviations between the two parameters were likely due to CaCO3 precipitation in ice and/or CO2 exchange with the air. Temperature and Doppler backscatter data suggest that whereas more TIC is exported from the ice at higher temperatures (slower ice formation rates), the vertical extent of that transport (e.g., to the bottom of the mixed layer) might be more efficient when ice is formed at colder temperatures (at a faster growth rate). For a more satisfactory assessment of the potential for TIC export to deeper water layers, the stratification and brine release processes should be examined more closely, ideally by the installation of additional salinity probes in the water and possibly the ice.

## DATA AVAILABILITY

The data used in this laboratory-based study are included in this manuscript and König (2017).

## AUTHOR CONTRIBUTIONS

DK designed and conducted the experiments with substantial guidance from LM and support from KS. SV provided, installed,

## REFERENCES

and processed the data of the Dopbeam sensors. DK conducted chemical analyses, processed all datasets, and drafted the manuscript. All authors contributed to manuscript editing and revision.

## FUNDING

This project was supported by the Fisheries and Oceans Canada Arctic Sciences Fund. Further financial support was provided by the Zeno Karl Schindler Foundation in the form of a master thesis grant, and by the ETH Zurich Financial Aid Office, which granted a personal travel allowance to DK.

## ACKNOWLEDGMENTS

We would like to thank the staff of the Institute of Ocean Sciences for their support, especially Kenny Scozzafava and Jasmine Wietzke for salinity analyses, Marty Davelaar for help with the TIC and alkalinity analyses, and Sarah Zimmermann and Jonathan Velarde for help with ice coring. Further thanks are due to Nicolas Gruber for his support during the master's thesis this paper is based on, and to him and his group, as well as the BEPSII community for all the fruitful discussions during the preparation of this manuscript.

## SUPPLEMENTARY MATERIAL

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

**Conflict of Interest Statement:** 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.

Copyright © 2018 König, Miller, Simpson and Vagle. This is an open-access article distributed under the terms of the Creative Commons Attribution License (CC BY). The use, distribution or reproduction in other forums is permitted, provided the original author(s) and the copyright owner(s) are credited and that the original publication in this journal is cited, in accordance with accepted academic practice. No use, distribution or reproduction is permitted which does not comply with these terms.

# Application of Sentinel-2 MSI in Arctic Research: Evaluating the Performance of Atmospheric Correction Approaches Over Arctic Sea Ice

#### Marcel König1 \*, Martin Hieronymi 2 and Natascha Oppelt 1

1 Earth Observation and Modelling, Department of Geography, Kiel University, Kiel, Germany, 2 Institute of Coastal Research, Helmholtz-Zentrum Geesthacht, Geesthacht, Germany

#### Edited by:

Benjamin Allen Lange, Freshwater Institute, Fisheries and Oceans Canada, Canada

#### Reviewed by:

Mallik Mahmud, University of Calgary, Canada Shujie Wang, Lamont Doherty Earth Observatory (LDEO), United States John Alec Casey, York University, Canada

> \*Correspondence: Marcel König<EMAIL_ADDRESS>
#### Specialty section:

This article was submitted to Cryospheric Sciences, a section of the journal Frontiers in Earth Science

Received: 27 July 2018 Accepted: 04 February 2019 Published: 22 February 2019

#### Citation:

König M, Hieronymi M and Oppelt N (2019) Application of Sentinel-2 MSI in Arctic Research: Evaluating the Performance of Atmospheric Correction Approaches Over Arctic Sea Ice. Front. Earth Sci. 7:22. doi: 10.3389/feart.2019.00022 Multispectral remote sensing may be a powerful tool for areal retrieval of biogeophysical parameters in the Arctic sea ice. The MultiSpectral Instrument on board the Sentinel-2 (S-2) satellites of the European Space Agency offers new possibilities for Arctic research; S-2A and S-2B provide 13 spectral bands between 443 and 2,202 nm and spatial resolutions between 10 and 60 m, which may enable the monitoring of large areas of Arctic sea ice. For an accurate retrieval of parameters such as surface albedo, the elimination of atmospheric influences in the data is essential. We therefore provide an evaluation of five currently available atmospheric correction processors for S-2 (ACOLITE, ATCOR, iCOR, Polymer, and Sen2Cor). We evaluate the results of the different processors using in situ spectral measurements of ice and snow and open water gathered north of Svalbard during RV Polarstern cruise PS106.1 in summer 2017. We used spectral shapes to assess performance for ice and snow surfaces. For open water, we additionally evaluated intensities. ACOLITE, ATCOR, and iCOR performed well over sea ice and Polymer generated the best results over open water. ATCOR, iCOR and Sen2Cor failed in the image-based retrieval of atmospheric parameters (aerosol optical thickness, water vapor). ACOLITE estimated AOT within the uncertainty range of AERONET measurements. Parameterization based on external data, therefore, was necessary to obtain reliable results. To illustrate consequences of processor selection on secondary products we computed average surface reflectance of six bands and normalized difference melt index (NDMI) on an image subset. Medians of average reflectance and NDMI range from 0.80–0.97 to 0.12–0.18 while medians for TOA are 0.75 and 0.06, respectively.

Keywords: remote sensing, atmospheric correction, sea ice, Arctic, Sentinel-2

## INTRODUCTION

The limited accessibility of Arctic sea ice makes satellite remote sensing an important tool for synoptic observations in this region. The retrieval of parameters, such as the spectral albedo of sea ice or snow and properties of melt ponds and leads, relies on spectral information; the retrieval of optical properties of target surfaces and waterbodies therefore is necessary. Optical sensors

such as the Advanced Visible High Resolution Radiometer (AVHRR; e.g., Huck et al. (2007)) and the MODerate Resolution Imaging Spectroradiometer (MODIS; e.g., Rösel et al. (2012), Tschudi et al. (2008)) have long been used in Arctic research (Pope et al., 2014; Nasonova et al., 2017) providing observations of different parameters such as sea ice extent, sea ice thickness and albedo. The Sentinel-3 satellites (S-3) from the European Space Agency's (ESA) Copernicus program, carrying the Ocean and Land Color Instrument (OLCI), continue this tradition (Donlon et al., 2012; Malenovský et al., 2012). Recently, Istomina and Heygster (2017) published a retrieval algorithm for sea ice albedo and melt pond fraction from S-3/OLCI observations. Despite a high temporal and spatial coverage of the abovementioned sensors (swath widths: AVHRR = 2,900 km, MODIS = 2,330 km, S-3/OLCI = 1,270 km), their coarse spatial resolution (AVHRR: 1.1–4 km, MODIS: 250– 1,000 m, S-3/OLCI: 300 m) impedes detailed observations of sea ice features such as melt ponds and ridges, which exist at spatial scales of meters to tens of meters (**Figure 1**). Studies based on optical sensors with a higher spatial resolution (≤30 m),

however, are rare. Markus et al. (2002) demonstrated the potential of Landsat 7 ETM+ (L-7; 30 m) for the classification of summertime sea ice surface conditions and retrieved the spatial distribution of ponded/unponded ice and open water (Markus et al., 2003), while Landy et al. (2014) used this sensor to observe sea ice development. Rösel (2013) underlined that high resolution satellite data enable determining melting features on Arctic sea ice and provide a basis for the comparison with coarse resolution satellite data. The ETM+ 8-bit sensor, however, showed saturation problems in the contrast-rich Arctic environment (Bindschadler et al., 2008; Rösel, 2013).

ESA's Sentinel-2 mission (S-2) offers new opportunities for optical sensors in the Landsat-like spatial domain. S-2 consists of the twin satellites S-2A and S-2B, launched in 2015 and 2017, respectively. Both satellites are equipped with the MultiSpectral Instrument (MSI), which provides 13 spectral bands with a 12 bit radiometric resolution in the wavelength region from 443 to 2,202 nm. Depending on the band, spatial resolution varies between 10, 20, and 60 m (European Space Agency, 2015, 2018a). **Figure 2** illustrates the spatial and spectral settings of S-2A. Compared to L-7, S-2 has a larger swath width (290 km) and geographical coverage reaches up to 83◦ north. The two S-2 satellites are shifted by 180◦ on the same orbit, which results in a return period of 5 days. In higher latitudes, overlapping swaths further increase the return period to up to one image per day (European Space Agency, 2015, 2018b; Copernicus EO Support, personal communication). Due to the terrestrial focus of the S-2 mission, geographical coverage in the Arctic Ocean is currently limited to areas around islands but areas of interest may be added to the Mission baseline "if sufficient resources are identified" (European Space Agency, 2015). Malenovský et al. (2012) mentioned the potential of S-2 for albedo retrieval, derivation of snow properties and mapping of polynyas and leads; nevertheless, they also highlighted the demand for "precise corrections for atmospheric propagation, topography and directional reflectance behavior" (Malenovský et al., 2012), pointing to the importance of atmospheric correction (AC) of optical remote sensing data. Zege et al. (2015) also emphasized the importance of AC for the retrieval of melt pond fraction and surface albedo.

Radiative transfer through the atmosphere is strongly influenced by Rayleigh and Mie scattering as well as aerosol and gas absorption. This combination is particularly challenging over water areas, where the atmospheric path radiance is typically >85% of the total signal in oceanic waters, >60% in sediment-rich waters, and >94% in very dark waters (IOCCG, 2010). Low sun-zenith angles (>70◦ ) and correspondingly long atmospheric paths further aggravate the processing of data at the beginning and the end of summer in high latitudes (IOCCG, 2015). An adequate AC is therefore essential for the derivation of physical surface properties and follow-up multi-temporal analyses. Another challenge resulting from atmospheric transfer and the scattering of radiation is the neighboring or adjacency effect, i.e., the scattering of light from neighboring surfaces into the sensors field of view resulting in information overlay (Sterckx et al., 2015b). Adjacency effects depend on the brightness contrast between a target pixel and its neighborhood (Richter et al., 2006). The Arctic sea ice is a mix of snow, ice, melt ponds and ocean water, with typical broadband albedos in the visible wavelength region (VIS) ranging from >0.90 for fresh snow to <0.05–0.1 for ocean water (Perovich et al., 1998). Moreover, while being negligible for wavelength regions >1.5µm, the adjacency effect increases with shorter wavelengths (Richter and Schläpfer, 2017), i.e., the wavelength region used for water remote sensing is heavily affected. Due to the characteristic spatial mixture and high contrasts of deep, clear ocean waters, melt ponds, bright ice and snow surfaces in the sea ice, we expect the adjacency effect to have a large impact on S-2 imagery of Arctic sea ice.

Several AC processors are available for S-2, and some papers already compared processors for lakes (e.g., Dörnhöfer et al., 2016; Martins et al., 2017). Only recently, Doxani et al. (2018) performed an inter-comparison of AC processors for different surface types in temperate and tropical zones. Yet, to our knowledge, no study exists for Arctic sea ice; therefore, the objective of this paper is to evaluate the performance of AC processors for S-2 over different surface types in the Arctic sea ice. For this, we processed S-2A data applying five AC processors (ACOLITE, ATCOR, iCOR, Polymer, and Sen2Cor) and evaluated the results using in situ data from open water and ice floe surface measurements acquired during an RV

FIGURE 3 | Location of S-2 granule (A) and ice floe (B) in the Arctic Ocean on June 10, 2017.

Polarstern cruise in June 2017. We further investigate the influence of different AC processors on the retrieval of apparent optical properties, i.e., surface reflectance and remote sensing reflectance, which are used to retrieve parameters such as the spectral albedo or NDMI.

## METHODS

## Study Area

#### Cruise Overview

During cruise PS106.1 RV Polarstern was moored to an ice floe north of Svalbard between June 3 and 16, 2017 (**Figure 3**). At the beginning of the measurement period, snow and white ice covered most of the floe. Some few ponds had already formed at pressure ridges but we did not observe any pond formation on the flat parts of the floe; bare ice was exposed only very occasionally. Neither ponds nor bare ice were present near the measurement sites. Melting, however, started during the measurement period. The specific surface area of snow samples from the top 1– 2 cm has been analyzed with an IceCube (A2 Photonic Sensors, France) at different locations in the area of interest. Snow melt results in a decrease of the snow's specific surface area. **Figure 4** illustrates that melt occurred between June 4 and 10.

#### Atmospheric Parameters

To provide input data for the different AC processors, we used sun photometer measurements conducted on board of RV Polarstern and published in the Maritime Aerosol Network (AERONET-MAN, Holben et al., 1998). For later analysis, we computed the aerosol optical thickness (AOT) at 550 nm from the sun photometer measurements as:

$$\beta = \frac{AOT\_{\lambda}}{\left(\lambda / 550\right)^{-\alpha}}\tag{1}$$

with β being turbidity, i.e., AOT at 550 nm, AOTλ being the mean AOT at λ (500 nm) and α being the mean Angstrom exponent between 440 and 870 nm. We further used data from the weather station on board RV Polarstern (Schmithüsen, 2018). **Table 1** summarizes minimum, maximum, mean and standard deviation of the atmospheric parameters used to atmospherically correct the S-2 data acquired on June 10, 2017 (details on S-2 data are provided in section Sentinel-2 Data).

#### TABLE 1 | Atmospheric parameters.

## In Situ Data

To validate the results of the AC processors, we performed radiometric measurements of snow and ice covered areas (section Snow and Ice Measurements on the Floe) as well as open water (section Open-Water Measurements). Simultaneously to the radiometric open water measurements, we took water samples for water constituent analysis (section Water Sampling). All measurements have been localized with Global Positioning System (GPS) measurements (section Positioning of Sampling Points).

#### Radiometric Measurements

#### **Snow and ice measurements on the floe**

For spectral measurements on the floe, we used an ASD LabSpec5000 spectrometer (Analytical Spectral Devices Inc., United States), covering the wavelength region from 350 to 2,500 nm with a 1 nm spectral sampling rate. The ASD was equipped with an optical fiber with a 23◦ field of view (FOV) mounted to a pistol grip attached to a pole to avoid influences of the polar suits during the measurements (**Figure 5A**). We used a Labsphere Spectralon reflectance standard with 95% diffuse reflectance (Labsphere Inc., United States) as white reference. In addition, we used a set of Ocean Optics STS-VIS spectrometers (Ocean Optics Inc., United States). The setup consists of two synchronized spectrometers; one pointing downwards with a 1◦ FOV and the other, equipped with a diffusor, pointing upwards to track planar incoming radiation. Both instruments cover a wavelength region from 340 to 820 nm with a spectral sampling rate <1 nm. The spectrometers have also been attached to the end of a pole to avoid influences of the polar suits on the measurements (**Figure 5B**). For the Ocean Optics Instruments, we used a Labsphere Spectralon 99% diffuse reflectance standard (Labsphere Inc., United States) as white reference. For both instrument setups, we took reference measurements at intervals <10 min and when illumination conditions changed. Before conducting a white reference measurement with the ASD, we measured the Spectralon panel to control whether illumination has changed since the last reference measurement. If we detected changes, we performed a visual quality check on all measurements since the previous white reference and deleted spectra with unusual shapes. We corrected the Ocean Optics measurements for irradiance changes using the upwards pointing sensor. To account for the spatial heterogeneity within one S-2 pixel, we either walked around during the measurement process or performed measurements at different sites that we chose randomly. **Table 2** summarizes measurement dates, instruments and illumination conditions.

#### **Open-water measurements**

In open-sea water of a polynya, we conducted radiometric measurements from the deck of the vessel using three hyperspectral (380–950 nm) RAMSES sensors (TriOS GmbH, Germany). The sensors were arranged to measure downwelling irradiance (Ed), upwelling radiance just above the water surface (Lu; sensor viewing angle of 40◦ with 90◦ azimuth angle to the sun), and corresponding sky radiance (Lsky). We

TABLE 2 | Illumination conditions and instruments used.

ASD, Analytical Spectral Devices spectrometer; OO, Ocean Optics spectrometer.

calculated remote sensing reflectance, Rrs, according to the following equation:

$$R\_{rs}\left(40^{\circ}, 90^{\circ}\right) = \frac{L\_{\mu}\left(40^{\circ}, 90^{\circ}\right) - \rho \, L\_{\text{sky}}\left(40^{\circ}, 90^{\circ}\right)}{E\_d} \tag{2}$$

The sea surface reflectance factor, ρ, depends on the sun-viewing geometry and wind-dependent roughness of the surface. Wind was light to moderate, but the limited fetch reduced roughness (mean square slope of the waves) compared to open sea. We applied the usually recommended surface reflactance factor of Mobley (1999) to determine the mean reflectance, in this case ρ = 0.0272. Moreover, we applied four different surface reflectance factors between 0.0176 and 0.03 (Mobley, 2015; Hieronymi, 2016) to account for potential effects of reduced roughness and polarization. This methodical variability together with temporal changes of the single measurements were taken into account to determine the standard deviation of the in situ measurements.

We resampled all field spectra to S-2A bands using the spectral response functions and central wavelengths (European Space Agency, 2018a,c).

#### Water Sampling

We took water samples from the polynya simultaneously to the radiometric measurements. Water samples were filtrated and optically analyzed in the ship's laboratory. We determined the spectral absorption properties using a PSICAM (Röttgers and Doerffer, 2007), a QFT-ICAM (Röttgers et al., 2016), and a Liquid Waveguide Capillary Cell (World Precision Instruments, United States). During the S-2A acquisition on June 10, the −1.7◦C cold seawater had a salinity of ∼34 PSU. The particle concentration was very low, yielding a total particulate absorption coefficient at 440 nm, ap(440) of approximately 0.013 m−1 . The colored dissolved organic matter (CDOM) absorption at 440 nm, acdom(440), was 0.038 m−1 . The spectral shape of the particulate absorption indicates only to the presence of organic phytoplankton particles and no sediments; there is no significant particulate and CDOM absorption in the near infrared (>700 nm). The chlorophyll-a concentration (Chl-a) from the photometric measurements was 0.19 (±0.1) mg m−3 ; additional measurements with an AlgaeTorch (bbe Moldaenke GmbH, Germany) yielded Chl between 0.2 and 0.3 mg m−3 . These results confirm that the sampled water was very clear during the measurement period.

#### Positioning of Sampling Points

During the measurement period (**Table 2**) the ice floe drifted several tens of kilometers. To relocate all field measurements in a satellite image, we performed a drift-correction. Three stationary GPS devices constantly tracked floe movement in a 10-s interval. Additionally, we determined the locations of all sampling points with GPS measurements (F5521gw, Ericsson, Sweden; Galaxy S7, Samsung, South Korea; eTrex 10, Garmin, United States) in the World Geodetic System 1984 reference system. To locate all measurements from different dates and times in a single satellite image, we first applied a Savitzky-Golay filter to smooth the positioning data of each of the stationary GPS devices; then we used a second order polynomial to interpolate points in a 1-s interval. In this way, we could allocate every field measurement to three stationary GPS positions. Using this data, we computed the distances from every field measurement position to each of the three stationary GPS devices at the respective time of data acquisition (tin situ). To find the measurement locations on the floe at the time of satellite overpass (tsatellite) we fit the distances to the three stationary GPS devices at tsatellite to the distances computed at tin situ using the method of least squares. Finally, we performed a visual quality control.

#### Sentinel-2 Data

On June 10, 2017, clear sky conditions allowed acquisition of S-2A data of the floe at 14:58:01 UTC. For this study, we used the S-2 L1C top of atmosphere (TOA) reflectance product (processing baseline: 02.05; European Space Agency, 2015). For AC evaluation, we used S-2 bands which are in the spectral range of the field instruments, i.e. S-2A bands 2–7 (490–783 nm) for sea ice surfaces and bands 2–7 plus band 8A (490–865 nm) for the open water surface. All processors provide a 10 m output, except ATCOR. Which only offers a 20 m output. To provide spatially comparable results, we therefore downsampled 10 m output data to a spatial resolution of 20 m by 2 × 2 pixel block averaging.

#### Regions of Interest

In general, a point-by-pixel comparison is necessary for an accurate validation of remote sensing data via in situ measurements. With a steadily drifting floe, however, we are unable to assess precise accuracies of GPS measurements and trilateration; we therefore dismissed a point-based validation. Instead, we defined a region of interest (ROI) for each measurement Section by computing a circular buffer with a 20 m radius around all GPS points belonging to one measurement section. We defined all pixels intersecting the respective buffer area to be one ROI. **Table 3** lists the respective number of pixels and field spectra as well as the measurement device used for each ROI.

Regarding the RAMSES measurements, we manually defined one ROI in the open water area next to RV Polarstern and excluded pixels with unusual shape presumably influenced by floating ice. **Figure 6** illustrates locations of the respective ROIs.

#### Atmospheric Correction Processors

For our study, we selected common, publicly available AC algorithms. The level of user-controlled parametrization depends

TABLE 3 | Number of pixels, number of in situ spectra and measurement device.

OO, Ocean Optics; ASD, Analytical Spectral Devices spectrometer for each ROI (Figure 3 shows the locations of ROIs).

on the respective processor and is described in the following. **Table 4** summarizes the most relevant parameters.

ACOLITE is developed at the Royal Belgian Institute of Natural Sciences and includes AC processors specifically designed for aquatic applications, e.g., it addresses diffuse sky reflectance at the water surface. We used the most recent version that applies a new dark spectrum fitting AC by default (Vanhellemont, submitted). The algorithm consists of five steps (Vanhellemont and Ruddick, 2018): (1) Correction of atmospheric gas transmittance and sky reflectance in the TOA signal; (2) Construction of a dark spectrum using the darkest pixels in each band; (3) Estimation of AOT at 550 nm for different aerosol models (continental, maritime) by linearly interpolating path reflectance from look up tables (LUTs) generated with the Second Simulation of a Satellite Signal in the Solar Spectrum–Vector (6SV) radiative transfer model (Vermote et al., 1997); (4) Retaining the band resulting in the lowest AOT for each aerosol model; (5) AC with band and model combination resulting in overall lowest AOT. This approach is presumably more appropriate for ice surfaces than the black short-wave infrared approach of previous versions (Vanhellemont and Ruddick, 2015). Information on precipitable water, atmospheric pressure and ozone concentration are obtained from ancillary data. Adjacency effects are to some extent addressed in the aerosol correction procedure, but are not entirely extinguished (Vanhellemont and Ruddick, 2018). We used the default parameterization but changed the WV default value to 1.14 cm according to the AERONET data. The output of ACOLITE may be images of different reflectances such as surface reflectance or Rrs with a 10 m spatial resolution, which was downsampled afterwards.

ATCOR (Atmospheric/Topographic Correction for Satellite Imagery; Richter and Schläpfer, 2017) uses pre-calculated LUTs computed with the Moderate-Resolution Atmospheric Radiance and Transmittance Model–version 5 (MODTRAN 5) (Berk et al., 2008) and image-based retrievals of atmospheric properties. The LUTs cover four aerosol models (rural, urban, maritime, and desert) and six water vapor column contents (0.4–5.0 cm). Varying ozone concentrations are addressed via an extra ozone database. The dense dark vegetation method is used to derive AOT and the aerosol type can be estimated by comparison of derived path radiance with MODTRAN standard aerosol types. WV is estimated using the Atmospheric Pre-corrected Differential Absorption algorithm. The resulting unit is "bottom of atmosphere" (BOA) reflectance (Richter and Schläpfer, 2017). The interface of ATCOR guides the user through the parametrization of the processor, which enables to address emerging problems. The small number of dark pixels in the image hampered an image-based retrieval of aerosol type and AOT. A WV map could be generated, but values near the vessel were high (∼1.8 cm) compared to the AERONET data. To parametrize the processor in accordance to the AERONET data, we therefore selected the pre-defined maritime aerosol type and a stable WV of 1.0 cm, which is the closest selectable default value, compared to the AERONET data. Visibility as a measure for AOT was set to 65 km. ATCOR considers adjacency effects as averaged reflectance of the neighborhood of a pixel (Richter and Schläpfer, 2017); its impact therefore depends on the difference in reflectance of a pixel and its neighborhood and decreases exponentially with increasing distance. Thus, for recurrent spatial patterns, adjacency effects remain similar at different spatial scales; consequently, the absolute range of the adjacency correction is uncritical. Therefore, we selected a typical range of 1 km (Richter and Schläpfer, 2017). The number of adjacency zones remained unchanged (1). We further applied a cirrus correction and selected the 20 m spatial resolution output.

In this study, we also used iCOR version 0.1 in ESA's Sentinel Application Platform (SNAP; v5.0). iCOR, previously known as OPERA (Sterckx et al., 2015a), is an AC processor for land and water as it accounts for the non-Lambertian reflectance of water surfaces. The iCOR workflow comprises four steps (De Keukelaere et al., 2018): (1) Classification of land/water pixels; (2) AOT retrieval over land following the approach in Guanter (2006), and extension to adjacent water pixels assuming a spatially homogeneous atmosphere; (3) adjacency correction; (4) AC using pre-calculated MODTRAN 5 LUTs based on a rural aerosol model (De Keukelaere et al., 2018). In the current iCOR SNAP version, WV is fixed to 2.0 cm. Over water the adjacency effect is corrected using the SIMilarity Environment Correction (SIMEC) approach (Sterckx et al., 2015b) while over land the user defines a fixed range (Sterckx et al., 2015a; VITO, 2017). We applied the adjacency correction and increased the default adjacency window for land surfaces to three pixels (tests with larger values introduced a reflectance peak at 705 nm). Further, we adjusted the default AOT value to 0.02 according to the AERONET data. Besides that, we used the default parameter setting. iCOR returns BOA reflectances in the native resolution of the respective band.

The polynomial-based atmospheric correction algorithm Polymer (Hygeos France; Steinmetz et al., 2011) is an AC processor for water bodies applicable to multiple sensors, including MSI and OLCI. Polymer is a spectral matching algorithm, which decouples the reflectance signal of the water body from atmospheric and water surface reflectance. It makes

FIGURE 6 | ROI locations on the investigated floe. Background: S-2A stretched pseudo-true color RGB image (RGB: 665, 560, and 490 nm, spatial resolution: 20 m) acquired at 2017-06-10 14:58:01 UTC.

use of the observation that the sum of atmospheric path reflectance and adjacency effects can be approximated by a polynomial consisting of three terms that address (1) nonspectral scattering or reflection, (2) fine aerosol scattering and (3) molecular scattering and adjacency effects, and does not require an aerosol model (IOCCG, 2015). One of the strengths of Polymer is the possibility to retrieve Rrs in presence of sun glint, which leads to a higher yield of evaluable pixels compared to other algorithms. Recently, the Polymer algorithm has been modified to improve results in high latitudes (Steinmetz and Ramon, 2018). In the processing chain, non-water surfaces are masked, i.e., Polymer has only been applied to ROI 10A. For this study, the Polymer 10 m outputs were projected and resampled to the Sentinel-2 grid using the Nearest Neighbor approach, and spatially downsampled to 20 m resolution by 2 × 2 pixel block averaging.

Sen2Cor is the Sentinel-2 Level-2A Prototype processor for land surfaces. For this paper, we used the stand-alone version of Sen2Cor 2.4.0, which also allows for an image-based retrieval of atmospheric parameters. Water vapor is derived with the Atmospheric Pre-corrected Differential Absorption algorithm applied to S-2 bands 8A and 9. Sen2Cor estimates AOT by means of the dense dark vegetation method using the relation of the reflectance in bands 4, 2, and 12. The image-derived atmospheric parameters are then associated with pre-computed LUTs generated by a libRadtran based radiative transfer model. The LUTs include two atmospheric models (mid-latitude summer and winter), two aerosol types (maritime, rural), six or four different WV column values and six ozone concentrations (Müller-Wilm, 2016). We chose the maritime aerosol type and selected the mid-latitude winter model atmosphere as suggested in Harris Geospatial Solutions Inc (2018). We activated the variable visibility option and changed the default value to 65 km.

#### Performance Measures

We performed the spectral measurements on the floe within 5 days before and after the satellite acquisition. Illumination conditions therefore differed for field measurements and the satellite observation (see section Radiometric Measurements) which may affect measurement results (Malinka et al., 2016). Moreover, due to the different spectrometer setups described in section Snow and Ice Measurements on the Floe, ASD and Ocean Optics show varying sensitivities to changing illumination conditions. The different FOVs result in varying sensitivities to changing sensor-surface geometries. An inter-comparison of ASD and Ocean Optics under field conditions revealed mean differences in surface reflectance of about 0.1; a comparison of in situ spectra and BOA reflectances by means of classic performance metrics such as bias or root mean squared error therefore is unreasonable.

To identify systematic differences between the respective processors, i.e., the average tendency of the BOA reflectance to be larger or smaller than TOA reflectance, we computed the percentage bias (PBIAS) as

$$PBIAS = 100 \cdot \frac{\sum\_{i=1}^{n} \left( R^{AC} (\lambda\_i) - R^{TOA} (\lambda\_i) \right)}{\sum\_{i=1}^{n} \left( R^{TOA} (\lambda\_i) \right)} \tag{3}$$

with R AC and R TOA being the bottom and top of atmosphere reflectance at band λi , respectively; n is the number of bands.

Compared to the absolute values, spectral shapes are more stable, especially at the beginning of the melt season and for the visible wavelength regions (Perovich, 1994). We therefore evaluated the spectral shapes of ice and snow measurements on the floe, computing the slope between two adjacent bands and then calculate the slopes' Mean Absolute Error (slope MAE):

$$slope\ MAE = \frac{1}{n} \sum\_{i=1}^{n} \left| \mathbf{S}^{AC}(\lambda\_i, \lambda\_{i+1}) - \mathbf{S}^{in\ situ}(\lambda\_i, \lambda\_{i+1}) \right| \tag{4}$$

with S AC(λi , λi+1) being the slope of the BOA reflectance outputs of the respective AC processor and S in situ(λi , λi+1) being the slope of the resampled in situ reflectance between band λi and band λi+1, respectively; n is the number of band pairs. We further calculated the coefficient of determination (r²) as the square of Pearson's correlation coefficient to evaluate linear correlations between the respective AC BOA output and the in situ measurement:

$$r^2 = \left(\frac{\sum\_{l=1}^{\text{st}} \left(R^{\text{in situ}}\left(\lambda\_l\right) - \overline{R^{\text{in situ}}}\right) \cdot \left(R^{\text{AC}}\left(\lambda\_l\right) - \overline{R^{\text{AC}}}\right)}{\sqrt{\sum\_{l=1}^{\text{st}} \left(R^{\text{in situ}}\left(\lambda\_l\right) - \overline{R^{\text{in situ}}}\right)^2} \cdot \sqrt{\sum\_{l=1}^{\text{st}} \left(R^{\text{AC}}\left(\lambda\_l\right) - \overline{R^{\text{AC}}}\right)^2}\right) \tag{5}$$

with R AC being the output of the respective AC processor and R in situ being the reflectance of the resampled in situ spectrum at band λi (band 2–7); n is the number of bands.

For deep polynya water, we performed measurements of Rrs from the deck of the vessel parallel to the satellite overpass under clear-sky conditions. No white referencing was necessary due to the measurement setup, which already integrates changes in illumination conditions. We therefore compared resulting reflectance intensities. ACOLITE calculates Rrs, while ATCOR, iCOR, SEN2COR and Polymer yield surface reflectances. The surface reflectances were transformed into Rrs following (Mobley et al., 2018):

$$R\_{rs} = \frac{R^{AC}}{\pi} \tag{6}$$

with R AC being the surface reflectance output of the respective AC processor. We then computed MAE to account for differences in reflectance intensity according to:

$$MAE = \frac{1}{n} \sum\_{i=1}^{n} \left| R\_{rs}^{AC}(\lambda\_i) - R\_{rs}^{in \,\,site}(\lambda\_i) \right| \tag{7}$$

with R AC rs being Rrs of the respective AC processor and R in situ rs being Rrs of the resampled in situ spectrum at band λi (band 2–8A); n is the number of bands. To evaluate correspondence in spectral shapes, we computed slope MAE and r², replacing surface reflectance with Rrs and expanding S-2A band selection by band 8A.

### RESULTS

Due to the different measurement designs described above, we discuss the results for open water as well as snow and ice surfaces separately.

#### Snow and Ice Surfaces

**Figures 7A,C,E** illustrates the mean shape and intensity of S2- A TOA reflectances, BOA reflectances of the respective AC

processors and resampled in situ spectra from snow and ice surfaces. Due to the low variance, we computed the mean spectra of all ROIs on the floe. In situ spectra show the typical shapes of snow and ice with maximum reflectances in the blue-green wavelength region and a decrease toward longer wavelengths. Maximum mean in situ reflectances at band 2 (490 nm) range from approximately 0.8–1.1 with high standard deviations presumably resulting from a high spatial variability on sub-pixel scale. Values are in line with field measurements from other studies (e.g., Perovich, 1998; Goyens et al., 2018). In comparison to the TOA spectrum, BOA spectra for all AC processors show increased reflectances over the entire spectrum. For each ROI, Sen2Cor produces the highest BOA reflectances, also indicated by high PBIAS values (**Table 5**). On average, Sen2Cor BOA reflectances are about 30 % higher than TOA reflectances. BOA reflectances of iCOR (mean PBIAS: 10.8%) and ATCOR (mean PBIAS: 11.7%) are generally lower than Sen2Cor values and resemble each other in intensity and shape. Yet, in comparison to iCOR, ATCOR BOA reflectances decrease slightly with shorter wavelengths. ACOLITE BOA reflectances are lowest for all ROIs and only 6.6% higher than TOA reflectances, on average. **Figure 7E** illustrates the slopes of the respective spectra showing largest differences for Sen2Cor (mean slope MAE: 1.17E-03); high slope MAE values in **Table 5** confirm this observation. For iCOR (mean slope MAE: 2.64E-04) and ATCOR (mean slope MAE: 2.29E-04), **Figure 7E** highlights differences in bands 2–4 and similarities in bands 5–7. Overall, ACOLITE fits the shape of the resampled in situ spectrum best, also indicated by low slope MAE values in **Table 5** (mean slope MAE: 1.92E-04). Coefficients of determination in **Table 5** confirm this trend and are highest for ACOLITE (mean r²: 0.9758), followed by ATCOR (mean r²: 0.9728) and iCOR (mean r²: 0.9702), while Sen2Cor shows the weakest linear relationship (mean r²: 0.9137), on average.

TABLE 5 | Performance measures of resampled mean in situ and BOA spectra for different ROIs on the ice floe.

## Open Water

**Figures 7B,D,F** illustrates the results for the open water measurements. It is obvious that ACOLITE and Sen2Cor

TABLE 6 | Performance measures of resampled mean in situ and Rrs spectra for respective AC processors in comparison with mean resampled in situ spectrum at ROI 10A.

intensities are one order of magnitude higher than the in situ values, also demonstrated by high MAE values in **Table 6**. High coefficients of determination (r²∼0.96), however, indicate a strong linear correlation between resampled in situ and Rrs spectra, although the spectral slope is different (slope MAE in **Table 6** and **Figure 7**). iCOR average Rrs is almost 10 times higher than average in situ data (MAE∼0.01), but spectral shape is more similar compared to ACOLITE and Sen2Cor indicated by a small slope MAE (2.17E-05) and high r² (∼0.97). Polymer retrieved Rrs are in the same order of magnitude as measured Rrs, illustrated by a small MAE of 0.0015, while ACOLITE, iCOR and Sen2Cor overestimate Rrs. The lowest slope MAE (1.54E-05) indicates that Polymer resembles the shape of the in situ spectrum most accurately. The comparably weak correlation (r²∼0.85) results from a drop of the spectrum in bands 5–8A (**Figure 7B**). The ATCOR output is also in the correct order of magnitude over the full spectrum and has the smallest MAE and a comparable low slope MAE. It has to be noted, however, that the spectral shape in the visible range differs, also indicated by a weak linear correlation (r²∼0.18). The slope MAE is higher than for Polymer, highlighting the difference in shape in visible bands 2 (490 nm) to 4 (665 nm).

### Application Examples

To illustrate the consequences of different ACs on parameter retrieval, we computed the average BOA reflectance from band 2 to 7 (490–782 nm) as a rough approximation to broadband albedo (**Figure 8**). In addition, we calculated NDMI, which can be used as a proxy for melt pond fraction (Istomina and Heygster, 2017) according to

$$NDMI = \frac{\lambda\_{Green} - \lambda\_{NIR}}{\lambda\_{Green} + \lambda\_{NIR}} \tag{8}$$

where λGreen is S-2 band 3 (560 nm) and λNIR is S-2 band 8A (865 nm) from the respective AC product. Theoretically, NDMI values range from −1 to +1. Due to the specific reflectance characteristics of the different sea ice surfaces in these wavelength regions the value range, however, is practically limited to values >0. Increasing values indicate an increasing melt pond fraction. **Figure 8** illustrates that the AC processor influences parameter retrieval and, in reverse, that the application rules the choice of AC processor. Note that median TOA reflectance of the floe surface is 0.75, and BOA reflectances are approximately 0.80 in the ACOLITE product and ∼0.84 in the ATCOR and

iCOR products, reaching 0.97 in the Sen2Cor product. In general, NDMI results show the same pattern, whereas median NDMI values are lowest for TOA (0.06), followed by ACOLITE (0.12), ATCOR (0.13), iCOR (0.14), and Sen2Cor (0.18). This corresponds to the spectral behavior illustrated in **Figure 7**.

## DISCUSSION

## Comparability of Field and Satellite Data

The comparison of field spectroscopy and satellite observations involves a number of difficulties, which particularly applies to remote regions such as the Arctic where ideal measurement opportunities and coincident satellite overpasses are rare. We conducted spectral measurements of the floe surface within a period of ±5 days of the satellite overpass under varying illumination conditions and possibly changing surface properties due to the onset of melt. We postulate, however, that an analysis of spectral shapes is reasonable as the anisotropic reflectance factor of snow at nadir is ∼1 (Perovich, 1994) and the shape of the spectral albedo is constant under different illumination angles (Malinka et al., 2016). For snow surfaces, Goyens et al. (2018) and Bourgeois et al. (2006) observed no significant relationship between variations in reflectance anisotropy with illumination conditions for sun zenith angles <65◦ . In addition, the effect of melt on the spectral shape of nadir reflectance of snow and ice is comparably low in the VIS (Perovich, 1994, 1998). Nevertheless, we are unable to exclude potential differences in spectral slope due to melt processes. Besides the temporal differences, it is challenging to represent the spatial variability within a 20 m × 20 m area measured by S2-A via field measurements, even when measured at different locations. Moreover, a re-location of the measurement positions on a drifting floe in combination with spatial uncertainties of both GPS devices and satellite image positioning exacerbates a pointby-pixel comparison. For this reason, we included surrounding pixels into the analysis. However, this method is vulnerable to extreme spectral features in a ROI that can have a great impact on the average spectrum. The standard deviations of the ROIs located on the floe, however, are small indicating low spectral variability (**Figure 7C**); we therefore assume the ROI approach to be a valid method for evaluating the AC

processors. ROI 10A over open water was defined manually and excluded pixel with unusual shape, presumably resulting from sub-pixel contamination due to floating ice. This preselection may substantially influence the outcome of the study, as pixel with expected shape are preselected. The pre-selection, however, avoids integration of extreme spectral features, which may have a major influence on the average spectrum due to sub-pixel contamination.

## Quality of In Situ Data

For some ROIs, field spectrometer data show reflectances >1.0, although both field spectrometer setups were calibrated with reflectance standards. The Ocean Optics setup also enabled the correction of illumination changes. For ASD measurements, we surveyed changes in illumination by measuring the reflectance standards before each subsequent calibration. We manually deleted ASD measurements affected by illumination changes, i.e., spectra with unusual shape. Successive spectra with stable shape showing reflectances >1.0, however, remained in the dataset. In addition to the differences regarding sensitivity to changes in illumination, both setups vary regarding their FOVs. The smaller FOV of the Ocean Optics makes it more vulnerable to twisting of the pole and anisotropic effects due to ice surface topography. An inter-comparison of spectrometers in the field revealed that Ocean Optics reflectance values are ∼0.1 higher than ASD reflectances, on average, while the spectral shape was very similar. Warren (1982) stated that cloud cover affects both the spectral distribution of irradiance and the effective incident zenith angle, which causes an increase in all-wave snow albedo of up to 11%. Painter and Dozier (2004) attributed reflectances >1.0 to anisotropic effects of snow, which depend on the sun zenith angle, the grain size and the viewing geometry. Anisotropic effects may also have influenced the measurements described by Goyens et al. (2018) and Perovich (1998), who report reflectances >1.0 at approximately 400 nm measured under clear sky conditions. Peltoniemi et al. (2005) included snow wetness as an additional parameter influencing the spectral albedo and anisotropic behavior of snow and ice surfaces. As we evaluate the spectra only regarding to their shape, however, results remain unaffected by values >1.0. Analysis of water samples in the deepwater polynia show no significant absorption or scattering of the water constituents in the near infrared spectral range (>700 nm). Accordingly, beyond the visible range, measured Rrs is smaller than 1.5 × 10−4 , i.e., near zero (black) in the near infrared. Measurements influenced by ice floating into the FOV of the RAMSES spectrometer were excluded from the evaluation. The AERONET AOT values (∼0.02) are within the range reported for typical Arctic background aerosol concentrations, with values between 0.02 and 0.08 (Tomasi et al., 2012; IOCCG, 2015).

## AC Processors

Differences between BOA or Rrs products may originate from the different processors, the integrated LUTs, different parameterization and image-based parameter retrieval regarding WV content, aerosol type and AOT, or different adjacency correction schemes, and the interplay of these components. Varying options for the user to modify the parameterization further complicate a direct comparison. Moreover, some of the processors are still under development and not yet entirely documented which hinders a detailed reasoning.

The LUTs applied in the AC processors have been generated with different radiative transfer models (MODTRAN 5, 6SV, libRadtran, and Successive Order of Scattering). These models differ regarding the implementation of physical equations, approximation of the radiative transfer equation and subroutines which can result in differences (e.g., Callieco and Dell'Acqua, 2011). Evaluating the impact of a certain radiative transfer model on AC, however, is beyond the scope of this study. The aerosol models covered by the pre-calculated LUTs cover common aerosol types. The options, however, differ between the AC processors (section Atmospheric Correction Processors). It is questionable how well these models represent Arctic aerosol components, e.g., mineral dust, sulfate aerosols and sea salt (at Barrow, Alaska; Tomasi et al., 2012; IOCCG, 2015). In late spring and summer, a phenomenon known as Arctic Haze results from anthropogenic aerosol accumulation in winter is common and occasionally aerosols from forest fires or volcanic eruptions are observed (IOCCG, 2015). To analyze the impact of the aerosol model we conducted a sensitivity analysis with different aerosol models in ATCOR. For the desert, maritime and rural aerosol models, differences in the VIS bands for the bright snow and ice surfaces were ∼2%, while the urban aerosol model resulted in an increase of ∼5–10% reflectance.

ACOLITE is designed for aquatic applications; the main differences compared to other AC processors are the AOT retrieval by the dark spectrum fitting algorithm and the consideration of skyglint. Its performance over ice surfaces, however, is good regarding the shape of BOA reflectance spectra. For sea ice surfaces differences in surface reflectance due to skyglint correction are <0.02. No water vapor and ozone values could be obtained from ancillary data, which might due to the developmental state of the processor; thus, the default values were used for AC and it remains unclear if the retrieval of water vapor and ozone in the Arctic from ancillary data is possible. Given the applied parameterization, the dark spectrum fitting algorithm results in the application of a continental aerosol model and an AOT of 0.0084 at 550 nm, i.e., comparably close to the measured value of 0.02 (AERONET-MAN reports uncertainties ≤0.02 per band; NASA Goddard Space Flight Center, 2017). We assume that the high Rrs values over water result from adjacency effects, which are not addressed by ACOLITE.

ATCOR performs well over ice surfaces while its performance over open water is poor as the spectral shapes of BOA and in situ data do not correspond. Due to the lack of dark pixels in the S2-A data set we defined AOT by visibility (65 km), illustrating that the dark vegetation method is unsuitable for AOT retrieval in this image. For WV, image-based retrieval resulted in an overestimation of WV (∼1.8 cm). We therefore defined a default value (1.0 cm) referring to the AERONET data. Few pixels (1.2– 1.4%) of bands 5–9 show negative reflectance values at open water areas, which may be caused by underestimated visibility, the aerosol type or calibration issues. Application of haze removal was not possible due to the small amount of clear pixels (<4%). Increasing the visibility parameter (up to 120 km) only had negligible consequences for the snow and ice surfaces and did not improve the spectral shape of open water pixels. This corresponds to Huck et al. (2007) who found that increasing visibility in a 6S forward simulation for AVHRR only affects open water surfaces. Changing the WV content from 1.8 to 1.0 cm had a significant influence, especially in S-2 bands that cover WV absorption features, i.e., band 8 (842 nm), band 9 (940 nm), and also band 5 (705 nm). For simulated Sea-viewing Wide Field-of-view Sensor (SeaWiFS) like spectra at the edge between sea ice and open water, Bélanger et al. (2007) found contaminations due to the adjacency effect up to 24 km off the ice edge resulting in unreliable results in derived chlorophyll concentrations and optical properties of open ocean water. Bulgarelli and Zibordi (2018) found that for S-2/MSI adjacency effects caused by highly reflecting surfaces (e.g., snow) are theoretically detectable at least 20 km off the coast in mid-latitude coastal environments. In contrast to these simulations, the area of interest is characterized by a small-scale mixture of snow, ice, melt ponds and open water, where spatial patterns and average reflectance are similar on different spatial scales. To check whether this statement applies for sea ice, we tested different ranges (1, 3, 5, and 10 km) and results confirmed that the influence is negligible for the spectral shape of the snow surface and intensity of open water pixels.

iCOR addresses the requirements for aquatic and terrestrial surfaces including a correction of adjacency effects. iCOR was unable to assess AOT from the image as the retrieval algorithm relies on vegetation and soil pixels (De Keukelaere et al., 2018), questioning the suitability of this approach for the Arctic sea ice in the absence of in situ AOT values. Even though, the fixed WV column height of 2.0 cm showed good results over ice, iCOR's unknown sensitivity to WV aggravates AC processor inter-comparison. The good performance over water compared to ACOLITE and Sen2Cor is probably attributed to the correction of adjacency effects using the SIMEC algorithm. For the application of SIMEC, Sterckx et al. (2015b) point out that aerosol information (AOT and aerosol model) should be either estimated using a land-based aerosol retrieval or from sun photometer data. While we used AERONET data to set the default AOT, we could not parameterize the aerosol type. Further, SIMEC is based on the assumption that reflectance differences between water and surrounding surfaces are strongest in the NIR, which is true for vegetated surfaces but not for ice, where differences are strongest in the VIS. We found that iCOR is very sensitive to changes in AOT. Increasing AOT from 0.02 to 0.1 resulted in an overcorrection and, consequently, a decrease of the spectrum matching the in situ measurements very well (slope MAE ∼8.49E-06).

Polymer performs best over water and produces reasonable results even in the high contrast Arctic sea ice matching findings of Steinmetz and Ramon (2018). It is robust to adjacency effects because the atmospheric model fits any spectrally smooth nonwater component (Steinmetz and Ramon, 2018).The manual definition of ROI 10A guaranteed that we only considered deep-water pixels, which may also be important for the good performance of Polymer. Nevertheless, we observe a high spectral variability in the VIS bands. Moreover, some negative values appear in the red-edge and NIR regions. Sub-pixel contamination due to floating ice should be corrected due to its spectrally-flat contribution to the reflectance signal.

Usually, Sen2Cor applies image-based retrievals to generate maps of AOT and WV. The lack of dark pixels in the S-2A image, however, prohibited an AOT retrieval based on the DDV method; thus, AOT was estimated from the provided visibility (65 km). Values for scene average WV and AOT at 550 nm are 1.51 cm and 0.132, respectively. Compared to the AERONET data, however, both values are very high. A further test with a fixed visibility of 120 km, corresponding to an AOT of ∼0.0780, as proposed by Pflug et al. (2016) as a good practice for clear air and low AOT over Antarctica, did not improve the results for BOA reflectance and Rrs significantly. The increased BOA reflectance in comparison to the other processors fit findings of Li et al. (2018) who report that surface reflectance derived with Sen2Cor is generally overestimated, in particular for bright pixels. The high values over water may be attributed to the lacking correction of skyglint and adjacency effects or the provided LUTs.

A general challenge in high latitudes is the low sun elevation throughout long periods of the year. Huck et al. (2007) found that the solar zenith angle has a major influence on the retrieved TOA reflectance of a nadir-viewing sensor. The low reflectance levels of water bodies further increase this effect, because atmospheric path length increases with decreasing sun elevation. Consequently, the amount of energy absorbed in the atmosphere increases while the amount that arrives at the Earth's surface decreases. At the same time, the amount of light that is specularly reflected at the water surface increases with decreasing sun elevation, while a lower amount is transmitted into the water body. As a consequence, the water body scatters less light toward the sensor (Hieronymi, 2013). Further, the influence of surface roughness on the transmission of light into the waterbody increases with sun zenith angle which challenge wind speed-dependent AC (Hieronymi, 2016). When the sun is close to the horizon, adjacent ice floes may also have a shading effect. The same applies for ice floes, where shadows of ridges increase with lower sun elevation and anisotropic effects increase for sun-zenith angles >65◦ (Bourgeois et al., 2006). Another issue regarding low sun elevation is the limitation of most processors to sun zenith angles ≤70◦ . The reason is that the LUTs are generated with radiative transfer models that base on the assumption of a plane parallel atmosphere rather than an actual spherical-shell atmosphere (IOCCG, 2015). Since in our study the sun zenith angle is 62.5◦ , we assume that this effect is negligible.

## Sentinel-2A Data

Typically, reflectance spectra of snow and ice (like in situ spectra) are characterized by a steady decrease from the blue domain toward longer wavelength regions. S-2A TOA spectra, however, show a steep decrease from band 2 to 3 (490 and 560 nm, respectively), followed by an increase from band 3 to 4 (560– 665 nm; **Figure 7A**). A comparison of BOA reflectance spectra with typical spectra found in literature (e.g., Perovich, 1994, 1998) indicates that ACOLITE, ATCOR and iCOR produce realistic results. As described in section Comparability of Field and Satellite Data, however, melt processes may have caused reflectance changes between in situ and satellite data. Thus, the evaluation of BOA spectral shape is challenging. In general, however, BOA shapes of ACOLITE, ATCOR, and iCOR are similar to the shape of the in situ spectra, i.e., the differences in slope MAE and r² are small. The shapes of Sen2Cor BOA reflectances do not resemble the typical shapes of snow and ice spectra as they show a steep decrease from the blue wavelength region toward longer wavelengths, interrupted by a local peak in band 5 at 705 nm. Furthermore, Sen2Cor is the only processor that produces reflectance values >1.0. Nadir reflectances of snow and white ice may exceed 1.0 in the blue part of the electromagnetic spectrum (e.g., Perovich, 1994; Goyens et al., 2018) but Sen2Cor BOA reflectances exceed 1.0 even in bands 3 and 4 (560 and 665 nm, respectively). This may result from anisotropic effects due to snow parameters or floe topography and matches observations of Li et al. (2018) who found that Sen2Cor overestimates BOA reflectance; ACOLITE, ATCOR, and iCOR BOA reflectance values, however, are <1.0 in all bands.

## Application Examples

As indicated in section Application Examples, the choice of AC processor influences parameter retrieval. The differences in BOA reflectance illustrated in **Figure 7** have a great impact on NDMI and surface albedo (**Figure 8**), which is roughly approximated here by the mean BOA reflectance of bands 2–7 (490–783 nm). The shifted peaks in the histograms refer to the offset-like behavior of the spectra from different AC processors illustrated in **Figure 7**. The example illustrates well that the application of a certain AC processor and resulting BOA reflectances affects the retrieval of geophysical parameters from S-2A imagery.

## CONCLUSION