Articles | Volume 13, issue 11
The Cryosphere, 13, 2869–2885, 2019
The Cryosphere, 13, 2869–2885, 2019

Research article 08 Nov 2019

Research article | 08 Nov 2019

Estimating the sea ice floe size distribution using satellite altimetry: theory, climatology, and model comparison

Estimating the sea ice floe size distribution using satellite altimetry: theory, climatology, and model comparison
Christopher Horvat1, Lettie A. Roach2,3, Rachel Tilling4,5, Cecilia M. Bitz2, Baylor Fox-Kemper1, Colin Guider6, Kaitlin Hill7, Andy Ridout8, and Andrew Shepherd9 Christopher Horvat et al.
  • 1Institute at Brown for Environment and Society, Brown University, Providence, RI, USA
  • 2Department of Atmospheric Sciences, University of Washington, Seattle, WA, USA
  • 3National Institute of Water and Atmospheric Research, Wellington, New Zealand
  • 4Cryospheric Sciences Laboratory, NASA Goddard Space Flight Center, Greenbelt, MD, USA
  • 5Earth System Science Interdisciplinary Center, University of Maryland, College Park, MD, USA
  • 6Department of Mathematics, University of North Carolina, Chapel Hill, NC, USA
  • 7School of Mathematics, University of Minnesota, Minneapolis, MN, USA
  • 8Centre for Polar Observation and Modelling, University College London, London, UK
  • 9Centre for Polar Observation and Modelling, University of Leeds, Leeds, UK

