Seasonal evolution of the sea ice floe size distribution in the Beaufort Sea from 2 decades of MODIS data

Arctic sea ice cover evolves seasonally from large plates separated by long, linear leads in the winter to a mosaic of smaller sea ice floes in the summer. The interplay between physical and thermodynamic mechanisms during this process ultimately creates the observed sea ice floe size distribution (FSD), which is an important metric for characterizing the sea ice cover and assessing model performance. Historically, the FSD has been studied at fixed locations over short periods, leaving a gap in our understanding of the spatial and temporal evolution of the FSD at large scales. Here, we present an automated framework for image segmentation, allowing the identification and labeling of individual ice floes in Moderate Resolution Imaging Spectroradiometer (MODIS) data. Using this algorithm, we automatically process and segment 4861 images, identifying more than 9.4 million floes over 23 years. The extracted characteristics of the floes – including area, perimeter, and orientation – evolve throughout the spring and summer in the Beaufort Sea. We find seasonal patterns of decreasing mean floe areas, increasing FSD power law slopes, and increasing variability in the floe orientation as the summer progresses.

1 Introduction

Arctic sea ice cover controls heat and moisture flux between the atmosphere and the ocean. It has an annual cycle characterized by the growth and melt of ice in which large, heterogeneous snow-covered winter ice floes fragment into an ensemble of smaller floes in the summer. The breakup of the ice cover and the resulting floe size distribution (FSD) are set by complex feedback loops involving physical and thermodynamic processes. For instance, given the higher perimeter-to-area ratios of small floes compared to those of larger floes, small floes experience proportionally more lateral melting. This effect is observed to be especially pronounced for floes smaller than 50 m in diameter (Steele, 1992; Horvat et al., 2016). Lateral melting shrinks floes, further raising their perimeter-to-area ratios and thus leading to a positive feedback cycle. Also, lateral melt creates density gradients that contribute to the non-homogeneous stratification of the upper-ocean mixed layer, enhancing mixing and eddy formation (Horvat et al., 2016). Long waves from summer Arctic storms fracture the ice pack (Asplin et al., 2012), an effect amplified by the retreating sea ice edge and lengthened open-water fetch enhancing wave energy. The FSD also yields information about how a sea ice field will respond to oceanic and atmospheric forcing. Since the drag coefficient between an ice floe and the ocean depends on the ice floe size, the FSD is related to ocean–atmosphere energy and momentum transfer (Birnbaum and Lüpkes, 2002).

Following the original conceptualization of the FSD in Rothrock and Thorndike (1984), numerous studies have documented FSDs in various regions throughout the Arctic and the Antarctic (see Stern et al., 2018b, for a comprehensive list of FSD studies). We focus on the Beaufort Sea, where FSDs have previously been determined from radar imagery (Holt and Martin, 2001; Hwang et al., 2017), high-resolution optical satellite imagery (Wang et al., 2016; Denton and Timmermans, 2022), and aerial photography (Rothrock and Thorndike, 1984; Perovich and Jones, 2014). While these studies have advanced our understanding of the FSD seasonal evolution (e.g., the effect of storms on floe breakup and the relationship between sea ice concentration (SIC) and FSDs), they are limited to small areas over short periods.

Figure 1 Study region. (a) True-color reflectance MODIS imagery of the Beaufort Sea on 22 May 2021: the study region is outlined in a green box. (b) Location density of observed floes in the study region. The 2-D histogram shows the total number of floes observed in each 1° × 1° box. (c) Distribution of observed floes by year.

Moderate Resolution Imaging Spectroradiometer (MODIS) data have been used previously in FSD studies (Toyota et al., 2016; Zhang et al., 2016; Stern et al., 2018a), but only for a short term of a few years at most. The Ice Floe Tracker (IFT) algorithm was developed (Lopez-Acosta et al., 2019) to segment MODIS images and track floes between consecutive images to evaluate sea ice–ocean interactions but does not prioritize the capture of the full FSD. We develop a new algorithm based on the work by Denton and Timmermans (2022) focusing on high identification rates of ice floes and retrieval of accurate geometric properties. With 23 years of available data in the spring-to-summer transition, the MODIS optical imagery provides the potential to study interannual and decadal changes in floe properties and analyze whether there is a change in the observed FSD in the Arctic. The paper is organized as follows. First we present our algorithm for image segmentation and floe identification in optical satellite imagery in a wide range of ice concentration and melt states. We then present and validate the extracted sea ice floe properties using higher-resolution data and existing datasets. We also discuss the algorithm and the data limitations and uncertainty in the segmentation algorithm output. Then, the floe identification algorithm is applied to thousands of MODIS images spanning from 2000 to 2022 in March through September in the Beaufort Sea. The seasonal evolution and interannual variability of the floe sizes are presented and discussed. We conclude the study with an overview of our findings and suggestions for future directions.

2 Study area and data

