Characterizing the size and shape of sea ice floes

Monitoring drift ice in the Arctic and Antarctic regions directly and by remote sensing is important for the study of climate, but a unified modeling framework is lacking. Hence, interpretation of the data, as well as the decision of what to measure, represent a challenge for different fields of science. To address this point, we analyzed, using statistical physics tools, satellite images of sea ice from four different locations in both the northern and southern hemispheres, and measured the size and the elongation of ice floes (floating pieces of ice). We find that (i) floe size follows a distribution that can be characterized with good approximation by a single length scale , which we discuss in the framework of stochastic fragmentation models, and (ii) the deviation of their shape from circularity is reproduced with remarkable precision by a geometric model of coalescence by freezing, based on random Voronoi tessellations, with a single free parameter expressing the shape disorder. Although the physical interpretations remain open, this advocates the parameters and as two independent indicators of the environment in the polar regions, which are easily accessible by remote sensing.

All images were segmented by the Potrace algorithm, and contour information for each detectable floe was gathered. More details about the satellite imagery and the data analysis can be found in the Methods section. We obtained approximately 27 000, 18 000, 16 000, and 7000 contours for the four images respectively. This data set is limited, but attempts to capture different locations and different seasons, primarily in the marginal ice zones. We regard it as a test case for our observations and methods. Larger sets would be needed for developing more realistic physical models.
First, we focus on the linear size of the floes. The segmentation procedure yields each contour ω as a closed polyline, identified by a set of ω n nodes ω ∈ i 2  with = , …, ω i n 1 . To extract a linear measurement of floe size, we consider the square root of the mean square distance of all contour points from the center of mass (an analog of the "radius of gyration" used in polymer physics). We weigh each point by Scientific RepoRts | 5:10226 | DOi: 10.1038/srep10226 its distance from the preceding one, in order to compensate for the nodes being spaced unevenly along the contour. This defines a floe "size" l as Note that measuring the linear size via the perimeter as π / p , which is more sensitive to the roughness of the contours, does not affect heavily the results presented hereafter.
The measurement of l for all floes in the four data sets gives the four corresponding size density distributions ( ) P l . Inspection of the curves indicates that a power-law regime is present in all data sets (with exponent ≈ − 2), followed by a smooth cutoff for large floes. We found that these curves collapse rather well onto a single curve by a simple rescaling of the two axes (Fig. 2). This suggests that each curve is characterized by a single length scale. The small-size drop-offs visible in the plot are due to the underestimation of floes of small size, close to the resolution of the images, and are therefore not universal. Unfortunately, this purely technical feature makes it impossible to follow an accurate approach to data collapse, in the spirit of 25 ; the rescaling has to be eyeballed. However, albeit of an empirical nature, the collapse is remarkable, considering that the data come from diverse locations and seasons, and involve broadly dissimilar scales.

