An improved sea ice detection algorithm using MODIS: application as a new European sea ice extent indicator

. The continued loss of sea ice in the Northern Hemisphere due to global warming poses a threat on biota and human activities, evidencing the necessity of efﬁcient sea ice monitoring tools. Aiming at the creation of an improved European sea ice extent indicator, the IceMap250 algorithm has been reworked to generate improved sea ice extent maps at 500 m resolution at nadir. Changes in the classiﬁcation approach and a new method to correct artefacts arising from the MODIS cloud mask allow the enlargement of the mapped area, the reduction of potential error sources and a qualitative improvement of the 5 resulting maps, while systematically achieving accuracies above 90 % . Monthly sea ice extent maps have been derived using a new synthesis method which acts as an additional error ﬁlter. Our results, covering the months of maximum (March) and minimum (September) sea ice extent during two decades (from 2000 to 2019), are a proof of the algorithm’s applicability as an indicator, illustrating the sea ice decline in the European regional seas. We observed no signiﬁcant trends in the Baltic (-2.75 ± 2.05 × 10 3 km 2 yr − 1 ) although, on the contrary, the European Arctic seas display clear negative trends both in March 10 (-27.98 ± 6.01 × 10 3 km 2 yr − 1 ) and September (-16.47 ± 5.66 × 10 3 km 2 yr − 1 ). Such trends indicate that the sea ice cover in March and September is shrinking at a rate of ∼ 9 % and ∼ 13 % per decade, respectively, even though the sea ice extent loss is comparatively ∼ 70 % greater in March. Therefore,


Introduction
The Arctic sea ice cover has been changing rapidly over the last decades, with its overall extent declining steadily since the first satellite observations in the late 1970s Comiso et al., 2008;Cavalieri and Parkinson, 2012;Massonnet et al., 2012;Meier et al., 2014), shrinking at a rate of about 10 % per decade in the last years (Comiso et al., 2008) and 20 reaching its historical minimum on September 2012 (NSIDC). Moreover, sea ice thickness has decreased as much as 65 % in the period extending from 1975 to 2012 (Lindsay and Schweiger, 2015). This massive loss of ice is unprecedented in the last few thousand years (Polyak et al., 2010). Although it is attributed both to climatic variability and to external forcing caused by an anthropogenic release of greenhouse gases Myhre et al., 2013), nowadays human influences have driven climate to exceed the bounds of natural variability (Karl and Trenberth, 2003). All projection models 25 agree that Arctic sea ice will continue shrinking and thinning, eventually leading to ice-free summers in the following decades (Massonnet et al., 2012;Stroeve et al., 2012;Collins et al., 2013;Notz and Stroeve, 2016) and even as soon as in the late 2030s (AMAP, 2017).
The dynamism of the sea ice and the role it plays on the climate, biota and on human activities makes necessary its moni-30 toring. Due to the difficulty of acquiring in situ observations, nowadays satellite imagery is the main tool to monitor sea ice at a global scale (Teleti and Luis, 2013). Several sea ice variables are continuously obtained and distributed by institutions such as the EUMETSAT Ocean and Sea Ice Satellite Application Facility (OSI-SAF) or the National Snow and Ice Data Center (NSIDC), commonly at resolutions between 10 km and 25 km. 35 In 2016, the European Environment Agency (EEA) published a sea ice extent indicator (EEA, 2016) aiming at the monitoring of sea ice trends both in the Arctic Ocean and the Baltic Sea. However, observations in both regions are not directly comparable, as sea ice extent in the Arctic was derived from the OSI-SAF passive microwave satellite data at 10 km resolution, while data for the Baltic came from several sources, including in situ observations and air temperature proxies. Therefore, in order to homogenize data acquisition in both regions, we tested the IceMap250 algorithm (Gignac et al., 2017), which produces 40 sea ice extent maps at 250 m thanks to a downscaling technique by Trishchenko et al. (2006). Testing revealed that IceMap250 may be severely affected by resolution-breaking artefacts found in the MODIS cloud mask, as happens with MODIS sea ice products MOD29 and MYD29. We also found that the mechanisms to avoid water and ice false positives are not optimal when one of those surfaces is absent. Additionally, Xiong et al. (2006) and Khlopenkov and Trishchenko (2008) argue that the band to band misregistration of MODIS Aqua may exceed the resolution achieved with the downscaling. Thus, as the usefulness of 45 the downscaling for sea ice detection was already demonstrated in Gignac et al. (2017), we keep the standard 500 m resolution to reduce processing time and to allow the use of both the instruments in Terra and Aqua. Therefore, the present work has two main objectives: 1) to develop an improved 500 m resolution sea ice detection algorithm (IceMap500) based on IceMap250, and 2) to prove the utility of IceMap500 as a new European sea ice extent indicator by 50 analysing the sea ice trends in the European regional seas from 2000 to 2019. The analysis covers the northernmost European sea regions defined by the European Union's Marine Strategy Framework Directive (MSFD) where sea ice might occur, and is restricted to the months when the maximum and minimum sea ice extent is reached in the Northern Hemisphere, that is, March and September, respectively.