The Beaufort Sea has experienced profound changes over recent decades. The end-of-summer sea ice area, measured as the September monthly average, is decreasing at a rate of 10 %–30 % per decade, resulting in substantially more solar heating in the upper ocean (Timmermans and Toole, 2023). There is a loss of older and thicker multiyear ice in the Beaufort Sea (Maslanik et al., 2007; Kwok and Cunningham, 2010), and the multiyear ice edge is receding (Galley et al., 2016). Ice loss in the Beaufort Sea has influenced the Beaufort Gyre circulation and freshwater (Timmermans and Toole, 2023). We focus this study in the Beaufort Sea on capturing the changes in the ice cover to further understand these important processes (Fig. 1a). We hypothesize that there is a quantifiable seasonal transition in the FSD from high-concentration winter sea ice to the melting-fractured sea ice of summer.

2.1 MODIS imagery

The MODIS instrument is on board two NASA satellites, Terra and Aqua, launched in 1999 and 2002, respectively. MODIS acquires data in 36 spectral bands spanning wavelengths from 0.4 to 14.4 µm with varying resolution. In this work we use the true-color reflectance imagery (Fig. 1a), a composite of Band 1 (red, 645 nm), Band 4 (green, 555 nm), and Band 3 (blue, 469 nm) (Vermote, 2015). This composite product is available at 250 m resolution. Full Arctic coverage imagery is provided four times daily outside of polar night (March through mid-October in the Beaufort Sea). The MODIS Cloud product (MOD06/MYD06, Platnick et al., 2016) reports calculated cloud properties for each pixel in the MODIS images and includes a cloud fraction derived from the infrared imagery. We use the cloud fraction parameter to mask out cloudy regions in the preprocessing step of the algorithm (Sect. 3.1).

2.2 Sentinel-2 imagery

The Sentinel-2 satellites A and B, launched in 2015 and 2017, carry the multispectral instrument, acquiring data in 13 spectral bands. Sentinel-2 imagery is captured approximately twice daily in the Arctic and is available up to 20 km off the coast. The Sentinel-2 Level-1C Top of Atmosphere reflectance product includes four bands that provide data at 10 m resolution: Band 2 (blue, centered at 492.3 nm), Band 3 (green, 558.9 nm), Band 4 (red, 664.9 nm), and Band 8 (near-infrared, 832.9 nm) (Drusch et al., 2012). In this work we use the high-resolution Sentinel-2 data to understand the limitations of the lower-resolution MODIS imagery.

3 Image segmentation methodology

Paget et al. (2001) describe an erosion–expansion algorithm that erodes the boundaries of floes to separate them, and it subsequently regrows them to their original shape. Denton and Timmermans (2022) and Stern et al. (2018a) introduce an iterative procedure that cycles through rounds of erosion and expansion, varying the amount of erosion to identify floes of different sizes. Building on the work in Denton and Timmermans (2022), we develop a new algorithm for the identification of ice floes that capture FSDs in the MODIS dataset. We rewrite the algorithm in Python and automate it to process thousands of images consecutively, introducing adaptive thresholds. The algorithm consists of four steps: preprocessing, ice–water discrimination, segmentation, and postprocessing (Fig. 2).

Figure 2 Image segmentation methodology. (a) Flow diagram showing the steps of the algorithm: preprocessing (pink), ice–water discrimination (orange), segmentation (green), and postprocessing (blue). (b) The original MODIS true-color reflectance imagery in the region of interest. The blue square shows the location of panels (d) and (e). (c) The MODIS true-color reflectance data with a mask (black) over the land and cloudy regions. (d) Ice–water discrimination results showing ice (white) and water or a mask (black). (e) Final segmented image showing individual floes in different colors (superimposed onto the original optical image). The example MODIS image shown here is from 10 June 2017.

3.1 Preprocessing

We use the MODIS true-color reflectance data as a starting point for our analysis. To ensure that we process only clear-sky ocean imagery, we mask the land and any cloud-covered areas. Atmospheric noise can blur the edges of the floes. Hence, we are conservative with cloud cover and use the existing MODIS cloud fraction data to eliminate areas that may be affected by clouds. We mask the areas with ≥95 % cloud coverage according to the MODIS cloud fraction product (Sect. 2.1, Fig. 2c). The presence of clouds in the Arctic in the summer is ubiquitous: on average 58 % of an image is covered in opaque clouds and unavailable for image processing (see the discussion in Sect. 4 and Fig. 5a, b).

3.2 Ice–water discrimination