Fragmentation theory for the size distribution
To rationalize the observed scaling behavior, we consider statistical models of fragmentation. We briefly review two classes of existing models; the first is a stochastic dynamics for the number of floes of a given size, the second is a deterministic geometric description of crack propagation in brittle fracture 26,27 . Then we introduce a third alternative model, and argue that the emergence of a single characteric length scale is compatible with all these scenarios. In the first class of models, the physical aspects of the fragmentation process are summarized by a rate π ( ) l at which a floe of size l breaks down into two smaller objects. At a mean-field level description, each floe is supposed to experience the average environment, i.e. spatial fluctuations are neglected. We consider a system in which the fraction of floes having linear size l at time t is ( ) P l t . The distribution P obeys the following rate equation 28 : the first term is due to floes of size l breaking down into smaller floes, the second term is due to larger floes of size λ generating fragments of size l (in the 2 possible ways), the factor λ / 1 has geometric origin (in its absence, the in-flow from larger floes would be rescaled by a factor proportional to their size); the last term ( ) s l is a generic source, which may represent for instance the generation of large floes by the fracture of ice fields.
When no source is present, ( ) = s l 0, the fragmentation dynamics (3) is known 29 to produce a scaling size distribution φ ( ) = ( / ( )) P l l s t t , where ( ) − / s t t : 1 2 is the decreasing typical floe size. If the breaking rate scales as the area of the floe, π ( ) ∼ l l 2 , then the scaling function, at large x, is φ ( ) = ( − ) − x x ax exp 2 2 , which reproduces the exponent −2 exposed in Fig. 2, and also predicts a universal exponential cutoff. An interpretation based on such a dynamical picture links the characteristic scale of floes to the time passed since the beginning of the fragmentation process, that is of the melting season. The data considered here were not sufficient to test these correlations and further investigations are needed. The second class of available models is that of brittle fracture. Fragment size distributions are obtained by considering the propagation of cracks in the material. Older theories 30 considered smooth cracks, while more recently it has been recognized that beyond a critical speed they can become unstable and branch 31 . Moreover branches attract each other and therefore merging of cracks has been included in the picture. A well-established theory 32,33 yields a fragment size distribution dominated by a hierarchical process, whereby a number of cracks stem from a main fissure and then merge, the longer-lived cracks giving rise to the larger fragments. In this description, again, the distribution of the linear floe sizes l is a power law of exponent −2 cut off by an exponential in l 2 . The upper cutoff represents the single length scale of the system. Additionally, we propose here a complementary formulation belonging to the first kind of models, which assumes that the empirical distributions could be the stationary states of processes happening at time scales much faster than the seasonal variations. The advantage of this description is that it provides a natural setting for the observation that the same scaling form agrees with empirical data taken in different periods of the year (compare e.g. data from the Svalbard area, taken in June, with those from the Kara Sea, taken in March).The steady state ( ) ∞ P l is obtained by setting the time derivative of ( ) P l t to zero in Eq. (3). Note that technically the presence of a source ( ) s l is necessary in order to have a steady state at all, since its absence would require an unbounded distribution ( ) ∞ P l . Setting ∂ ( ) = P l 0 t t in (3) and taking the derivative with respect to l yields an ordinary differential equation whose solution is The rescalings on the x-axis were adjusted by eye from the cutoff points, and the normalization constants corrected accordingly. Circles, squares, triangles, and rhombi correspond to data sets 1, 2, 3, and 4, respectively. The dashed line shows a power law of exponent −2 as a reference. The dropoff at small sizes is due to the finite resolution of the images; data in this region are discarded in Fig. 3.
We consider the case where the source generates floes of a fixed size l with a rate ν, i.e. νδ ( ) = ( − ) . s l l l This includes the case of no source (ν = 0). The solution, away from = l l, : the size dependence of the rate π ( ) l completely determines ( ) ∞ P l . In order to proceed, one has to specify π ( ) l . We present here a generic argument for the scaling of π ( ) l , which is suggested by the yearly process of sea-ice refreezing 5,18 whence water interfaces between different floes are rejoined. We make the simplifying assumption that the only relevant scale is that of the coalesced fragments, δl. Supposing that a floe of linear size l, resulting from the coalescence of δ ( / ) l l 2 smaller floes, breaks along the junction lines between them, then it can fracture in roughly different ways. If these are admitted to be independent then the breaking rate will be π π where π 0 is a constant. The form of the breakage rate resulting from this argument can be used in Eq.
(4) to obtain a steady-state solution for the floe-size distribution. In particular, the steady state solution in the absence of source is where ( ) Z l min is a normalization constant. The divergence in = l 0 is not integrable, therefore a cutoff > l 0 min needs to exist. Physically, this cutoff can represent an elementary scale under which fragments cease to divide (this happens for instance in dust aerosols 11 ), or under which melting becomes dominant. In the analysis of the satellite images considered here, the lower cutoff is not physical, but corresponds to the resolution of the images (pixel size), as can be seen in Fig. 1, where tiny floes are visible that are not revealed by the algorithm. The model-that is π ( ) l , and therefore ( ) ∞ P l -is then characterized essentially by the single scale δl (since the lower cutoff only affects the normalization).
It should be stressed that all models are necessarily simplified descriptions and neglect some of the phenomena caused by the very complex interaction between sea ice and the environment 34 . A most important one is melting due to the seasonal temperature variations. Lateral melting 35 would be represented by a term ∂ ( ) P l l in (3), thus possibly breaking the scaling form of the stationary state. The regularity observed can be then interpreted as a measure of the marginality of lateral ice melting as opposed to fracturing, in the regimes considered. Note that melting is supposed to become more relevant at smaller scales; our data suggest that higher-resolution images are necessary to this aim. Indeed, aerial photography suggests that a different behavior sets in at small sizes 36,37 Moreover, the source term could introduce deviations from the source-free steady state, as expressed by (4), which could be detected by studying the size distribution in the marginal ice zone close to the fracturing pack ice 22 . Additionally, these "mean-field" descriptions disregard floe shape, as well as the isotropy of fissures 20,38 .
Summing up, the three models described above are compatible with the following scaling form for the fragment size distribution, where the scaling variable is λ δ = / l l. This expression is in very satisfactory agreement with the empirical data. Fig. 3 compares the theoretical prediction (5) with the empirical data from our four data sets. Some previous studies claimed power-law regimes with exponents possibly deviating from −2 36,39 However, there is some debate around the numerical solidity, the interpretation, and the practical significance of this variability 37 . It is possible that some of the deviations from the scaling form (6) found in other studies came from the superposition of data having two or more length scales.