Study area
This work focuses on the European regional seas established by the MSFD. As sea ice only occurs in the northernmost oceanic sea regions or in enclosed, low-salinity water bodies such as the Baltic Sea, spatial coverage has been significantly reduced to avoid the processing of uninformative data. Therefore, the final indicator extends over the sea regions belonging to the Arctic, North-East Atlantic Ocean and the Baltic Sea, as is shown in Fig. 1, covering an area roughly 4 ×10 6 km 2 . With the inclusion 60 of a 400 km buffer to coherently join all the target regional seas in a single study region, the totality of the processed area ascends up to approximately 8 ×10 6 km 2 . Figure 1. Northern European regional seas, as defined by the MSDF: 1) Iceland Sea, 2) Norwegian Sea, 3) Barents Sea, 4) White Sea, 5) Baltic Sea, 6) Greater North Sea, and 7) Celtic Seas. In medium blue are shown the target sea regions, whereas in light blue is represented the generated buffer, whose external limit corresponds to the total processed area.
Oceanic sea ice in the Northern Hemisphere has both a perennial and a seasonal fraction. Typically, maximum and minimum sea ice extent are reached in March and September (Stroeve et al., 2008), respectively, with the perennial fraction being mostly enclosed in the Arctic basin (Comiso, 2009). According to NSIDC's Sea Ice Index (Fetterer et al., 2017), sea ice is present during the Arctic winter months in some of the European sea regions (i.e. the Barents Sea, the White Sea, and the northernmost areas of the Norwegian Sea). As sea ice is also found along the eastern coast of Greenland, it may occasionally reach the Iceland Sea or the waters surrounding the Jan Mayen island.

70
The ice cover in the Baltic Sea, which has no perennial fraction, can be highly variable due to the milder climate, often resulting in different freezing and melting periods during the same winter (Granskog et al., 2006). The sea ice season usually lasts for six to eight months, starting in October or November in the Bothnian Bay and the Gulf of Finland. Maximum extent is normally reached in March (Haapala et al., 2015).

Selected data
Due to its balance between temporal and spatial coverage, we use MODIS visible and infrared imagery to generate sea ice extent maps at 500 m resolution at nadir. MODIS is on board of NASA's sun-synchronous satellites Terra and Aqua, launched in 1999 and 2002, respectively. It acquires data in 36 spectral bands, ranging from the visible spectrum to the thermal infrared.
Spatial resolution at nadir varies from 250 m (bands 1 and 2) to 500 m (bands 3-7) and 1 km (bands 8-36). MODIS has a large 80 swath width of 2330 km, allowing a revisit time of 1 to 2 days. Although MODIS is severely affected by weather and lighting conditions, its resolution is much higher than that of passive microwave sensors: widely used microwave radiometers such as SSM/I-SSMIS provide data at 25 km cell size. Moreover, its swath width is greater than that of the synthetic-aperture radar and other sensors operating in the visible-infrared spectrum such as those carried by the Landsat series and Sentinel-2, which nonetheless acquire data at even higher spatial resolutions (30 to 10 m). Thus, we use the data shown in Table 1  The algorithm uses TOA radiance as input data, which is converted to TOA reflectance as in previous sea ice detection works (Hall et al., 2001;Gignac et al., 2017). Therefore, the threshold values used in the classification are higher than if surface re-90 flectance was used due to the contribution of the atmosphere. Although TOA data does not reflect the physical properties of sea ice and water, it avoids extensive processing due to atmospheric correction, facilitates the algorithm's replicability and ensures the consistency of the dataset. Note that the objective of the algorithm is to map sea ice presence rather than using reflectance as a proxy to get other physical variables such as sea ice concentration. In addition, it must be considered that the use of data at higher resolutions than MODIS would cause the processing to be computationally very demanding, especially if 95 covering large areas, as is the case of the present study. It would also render data with sparse spatial and temporal coverage as higher-resolution sensors typically have smaller swath sizes and longer revisit times. This would be especially problematic for areas located at mid-latitudes, although this effect is minimized at the poles.