The images cover a large area with a range of ice conditions, i.e., from the ice edge into the pack ice region, and we cannot apply a simple threshold pixel value to distinguish the bright ice from the dark open-water pixels (as was done by Denton and Timmermans, 2022). The thin leads that separate the large floes in the spring may be covered by a thin layer of ice increasing the brightness of the lead beyond the typical brightness values of the open ocean. Similarly, the ubiquitous low-lying fog in the Arctic summer may locally brighten the appearance of open water. Considering the variable pixel values of open water, we apply an adaptive threshold to determine local values for the ice–water discrimination. At this point, the land and cloud pixels are masked and are not considered in this step. We apply this method to the red-channel band of the MODIS imagery, which exhibits the highest contrast in pixel values. The dynamic threshold value is the weighted mean for the 399×399-pixel neighborhood (∼100 km × ∼100 km). All pixels with brightness greater (less) than the threshold are identified as ice (open water). In this way, we are able to account for varying brightness levels of open water given different ice concentrations and atmospheric conditions.

3.3 Image segmentation

The image segmentation step follows an erosion–expansion routine similar to that described by Paget et al. (2001) and subsequent studies segmenting airborne and satellite imagery of sea ice floes (e.g., Denton and Timmermans, 2022; Steer et al., 2008). The input to the segmentation routine is the binary classified image created in the previous step (Sect. 3.2), where 1 represents sea ice and 0 represents open water. The morphological erosion operation is applied to the binary image, removing pixels on the object boundaries using a diamond-shaped structuring element with a radius of 1 px. The binary image is eroded a total of eight times in order to ensure that floes are separated. This extensive erosion removes small floes from the image. At this point, the remaining distinct floes are tagged and then regrown (dilated) to their original state. Any identified floes touching the image border, the land mask, or the cloud mask are removed from the binary image. The remaining tagged floes are saved in the floe library and then removed from the image, i.e., changed to 0, and the next iteration round begins for identification of smaller floes. This erosion–tag–expansion process is repeated, with fewer erosions each time, allowing for subsequent smaller floes to be identified with each iteration. In this way, sea ice floes of varying sizes are separated and identified. At the end of the image segmentation routine, an image with each unique object labeled is produced. Over the thousands of processed images, on average, 26 % of the classified sea ice area is identified as individual floes, with the remaining sections consisting of ice filaments, brash ice, or pieces of ice smaller than the minimum detectable floe size. Here, we processed the years 2000–2022 from days of year 60–274 (approximately 1 March–30 September, encompassing the time of year when light levels are sufficiently high for optical imagery) and segmented 4861 images and days, thereby identifying 9 448 563 floes.

3.4 Postprocessing

The following geometrical parameters are calculated for each identified floe: centroid position, floe orientation, area, perimeter, major and minor axes, circularity, and intensity mean. The major (minor) axis is the length of the major (minor) axis of an ellipse with the same normalized second central moments as the identified floe shape. The intensity mean is the average red-channel value for the area of the floe. Floe orientation is defined as the angle of the major axis of the floe from polar stereographic north. Orientation values range from $-\frac{\mathit{\pi }}{\mathrm{2}}$ to $\frac{\mathit{\pi }}{\mathrm{2}}$. The circular standard deviation of the orientation is calculated to represent the variability of the floe orientation, with low values representing the floe alignment. We also calculate the circularity of a floe:

$$\begin{array}{}\text{(1)}& C={\displaystyle \frac{\mathrm{4}\mathit{\pi }A}{P^{\mathrm{2}}}},\end{array}$$

where A is the floe area, P is the floe perimeter, and a circle has maximum circularity C=1.

Finally, we examine the properties of each identified object to ensure that it is a floe. We examine the red-channel pixel-intensity mean of the floes and discard floes with a value of less than 150, which is an empirically determined value. Image pixels on ice floes have high red-channel values, and objects with a low-intensity mean may be incorrectly identified as floes. Rather, low-intensity values may indicate clusters of brash ice, for example. This quality assurance step results in the elimination of 86 501 floes (<1 % of the total number of floes).

Figure 3 Validation of MODIS image segmentation with higher-resolution (10 m) Sentinel-2 imagery. Each row shows spatially coincident MODIS (first column) and Sentinel-2 (second column) imagery captured on the same day. The segmentation algorithm is applied to each of these images, and the individual identified floes are colored in the image. The third column shows the probability density function (PDF) of the FSD for the Sentinel-2 (blue) and MODIS (red) images, with the best fit power law shown as the dashed line. The three rows are for imagery captured on 4 September 2019 (40 % SIC), 12 June 2020 (70 % SIC), and 14 May 2021 (98 % SIC).

Figure 4 Area comparison of identifications in MODIS and Sentinel-2 imagery (Fig. 3). The 82 floes identified in both datasets are shown, with the dot color corresponding to the date of the images as indicated in the legend.

3.5 Floe size distribution

The FSD contributes to the characterization of the ice floe field by providing a quantitative description of the ice floe area statistics. Taken together with the floe geometry properties outlined in Sect. 3.4, the FSD and other parameters commonly used to describe floe fields, such as the SIC and average ice thickness, allow us to study the physical processes that shape the structure and evolution of sea ice. We utilize the powerlaw Python package (Alstott et al., 2014) based on the maximum likelihood estimation power-law-fitting methods described by Clauset et al. (2009) and Klaus et al. (2011). The noncumulative power law is described by