Correspondence: Christopher Horvat (


In sea-ice-covered areas, the sea ice floe size distribution (FSD) plays an important role in many processes affecting the coupled sea–ice–ocean–atmosphere system. Observations of the FSD are sparse – traditionally taken via a painstaking analysis of ice surface photography – and the seasonal and inter-annual evolution of floe size regionally and globally is largely unknown. Frequently, measured FSDs are assessed using a single number, the scaling exponent of the closest power-law fit to the observed floe size data, although in the absence of adequate datasets there have been limited tests of this “power-law hypothesis”. Here we derive and explain a mathematical technique for deriving statistics of the sea ice FSD from polar-orbiting altimeters, satellites with sub-daily return times to polar regions with high along-track resolutions. Applied to the CryoSat-2 radar altimetric record, covering the period from 2010 to 2018, and incorporating 11 million individual floe samples, we produce the first pan-Arctic climatology and seasonal cycle of sea ice floe size statistics. We then perform the first pan-Arctic test of the power-law hypothesis, finding limited support in the range of floe sizes typically analyzed in photographic observational studies. We compare the seasonal variability in observed floe size to fully coupled climate model simulations including a prognostic floe size and thickness distribution and coupled wave model, finding good agreement in regions where modeled ocean surface waves cause sea ice fracture.

1 Introduction

Earth's polar oceans are covered with sea ice: a thin, heterogeneous interface that plays an important role in the coupling between ocean and atmosphere. Sea ice is a collection of many individual pieces, called floes, which may be characterized in terms of a horizontal length scale, their “size”. On the large scales relevant to global climate modeling, the statistical variability of floe size is described using the floe size distribution (FSD; Rothrock and Thorndike1984).

The FSD is an important property of the sea ice cover that influences the multiscale temporal and geographic variability of sea ice, akin to the grain size in sedimentology or particle size distribution in atmospheric chemistry. The scale of individual floes plays a role in many sea-ice-related processes: sea ice melt rate (Steele1992; Horvat et al.2016; Horvat and Tziperman2018), the evolution of the oceanic mixed layer (Manucharyan and Thompson2017), atmospheric boundary layer exchange (Birnbaum and Lüpkes2002; Lüpkes and Birnbaum2005; Tsamados et al.2014), the sea ice response to applied stress (Feltham2008; Wilchinsky and Feltham2011), and the propagation of waves into the ice (Squire et al.1995; Squire2007; Smith and Thomson2016). The importance of the sea ice FSD has led to the development of diagnostic FSD models of varying complexity (Williams et al.2013; Zhang et al.2015; Bateson et al.2019) and a prognostic floe size and thickness distribution (FSTD) scheme (Horvat and Tziperman2015; Roach et al.2018a).

Despite the potential relevance of sea ice floe size to polar climate evolution, there remain no climate-scale assessments of average floe size or the FSD. The observational record of floe statistics derives from visual imagery localized in space and time (i.e., Rothrock and Thorndike1984; Toyota et al.2006, 2011; Steer et al.2008) or from repeat measurements in the same region over multiple months (Hwang et al.2017; Stern et al.2018a), which may subsequently be used to compile a seasonal cycle of the FSD (Perovich and Jones2014; Stern et al.2018a). FSD measurements are obtained by identifying individual floes within a two-dimensional image of the sea ice surface. Because floe sizes span several orders of magnitude, accurate representations of the FSD – even in relatively small geographical domains and in perfect lighting and surface conditions – require high resolution and high observational coverage. Nearly all measurements of the FSD have been made in accordance with a “power-law” scaling hypothesis commonly used to describe multiscale systems (Mandelbrot and Wheeler1983), in which the resulting FSD is fit to a straight line in logarithmic coordinates, whose slope, α, is reported as an intrinsic property of the floe mosaic. There is large uncertainty in these scaling coefficients, the range they apply over, and their applicability and origin (Herman2011; Horvat and Tziperman2017; Herman et al.2018; Stern et al.2018b). Improvements in the quality and quantity of available FSD data are needed before arriving at consensus-derived FSD statistics to guide and assess model performance.

Here we outline a method that exploits satellite radar altimetry to construct the FSD and its moments across polar regions with sub-kilometer spatial resolution, sub-daily temporal resolution, and spanning multiple orders of magnitude in size. Altimeters, like the ones carried on the Envisat, ICESat, CryoSat-2, and ICESat-2 satellites, make repeated, frequent passes over polar oceans, and substantial efforts have been made to process the satellite returns to discriminate between open water, floes, and leads. The altimetric returns have found many uses, including reconstructing the sea ice thickness field (Laxon et al.2013; Tilling et al.2016, 2018) and ocean surface circulation under sea ice (Peacock and Laxon2004; Armitage et al.2018). Fields inferred from altimetry have led to advances in understanding polar systems: from forecast and climate prediction (Day et al.2014) to model validation (Schröder et al.2018; Allard et al.2018) to climate change studies (Laxon et al.2003; Kwok2018), and have been evaluated and validated using field campaign data (Skourup et al.2017; Sandberg Sorensen et al.2018; Tilling et al.2018).

One-dimensional measurements of sea ice properties, like along-track altimetric measurements of ice open water, have long been sought to describe the two-dimensional ice surface. Rothrock and Thorndike (1984) originally described a method for reconstructing the sea ice floe size distribution in a region using straight-line measurements over the geometry of floes. Lindsay and Rothrock (1995) later compiled the statistics of lead and ice spacings in two-dimensional imagery. Other work has taken place to derive and understand the width distribution of individual leads in visual imagery and altimetry (Wadhams et al.1988; Key and Peckham1991; Key1993; Wernecke and Kaleschke2015), which can be used to estimate heat fluxes and turbulent transfer between the ocean and atmosphere. To date, however, these studies have not been designed to facilitate a comparison with model data, nor have altimetric studies been used to compile floe size statistics. These objectives are the focus of this work.

We outline the mathematical theory that allows for comparison of altimetric datasets and the FSD in Sect. 2. In Sect. 3 we apply this method to a new dataset of segmented CryoSat-2 sea-ice-type data from 2010 to 2018. Using these data we produce the first climatological maps of mean sea ice floe size and fragmentation for the Arctic Ocean. We then test the power-law hypothesis, finding limited support for power-law scaling across most of the dataset in Sect. 4. One of the key aims of the paper is to develop floe size distribution measurements that are useful for model validation and calibration. In Sect. 5, we show a proof of concept, demonstrating how altimetric data can be used to constrain and evaluate new models of the FSD, comparing the CryoSat-2 FSD data to a climate model simulation with a prognostic FSTD model. We conclude in Sect. 6.

2 Floe chords and the floe size distribution

For an individual pass over sea ice by a polar-orbiting satellite altimeter, return waveforms along the satellite orbit track are assigned a surface type depending on the waveform shape and coincident sea ice concentration (Tilling et al.2018). A “floe chord” of length D is a continuous series of points identified as sea ice, covering a geographic distance D (Tilling et al.2019a, b). Define a floe's size, r, as its “effective radius” – the square root of the floe's area divided by π (Rothrock and Thorndike1984; Horvat and Tziperman2015) We use radius instead of diameter, as appears in some other observational studies, for comparison with model output in Sect. 5. Because the satellite path is at an unknown angle with respect to the (also unknown) floe geometry, any individual floe chord measurement is not a floe size measurement. Converting between suitably processed altimetric floe chord measurements and floe size statistics is therefore the subject of this section. Details on the processing of the CryoSat-2 waveform, used to produce a dataset of floe chords spanning the period 2010–2018, are outlined in Sect. 3 and Tilling et al. (2019b).

For a domain of horizontal area A, and over a period of time ΔT that corresponds to several repeat satellite passes, we bin the set of recorded floe chords to form a probability distribution S(D), which we term the “floe chord distribution” (FCD), where S(D)dD is equal to the number fraction of floe chords in A over ΔT with length between D and D+dD, and is normalized to 1. To collapse all measured chords onto a single independent scalar coordinate (D), we follow the example of turbulence statistics (Batchelor1953) and assume that the floe chord distribution data are homogeneous, isotropic, and stationary within the region and time data are collected. In the same region, we define the (non-cumulative) number FSD P(r), where P(r)dr is the fractional number of floes with a size between r and r+dr in A, and is also normalized to 1. The FSD inherits the assumptions of homogeneity, isotropy, and stationarity from the FCD. Our objective is to relate the FCD, S(D), or quantities derived from the FCD, to the statistics of the FSD, P(r).

Bayes' theorem relates S(D) and P(r) through conditional probabilities,

(1) F ( r ; D ) S ( D ) = F ̃ ( D ; r ) P ( r ) .

The conditional probability F(r;D) relates given chord lengths to the floe size distribution that could generate them: F(r;D)dr is the probability that floes with size in the range from r to r+dR were sampled given a chord of length D. The conditional probability F̃(D;r) relates given floe sizes to the chord length distribution they generate: F̃(D;r)dD is the probability of measuring a floe chord of length from D to D+dD given that a floe of size r was measured.

This second probability distribution F̃(D;r) can be derived from first principles under a single assumption: that the chord length distribution that would be sampled from a set of floes of size r is independent of r (equivalently, the floe shape distribution is scale-invariant). Formally, this requirement is

(2) F ̃ ( D ; r ) d D = G ( ξ ) d ξ ,

where G(ξ)=G(D2r) is an unknown function that integrates to 1 over the interval from ξ=0 to 1. Under this assumption, the distribution of possible chord lengths measured from floes of size r has the same functional form independent of r. The probability distribution F(D;r) may be derived by considering the geometric relationship between straight-line satellite passes and the geometry of the floes they pass over. Individual floe shapes are highly variable: making an assumption about the distribution of floe shapes may introduce biases in the statistics derived from the FCD. Yet as we prove in Appendix A, the ability to derive FSD statistics from the FCD does not depend on the precise form of F̃(D;r) so long as the homogeneous, isotropic, stationary, and scale-invariance assumptions are retained, and the evaluation of power-law scaling is in fact independent of F̃(D;r).

To proceed and arrive at a concrete (although not general) realization of these functions, we will assume all floes are perfect circles. In assessments of the relationship between major and minor axes of individual floes, the “roundness” parameter for a floe is typically within 15 % of 1 (Rothrock and Thorndike1984; Toyota et al.2011; Perovich and Jones2014; Gherardi and Lagomarsino2015; Alberello et al.2019), suggesting that this circular assumption, while simplistic, is broadly appropriate. Nevertheless, it will likely be necessary to amend the analysis below in the future to account for more realistic shape distributions and geometries (e.g., diamonds; Wilchinsky and Feltham2006), regional differences in floe shape properties (such as in regions where shear stress determines fracture patterns and floe shapes; Schulson and Hibler1991), or to evaluate the sensitivity of the results that follow to the assumed shape distribution. Solving for F̃(D;r) is a geometric problem that relates the possible measured chord lengths to the underlying floe size, and we solve this explicitly for circular floes here. Similar geometric problems have been identified and solved in other fields (e.g., Pons et al.2006; Nere et al.2007), and we therefore leave refinement of F̃(D;r) to future work.

Figure 1Relating a floe chord to floe size for a circular floe. A satellite track (dashed black line) passes over a floe of radius r (solid black line). The track records a series of echoes of length D, which is the length of a chord (red line) identified by its interior angle, θ.


Consider the special case that all floes are perfect circles, illustrated in Fig. 1. Because there is no correlation between the statistics of local sea ice deformation and predetermined satellite tracks, an individual recorded floe chord, D, originating from a floe of radius r, was obtained from a satellite trajectory that crosses the floe at a random interior angle θ; thus the distribution of θ is uniform. Because of rotational symmetry, we need only consider θ[0,π), sampled according to a probability distribution T(θ;r)=π-1. The length D is thus a chord of this circular floe, with D=2rsinθ/2. Accordingly,

(3) F ̃ ( D ; r ) = T ( θ ; r ) θ D = 2 π 1 ( 2 r ) 2 - D 2 r > D / 2 , 0 otherwise ,

which is a probability function that meets the above criterion (Eq. 2).

The nth moment of the floe chord distribution S(D) is defined

(4) D n 0 D n S ( D ) d D = 0 d r P ( r ) 0 D n F ̃ ( D ; r ) d D .

For any function F̃(D;r) satisfying the scale invariance above, the right-hand side may be expressed in terms of moments of P(r) (see Appendix A). For circular floes, using Eq. (3),

(5) D n = 0 d r P ( r ) 0 2 r 2 π D n ( 2 r ) 2 - D 2 d D = 0 d r P ( r ) 2 n + 1 π r n 0 π 2 sin ( x ) n d x = A n r n ,

where D2rξ=sin(x), rn is the nth moment of P(r), and the coefficient An is


where B is the beta function. For n=0, 1, 2, or 3, then An is 1, 4π, 2, or 323π. Two important FSD-derived quantities are derived from ratios of FSD moments, and therefore can be obtained from the FCD directly: the “representative radius” (Horvat and Tziperman2017; Roach et al.2018a),

(6) r 0 r 3 P ( r ) d r 0 r 2 P ( r ) d r = r 3 r 2 = 3 π 16 D 3 D 2 ,

and the floe perimeter per ice area, a measure of sea ice fragmentation,

(7) P 0 r P ( r ) d r 0 r 2 P ( r ) d r = π 2 D 1 D 2 .

These derived quantities are useful because they require no further information about the sea ice (such as its concentration) to compare against modeled FSDs. However, both r and 𝒫 can represent only those floes whose size is larger than rmin=Dmin/2, the smallest possible floe size sampled. For perfect power-law distributions beginning at a scale of rmin or before, both metrics are functions of rmin. However, for the real FCDs measured here, a maximum floe size exists, and a power-law scaling is not found approaching rmin, so the use of such metrics is justified (see Sect. 4). Because of the finite sampling resolution of the altimeter, chords that would originate from floes with a diameter near the sampling resolution may not be observed, and thus Dn〉≤Anrn. We explore this uncertainty in Appendix B. For a known floe size distribution, the error decreases exponentially as a function of the distributional moment being considered, though it can be large (20 % or more) in pathological cases. For distributional tails characterized by observed scaling exponents (Stern et al.2018b), and for moments considered here, this uncertainty can be determined systematically and vanishes for measurement spacings smaller than the radius of the most common floe size. This resolution error does not affect the analysis of the power-law hypothesis, as that analysis is focused on the distributional tail. However, because 𝒫 is proportional to a negative moment of the FCD, it is sensitive to changes in the number of small chord lengths. Because of the measurement uncertainty for smaller chord lengths we will focus instead on r, which is a positive moment of the FCD.

2.1 Evaluating the floe size power-law hypothesis with floe chord data

Suppose the FSD P(r) has a power-law tail that begins at some specified value r1. Then for r>r1, P(r)P(r;α,C)=Cr-α, for an unknown coefficient C and power-law slope α. Integrating Eq. (1) over all r,

(8) S ( D ) = 0 F ̃ ( D ; r ) P ( r ) d r ,

where the integral of the left-hand side of Eq. (1) is equal to S(D) as F(r;D)dr=1. Under the assumption of Eq. (2), if P is a power law, so is S(D) (Appendix A). For circular floes,

(9) S ( D ) = 2 C π r 1 r - α ( 2 r ) 2 - D 2 d r .

Because of the sampling resolution of the altimeter there is a minimum resolved chord scale Dmin. If DminD*2r1, there is an explicit solution for S(D), a power-law distribution over the range (D*,),

(10) S ( D ) = C B 1 2 , α 2 2 α - 1 π D - α C α D - α ,

where B is the beta function. The coefficient Cα is a multiplicative factor independent of size, and the power-law exponent for a FCD is the same as the exponent for FSD, where the two are related by Eq. (1).

Moments of a power-law tail can be evaluated explicitly (for α>n+1),

(11) r n = C r 1 r n - α d r = C r 1 n + 1 - α n + 1 - α .

Then for both the FCD and FSD, the ratio of two moments is independent of the unknown coefficient C, i.e.,

(12) R n , ϵ D n + ϵ - 1 D n - 1 = D min ϵ n - α n + ϵ - α ,

valid for n+ϵ<α. The power-law coefficient can be obtained for any n, ϵ as

(13) α n , ϵ = n + ϵ R n , ϵ R n , ϵ - D min ϵ =  constant .

In the analysis below we will arbitrarily select only n=0.5,ϵ=1 for comparison (for scaling coefficients α>1.5, the bulk of reported power-law coefficients are in this range; Stern et al.2018b). Because the observations will not be perfect power-law distributions, we will use α0.5,1α* as an estimator. A second estimate of the power-law scaling coefficient, α^, is computed via the maximum likelihood estimator (Muniruzzaman1957; Clauset et al.2009; Virkar and Clauset2014) (details in Appendix C) as

(14) α ^ = 1 + N i = 1 N ln D i D min ,

where N is the number of chords. If the power-law hypothesis holds, then the two estimates of α agree, although the agreement of α^ and αn,ϵ is not sufficient to confirm the power-law hypothesis. In the Supplement Sect. S1, we compare these two estimates when they are evaluated against synthetic datasets drawn from a true power-law distribution. The two agree even when the size of the data is relatively small (N<25). While in practice Eq. (13) is easy to apply, it only holds when αn,ϵ>n+1, and unlike the method of Clauset et al. (2009), it does not allow for a robust statistical analysis of the power-law fit, and should only be used when the data are assumed to follow a power-law already.

Figure 2Constructing a FCD from altimetry. (a) Base 10 logarithm of the number of floe chords identified and binned into the CESM grid, across all CryoSat returns in the Arctic from 2010 to 2018. Black line is a single satellite track on 21 January 2014. (b) Subsection of the track centered on the blue dot in (a). Blue line is freeboard of sea ice in radar echoes defined as “floes” along the track. Black lines are chords identified from the freeboard retrieval. (c) Total number of chords measured in each month in the Arctic. Plot is centered on 1 January. (d) FCD for the satellite track depicted in (a). Black marks on the x axis are the logarithmically spaced chord length bins.

3 Climatology and trends in floe properties derived from CryoSat-2 altimetry

We apply the analytic technique described in Sect. 2 to a floe chord dataset constructed from the CryoSat-2 radar altimeter processed by the Center for Polar Observation and Modelling (CPOM) over the period from October 2010 to present (CPOM data products are available at, last access: 9 January 2019). CryoSat-2 radar echo returns are defined as “lead”, “floe”, “open ocean”, or “ambiguous” according to waveform shape and sea ice concentration (Tilling et al.2016, 2018), at an approximately constant along-track spacing Dmin=300 m. Floe chords are defined as a continuous sequence of one or more “floe echoes”, with a gap of one ambiguous echo permitted within a floe sequence to allow for anomalous returns. A chord length is taken from the midpoint of the first to the midpoint of the last radar echo. Individual chord lengths can be underestimated when continuous floes are separated artificially by producing two or more ambiguous echoes in sequence, or when highly reflective leads dominate the waveform return close to the floe edge and cause measurement dropout (Tilling et al.2019b). Lead contamination, or “snagging” (Armitage and Davidson2014), is more likely when the altimeter cuts off a small section of a floe, i.e., for small values of θ. Overestimates of chord length can also occur when ice floes are in close contact with neighboring floes. Therefore, floe chord lengths should be considered a satellite-derived product, not a true measurement of floe size. The minimum chord length retrieval Dmin is limited to the CryoSat-2 footprint (∼300 m along-track) (see the discussion in Appendix B). However, surface discrimination via altimetry is highly accurate in months without melt ponds (Peacock and Laxon2004; Guerreiro et al.2017; Quartly et al.2019), giving confidence that floe echoes represent a coherent length of ice. More details on the details of chord identification may be found in Tilling et al. (2019b). Indeed, these raw floe chord data have been used successfully to reduce biases in altimeter-observed satellite sea ice thickness estimates from altimeters with different footprint sizes (Tilling et al.2019b). Here we analyze the sea ice floe size distribution using that floe chord product.

Figure 2 shows an example of floe chord data for a single CryoSat-2 track over the Arctic on 21 January 2014. Freeboard values for echoes discriminated as floe are plotted in Fig. 2b as a function of the along-track distance in kilometers, and correspond to the blue circle in Fig. 2a. Floe chords are identified as black segments in Fig. 2b. The histogram of all 741 identified chords for this single satellite pass is shown in log-log space in Fig. 2d.

The full CryoSat-2 dataset examined here spans the time period from October 2010 to November 2018, and floe chords measured using the above technique are binned into the CICE sea ice model's two-dimensional sea ice grid for each month and year to facilitate comparison with model products. This implies that we invoke the principles of isotropy, homogeneity, and stationarity of the FCD, required to produce such a distribution, on the length scale of the CICE model grid (O(25 km)) and timescale of a month. For every grid cell i, month m, and year y, we have a vector of floe chords {Di,m,y} from which we build a FCD. The base 10 logarithm of the total number of floe chords recorded in each grid cell per month is shown in Fig. 2a. Because the satellite passes are densest near the pole, the measurement density is highest near the pole as well. Figure 2c shows the number of Arctic measurements in each month. Sea ice type from CryoSat-2 is not available during summer months, as melt ponds make it difficult to discriminate between leads and ponded floe surfaces, and we do not include measurements from May to September. Across the entire set of satellite tracks included here, 11 million chord lengths are recorded in the Arctic.

Figure 3a shows the seasonal cycle of Arctic representative radius over the CryoSat-2 period obtained by applying Eq. (6) to the binned CryoSat-2 floe chord product. Individual years are plotted as thin lines, and the climatological average is shown in red. Details on how temporal and spatial average statistics are computed are included in Appendix D. During the months of October–December, the climatological representative radius is roughly 35 % larger (7.06 km vs. 5.18 km) than during February–April. This seasonal cycle is broadly consistent across years. A possible interpretation of this seasonal cycle is that large first-year ice pans form in October and are later fractured into smaller floes throughout the winter months. This concept is supported by observations that large-scale fracturing of sea ice in the Beaufort Sea is dominated by coastal processes and therefore can only occur once sea ice freezes to the coast in midwinter (Richter-Menge2002), although such an interpretation is speculative and must be evaluated further as this method is refined. Figure 3b shows annual-average representative radius in red for each full year from 2011 to 2017, with thin lines corresponding to the individual months within that year. Seasonal variability is significantly larger than inter-annual variability. There is no statistically significant linear trend at the p=0.05 level.

The geographic variability of representative radius over the “early winter” (October–December) and “late winter” (February–April) periods is shown in Fig. 3c–d, for all grid areas. We display only those areas with at least 25 recorded floe lengths in each month during the averaging period. In Sect. S2 and Fig. S1, we examine the sensitivity of bulk FSD statistics to this threshold, finding similar seasonal cycles and climatologies. The largest representative radii in the Arctic lie in the interior Arctic near the pole, with a tongue of large floes that extends along the Canadian Arctic in late winter. There is a notable increase in representative radius with latitude. In Fig. S2, we show that this relationship cannot be explained as a result of the increasing density of measurements near the pole and may therefore be a geophysical signal. The smallest representative radii (below 1 km) lie in the Bering Strait and the Russian Arctic in early winter and in the Laptev Sea in late winter. The difference in representative radius between fall and spring is accounted for by the reduction of floe sizes in regions near the Arctic interior (see Fig. 6).

Figure 3(a, b) Temporal and geographic variability of Arctic representative radius. (a) Climatology of Arctic-average representative radius in units of kilometers (red line). Thin lines are individual CryoSat-2 years. (b) Annual-average Arctic representative radius (red line). Thin lines are the average in individual months. (c) Climatological representative radius in the months October–December. (d) Same as (c) but for February–April.

4 Evaluating the power-law hypothesis using floe size statistics derived from CryoSat-2

Given a collection of chord lengths, we would like to examine whether it is distributed according to a power law. Under the assumptions of Sect. 2, the scaling behavior of the FSD is the same as of the FCD (see Appendix A). We use the statistical methodology outlined in Clauset et al. (2007, 2009) and Virkar and Clauset (2014) (which we term the MLE method) to evaluate shape parameters of the most likely power-law fit and to test its plausibility. This method has been used to evaluate power-law behavior in a recent FSD model (Horvat and Tziperman2017) and observational studies (Hwang et al.2017; Stern et al.2018b) and proceeds as follows.

  1. Lower-truncate the FCD. First identify a minimum chord scale, D*, above which we hypothesize a power-law tail, and analyze only those floe chord measurements. We either (a) choose D* to be 900 m (to reduce the impact of small-size sampling errors discussed in Sect. 2) or (b) use the scheme described in Clauset et al. (2007) to evaluate the most likely value of D* for a power-law tail. The length of this lower-truncated distribution is N. In the descriptions that follow, we use the subscript “all” to describe case (a) and “tail” to describe case (b).

  2. Compute power-law scaling estimates and parameter uncertainty. We obtain two estimates of the FCD scaling estimate, either computing α* via Eq. (13) or computing α^, and uncertainty estimates in both α^ and D* via the MLE method (Eq. 14). That the two estimates of α agree is a necessary condition for the FCD (and thus FSD) to be power-law distributed.

  3. Examine the plausibility of the power-law fit. We generate M FCDs of size N (the same number of synthetic chords as observed chords), with each synthetic FCD drawn from the hypothesized power-law distribution P(α^,D*). For each of these synthetic FCDs, we compute the Kolmogorov–Smirnov distance between it and the hypothesized power-law model that generated it, P(α^,D*). We also compute the distance between the observed FCD and P(α^,D*). A p value, p, is equal to the fraction of those M synthetic FCDs that are “further away” from the hypothesized power-law model than is the observed FCD. We use M=10 000, which permits computation of p within 0.005 (Clauset et al.2009), and rule out the power-law hypothesis under the condition p<0.1 (Virkar and Clauset2014).

We note that a power law describes the scaling of a distribution's tail. Previous observational studies have discussed “double power laws” (i.e., Toyota et al.2011), i.e., two power-law distributions of a different exponent joined at a specified scale. The methods employed here would capably capture the large-size power-law scaling but not the small-scale scaling. Such double power laws are necessarily scale-variant, and require at least three parameters to describe. The conceptual and mathematical simplicity of the power-law hypothesis does not apply in such a case, and we do not consider them here.

Figure 4Examining the power-law hypothesis. (a) Histogram of all chord lengths recorded in the Arctic for the months November–April (black). Bin centers are indicated by hashes and are logarithmically spaced. The blue line is power-law fit to all observed sizes according to Eq. (14). The red line is power-law fit to the tail. Dashed red lines are fit lines using the ±1 standard deviation values of α^. The red vertical line is the most likely beginning of the power-law tail, D*, with shaded region ±1 standard deviation in D*. (b) Same as (a), but for measurements in April. (c) Maximum likelihood estimate of the beginning of the power-law tail, D* (km) for all measurements at each geographic location over the observational period. Only locations with N>1000 are plotted. (d) Maximum likelihood estimate of power-law tail exponent, α^tail, for the same points. Colored values have more than 200 chord lengths in the tail and p>0.1. Zero values are those locations plotted in (c) but where either p<0.1 or there are fewer than 200 measurements in the tail. (e) Number of chord lengths in the tail (above D*) at each location.

The MLE method is a rigorous test of the power-law hypothesis that eliminates potential human bias when interpreting observational data. To illustrate why this is important, we first consider the entire set of 11 million chord lengths recorded in the Arctic in all months (October–April), spanning a length range from 300 m to 100 km. The histogram of these floe chords is the black line in Fig. 4a (hashes on black line are the logarithmically spaced bin centers). Beginning from D*=Dmin=900 m, α^all=1.97 (blue line) and αall*=2.05 (not shown). The observations are further away from synthetic data drawn from P(α^all,D*) in each of the M=1000 random draws (pall=0/1000) and we reject the power-law hypothesis for these measurements. We note that if the resolution bias explored in Appendix B proves to be larger than expected, the underrepresentation of small floe lengths may affect the analysis of the full distribution.

Examining the tail of the distribution in Fig. 4a, the maximum likelihood estimate of D* is ≈15.0 km (red vertical line, vertical shaded region is the range of uncertainty for D*), above which there are ∼40 000 chord length measurements between 24.7 and 99 km (0.4 % of the dataset). On the truncated FCD, α^tail=4.65 (red line, dashed lines are uncertainty ranges for α^tail) and αtail*=4.67 (not shown), similar to the large-scale roll-off reported in observations (Toyota et al.2016). Even when restricted to the FCD tail, ptail=0/1000.

Finding no statistical basis for a power-law fit to the tail in Fig. 4a underscores the challenge in using the human eye to observe power-law scaling. While the black and red lines in Fig. 4a appear similar across much of the range of sizes above 24.7 km, examining the misfit between the power-law estimates and the data shows that the two curves in fact differ significantly across the entire fit range. A misfit error can be defined as

(15) E = P ( x i , α ^ tail , D * ) - P ( x i ) P ( x i , α ^ tail , D * ) ,

where xi represents the bin locations, angle brackets denote an average over the relevant bins, and P(xi) represents the observed histogram values. Over the range from 24.7 to 100 km, the misfit error is 33 %. The visual agreement, misfit error, and apparent slope and shape of the distribution depend sensitively on the bin spacing and the logarithmic plotting.

Sea ice parameterizations that assume a power-law distribution may significantly bias sea ice statistics. The imposition of any fixed distributional shape, when FSD dynamics are scale-variant, leads to implicit nonlocal redistribution of sea ice between floe size categories (Horvat and Tziperman2017). To see this in practice we compare the difference in Arctic-wide representative radius, r, which is used in parameterizations of wave attenuation and ice thermodynamics, between the most-likely power-law fit to the data and the “true” value obtained via Eq. (6). The observations yield r=10.2 km, versus 34.5 km for the power-law fit. Examining only the tail of the distribution (chord lengths above 24.7 km) yields better agreement: 23.7 km for the observations and 24.4 km for the fit line. Yet this tail constitutes just 1 % of all measured chord lengths, corresponding to just 18 % of total ice area and 4.5 % of the perimeter per square meter (Eq. 7).

Segmenting the chord length data into individual months in the Arctic, there are none where pall>0. Examining only the tail of each month's distribution, ptail<0.1 in all months. Only in April is there a nonzero ptail=0.04, for which the analysis of Fig. 4a is repeated as Fig. 4b. In April, α^all=1.99, α^tail=5.70, and D*=30.7 km. The tail consists of 1618 measured chord lengths up to 97.5 km, accounting for 8 % of the total floe area and 1.4 % of the perimeter per square meter. The misfit error between the April FCD tail and P(α^tail,D*) is 76 %. Accumulating all measured chord lengths from October to May into the CESM model grid, we find zero locations that support a power-law distribution across the range of measurements (i.e., pall>0.1). For grid areas with N>1000, we show the value of D* computed using the local FCD in Fig. 4c. Values of D* range from 2 km along the Russian Arctic to more than 10 km near the North Pole.

Figure 5Top row: Temporal variability of power-law fits to Arctic FCDs. (a) Estimate of the most likely power-law scaling coefficient for all recorded floe chords as a function of month over all years, calculated from the MLE method Eq. (14) (red lines) or Eq. (13) (blue lines). Thick lines are climatological averages, and thin lines are individual years. Plot is centered on 1 January. (b) Like (a), but plotted for individual years over all months. Thick lines are the average over months plotted in (b) and thin lines are individual months in each year. (c–d) Same as (a–b), but for the distributional tail starting from D* computed using the MLE method. “Arctic” refers to points above 60 N.


While most of the Arctic has at least 1000 total measurements across all years, FCD tails (D>D*) are not as well-sampled. We investigate these tails including regions with at least 200 measured floe chords larger than D*. The percentage of geographic areas with at least 1000 total measurements that have a tail with at least 200 measurements is 44 %; on average D* is 5.4 km for these regions. For most of these regions we can not rule out a power-law tail. For the subset of regions with 1000 total measurements, 200 measurements in the tail, and where the power-law hypothesis cannot be ruled out, the average D* is 6.5 km and average α^tail is 3.34, within the typical range of Arctic FSD measurements (Stern et al.2018b). In Fig. 4d we show the values of α^tail at these locations. Colored cells are those with p>0.1 and a tail with at least 200 measurements. In Fig. 4e we show the base 10 logarithm of the MLE tail for all geographic locations. Those regions for which a power law cannot be ruled out are generally those with the largest floes and the highest sampling, clustered near the central Arctic. The weakest support for a power-law tail is in the Chukchi and Beaufort seas, where power-law floe size distributions have often been reported. We note that our choice of tail length plays an important role in whether the power-law hypothesis is rejected in the tail across the Arctic. For example, the fraction of Arctic regions with at least 1000 total measurements, a tail of at least 100, 200, and 400 measurements, and that does not reject the power-law hypothesis is 72 %, 52 %, and 15 %, respectively. The better-sampled the FCD–FSD, the more likely the power-law hypothesis is rejected.

Scaling coefficients can provide useful information about the distributional shape. In Fig. 5a–d we show the seasonal and inter-annual variability of power-law estimates in the Arctic. Figure 5a plots the climatology of the power-law scaling estimates when including all measured chord lengths in dark red (using Eq. 13) or blue (using 14). Individual years are thin red or blue lines. The two estimates disagree. Because agreement between the two estimates is necessary for the power-law hypothesis to be true (see Sect. 4, Sect. S1), this alone is sufficient to rule it out. There is a seasonal cycle in the power-law fitting to the full distribution, with αall increasing (steepening) from September to January and remaining flat until April and no significant linear trend at the p=0.05 level for the annual-average value of αall. Figure 5c–d repeat this analysis on the tail of those monthly distributions. In this case, the two estimates agree well. There is a different seasonal cycle in the steepness of the distributional tail: shallowest in early winter and steeper in late winter. This indicates that the changes across the winter months may be due to a reduction of the largest floes and a steepening of the distributional tail, although there is significant inter-annual variability among these estimates. A similar seasonal cycle to that found in Fig. 6a, c, with an FSD that steepens from September to April, was found in image analysis of floes in the Beaufort and Chukchi seas (Stern et al.2018b), with α≈2.5, although the distribution steepened monotonically over that period. There is no significant linear trend at the p=0.05 level in the annual-averaged FSD tail slope (Fig. 5d).

5 An example model–observation comparison of floe size variability

A principal aim of this work is to allow model–data comparisons and facilitate testing rapidly developing FSD–FSTD models. Here we demonstrate how such a comparison can be made and provide useful information to modelers, even in the presence of the high uncertainties in this nascent FSD reconstruction technique. With the gridded data provided above, we may now directly compare development-stage sea ice models that incorporate FSD effects to observations. To do so, we use the Roach et al. (2018a) prognostic model for the FSD–FSTD, based on the Horvat and Tziperman (2015) theoretical FSTD framework, implemented into the CICE 5.1.2 (Hunke et al.2015) sea ice model. The FSTD is a sea ice state variable, subject to interaction of five key physical processes: lateral growth, lateral melt, fracture by ocean surface waves, welding of floes in freezing conditions, and wave-dependent new ice growth (Horvat and Tziperman2015, 2017; Roach et al.2018a, b). Previously published model runs (Roach et al.2018a) focused on the impact of the FSD on lateral melt, which is largely driven by small floes (Steele1992), and so floe sizes above 1 km were not considered. As a larger range of scales is resolved in the CryoSat-2 observational product, we conducted a model run that extended the floe size categories to scales larger than 1 km, using 24 logarithmically spaced floe size categories from 0.5 m to 33 km.

Figure 6Geographic and climatological comparison of modeled and observed representative radii. (a–b) Average representative radius from November to December in (a) the CryoSat-2 observational dataset and (b) the FSTD model. Grey shaded regions in (b) are the interior of contours in (a), which represent “pack ice” unaffected by waves in the model simulations. (c) Climatology of Arctic-average representative radius in units of kilometer for the MIZ in observations (red) and the model (blue). Green line is the annual average for the “pack”, the excluded regions in (b). Thin lines are averages in individual years from 2011 to 2016 in the MIZ. (d–e) Same as (a–b), but for the months of February–April. (f) Annual-average Arctic representative radius for wave-affected regions in MIZ observations (red), MIZ model (blue), and pack ice observations. Thin lines are the average in individual months in the MIZ observations.

This FSTD model simulation is coupled to a slab ocean model and the WAVEWATCH III ocean surface wave model (Tolman2009), forced by the JRA-55 atmospheric reanalysis (JRA-552013) over the period from 2000 to 2016. These wave-coupled runs are branched at year 2000 from a stand-alone sea ice run from 1975 to 2000, spun up using repeated 1975 atmospheric forcing. Additional model physics beyond those processes outlined in Roach et al. (2018a) have been added to determine the initial size of newly formed sea ice floes as a function of the ocean surface wave field. Details on this new parameterization, model initialization, and spin-up, are described in Roach et al. (2019). Recalling the finite measurement resolution of the CryoSat-2 dataset, the modeled representative radius is calculated only including floe size categories from 150 m (as floe sizes are radii, this corresponds to a radius equal to the minimum chord length) and larger. We include all FCD measurements here (chord lengths above 300 m) to make the broadest comparison, but note that the potential underrepresentation of floes with diameters near the sampling resolution may lead to inaccurate values of r in regions mainly consisting of such floes.

Figure 6a–b, d–e compare modeled and observed climatologies of Arctic representative radius (for floe diameters 300 m and larger) averaged over 2011–2016 and the months of October–December (a, b) and February–May (c, d). Geographic variability of representative radius is broadly similar between model and observation: the largest floes lie in the Arctic interior, with regions of smaller floes in the straits and continental margins. Across the interior Arctic, simulated representative radii are significantly larger than are found in the observations, as the Roach et al. (2018a) FSTD model does not include processes that break up large floes in the absence of ocean surface waves. To compare seasonality between model and observations, we compare only those regions that experience wave fracture in the model runs, areas we collectively term the marginal ice zone (MIZ). The MIZ is defined by excluding categories that do not experience wave fracture in a given month (see Appendix D), shown as the contoured regions in Fig. 6a–b, d–f and greyed out in Fig. 6b, e. All excluded “pack ice” regions have modeled representative radii greater than 18 km. The MIZ region accounts for 37 % of grid areas with at least 25 chord measurements in months from October to December and 35 % of such areas for the period February–March. Note that the month of October is absent from these plots because no well-sampled regions are classified as MIZ across all model years according to the criteria outlined in Appendix D.

Figure 6c compares the observed (red) and modeled (blue) Arctic-average representative radii for the MIZ over the period 2011–2016 as in Fig. 3a. The seasonal cycle of representative radius in the MIZ is different in the observations (red line; thin orange lines are individual months) than when all geographic regions are included (Fig. 3a). The seasonal cycle of representative radius in the pack ice region (i.e., not the MIZ) is shown as a green line in Fig. 6c. In the MIZ, average representative radii are smaller (on average 4.17 km vs. 6.49 km in the pack ice region). In contrast to the seasonal variation across all geographic regions (Fig. 3a) as well as in the pack ice, floes are larger in February–April than in November–December (5.40 km vs. 3.15 km). In both the MIZ and pack ice regions, however, average representative radius is similar in late winter. The largest difference between the two regions is from November to December, where representative radii are more than twice as large in the pack ice than the MIZ.

Figure 6f shows the annual average representative radius in the MIZ (red), pack ice (green), and modeled MIZ regions (blue). Modeled MIZ representative radii have a similar magnitude compared to the MIZ observations, though these regions have smaller floes than the interior. To address the scale mismatch between the too-high modeled floe sizes and observed representative radii in the interior Arctic, as well as the strong and different seasonal cycle in representative radius in both regions, modeling efforts must include additional mechanisms for reducing floe size in the Arctic interior away from waves, such as mechanical fragmentation (Toyota et al.2006; Rynders et al.2016) or ridge dynamics (Roberts et al.2019), to obtain realistic representative radii across the entire Arctic, as these processes are not present in the model used to make this comparison.

6 Conclusions

Here we developed and demonstrated a method for deriving the statistics of the sea ice FSD from satellite radar altimeter measurements of chord length. This method provides the first pan-Arctic accounting of climate-relevant quantities derived from the FSD, permits testing of existing scaling laws previously used to characterize distributions of floe size, and allows for gridded comparisons between FSD models and observations. Using this new technique we produced climatological, annual-average, and geographic mean moments of the Arctic FSD across a range of resolved length scales from 300 m to 100 km.

With the combination of satellite altimetry and mathematical theory, we were able to rigorously examine the power-law hypothesis related to the FSD under simple assumptions about the underlying floe chord data and the fidelity of CryoSat-2 satellite retrievals. Segmenting measurements by geographic location, by month, and by year, we find limited statistical basis for a power-law scaling beginning below about 6.5 km. In a limited number of geographic locations, we find the observational data cannot rule out power-law scaling, except for typical sizes above about 6.5 km. Assuming a power-law floe size distribution can bias sea ice model output and conceptual understanding, the geographic variability and lack of consistent multi-scale behavior reinforces the need for sea ice models to account for floe-scale processes rather than diagnose a distributional shape.

Observations that span the polar regions and different years and seasons are valuable for future refinement of process-based models of the FSD. In Sect. 5, we demonstrated how such model–observation comparisons can be made and can provide useful insights for model developers. At present, some general features of floe size evolution (in particular the magnitude and seasonal cycle of the representative radius) are broadly similar between model and observation in the marginal ice zone. Yet there is a significant scale mismatch in the interior Arctic between the presented simulations and this observational product because of missing fragmentation physics in the absence of ocean surface waves. Floe size modeling efforts have focused on marginal ice zone processes (Horvat and Tziperman2015; Zhang et al.2015), and particularly floe sizes below about 1 km because these small floes play an important role in sea ice thermodynamics for floe sizes. The CryoSat-2 observations, however, are best suited to resolving floe chords of several hundred meters and above. New satellite altimeters like ICESat-2 have the potential to increase the chord length resolution to scales of 20–100 m and provide insight at smaller scales.

We emphasize strongly that refinement may be necessary to apply this method for operational purposes, trend analysis, and further model validation. This paper has focused on the framework for making altimetric measurements of the FSD and comparison to model output, but the obtained chord lengths and distributions have not been carefully validated against other observational methods, and this will be necessary before further application of this method. Before doing so, we have tried to outline the most significant uncertainties in the method. The typical assumptions of homogeneity, isotropy, and stationarity are invoked here at the length scale of the CICE model grid (O(25) km on each side) and timescale of 1 month. These statistical assumptions may not be satisfied if, for example, the number of measurements in a given region in 1 month is insufficient to sample the known anisotropy of the sea ice floe field, and additional passes change the mean chord length significantly (see Sect. S2 and Fig. S1). While we found little evidence for power-law scaling throughout most areas of the Arctic, this may be sensitive to the geographic (here the CICE model grid of approximately 25 km × 25 km) and temporal (here all measurements from 2010 to 2018) windows we use to collect and evaluate chord length measurements for a power law. The assumption of scale-invariant sampling, observational uncertainty because of the finite sampling resolution, analysis of ambiguous returns, and the accuracy of retrievals in regions of thin sea ice may also affect the inferred size of sea ice floes. This in turn may affect the climatologies described in this study.

While processed CryoSat-2 data have been validated against both visual imagery and ground-based observations, they were not designed with this application in mind – additional quality control may be necessary for climate studies of changing floe properties. The positive comparison between model and observation in Sect. 5 could also be due to a compensation between these measurement uncertainties and will need to be re-examined in future validation work. Yet observational uncertainties regarding, for example, the floe shape distribution can be roughly estimated at the order of the error in effective radius obtained for circular floes (r=A/π) or a square (r=A/4), with a relative error of 25 %. Constraining model results beyond this scale of error will require further refinement. However, as shown in Fig. 6, at present the model–data mismatch in the interior Arctic can exceed a factor of 3. Even with expected levels of error in the present derived FCD–FSD product, some constraints on the model can be considered at present with this method. A future comparison of results from the Ice-Sat2 and CryoSat-2 altimeters will provide insights into the relevance of measurement and statistical uncertainties, as will comparison of altimetrically derived floe chord measurements with visual imagery.

Even accounting for important caveats that arise from making satellite measurements, remotely sensing the sea ice FSD from altimeters at sub-daily resolutions can provide a significant increase in data for comparison and analysis of new sea ice models that parameterize the FSD. Previously the difficulty of making measurements of the FSD at relevant spatial and temporal scales has inhibited the widespread adoption of such floe-sensitive sea ice models. Understanding sea ice variability at the floe scale is also an important aspect of sea ice forecasting, and the ability to remotely assess the sea ice FSD at near-real time will allow for further improvement of operational forecasting networks.

Data availability

CPOM sea ice data, including raw floe length data, are available through the CPOM data portal at (last access: 1 November 2019). The processed FCD–FSD statistics are available at (last access: 1 November 2019). The Roach et al. (2018a) FSTD model is publicly developed and available at (last access: 1 November 2019).

Appendix A: Proof that the FCD and FSD have the same statistical properties

For generic probability distributions S(D) and P(r), and a probability function F̃(D;r), via Eq. (4) we have the relationship

(A1) D n = 0 d r P ( r ) 0 r D n F ̃ ( D ; r ) d D ,

where we restrict the upper bounds on the second integral because F̃(D;r) is zero for D>r. Under the scale-invariant sampling assumption F̃(D;r)dD=G(ξ)dξ, where ξ=D2r for D<2r (ξ<1). Therefore,


where An is the nth moment of G(ξ), a constant that depends on the functional form of G. For any such probability function (for example that derived in Sect. 2 for circular floes), the moments of the FSD and the moments of the FCD are proportional. Most of the hypothetical statistical distributions we would consider (for example, power laws) can be fully determined in terms of their moments, and thus the relationship between moments of the FSD and FCD is typically sufficient to reconstruct the underlying FSD.

Supposing P(r) were a power-law function, converting Eq. () to an integral over ξ from 0 to 1, we have

(A5) S ( D ) = 0 F ̃ ( D ; r ) P ( r ) d r = 0 1 P ( D / 2 ξ ) G ( ξ ) ξ d ξ .

For a power-law function, P(D/(2ξ))Dξ-α and

(A6) S ( D ) D - α 0 1 ξ α - 1 G ( ξ ) d ξ = A α - 1 D - α .

From Eqs. (A4) and (A6), and under the assumptions of Sect. 2, all moments of the FSD and FCD are related by a computable function of the moment only, and power-law FSDs are derived from power-law FCDs with the same scaling law. While the proportionality of moments and Eq. (A6) prove that an observed power-law FCD must reflect an underlying power-law FSD, the same analysis used to arrive at Eq. (A6) can be repeated to find P(r) given a power-law-distributed S(D) as well.

Appendix B: Bounds on the relationship between chord length and floe size moments

The real altimetric data product has a finite sampling resolution Dmin, which can bias the computed FSD moments and power-law decay profile. For example, applied to real data with a finite sampling resolution, the integrals in Eqs. (4) to (5) are taken beginning at the minimum observed chord lengths Dmin and floe sizes rmin=Dmin/2. Moments of the distributions S and P reflect only statistics for floes larger than Dmin and rmin, respectively. All other aspects of this derivation remain the same, as F̃(D;r) is zero for any r<D/2. However, the relationship expressed in Eq. (4) becomes


where Y(r)sin-1(Dmin2r), Sn(y)=0ysinn(x)dx, and E is the error in relating the nth moments of S(D) and P(r). Since P(r) is unknown, E cannot be computed a priori. The function Sn(Y(r)) expresses the percentage of chords formed from floes of size r that would be smaller than Dmin, although it is not readily expressed as a function of n. The most pathological distribution is when P(r) is a delta function at rmin, P(r)=δ(r-rmin), Y(rmin)=π/2, and E=1 as no chord lengths would be measured.

We can compute the error function for any delta function distribution as

(B4) E ( δ ( r - r * ) ; n ) = S n ( Y ( r * ) ) S n ( π 2 ) ,

and the misfit is the proportion of the integral of sin n(x) between 0 and Y(r*). Because sin (x) is monotonically increasing from x=0 to π∕2, the integral of Sn is bounded above:

(B5) S n ( Y ( r * ) ) Y ( r * ) sin n ( Y ( r * ) ) = Y ( r * ) D min 2 r * n ,

and the misfit error is bounded above by

(B6) E ( δ ( r - r * ; n ) ) D min 2 r * n Y ( r * ) B ( n + 1 2 , 1 2 ) .

The reciprocal of B is equal to π at n=0 and decreases sub-linearly, and so away from rmin the error term decays exponentially with n and is small even for nearly pathological distributions (for n=1, r*=Dmin, for example, Eπ/2414%. Knowing the distribution of errors behaves in this way allows us to establish upper bounds by integrating P as a sum of δ functions.

We note that increasing resolution of floe chords will result in tighter bounds on this error. When Y(r)*1, which occurs when r*Dmin2sin(1)0.59Dmin, we can exploit a tighter bound using the fact that sin n(x)≤xn,

(B7) S n ( Y ( r * ) ) Y ( r * ) n + 1 n + 1 Y ( r * ) D min 2 r * n .

Using the same example as above (n=1, r*=Dmin) bounds the error Eπ2/1447%. A real-world distribution of floe sizes must have a peak value above zero; thus by increasing the sampling resolution (say, for example, to near the size of pancakes, i.e., Dmin≈20 m or less, approached by the ICESat-2 altimeter), this bound takes over and errors are reduced substantially.

We can explicitly solve Eq. (B3) for distributions with power-law tails. These distributions are peaked at the minimum floe size, and so will have high moment error. For power laws with α=-1,-2,-3, or −4, E(P(r;α,rmin),1) is 1, 4, 16, or 25 %. For n=2, E(P(r;α,rmin),2) is .003, .04, 2, or 9.6 %: the increase in error with decreasing α is because sharper power-law slopes concentrate most of the distribution towards the smallest scale.

Appendix C: Maximum likelihood estimation for chord length distributions

Given a set of floe chords {D}i and an estimate of the beginning of a power-law tail D*, we would like to find the most likely power-law floe size distribution P(r;α,rmin) that generated them. As discussed in Appendix A, moments of the FSD and FCD are related by a multiplicative factor, and the distributions themselves will share the same power-law exponent. Thus we may test the power-law hypothesis directly on the FCD S(D). The power-law hypothesis means that S(D) is of the form

(C1) S ( D ) = ( α - 1 ) D * D D * - α .

Following Muniruzzaman (1957) and Clauset et al. (2009) (see also the derivation in Stern et al.2018a), we compute the log-likelihood of the observations for a given α (Eq. 10),


As the natural log is monotonically increasing in its argument, to find the most likely α, denoted α^, we take the derivative with respect to α and set to zero,

(C4) 1 α - 1 + ln ( D * ) = 1 N i = 1 N ln D i D * ,

which resolves as a solution for the most likely α:

(C5) α ^ = 1 + N i = 1 N ln D i D * .

The above analysis concerns the most likely α that explains the FCD. We may ask a separate question: what is the most likely α, which we define as αP, that would explain the FSD, given the explicit relationship that can be derived between S(D) and a power-law-distributed P(r) examined in Eq. (10).

(C6) S ( D ) = C B 1 2 , α 2 2 α - 1 π D - α ,

where C is unknown. Repeating the above analysis,

(C7) L ln i = 1 N S ( D i ) = ln C N B 1 2 , α P 2 N 2 ( α P - 1 ) π N i = 1 N D i - α P

(C8) = N ln C + N ln B 1 2 , α P 2 + N ( α P - 1 ) ln 2 - N ln π - α P i N ln D i .

Next we take the derivative of with respect to αP and setting to zero. We use the fact that B(x,y)=B(y,x) and B(x,y)x=B(x,y)ψ(x)-ψ(x+y), where ψ is the digamma function, to find

(C9) ln B 1 2 , α P 2 α P = 1 2 ψ α P 2 - ψ α P + 1 2 .

The maximum likelihood αP is the solution to the transcendental equation,

(C10) 1 2 ψ α P 2 - ψ α P + 1 2 + ln 2 = 1 N i = 1 N ln D i ,

which is an alternative method for obtaining the FSD scaling.

Appendix D: Averaging and segmenting FSD statistics

Due to limitations in the number of floe chords recorded at any particular location over time, we do not include all geographic locations when computing hemispheric means. Averaging is performed by including only geographic regions where there are at least 25 recorded floe chords. The area being averaged over is thus not fixed in time. For seasonal cycle plots, we only include months which have enough measurements for all fully sampled CryoSat-2 years (2011–2018). For annual averages, we include only those years where all CryoSat-2 months (excluding June–September) have enough measurements.

When masking additional regions to perform the model–observation comparisons in Fig. 6, we note that because the Roach et al. (2018a) model does not include processes that fragment larger floes into smaller floes in the absence of ocean surface waves, regions in the interior Arctic without wave activity have nearly all sea ice area belonging to the highest floe size categories. Nearly all regions where wave fracture is an active process also have representative radii below about 10 km (Roach et al.2019). We define regions that do not experience wave fracture as those with an abnormally high simulated representative radius, which we choose to be the 22nd floe size category (r=18.6 km) or above. The mask and comparisons in Fig. 6 are made by excluding all such areas.


The supplement related to this article is available online at:

Author contributions

CH derived the mathematical theory and wrote the paper. LR built and performed the climate model simulation. RT, AR, and AS provided and interpreted the CryoSat-2 data. KH, CG, CB, and BK contributed to the study design. All authors have participated in paper preparation.

Competing interests

The authors declare no competing interests.


CH was supported by the NOAA Climate and Global Change Postdoctoral Fellowship Program, administered by UCAR's Cooperative Programs for the Advancement of Earth System Science (CPAESS), sponsored in part through cooperative agreement number NA16NWS4620043, years 2017–2021, with the National Oceanic and Atmospheric Administration (NOAA) and the U.S. Department of Commerce (DOC). CH, CG, and KH thank the American Mathematical Society for their support through the Mathematics Research Community “Differential Equations, Probability, and Sea Ice”, funded by NSF grants 1321794 and 1641020. LR was funded via Marsden contract VUW‐1408 and the New Zealand Deep South National Science Challenge, MBIE contract number C01X1445. CMB was supported by the National Science Foundation grant PLR-1643431. BFK was supported by ONR grant N00014-17-1-2963 and NSF grant 1350795. RT, AR, and AS were supported by the UK NERC Centre for Polar Observation and Modelling and the European Space Agency.

Financial support

This research has been supported by the National Oceanic and Atmospheric Administration, Climate Program Office (grant no. NA16NWS4620043), the National Science Foundation, Division of Mathematical Sciences (grant nos. 1321794, 1641020, 1350795, and 1643431), the Office of Naval Research (grant no. N00014-17-1-2963), the Marsden Fund (grant no. VUW‐1408), and the New Zealand Deep South National Challenge (grant no. C01X1445).

Review statement

This paper was edited by Jennifer Hutchings and reviewed by Thomas Armitage and one anonymous referee.


Alberello, A., Onorato, M., Bennetts, L., Vichi, M., Eayrs, C., MacHutchon, K., and Toffoli, A.: Brief communication: Pancake ice floe size distribution during the winter expansion of the Antarctic marginal ice zone, The Cryosphere, 13, 41–48,, 2019. a

Allard, R. A., Farrell, S. L., Hebert, D. A., Johnston, W. F., Li, L., Kurtz, N. T., Phelps, M. W., Posey, P. G., Tilling, R., Ridout, A., and Wallcraft, A. J.: Utilizing CryoSat-2 sea ice thickness to initialize a coupled ice-ocean modeling system, J. Geophys. Res.-Oceans, 62, 1265–1280,, 2018. a

Armitage, T. W. K. and Davidson, M. W. J.: Using the Interferometric Capabilities of the ESA CryoSat-2 Mission to Improve the Accuracy of Sea Ice Freeboard Retrievals, IEEE Trans. Geosci. Remote Sens., 52, 529–536,, 2014. a

Armitage, T. W. K., Bacon, S., and Kwok, R.: Arctic Sea Level and Surface Circulation Response to the Arctic Oscillation, Geophys. Res. Lett., 45, 6576–6584,, 2018. a

Batchelor, G. K.: The theory of homogeneous turbulence, Cambridge University Press, 1953. a

Bateson, A. W., Feltham, D. L., Schröder, D., Hosekova, L., Ridley, J. K., and Aksenov, Y.: Impact of floe size distribution on seasonal fragmentation and melt of Arctic sea ice, The Cryosphere Discuss.,, in review, 2019. a

Birnbaum, G. and Lüpkes, C.: A new parameterization of surface drag in the marginal sea ice zone, Tellus A, 54, 107–123,, 2002. a

Clauset, A., Young, M., and Gleditsch, K. S.: On the Frequency of Severe Terrorist Events, J. Conf. Resolut., 51, 58–87, 2007. a, b

Clauset, A., Shalizi, C. R., Newman, M. E. J., Rohilla Shalizi, C., and J Newman, M. E.: Power-Law Distributions in Empirical Data, SIAM Rev., 51, 661–703,, 2009. a, b, c, d, e

Day, J. J., Hawkins, E., and Tietsche, S.: Will Arctic sea ice thickness initialization improve seasonal forecast skill?, Geophys. Res. Lett., 41, 7566–7575,, 2014. a

Feltham, D. L.: Sea Ice Rheology, Ann. Rev. Fluid Mechan., 40, 91–112,, 2008. a

Gherardi, M. and Lagomarsino, M. C.: Characterizing the size and shape of sea ice floes, Sci. Rep., 5, 10226,, 2015. a

Guerreiro, K., Fleury, S., Zakharova, E., Kouraev, A., Rémy, F., and Maisongrande, P.: Comparison of CryoSat-2 and ENVISAT radar freeboard over Arctic sea ice: toward an improved Envisat freeboard retrieval, The Cryosphere, 11, 2059–2073,, 2017. a

Herman, A.: Molecular-dynamics simulation of clustering processes in sea-ice floes, Phys. Rev. E, 84, 1–11,, 2011. a

Herman, A., Evers, K.-U., and Reimer, N.: Floe-size distributions in laboratory ice broken by waves, The Cryosphere, 12, 685–699,, 2018 a

Horvat, C. and Tziperman, E.: A prognostic model of the sea-ice floe size and thickness distribution, The Cryosphere, 9, 2119–2134,, 2015. a, b, c, d, e

Horvat, C. and Tziperman, E.: The evolution of scaling laws in the sea ice floe size distribution, J. Geophys. Res.-Oceans, 122, 7630–7650,, 2017. a, b, c, d, e

Horvat, C. and Tziperman, E.: Understanding Melting due to Ocean Eddy Heat Fluxes at the Edge of Sea-Ice Floes, Geophys. Res. Lett., 45, 9721–9730,, 2018. a

Horvat, C., Tziperman, E., and Campin, J.-M.: Interaction of sea ice floe size, ocean eddies, and sea ice melting, Geophys. Res. Lett., 43, 8083–8090,, 2016. a

Hunke, E. C., Lipscomb, W. H., Turner, A. K., Jeffery, N., and Elliott, S.: CICE : the Los Alamos Sea Ice Model Documentation and Software User's Manual Version 5.1 LA-CC-06-012, Tech. rep., Los Alamos National Laboratory, available at: (last access: 1 November 2019), 2015. a

Hwang, B., Wilkinson, J., Maksym, E., Graber, H. C., Schweiger, A., Horvat, C., Perovich, D. K., Arntsen, A. E., Stanton, T. P., Ren, J., Wadhams, P., Maksym, T., Graber, H. C., Schweiger, A., Horvat, C., Perovich, D. K., Arntsen, A. E., Timothy, P., Ren, J., and Wadhams, P.: Winter-to-summer transition of Arctic sea ice breakup and floe size distribution in the Beaufort Sea, Elementa: Sci. Anthrop., 5, 40,, 2017. a, b

JRA-55: JRA-55: Japanese 55-year Reanalysis, Daily 3-Hourly and 6-Hourly Data, Research Data Archive at the National Center for Atmospheric Research, Computational and Information Systems Laboratory, 2013. a

Key, J.: Estimating the area fraction of geophysical fields from measurements along a transect, IEEE Tran. Geosci. Remote Sens., 31, 1099–1102,, 1993. a

Key, J. and Peckham, S.: Probable errors in width distributions of sea ice leads measured along a transect, J. Geophys. Res., 96, 18417,, 1991. a

Kwok, R.: Arctic sea ice thickness, volume, and multiyear ice coverage: losses and coupled variability (1958–2018), Environ. Res. Lett., 13, 105005,, 2018. a

Laxon, S., Peacock, N., and Smith, D.: High interannual variability of sea ice thickness in the Arctic region, Nature, 425, 947–950,, 2003. a

Laxon, S. W., Giles, K. A., Ridout, A. L., Wingham, D. J., Willatt, R., Cullen, R., Kwok, R., Schweiger, A., Zhang, J., Haas, C., Hendricks, S., Krishfield, R., Kurtz, N., Farrell, S., and Davidson, M.: CryoSat-2 estimates of Arctic sea ice thickness and volume, Geophys. Res. Lett., 40, 732–737,, 2013. a

Lindsay, R. W. and Rothrock, D. A.: Arctic sea ice leads from advanced very high resolution radiometer images, J. Geophys. Res., 100, 4533,, 1995. a

Lüpkes, C. and Birnbaum, G.: Surface drag in the Arctic marginal sea-ice zone: A comparison of different parameterisation concepts, Bound.-Lay. Meteorol., 117, 179–211,, 2005. a

Mandelbrot, B. B. and Wheeler, J. A.: The Fractal Geometry of Nature, Vol. 51, W. H. Freeman,, 1983. a

Manucharyan, G. E. and Thompson, A. F.: Submesoscale Sea Ice-Ocean Interactions in Marginal Ice Zones, J. Geophys. Res.-Oceans, 122, 9455–9475,, 2017. a

Muniruzzaman, A. N. M.: On Measures of Location and Dispersion and Tests of Hypotheses in a Pareto Population, Calc. Stat. Assoc. Bull., 7, 115–123,, 1957. a, b

Nere, N. K., Ramkrishna, D., Parker, B. E., Bell, W. V., and Mohan, P.: Transformation of the Chord-Length Distributions to Size Distributions for Nonspherical Particles with Orientation Bias †, Indust. Eng. Chem. Res., 46, 3041–3047,, 2007. a

Peacock, N. R. and Laxon, S. W.: Sea surface height determination in the Arctic Ocean from ERS altimetry, J. Geophys. Res.-Oceans, 109, C07001,, 2004. a, b

Perovich, D. K. and Jones, K. F.: The seasonal evolution of sea ice floe size distribution, J. Geophys. Res.-Oceans, 119, 8767–8777,, 2014. a, b

Pons, M.-N., Milferstedt, K., and Morgenroth, E.: Modeling of chord length distributions, Chem. Eng. Sci., 61, 3962–3973,, 2006. a

Quartly, G. D., Rinne, E., Passaro, M., Andersen, O. B., Dinardo, S., Fleury, S., Guillot, A., Hendricks, S., Kurekin, A. A., Müller, F. L., Ricker, R., Skourup, H., and Tsamados, M.: Retrieving Sea Level and Freeboard in the Arctic: A Review of Current Radar Altimetry Methodologies and Future Perspectives, Remote Sens., 11, 881,, 2019. a

Richter-Menge, J. A.: Relating arctic pack ice stress and deformation under winter conditions, J. Geophys. Res., 107, 8040,, 2002. a

Roach, L. A., Horvat, C., Dean, S. M., and Bitz, C. M.: An emergent sea ice floe size distribution in a global coupled ocean–sea ice model, J. Geophys. Res.-Oceans, 123, 4322–4337,, 2018a. a, b, c, d, e, f, g, h, i

Roach, L. A., Smith, M. M., and Dean, S. M.: Quantifying growth of pancake sea ice floes using images from drifting buoys, J. Geophys. Res.-Oceans, 123, 2851–2866,, 2018b. a

Roach, L., Bitz, C., Horvat, C., and Dean, S.: Advances in modelling interactions between sea ice and ocean surface waves, Journal of Advances in Modeling Earth Systems, 1–14, in review, 2019. a, b

Roberts, A. F., Hunke, E. C., Kamal, S. M., Lipscomb, W. H., Horvat, C., and Maslowski, W.: A Variational Method for Sea Ice Ridging in Earth System Models, J. Adv. Model. Earth Syst., p. 2018MS001395,, 2019. a

Rothrock, D. A. and Thorndike, A. S.: Measuring the Sea Ice Floe Size Distribution, J. Geophys. Res., 89, 6477–6486,, 1984. a, b, c, d, e

Rynders, S., Aksenov, Y., Feltham, D., Nurser, G., and Naveira Garabato, A.: Modelling MIZ dynamics in a global model, in: EGU General Assembly Conference Abstracts, p. 1004, 2016. a

Sandberg Sorensen, L., Simonsen, S., Langley, K., Gray, L., Helm, V., Nilsson, J., Stenseng, L., Skourup, H., Forsberg, R., and Davidson, M.: Validation of CryoSat-2 SARIn Data over Austfonna Ice Cap Using Airborne Laser Scanner Measurements, Remote Sens., 10, 1354,, 2018. a

Schröder, D., Feltham, D. L., Tsamados, M., Ridout, A., and Tilling, R.: New insight from CryoSat-2 sea ice thickness for sea ice modelling, The Cryosphere, 13, 125–139,, 2019. a

Schulson, E. M. and Hibler, W. D.: The fracture of ice on scales large and small: Arctic leads and wing cracks, J. Glaciol., 37, 319–322,, 1991. a

Skourup, H., Simonsen, S. B., Sorensen, L., Bella, A., Forsberg, R., Hvidegaard, S., and Helm, V.: ESA CryoVEx/EU ICE-ARC 2016-Airborne Field Campaign with ASIRAS Radar and Laser Scanner over Austfonna, Fram Strait and the Wandel Sea, Tech. rep., National Space Institute, Technical University of Denmark, 2017. a

Smith, M. and Thomson, J.: Scaling observations of surface waves in the Beaufort Sea, Elementa: Science of the Anthropocene, 4, 000097,, 2016. a

Squire, V. A.: Of ocean waves and sea-ice revisited, Cold Reg. Sci. Technol., 49, 110–133,, 2007. a

Squire, V. A., Dugan, J. P., Wadhams, P., Rottier, P. J., and Liu, A. K.: Of Ocean Waves and Sea Ice, Annu. Rev. Fluid Mechan., 27, 115–168,, 1995. a

Steele, M.: Sea ice melting and floe geometry in a simple ice-ocean model, J. Geophys. Res.-Oceans, 97, 17729–17738,, 1992. a, b

Steer, A., Worby, A., and Heil, P.: Observed changes in sea-ice floe size distribution during early summer in the western Weddell Sea, Deep Sea Res. Part II, 55, 933–942,, 2008. a

Stern, H. L., Schweiger, A. J., Stark, M., Zhang, J., Steele, M., and Hwang, B.: Seasonal evolution of the sea-ice floe size distribution in the Beaufort and Chukchi seas, Elem. Sci. Anth., 6, 48,, 2018a. a, b, c

Stern, H. L., Schweiger, A. J., Zhang, J., and Steele, M.: On reconciling disparate studies of the sea-ice floe size distribution, Elem. Sci. Anth., 6, 49,, 2018b. a, b, c, d, e, f

Tilling, R., Ridout, A., and Shepherd, A.: Near-real-time Arctic sea ice thickness and volume from CryoSat-2, The Cryosphere, 10, 2003–2012,, 2016. a, b

Tilling, R., Ridout, A., and Shepherd, A.: Estimating Arctic sea ice thickness and volume using CryoSat-2 radar altimeter data, Adv. Space Res.,, 2018. a, b, c, d

Tilling, R., Ridout, A., and Shepherd, A.: Assessing the impact of lead and floe sampling on Arctic sea ice thickness estimates from Envisat and CryoSat‐2, J. Geophys. Res.-Oceans, 124,, 2019a. a

Tilling, R., Ridout, A., and Shepherd, A.: Assessing the impact of lead and floe sampling on Arctic sea ice thickness estimates from Envisat and CryoSat‐2, J. Geophys. Res.-Oceans, 124,, 2019b. a, b, c, d, e

Tolman, H. L. H. L.: User manual and system documentation of WAVEWATCH III TM version 3.14, Technical Note 276, 276, 194, available at:\%0Apapers3://publication/uuid/E1C39B58-DBCB-4F8A-ADCD-F2D2210DDC46 (last access: 1 November 2019), 2009. a

Toyota, T., Takatsuji, S., and Nakayama, M.: Characteristics of sea ice floe size distribution in the seasonal ice zone, Geophys. Res. Lett., 33, L02616,, 2006. a, b

Toyota, T., Haas, C., and Tamura, T.: Size distribution and shape properties of relatively small sea-ice floes in the Antarctic marginal ice zone in late winter, Deep-Sea Res. Part II, 58, 1182–1193,, 2011. a, b, c

Toyota, T., Kohout, A., and Fraser, A. D.: Formation processes of sea ice floe size distribution in the interior pack and its relationship to the marginal ice zone off East Antarctica, Deep Sea Res. Pt II, 131, 28–40,, 2016. a

Tsamados, M., Feltham, D. L., Schroeder, D., Flocco, D., Farrell, S. L., Kurtz, N., Laxon, S. W., and Bacon, S.: Impact of Variable Atmospheric and Oceanic Form Drag on Simulations of Arctic Sea Ice*, J. Phys. Oceanogr., 44, 1329–1353,, 2014. a

Virkar, Y. and Clauset, A.: Power-law distributions in binned empirical data, The Ann. Appl. Stat., 8, 89–119,, 2014. a, b, c

Wadhams, P., Squire, V. a., Goodman, D. J., Cowan, A. M., and Moore, S. C.: The attenuation rates of ocean waves in the marginal ice zone, J. Geophys. Res., 93, 6799,, 1988. a

Wernecke, A. and Kaleschke, L.: Lead detection in Arctic sea ice from CryoSat-2: quality assessment, lead area fraction and width distribution, The Cryosphere, 9, 1955–1968,, 2015. a

Wilchinsky, A. V. and Feltham, D. L.: Modelling the rheology of sea ice as a collection of diamond-shaped floes, J. Non-Newton. Fluid Mechan., 138, 22–32,, 2006. a

Wilchinsky, A. V. and Feltham, D. L.: Modeling Coulombic failure of sea ice with leads, J. Geophys. Res.-Oceans, 116, C08040,, 2011. a

Williams, T. D., Bennetts, L. G., Squire, V. A., Dumont, D., and Bertino, L.: Wave-ice interactions in the marginal ice zone. Part 1: Theoretical foundations, Ocean Modell., 71, 81–91,, 2013. a

Zhang, J., Schweiger, A., Steele, M., and Stern, H.: Sea ice floe size distribution in the marginal ice zone: Theory and numerical experiments, J. Geophys. Res.-Oceans, 120, 3484–3498,, 2015. a, b

Short summary
Changes in the floe size distribution (FSD) are important for sea ice evolution but to date largely unobserved and unknown. Climate models, forecast centres, ship captains, and logistic specialists cannot currently obtain statistical information about sea ice floe size on demand. We develop a new method to observe the FSD at global scales and high temporal and spatial resolution. With refinement, this method can provide crucial information for polar ship routing and real-time forecasting.