100
The IceMap250 algorithm relies on a hybrid cloud masking approach designed to minimize error while maximizing the mapped area, using the MODIS MOD35 cloud mask and an additional visibility (VIS) mask, both at a 1 km resolution. Threshold tests based on the Normalized Difference Snow and Ice Index 2 (NDSII-2) (Keshri et al., 2009) and the TOA reflectance at 545-565 nm are used to classify sea ice and water in the masked datasets. However, this classification process faces some challenging potential errors. One of the most notable classification errors arises from the NDSII-2 test, which uses the Jenks 105 natural breaks optimization (Jenks, 1967) to split pixels in two groups, regardless of the surface classes present in a scene.
When batch processing MODIS data it may be likely to run into scenes lacking either ocean water or sea ice and, consequently, the Jenks optimization splits pixels into both surface classes erroneously. Clouds which are undetected by the MOD35 cloud mask algorithm (Ackerman et al., 2010) and sun glint over ocean water may also be common error sources due to the similar reflectance characteristics to sea ice. Additionally, as already discussed in Gignac et al. (2017), bare ice and melt ponds may 110 also fail the classification tests due to the similarity with ocean water.
Nevertheless, the most important issue concerning the quality of the classification arises from the MOD35 cloud mask. It features resolution-breaking square artefacts of 25 km side length along the ice edge ( Fig. 2) that originate in the setting of the snow/ice background flag during the production process of the mask (Riggs and Hall, 2015). The source of the artefacts 115 is NSIDC's Near-real-time Ice and Snow Extent (NISE) product, based on SSM/I-SSMIS passive microwave data at 25 km resolution, that is used to determine the flag's state. Such artefacts propagate from one product to another, and can also be seen in MODIS sea ice products MOD29 and MYD29. They can occupy extensive areas in some scenes, causing the loss of many cloud-free classifiable pixels and preventing the detection of the ice edge. To mitigate those potential classification errors, IceMap500 features changes in the data masking and the classification rules, additional threshold tests, a smaller artefact correction algorithm and a new monthly map synthesis approach (see the structure in Fig. 3).

The masking
IceMap500 uses the same hybrid cloud masking approach as IceMap250. Nevertheless, the MOD35 mask includes additional 125 constraints so not only cloud cover is considered, but also the lighting conditions, sun glint and the presence of land. This information is contained within the MODIS product MOD35_L2, which provides multiple quality assessment flags, as is summarized in the product's user's guide (Strabala, 2004). We use the following flag states: 1. Unobstructed FOV, selecting only pixels identified as confident clear. This flag is the cloud mask already used in IceMap250.

130
2. Day/Night, selecting only pixels identified as day. This flag is of special importance during the winter months, when the polar twilight zone reaches the lowest latitudes and, therefore, the available daytime area becomes scarcer.
3. Sun glint, selecting only pixels identified as no sun glint. This way, areas with sun glint caused by the reflection angle of the sun being between 0 • and 36 • are discarded. It is important to emphasize that other sun glint sources are not considered (Ackerman et al., 2010).  4. Land/Water, selecting only pixels identified as water. Land masking is crucial to ensure the quality of the resulting classification because, as already pointed out in Gignac et al. (2017), an incorrect masking may generate sea ice false positives due to the reflectance contrast of land with water.
The VIS mask is used and calculated as in IceMap250. This mask is intended to identify areas where visibility is sufficient to perform a classification, for the sole goal of detecting open water. It uses the normalized difference between the MODIS 140 thermal bands 20 and 32 as in Eq. (1).
Threshold values are obtained by calculating the standard score of R (B20/B32) , as seen in Eq.
(2), where µ and σ are the mean and standard deviation of R (B20/B32) , respectively. Pixels where VIS < 0.5 are tagged as having enough visibility. The masking produces the MOD35 and the VIS datasets, which are later classified separately. Note that while masking is done at 1 km resolution, the classification uses data at 500 m, so sea ice and water are mapped at 500 m within the mask limits.