Characterization of floe shapes
The previous section has shown that the characteristic length δl may be subject to very different interpretations, and even the question on whether the scenario is dynamic or stationary seems left open. In this section we propose an additional observable, based on the shape of the floes, which is, as we will show, independent of their size.
When considering contour data, the gyration radius used in the definition of l in Eq. (1) is only one of several scalar quantities that can be constructed starting from a tensorial object, called gyration tensor. The gyration tensor, in the simple case of a set of point particles with unit mass, is exactly the inertia tensor, describing the rotational degrees of freedom of the system (the points are supposed to have fixed relative positions). It provides a compact description of some properties related to the shape of contours and lattice walks, and is fruitfully employed for instance in polymer physics [40][41][42][43] Here we use a slightly modified version, which takes into account the different "masses" corresponding to the different step lengths between the contour points ω i (see also the definition of l above). In a fixed coordinate system ( , ) x y , the gyration tensor is a symmetric matrix αβ Q , where α and β can take values corresponding to the two coordinates , x y, defined as c c the center of mass ω c is defined as in (2), and p is the perimeter. In this form, Q is (an approximation of) the inertia tensor of a uniform mass distribution lying on the perimeter of the floe, and is therefore a proxy for symmetry properties related to its shape. The symmetric matrix Q has two real eigenvalues, q 1 and q 2 , representing the moments of inertia with respect to the two principal axes of rotation. Notice that the square of the gyration radius is / = + = l q q Q 4 T r which expresses the elongation of the object. Specifically, it is the inverse of the ratio between the area of the inertia ellipsoid of the floe (which approximates its surface) and that of the smallest circle that contains it. For a spherically symmetrical floe (and for all regular polygons) it takes its minimal value = r 1 2 , while = ∞ r 2 for a segment. Its square root r is the aspect ratio, i.e., the ratio between the lengths of the longer and the shorter axis.
We have measured r 2 for all contours in the four data sets. The empirical distribution functions (shown in Fig. 4) are peaked around = . ∼ . To interpret and characterize these data, we introduce a simple stochastic geometric model, whose only parameter is linked to the distribution of seed locations. Specifically, the seeds are placed randomly following what is known as a "simple point process" in probability; the actual implementation is described below. To fix the ideas, this model may be interpreted as describing the refreezing of floes during winter as the coalescence of radially growing ice floes started from randomly-located seeds, but its formulation is purely geometric and one cannot exclude that it may be derived from alternative physical interpretations. Once the seeds are fixed, the shapes of the resulting floes are obtained as the cells of the Voronoi tessellation generated by the seeds 44 . A Voronoi tessellation is a partitioning of the plane into regions (cells), based on the distances from a set of special points (seeds). Each point in the plane (minus a set of zero measure) is associated to its closest seed. A cell is the set of all points associated with a given seed. In other words, floes are the domains of proximity of the freezing seeds.
The simplest point process is arguably the Poisson point process, whose realizations are sets of random independent points. We employ the following generalization, aimed at introducing repulsive correlation between points in an easily controllable way. Refer to Fig. 5. Fix an integer n (n 2 will be the number of seeds) and consider a rectangular region  where the noise terms ξ i x and ξ i y are independent identically-distributed Gaussian random variables with zero mean and unit variance. The quantity σ then expresses the departure from perfect order. The perfect lattice is realized for σ = 0, in which case the Voronoi cells are all regular hexagons with shape factor = r 1 2 ; the opposite limit σ → ∞ recovers the Poisson point process. Simulating the model for several values of σ allowed us to choose, for each data set, the value of σ that best fitted the empirical distribution of shape factors. Fig. 6 displays the theoretical distributions obtained by simulation, and shows that, as expected, σ → ∞ converges to the random Poisson-Voronoi case, while a delta-shaped distribution centered at = r 1 2 is approached for σ → 0. The agreement between this geometric model and the empirical data is impressive (Fig. 4 and inset of Fig. 6), suggesting that the phenomenology encoded in the elongation of ice-floe shapes can be completely characterized by a single parameter, measured by the average shape factor, or, equivalently, by σ. We therefore propose σ as a measure of the "characteristic disorder" in the distribution of the ice floe seeds.
The solidity of such a measure is supported by the following observations. (i) Results are not sensitive to the number of seeds (in the simulations reported in the figures, we sampled 10 3 realizations with ≈ n 10 2 5 points for each σ, but the results are indistinguishable from those at ≈ n 100 2 ). (ii) Choosing a specific lattice (for instance square or hexagonal) and probability distribution of ξ i does not affect the agreement with the empirical data, apart from slight readjustments of the fitted values of σ. (iii) Surprisingly, and importantly, the shape factor of a single floe is uncorrelated to its size, as shown by the scatterplots in Fig. 7 for two data sets with very different characteristic lengths (the other two are similar). This last observation indicates that the shape properties are decoupled from the size distribution, and justifies the description of the floe shaping mechanism as a Voronoi tessellation for a wide range of length scales.