$$\begin{array}{}\text{(2)}& p\left(x\right)=c{x}^{-\mathit{\alpha }},\end{array}$$

where p(x) is the probability of a given instance of x (the chosen geometric property of the sea ice floe – we use floe area), c is a normalization constant ensuring that the function integrates to 1, and α is the fitted parameter and slope of the power law distribution. We specify minimum and maximum floe size values, x_min=5 km² and x_max=300 km², based on the range of observed floe sizes and the goodness of fit of a power law distribution to this range. Approximately 97 % of the floes fall in this size range. Setting x_min and x_max also allows for a finite integration of the power law. With given x_min and x_max values, the constant (c) is given by (see Appendix A)

$$\begin{array}{}\text{(3)}& c={\displaystyle \frac{\mathrm{1}-\mathit{\alpha }}{x_{\max }^{\mathrm{1}-\mathit{\alpha }}-x_{\min }^{\mathrm{1}-\mathit{\alpha }}}}.\end{array}$$

The standard error of α is defined by Clauset et al. (2009) as

$$\begin{array}{}\text{(4)}& \mathit{\sigma }={\displaystyle \frac{\mathit{\alpha }-\mathrm{1}}{\sqrt{n}}}+O(\mathrm{1}/n),\end{array}$$

where n is the sample size and the higher-order correction is positive. For these equations, α must be greater than 1. We calculate the FSD power law fit for different datasets to understand how the ice cover evolves, seasonally and annually.

3.6 Validation

We evaluate the consistency of the algorithm by examining the floes extracted from Aqua and Terra images on the same day. The images cover a large region and the acquisition times of Aqua and Terra are approximately 2 h apart, so we can assume that the images cover the same expanse of ice in same-day acquisitions. We can expect a similar power law distribution of identified floe sizes on the same day, albeit from different satellites, given that they carry the same MODIS instrument. We randomly selected 100 d from the 23-year collection of images to examine both Aqua and Terra images. We segment both images, calculate floe properties, and match floes based on the centroid location. We find a correlation value of 0.99 for the matched floe areas. For each pair of images where each image has at least 50 000 km² of identified floes, we fit a power law to the floe area distribution (Sect. 3.5) and determine the slope (α). We find an absolute mean difference in the Aqua and Terra α values of 0.009 and a standard deviation of 0.006. This suggests strong agreement between the floes identified in the Aqua and Terra satellite imagery, confirming that the segmentation algorithm is consistently identifying floes.

To understand and quantify the limitations of the moderate-resolution imagery used by the segmentation algorithm, we apply the algorithm to higher-resolution 10 m Sentinel-2 imagery. We examine spatially coincident MODIS imagery and Sentinel-2 imagery observed on the same day (Fig. 3). The extents of the two images are matched, so that the same area is analyzed. We evaluate imagery in a range of sea ice conditions: low SIC (40 %) seen at the end of the summer (Fig. 3a–c), ice in the marginal ice zone in summer (70 % SIC, Fig. 3d–f), and high SIC (98 %) at the beginning of the melt season (Fig. 3g–i). The SIC was determined from the ice–water discrimination step in the algorithm. We pair floes identified in the coincident imagery based on the centroid location and examine the corresponding floe properties. The areas of the 82 matching floes agree well, with a squared correlation of 0.99 and an absolute mean area difference of 0.18 km² (Fig. 4). The 82 floes form 25 % of the Sentinel-2 floes and 21 % of the MODIS floes in the observable floe area range for MODIS data (x_min=5 km², x_max=300 km²). Despite a low percentage of matching floes, the identified floes in both sets of images are a good representation of the floe areas in each image.

We fit a power law to the FSD for all floes within the floe range of each of the images (Fig. 3c, f, i). The FSDs for each pair of Sentinel-2 and MODIS images agree well. The difference in α values ranges from 0.06 to 0.25, with 1σ confidence intervals overlapping between Sentinel-2 and MODIS distributions for the low- and high-SIC instances (Fig. 3c and i), and there is a 2σ overlap for the medium-SIC example (Fig. 3f). Note that the identification of floes in the Sentinel-2 imagery was not limited to the MODIS range of floe sizes, and approximately 82 % of the floes identified in the Sentinel-2 imagery are less than the determined x_min value for the MODIS imagery. In the overlapping range of floe areas (x_min=5 km², x_max=300 km²), 14 % more floes are identified in the MODIS imagery compared to the Sentinel-2 imagery. The higher resolution of Sentinel-2 allows for identification of smaller floes, but the segmentation algorithm performs similarly to the lower-resolution imagery for larger floes. In the high-SIC scenario (Fig. 3g–i), the narrow leads seen in both images make the separation of floes more challenging. The segmented floes seen in panels (g) and (h) look different, but the agreement of the FSD α values indicates that the floes identified in the images are a good representation of the ice cover. We note that, in high-concentration regions, there is more uncertainty due to the difficulty in identifying individual floes separated by narrow leads. Despite the resolution limitations of the MODIS imagery, our validation shows that the algorithm applied to the MODIS imagery samples the floe sizes sufficiently to produce an accurate FSD α value and that the floes identified in both images have highly correlated floe areas.