The classification tests
The original thresholding method used in IceMap250 classifies as sea ice all pixels that pass any of the following two threshold tests: 1. NDSII-2 threshold test (t ndsii2 ). The threshold value k is determined by slicing the NDSII-2, shown in Eq. (3), into two classes with the Jenks natural breaks optimization (Jenks, 1967), which maximizes inter-class variance and minimizes intra-class variance. Pixels in the first group (i.e. below k) are classified as sea ice.
2. Green TOA reflectance threshold test (t b4 ). A pixel is tagged as sea ice if its reflectance is 17 % at 545-565 nm (band 155 4). This threshold is based on the contrast in reflectance between ice and water at visible wavelengths as suggested by Riggs et al. (1999) and validated by Gignac et al. (2017).
However, in IceMap500 a new threshold test is introduced: 3. Mid-range infrared temperature test (t b20 ). This new threshold is based on the Sea Surface Temperature (SST) using band 20 (3.660-3.840 µm). It is always used in conjunction with t b4 , but only during the MOD35 dataset classification.

160
Therefore, sea ice is classified only when both B4 17 % and SST < 1 • C. The goal of t b20 is to reduce potential sea ice false positives due to sun glint, as not all sources are considered in the MOD35 mask (see subsubsection 2.3.1).
This threshold intends to include melt ponds, leads, and water close to the ice edge to prevent breaking the 500 m resolution. The SST test relies on a simple atmospheric correction described in the MODIS SST algorithm theoretical basis document (Brown and Minnett, 1999) for mid-range infrared SST derivation, as in Eq. 4: where T B20 is the brightness temperature of MODIS band 20. Mid-range infrared has been selected instead of thermal infrared because the atmospheric correction is straightforward and may be affected by reflected solar radiation, making easier the exclusion of sun glint as a result of the temperature increase. The 1 • C threshold is selected so melt ponds and pixels around the ice edge are included (see Zhang et al. (2017) for a detailed discussion on the temperature of melt 170 ponds), while still leaving out most water in the study area susceptible of being affected by sun glint (refer, for instance, to global SST products by the NOAA).
In addition, in IceMap500 a more restrictive classification approach is adopted to compensate the output of t ndsii2 in scenes with a single surface class, as the Jenks optimization will still split data in two groups. The classification rules are datasetdependent. Nevertheless, due to the merging of the MOD35 and VIS maps, changes in a single dataset classification ultimately 175 affect the whole outcome. The IceMap500 classification rules are shown in Table 2: sea ice is only mapped in the MOD35 dataset when there is consensus between the tests, while in the VIS dataset it is mapped whenever t b4 is positive. A downside of this method is that it may leave some melt ponds as NoData, since in the most advanced melting states they tend to show NDSII-2 values similar to water (Gignac et al., 2017). Once the MOD35 map is created, an additional set of tests is introduced to attenuate the effects of the NISE artefacts or blocks present in the MOD35 mask, which propagate to the MOD35 classification and ultimately to the monthly extent maps, that may be extensively affected. Although the inclusion of this correction increases the chances of classification errors, it greatly improves the quality of the maps and increases the classified area. It is intended to affect only cloud-free areas set as NoData that are close enough to sea ice. To avoid error amplification, sea ice clusters below 100 pixels are deleted before the correction:

185
if those clusters are found far from the ice edge it is likely that they originate from sun glint or unmasked clouds, while those found close to large clusters of sea ice will be classified again as such. The MOD35 correction includes five tests, as illustrated in Fig. 4. 1. NoData test. NoData pixels pass the test, while classified areas remain the same. All pixels set as NoData during the MOD35 classification also undergo the tests, and may be finally labelled as sea ice or water. 2. Euclidean distance test. NoData pixels found at 35 km or closer to a cluster of sea ice pass the test; those found above this threshold are left as NoData. This distance is roughly equal to the diagonal of NISE's 25 km artefacts, and is used to reduce the chances of misclassifying clouds as sea ice by setting a buffer along the ice edge.
3. Band 7 TOA reflectance test (t b7 ). Pixels below 3.5 % TOA reflectance at 2.105-2.155 µm pass the test, otherwise they are left as NoData. This threshold is based on the low reflectance that water, snow, and ice display around 2 µm: Finally, the MOD35 map and the result of the MOD35 block correction are joined together in a single map, which is later combined with the VIS map according to the compositing rules in Table 3.

Monthly map synthesis and extent derivation
The corrected MOD35 and VIS maps created for each scene are combined to take advantage of the strengths of both the 210 MOD35 and the VIS classification methods, following the criteria seen in Table 3. The extensive cloud cover found in most scenes and the restrictiveness of the classification implies little area is finally mapped, although the new correction reduces the impact of the cloud mask. In any case, many scenes are required to map large expanses of the sea ice cover. In IceMap250 weekly maps are derived using a majority filter, with every pixel classified as sea ice assumed to be equally reliable. Here, a new monthly map synthesis method is proposed based on the number of coincident sea ice classifications achieved in each 215 pixel, meaning that pixels classified as sea ice in a large number of scenes will have higher reliability. The synthesis maps are generated by calculating the sum of composite maps where ice = 1 and water = 0, and later normalizing the results according to the maximum number of coincident sea ice observations achieved. The output provides information about where is sea ice more likely to be found according to the processed MODIS scenes, thus we appropriately refer to the resulting maps as sea ice presence likelihood maps (Fig. 6).

220
Pixels below a selected threshold value in the sea ice presence likelihood maps can be discarded to get rid of the least reliable observations, acting as an additional post-classification filter. In our case, the selected threshold is 10 %, which represents a balanced compromise between error filtering and area mapped. This synthesis approach generates sea ice extent maps, as the constant motion of the ice tends to hide the presence of features such as leads, cracks, polynyas and ice floes. Finally, the 225 euclidean distance from both sea ice and water is calculated, and is later used to fill NoData gaps by setting as sea ice those pixels closer to sea ice than to water, smoothing the ice edge and generating a continuous sea ice extent map for the given month. 3 Results All three cases display negative trends, indicating a shrinking of the sea ice cover. Table 4 shows numerically the decrease in sea ice extent. According to the calculated trends, in the European Arctic the sea ice decline is ∼70 % faster in March than in September. Although September's extent is comparatively smaller, the standard error of the trends is similar in both months (∼ 6×10 3 km 2 yr −1 ), indicating September displays a higher variability, as evidenced by the lower R 2 value. Nevertheless, 240 both trends have been found to be statistically significant, assuming the null hypothesis (H 0 ) that the slope of the trend lines is zero. In particular, the March trend displays a very low p-value, indicating a significance level of ∼ 99.98 %.

Sea ice trends
While the Baltic Sea trend line in Fig. 7 clearly shows a negative tendency, the monthly results also display large variability.
This causes R 2 to be very low and the standard error of the slope to be almost equal to the slope itself. Moreover, this trend 245 is not statistically significant, so H 0 can not be rejected. Therefore, the observed trend may reflect a real negative tendency masked by high natural variability or may simply result from stochastic sea ice extent measures independent of time.

Accuracy assessment
We randomly selected eight years to perform the quality assessment, from which a total number of 32 scenes have been used, 250 that is, two scenes per month to allow comparison. As a prerequisite, each validation scene must have both sea ice and water Table 4. Numerical results of the Arctic and Baltic sea ice trends and the standard error of the slope, along with two goodness of fit estimators: the coefficient of determination and the p-value. P-values were obtained from two-tailed tests assuming 16 degrees of freedom and null hypothesis that there is no correlation between the two variables, i.e. that the slope of the trend line is zero.

Month
Trend R 2 p-value (10 3 km 2 yr −1 ) Arctic  1966 pixels, otherwise it is discarded. Validation has been carried out with confusion matrices by generating 1500 random points per scene over the classified areas. Those points have been manually tagged as either sea ice, water, or clouds, with the help of the corresponding RGB composite. Although no clouds are mapped in the algorithm, points found over clouds opaque enough to avoid the identification of the Earth's surface add to the total sea ice commission error. As already noted in Gignac et al.