Discussion and Conclusions
The sea ice system is complex, and several processes are at play, influencing each other. Further data analysis and modeling may attempt to assess, for instance, the relevance of melt pond formation on the thinning surface, which enhances fracturing, or of the development of pressure ridges. Our data are relatively limited, and leave the physical interpretations open, but the regularities observed in the floe shapes express important constraints that any realistic theory should satisfy. The main result of this work is the introduction of two novel independent scalar observables, or "order parameters" for ice floes, summarizing the properties of the distribution of sizes and of shape anisotropies respectively. Regarding floe size, data and different models of fracture suggest a scenario where the spectrum of the different forces at play is reducible, to a good approximation, to a single dominant scale, δl. This parameter contains all the information concerning the complex environment, including the ocean, the winds, and the This shows that independent measurements linking floe size to environmental history, conditions, and perturbations are needed to gain most physical insight. Different models assign different interpretations to the length scale δl, which is the only distinguishing feature of our four data sets as long as floe size is concerned. Note that the assumption of stationarity corresponds to considering long times; this is realized in practice if the dynamics of fracture is faster than the seasonal variation of parameters. The scale δl is probably affected by the annual changes in temperature and by other environmental conditions, and it seems reasonable to expect it to vary on time scales longer than the one driving equation (3). A systematic study of the variations of δl in a given location could give important clues on the changing climate and environment. Possibly, the characteristic length could provide an indicator of ice thickness 18 , which is an important parameter in the interactions between ice and ocean waves 22 . While several fracture theory and fragmentation scenarios may account for the fact that floe sizes have ubiquitous  ) and Poisson-Voronoi random tessellation. Curves "a" to "f " are obtained by the model with increasing values of the shape disorder σ; the curve "g" is obtained by the Poisson-Voronoi model. The inset shows how the average shape factor (corresponding to the curves in the main panel and to empirical data) is related to the inverse of the shape disorder σ / 1 (red line).
Scientific RepoRts | 5:10226 | DOi: 10.1038/srep10226 characteristics, we cannot exclude that other physical processes are relevant to establish this. For example, the wave field on the surface of the ocean, might be determinant for setting the single length scale controlling the floe size distribution. Additionally, the emergent characteristic length is likely affected by multiple important factors such as season (growth vs melting), location (marginal vs consolidated) and ice age and type (saline recently formed vs fresher older ice).
Considering shape elongation, we find that, as in the case of floe size, a single parameter fully characterizes the distribution in the whole available range of variability. However, differently from size, we found no universal scaling function accounting for the shape anisotropy data. On the other hand, we showed how a simple geometric model of ice accretion and coalescence (related to random Voronoi tessellations) realizes a one-parameter family of distributions fully capturing the shape fluctuations of ice floes, condensing once again the statistics of the elongations to a single dimensionless parameter, the shape disorder σ. The remarkable agreement between model and data validates the interpretation of σ as a physically relevant observable. The precise identification of the physical processes responsible for shaping the observed distributions is a difficult task, and remains an open question. Floe elongation distributions may be sensitive to the specific conditions mentioned above (season, location, ice age and type), and future systematic data-analysis studies may be able to assess these features.
It would be interesting to study the shape of ice fragments in other systems, both geological (e.g., the CO 2 ice layer on Mars 45 ) and in the laboratory 16 . The parameter σ is possibly related to the packing properties of ice just before the freezing season 3 . In line with the stationary model for floe size, the random Voronoi tessellation model suggests that shape anisotropies might be the product of fracturing mechanisms facilitating separation along the junctions between coalesced floes. Plausibly, more regular lattices in the model (lower values of σ) correspond to closer packing, and thus to stronger interaction between the floes. Curiously, random Voronoi tessellations are used as a phenomenological protocol for generating realistically-looking fracture patterns in computer graphics 46 .
These results could have implications for the rheology of sea ice (a description of ice in the marginal ice zone as a non-Newtonian fluid has been attempted by some authors 47 ). A potential practical application of such results could be the reconstruction of the number of floes of a given scale starting from incomplete measurements or measurements at different scales. For instance, consider the situation where one has access to satellite data at low resolution (large l min ), which include only large floes, beyond the crossover scale. Then the full curve can be reconstructed, by fits against the scaling form, and the number of smaller floes can be indirectly evaluated.