Figure 5 Image statistics. Panels (a) and (b) show cloud fraction statistics from the cloud mask. Panels (c) and (d) show the total ice area identified from the ice–water discrimination (orange) and the total area of the identified ice floes (blue). The gray line shows the number of identified ice floes. Panels (a) and (c) show statistics over the season, and panels (b) and (d) show statistics over the years.

Finally, we apply a bootstrapping approach to quantify uncertainty in the derived parameters. We begin by randomly selecting 100 segmented MODIS images (corresponding to 100 different days), where each image has at least 50 000 km² of identified floes. We then create 1000 bootstrap datasets of floes for each of the 100 images, where each bootstrapped sample has the same number of identified floes (N) as the original image. This is done by randomly selecting N floes from an image, such that after each floe is selected it is returned to the original image (i.e., sampling with replacement). We then calculate the power law distributions of the bootstrap datasets and the standard deviation of the α values over the 1000 bootstrap datasets for each of the 100 images. The standard deviation of α generated from the bootstrapping of an image is on average 0.024, ranging from 0.007 to 0.08. The standard deviation increases as the α value increases, indicating more uncertainty.

3.7 Data and algorithm limitations

We note that there are limitations to our analysis due to the moderate resolution of the imagery. The separation of floes requires openings between floes (such as leads) that are at least the size of an image pixel (250 m), and typically multiple pixels are required to fully resolve a lead. This results in multiple floes being considered a single floe (examples can be seen in Fig. 3g). This phenomenon is especially prevalent in the early spring, when floes are still tightly packed and have not experienced lateral melt and floe divergence. In addition, larger floes are more likely to intersect the image border. Floes that intersect the image border are removed from the identified floes as the properties are not correct, and thus large floes are preferentially eliminated. For these reasons, we limit the range of the power law fit and do not draw conclusions about the power law distribution of floes greater than 300 km². Nonetheless, the error associated with the FSD is largest when the SIC is high in the early spring.

We chose to analyze MODIS imagery due to its long record and consistent coverage, but other higher-resolution imagery is required to examine the FSD for floes smaller than 5 km². Other studies have examined FSD in synthetic aperture radar (SAR) imagery, which eliminates the need for clear-sky conditions as SAR is not sensitive to clouds. However, SAR data are not as widely available for long-term applications and have other limitations and complications, such as speckle, granular noise, and ambiguous returns when meltwater is present on the ice surface.

Figure 6 The 1 March through 30 September evolution of the floe properties. The properties were aggregated into a 10 d running window and sampled every 5 d. (a) Mean floe area for the observed (solid) and analytical (dashed) values (see Appendix A2). The 75th percentile of the floe area is also shown (dark blue). (b) Median floe area for the observed (solid) and analytical (dashed) values (see Appendix A3). (c) Power law slope (α). (d) Standard deviation of the orientation of floes for all floes.

4 Spring-to-summer transition of floe characteristics

After processing 4861 images, we analyzed the basic statistics of the classification and segmentation routine. Over 9.4 million floes were identified, and the distribution of their locations is shown in Fig. 1b. Note that these are not 9.4 unique floes as each image is taken as an independent observation. The MODIS cloud fraction in the study region is consistent with the Arctic-wide pattern of the highest cloud fraction in the summer and fall (Fig. 5) (Schweiger, 2004). Although this study is not focused on atmospheric trends, we do not find a significant long-term trend in the cloud fraction over our study region and period (Fig. 5b). We note that our analysis of cloud-free areas only is limited due to a high percentage of cloud cover obscuring the sea ice in optical imagery. As expected, the total classified sea ice area and the total identified floe area in the imagery decrease throughout the summer as the ice melts (Fig. 5c, orange and blue). The largest difference between the ice area (orange) and the identified floe area (blue) is in the spring. There are a number of reasons why the majority of the ice cover is not able to be segmented into floes in the spring. The floes in the spring are larger (Fig. 6a) and therefore are more likely to intersect the border of an image and be eliminated to ensure that only floes fully captured by the imagery are used. Also, the ice floes are tightly packed in the spring and the MODIS imagery does not resolve small leads (Fig. 3h). Thus, floes cannot be separated as well compared to in the summer, when floes are separated by larger areas of open water (Fig. 3e). The number of identified floes (Fig. 5c, gray) is not correlated with the total area of observed floes (Fig. 5c, orange). More floes are identified as the ice separates in the summer, but as the ice melts and advects out of the study area, there is less ice and fewer floes are identified. There are no significant trends in the ice and floe areas or the number of floes identified over the 23 years of observations (Figs. 1c and 5d).