255
(2017), this method requires the scenes to be validated by the same analyst in order to maintain its coherence.
Accuracy assessment results have been summarized in Table 5. All scenes achieved overall accuracies above 90 %, resulting in an average accuracy of 95.96 %. The average kappa coefficient of 0.853 indicates a strong agreement between classification and ground truth, despite being affected by scenes with few water validation points, causing the kappa coefficient to drop due 260 to the disproportion between classes. Individually, only 5 out of 32 computed kappa coefficients are found below the 0.800 value, while 10 are found between 0.800-0.900 and 17 above 0.900, indicating very strong agreement. The primary source of error affecting the classification is sea ice commission, with its mean value alone adding up to 7.33 %, that is, more than sea ice omission, water commission, and water omission combined.

265
By analysing separately both months, mean accuracy is found to be higher in March than September, differing by 1.97 %.
Accuracy results in September are also slightly more variable, showing a σ of 2.78 % versus 2.48 % in March. On the contrary, as due to the extensive sea ice cover March scenes are especially prone to almost lack water and, therefore, water validation points, very low kappa values are eventually obtained resulting in a lower mean kappa coefficient in March than in September.  being 2.74 % in March and 3.30 % in September, while water commission is 2.5 % and 1.87 %, respectively. Thus, globally, the major error contribution is due to the misclassification of clouds as sea ice, especially in September, while misclassification of sea ice as water and water as sea ice remain lower in the first case and minimal in the latter.
According to Chan and Comiso (2013), the MOD35 cloud mask tends to underestimate the cloud cover over sea ice, whereas 280 over open water it is overestimated but closer to reality. Indeed, most sea ice commission error in our validation is due to the misclassification of clouds as sea ice within the limits of the sea ice cover; in fact, despite the cloud fraction being much larger over open ocean than over sea ice, in the first case sea ice commission errors are uncommon. Some of the clouds that are commonly left undetected by the MOD35 cloud mask include low-level (top below 2 km), high-level (top above 6 km), and thin clouds less than 2 km thick (Chan and Comiso, 2013). Additionally, our validation showed some cloud shadows cast over 285 cloudy areas may sometimes be classified as clear. The rise of sea ice commission error during September may be explained by the fact that, as shown by Chan and Comiso (2013), late summer in the Arctic is considerably cloudier than winter, as lower sea ice concentration relates to a larger cloud fraction.
Since sun glint issues have been mostly solved, as evidenced by the minimal impact of water omission error, and most sea The Sea Ice Index (Fetterer et al., 2017) is a widely used global sea ice extent and concentration product distributed by the NSIDC, which is derived from satellite passive microwave data at 25 km spatial resolution. It covers from 1978 to the present, being updated on a daily basis, and provides monthly median sea ice extent maps. In spite of the difference in spatial resolution between the Sea Ice Index and our results, measuring the agreement between both datasets acts as an estimator of the quality and consistency of the algorithm's monthly composites. Agreement has been calculated as the coincident sea ice area fraction 300 between both datasets, as compared to the total sea ice extent including coincident and non-coincident area.  Mean agreement in March is 89.46 % with a standard deviation of 1.08 %, whereas in September mean agreement is lower, 85.53 %, and displays higher variability, with a standard deviation of 3.07 %. Only in a single case does the agreement fall 305 below 80 %, corresponding to September 2013: this particular case will be further discussed in section 4.