Data and methods
Data sets. We give here more details about the satellite imagery used. Each data set is composed of a single satellite image. Image 1 (Montagu Island) was taken by the QuickBird 2 satellite on October 16, 2003; catalog ID: 1010010002631B00 48 ; we used the top right quadrant of the image covering an area of ~100 km 2 , available in high resolution through Google Maps. Image 2 (Hopen Island) was taken by the GeoEye 1 satellite on June 13, 2009; catalog ID: 1050410001E37000 48 ; We used a small portion of around × 10 10 km 2 in the north-east sector of the image, available in high resolution through the Google Earth application. Spectral bands and exact resolutions for images 1 and 2 could not be determined; if one assumes, as is likely, that Google services use multispectral (red+ green+ blue) images then the pixel resolutions are .  Thematic Mapper Plus" instrument; bands 1,2, and 3, corresponding to wavelengths in the red, green, and blue visible light respectively, were merged; the pixel resolution is 30 m. Both images have an approximate spatial extent of × 170 180 km 2 . Only the quadrants with less than % 10 cloud coverage were selected.
Segmentation method. All images were first made monochromatic by application of a threshold on pixel intensity. This threshold is gauged by sight, but we tested that its precise value does not have a relevant impact on the results, as long as it lies between the average ice value and the average background value. Then the contours of single ice floes were extracted (Fig. 1) using the Potrace algorithm 51 . The main steps of the algorithm are as follows. (I) The bitmap is decomposed into paths, separating black and white regions. Square lattices (as the pixels in an image) present an ambiguity in the definition of clusters, when two contours meet perpendicularly at a vertex; we used the prescription of connecting preferentially black components in these cases. (II) After despeckling (discarding paths enclosing only 1 or 2 pixels), paths on the square lattice are converted into polygons, following a parameter-free optimization phase (details are in the documentation 51 ). The full algorithm includes a further step-based on aesthetic principles-aimed at detecting sharp corners and smoothing out the others, which we skipped.
Remarks on the cutoffs. Figure 1 shows a non-negligible area of unsegmented ice, mainly due to two components. (I) Very small floes, below the despeckling threshold; this is due to the pixel resolution, as discussed above, and is exposed by the quick drop-offs of the curves in Fig. 3. (II) Darker regions, below the intensity threshold; this is a more serious limitation, as it reduces systematically the estimate of the distributions at low sizes. It is probably responsible for the slight deviations from power law that are visible in some of the curves in Fig. 3.
While pixel size imposes a cutoff on the minimum floe-size detectable, image size has an influence on the statistics of large objects. Underestimation of the number of large floes can be quantified 52 ; for instance, a circular floe of 5 km in diameter at least partially covered by a Landsat × 170 180 km 2 snapshot has a % 10 probability of being segmented incompletely. Notice that such a small change would be undetectable in the logarithmic plots in Fig. 2 and Fig. 3. We manually removed all partially covered floes from the data sets.