Figure 7 Observed floe properties compared to sea ice concentration (SIC). (a) Standard deviation of the orientation of floes on a given day compared to the image's SIC. (b) Mean floe area compared to SIC on a given day. The color bar indicates the total image count.

4.1 Floe area

In the spring, the Beaufort Sea has a high ice concentration consisting of large floes with rectilinear fractures. During the summer, as the ice edge recedes from the Alaskan and Canadian shorelines, the ice cover transitions to a dynamic collection of randomly oriented floes among brash ice. We analyze floe characteristics from 1 March through 30 September in the MODIS imagery and observe patterns corresponding to this transition. The mean floe area first increases from 32 km² on 6 March to the seasonal high of 36 km² on 10 April (Fig. 6a, blue). The mean floe area then decreases throughout the summer, plateauing in August around 21 km² and then increasing to reach a mean area of 23 km² by the end of September. The mean floe area closely follows the 75th percentile floe area (Fig. 6a, dark blue). The median floe area is consistently lower than the mean floe area (Fig. 6b, orange), as the distribution is positively skewed with many small floes. Because of the high skewness of the distribution, the standard deviation is large (38.9 km² on average) as it is influenced by floes much larger than the mean floe area. The median floe size exhibits a similar pattern to that of the mean, with a maximum median value of 16 km² on 10 April and decreasing to 11.5 km² on 30 September. We also see the largest floes in areas of the highest SICs (Fig. 7b). Large ice floes and high SICs exist in the early spring. Floes that exist in areas of low ice concentration are more likely to experience the effect of waves and may break up into smaller floes due to wave fracture (Squire et al., 1995).

We note that the magnitudes of the median and mean floe area values are sensitive to changes in the x_min and x_max values. However, the pattern of evolving floe sizes is related to the changing sea ice cover. The mean (expected value) and median (where the cumulative distribution function is 0.5) of the fitted power law can be expressed analytically (see Appendix A2 and A3). These values, with x_min and x_max set, are functions of the evolving α. We show the analytical mean and median as dot-dashed lines in Fig. 6a–b, respectively, which indicate that the analytical solutions are consistently less than the observed mean and median floe areas. This is because the power law fit is dominated by small floes, as there are more of them. Thus, there is greater error in the fit associated with the large observed floes.

Figure 8 α values of a power-law-fitted FSD for each month of analysis. Each strip shows the α value for all floes identified in that month over the 23-year analysis. The dotted line is the line of best fit, the slope of which is printed in the bottom-right-hand corner of each strip.

4.2 Power law fit to the floe size distribution

We fit a power law to the collection of floes grouped in a 10 d running window. The value of α (slope of the power law distribution) of all identified floes is 1.85 but ranges from 1.74 to 2.0 throughout the spring-to-summer period (Fig. 6c). The standard deviations and the confidence intervals of all the power law fits are less than 0.01. The power law α is inversely correlated with the mean floe area. The slope value decreases from the beginning of March to the minimum on 9 April, then increases throughout the summer to a maximum on 12 August, and slowly begins to decrease again for the remainder of the summer. The inverse relationship with the mean floe area is expected (see Appendix A2), as a higher frequency of small floes will decrease the mean floe area and increase the slope of the FSD fit. We find the magnitude and seasonal trends of α to be consistent with previous studies.

Stern et al. (2018a) examined 116 MODIS images in 2013 and 2014 and fit a power law to floes ranging from 2 to 30 km in size, which is approximately equivalent to 2.64 to 594 km² (Rothrock and Thorndike, 1984). To compare results, we must account for the fact that we have fit a power law to the floe area, while Stern et al. (2018a) found a power law slope (α_n) using the mean caliper diameter. This value is related to our reported α as (see, e.g., Denton and Timmermans, 2022)

$$\begin{array}{}\text{(5)}& {\mathit{\alpha }}_n=\mathrm{2}\mathit{\alpha }-\mathrm{1}.\end{array}$$

Stern et al. (2018a) found that α_n (α) was approximately 2.0 (1.5) in May, increased to about 2.9 (1.95) in July, and then decreased to about 2.2 (1.6) by October. We find comparable slopes and observe a similar evolution with an α value of about 1.75 in May, a maximum α in August of 2.0, and then a slight decrease in September (Fig. 6c). The small differences in values are likely due to the locations and times of the observed floes. Denton and Timmermans (2022) analyzed smaller floes (50 m² to 5 km²) in the Canada Basin and found a seasonal trend in the FSD, with α ranging from 1.65 to 2.03.