Discussion
Sea ice trends obtained from our monthly extent maps in the European Arctic are consistent with previous observations and both are statistically significant. In Cavalieri and Parkinson (2012), sea ice trends are shown by region: our study area approximately matches what the authors call Greenland Sea and Kara and Barents Seas. Data from 1979 to 2010 reveals in the Greenland Sea a trend of -9.5±1.9×10 3 km 2 yr −1 in winter and -4.8±1.6×10 3 km 2 yr −1 in summer, showing a larger sea ice loss during winter as in the present study. Trends in the Kara and Barents Seas are similar in winter and summer, being -12.2±2.4×10 3 km 2 yr −1 and -13.8±2.8×10 3 km 2 yr −1 respectively. However, our regional trends may not reflect the overall sea ice extent tendencies. This may be exemplified by the fact that in the Northern Hemisphere as a whole the sea ice loss is more pronounced in summer than in winter, and that in our case the minimum sea ice extent does not correspond to September 315 2012 (EEA, 2017; Cavalieri and Parkinson, 2012;. In the case of the Baltic Sea, no statistically significant trend can be inferred due to high interannual variability and the limited lifespan of MODIS. This, however, does not imply that H 0 (i.e. that the Baltic ice cover is stable) is true: previous research (Jevrejeva et al., 2004) based on data from coastal observatories covering years 1900 to 2000 reveals a significant decreasing 320 trend in sea ice occurrence probability in the southern Baltic Sea, while in the northern half ice occurs every winter. Moreover, it shows a shortening of the sea ice season and an advance in the date of break-up, especially in the northern areas. More recent analyses (Vihma and Haapala, 2009;Haapala et al., 2015) also indicate that over the last century the sea ice season has shortened and the occurrence of severe winters has fallen. According to EEA (2016), Baltic sea ice extent trends are affected by large interannual variability caused by the North Atlantic Oscillation that prevents them from being statistically significant.

325
The low water omission error obtained in the quality assessment reflects that most sun glint issues have been solved, both by the sun glint mask provided in MOD35_L2 and the more restrictive classification approach. Nonetheless, while sea ice omission and water commission are still low, they play a much more important role on the overall accuracy. The major source of error, according to the validation, are clouds not detected by the MOD35 mask. Additional thresholds could be introduced 330 to reduce unmasked cloud cover as much as possible, at the expense of increasing the running time of the algorithm which is already enlarged by the MOD35 block correction. As a result of its application, the area loss caused by the adopted restrictive classification approach is counteracted, as evidenced by the sea ice presence likelihood maps (see previous Fig. 6). Nonetheless, the potential presence of NoData gaps in the likelihood maps is an additional factor increasing uncertainty in our monthly sea ice extent derivation. Although those gaps are filled according to the minimum euclidean distance to sea ice or water, its 335 classification is not based on real observations and therefore uncertainty increases.
The resulting monthly sea ice extent maps show an agreement with NSIDC's Sea Ice Index almost always above 80 %, being higher in March than in September. Due to the difference in spatial resolution, agreements close to 100 % are not possible: aside from the position of the sea ice edge, this difference also affects the coastline, increasing error if sea ice is present.

340
Moreover, some fjords along the coast of Greenland which are permanently covered by glaciers are tagged as land by the Sea Ice Index, while the MOD35 mask includes them as ocean, thus being ultimately classified as sea ice by our algorithm. As the sea ice cover during September is considerably smaller and is mostly found along the coast of Greenland, non-coincident sea ice between both products due to the coastline discrepancy is proportionally larger in September, contributing to the lower agreement values. Disagreement also arises from the detection of fragmented sea ice and ice floes, which are frequent during the Arctic summer: due to the 25 km resolution of the Sea Ice Index, some of those areas may not exceed the 15 % sea ice concentration threshold used to determine extent and thus are tagged as water. September 2013, which displays the lowest agreement value, is an example of such behaviour (see Fig. 8).

350
IceMap500 has been shown to produce high quality sea ice extent maps by systematically achieving accuracies above 90 %.
Quality assessment revealed the most common error is sea ice commission caused by unmasked clouds, manifesting the key role that the MOD35 cloud mask plays on the overall accuracy of the algorithm. The addition of the NISE artefact correction substantially improves the delineation of the ice edge, preventing the propagation of such artefacts, and increases the area mapped, which is of capital importance when deriving daily and monthly maps due to the restrictiveness of the classification 355 and the weather dependence of MODIS visible and infrared data. However, although it has not been specifically designed to work in a single study area, its application in other regions has not been assessed and may yield different accuracies. High agreement between our monthly sea ice extent maps and NSIDC's Sea Ice Index prove the consistency of the monthly synthesis method and further exemplifies the overall good performance of the algorithm. Data produced by IceMap500 has proved useful to evaluate sea ice extent trends in the European Arctic and the Baltic Sea, exemplifying one of the potential applications it 360 may be used for. Significant negative trends have been observed both in March and September in the Arctic, while the Baltic