The FSD trend is also consistent with previous studies that have taken a Lagrangian approach, tracking the same ice throughout the summer. During the Surface Heat Budget of the Arctic (SHEBA) campaign in 1998 in the Beaufort and Chukchi seas, aerial photography was collected in the proximity of the ship. Perovich and Jones (2014) found that α values increased throughout the summer, reached a maximum on 10 August, and subsequently decreased into September as small floes froze and fused into larger floes (Perovich and Jones, 2014). Hwang et al. (2017) observed the FSD from satellite SAR imagery tracking four buoys in the Beaufort Sea in 2014 and capturing images of the same ice, finding an increase in FSD α values from July through August with enhanced floe breakup linked to wind events.

We examine the FSD power law slope over the 23-year MODIS record to quantify the interannual variability and discern any decadal trends. As the FSD exhibits seasonal variability, we look at the monthly FSD power law slopes and compare them over the years. We see no significant trends in the monthly FSD slopes over the 23-year period (Fig. 8). The error in all the calculated α values is <0.01, as determined by 1000 bootstrap samples to calculate a 95 % confidence interval and the standard error calculated as in Eq. (4). The month of September exhibits the most variability, with α values ranging from 1.85 to 2.17 (Figs. 6c and 8). One may expect a trend towards larger α values (steeper FSD power law slopes) as the Beaufort Sea ice cover exhibits earlier retreat (Fetterer et al., 2017) and a transition to first-year ice that is more susceptible to fracture (Galley et al., 2016). However, this is not exhibited in our data. This lack of a significant decadal trend on the FSD slope may be due to the large study region that simultaneously contains pack ice and open water, or, because MODIS cannot resolve floes smaller than 5 km² in area, there may be many more small floes in recent years that are not identified in the MODIS data.

4.3 Floe orientation

In the late winter and spring, the Beaufort Sea has a high SIC. As external wind stresses are applied to the ice pack, the ice pack experiences strain, leading to a fracture in a preferential direction that depends on the orientation of the force relative to the coast (Lewis and Hutchings, 2019; Jewell et al., 2023). When we examine the newly fractured floes during the spring, we find low variability in their orientation within the 10 d window (Fig. 6d). This effect is especially noticeable in areas of high ice concentration, where the ice movement and readjustment to external forces are limited by the surrounding ice and thus lower standard deviation of orientation with a greater SIC (Fig. 7a). The small fractures, or cracks, grow into leads that can be seen in the MODIS imagery (e.g., Fig. 3g). As the summer progresses, SIC decreases and ice floes break up and disperse, resulting in a decreasing average floe area (Fig. 6a) and an increasing standard deviation of floe orientations with lower sea ice concentrations (Figs. 6d and 7). We also examine the floes with the highest eccentricity (>75th percentile), which is the ratio of the major axis to the minor axis. This subset of floes exhibits a similar but less extreme trend of an increasing standard deviation of orientation throughout the summer. The orientations of the floes are an indication of the stresses that have caused ice breakup and can provide insight into the structural properties of the ice pack to inform how it may respond to future stresses. When ice floes can rotate with minimal to no interaction with other floes, their rotation rates can be related to ocean vorticity (Manucharyan et al., 2022).

5 Conclusions

The algorithm developed in this study establishes the ability to derive meaningful floe size information from the longest daily global satellite observation record of Earth to date. The MODIS dataset has Arctic-wide coverage and spans 23 years. The expansive dataset combined with the modified segmentation algorithm presented here allows for continued study of ice conditions and characteristics in the spring and summer. To validate the feasibility of using moderate-resolution imagery, the identified floes are validated visually and compared with floes identified in higher-resolution Sentinel-2 imagery. We demonstrate that, in the specified range of floes with areas greater than 5 km², the segmentation algorithm performs equally well on images of moderate resolution (MODIS) and high resolution (Sentinel-2). Although the algorithm is able to identify smaller floes in the Sentinel-2 imagery, the retrieved floe sizes and the FSD agree well between the Sentinel-2 and MODIS imagery. This allows us to confidently apply this systematic method to analyze thousands of images over 23 years covering a wide range of ice conditions.

We examine the seasonal evolution of the FSD and aspects of the floe geometry in the Beaufort Sea from mid-April through early August. We find a decrease in the mean and median floe size, increasing α (steepening power law slope), and an increase in the variability in the orientation of the floes. As the sea ice cover appears and behaves significantly differently depending on the time of year, it is essential to use floe characteristics from the specific time period of interest when evaluating or tuning models. While no significant decadal trends were observed in the monthly FSD over the 23-year period, future work considering smaller spatial domains may be necessary to investigate interannual and decadal variability in detail.

Expanding the analysis to new areas of the marginal ice zone may show regional differences in the floe characteristics and timing of floe breakup in the summer. Combining this new information with existing satellite measurements (e.g., ice drift, ice type, or ice thickness) can provide further insights into the behavior of the ice pack. We chose to analyze MODIS imagery due to its long record and consistent coverage, but other higher-resolution imagery is required to analyze how the FSD holds for small floes. It may also be of interest to examine the floe evolution throughout the winter, which would require active sensors that can produce images when there is no sun illumination available. Further work will include incorporating this image segmentation algorithm into the pipeline of the IFT algorithm (Lopez-Acosta et al., 2019). This routine tracks floes with similar characteristics between consecutive MODIS images and can thus determine ice velocities and rotation rates, inferring ocean dynamics in regions that are otherwise under-observed. With a new segmentation algorithm able to identify floes in a wider range of sizes, we can expand the IFT output to also include the FSD and uncover more information about the underlying ocean.

Appendix A: Mathematical derivations for a power law floe size distribution

A1 Calculation of the power law constant

To solve for c in the power law distribution (Eq. 2), we integrate the power law from x_min to x_max, noting that the sum of all probabilities in the range is 1:

$$\begin{array}{}\text{(A1)}& \underset{x_{\min }}{\overset{x_{\max }}{\int }}c{x}^{-\mathit{\alpha }}\mathrm{d}x={\displaystyle \frac{-c{x}^{\mathrm{1}-\mathit{\alpha }}}{\mathit{\alpha }-\mathrm{1}}}{|}_{x_{\min }}^{x_{\max }}=-{\displaystyle \frac{c{x}_{\max }^{\mathrm{1}-\mathit{\alpha }}}{\mathit{\alpha }-\mathrm{1}}}+{\displaystyle \frac{c{x}_{\min }^{\mathrm{1}-\mathit{\alpha }}}{\mathit{\alpha }-\mathrm{1}}}=\mathrm{1},\end{array}$$

where α>1.

This yields

$$\begin{array}{}\text{(A2)}& c={\displaystyle \frac{\mathrm{1}-\mathit{\alpha }}{x_{\max }^{\mathrm{1}-\mathit{\alpha }}-x_{\min }^{\mathrm{1}-\mathit{\alpha }}}}.\end{array}$$

A2 Analytical mean

The expected value of a distribution is

$$\begin{array}{}\text{(A3)}& E=\underset{x_{\min }}{\overset{x_{\max }}{\int }}xp\left(x\right)\mathrm{d}x.\end{array}$$

With Eq. (2), this is

$$\begin{array}{}\text{(A4)}& E=\underset{x_{\min }}{\overset{x_{\max }}{\int }}c{x}^{\mathrm{1}-\mathit{\alpha }}\mathrm{d}x,\end{array}$$

where c is given by Eq. (A2).

Integrating yields

$$\begin{array}{}\text{(A5)}& E={\displaystyle \frac{c}{\mathrm{2}-\mathit{\alpha }}}\left(x_{\max }^{\mathrm{2}-\mathit{\alpha }}-x_{\min }^{\mathrm{2}-\mathit{\alpha }}\right),\end{array}$$

where E depends only on α because x_min and x_max are set as 5 and 300 km², respectively. The full equation, substituting c from Eq. (A2), is

$$\begin{array}{}\text{(A6)}& E={\displaystyle \frac{\mathrm{1}-\mathit{\alpha }}{\mathrm{2}-\mathit{\alpha }}}{\displaystyle \frac{x_{\max }^{\mathrm{2}-\mathit{\alpha }}-x_{\min }^{\mathrm{2}-\mathit{\alpha }}}{x_{\max }^{\mathrm{1}-\mathit{\alpha }}-x_{\min }^{\mathrm{1}-\mathit{\alpha }}}}\end{array}$$

for α ≠ 2.

A3 Analytical median

The median of a distribution x_o is found where the cumulative distribution is equal to 0.5. That is,

$$\begin{array}{}\text{(A7)}& \underset{x_{\min }}{\overset{x_o}{\int }}p\left(x\right)\mathrm{d}x={\displaystyle \frac{c}{\mathrm{1}-\mathit{\alpha }}}{x}^{\mathrm{1}-\mathit{\alpha }}{|}_{x_{\min }}^{x_o}=\mathrm{0.5}.\end{array}$$

Solving for x_o yields

$$\begin{array}{}\text{(A8)}& x_o={\left[x_{\min }^{\mathrm{1}-\mathit{\alpha }}+{\displaystyle \frac{\mathrm{1}-\mathit{\alpha }}{\mathrm{2}c}}\right]}^{\frac{\mathrm{1}}{\mathrm{1}-\mathit{\alpha }}},\end{array}$$

which, with Eq. (A2), yields

$$\begin{array}{}\text{(A9)}& x_o={\left[{\displaystyle \frac{\mathrm{1}}{\mathrm{2}}}(x_{\min }^{\mathrm{1}-\mathit{\alpha }}+x_{\max }^{\mathrm{1}-\mathit{\alpha }})\right]}^{\frac{\mathrm{1}}{\mathrm{1}-\mathit{\alpha }}}.\end{array}$$