To better simulate the seasonal evolution of sea ice in the Arctic, with particular attention to the marginal ice zone, a sea ice model of the distribution of ice thickness, floe size, and enthalpy was implemented into the Pan-arctic Ice–Ocean Modeling and Assimilation System (PIOMAS). Theories on floe size distribution (FSD) and ice thickness distribution (ITD) were coupled in order to explicitly simulate multicategory FSD and ITD distributions simultaneously. The expanded PIOMAS was then used to estimate the seasonal evolution of the Arctic FSD in 2014 when FSD observations are available for model calibration and validation. Results indicate that the simulated FSD, commonly described equivalently as cumulative floe number distribution (CFND), generally follows a power law across space and time and agrees with the CFND observations derived from TerraSAR-X satellite images. The simulated power-law exponents also correlate with those derived using MODIS images, with a low mean bias of –2%. In the marginal ice zone, the modeled CFND shows a large number of small floes in winter because of stronger winds acting on thin, weak first-year ice in the ice edge region. In mid-spring and summer, the CFND resembles an upper truncated power law, with the largest floes mostly broken into smaller ones; however, the number of small floes is lower than in winter because floes of small sizes or first-year ice are easily melted away. In the ice pack interior there are fewer floes in late fall and winter than in summer because many of the floes are “welded” together into larger floes in freezing conditions, leading to a relatively flat CFND with low power-law exponents. The simulated mean floe size averaged over all ice-covered areas shows a clear annual cycle, large in winter and smaller in summer. However, there is no obvious annual cycle of mean floe size averaged over the marginal ice zone. The incorporation of FSD into PIOMAS results in reduced ice thickness, mainly in the marginal ice zone, which improves the simulation of ice extent and yields an earlier ice retreat.

Significant changes have occurred in Arctic sea ice, whether in the ice pack interior (IPI) or in the marginal ice zone (MIZ) (e.g., Comiso, 2012; Strong and Rigor, 2013; Meier et al., 2014; Kwok and Cunningham, 2015; Lindsay and Schweiger, 2015). The MIZ is a transition region from open water to pack ice and a research focus of the recent Office of Naval Research (ONR) MIZ Program. The MIZ (defined hereafter as the areas with ice concentration between 0.15 and 0.80 following Strong and Rigor, 2013) lies in the subarctic seas in winter and transitions into the interior of the Arctic Basin in summer. Ice in the MIZ in these different regions breaks up more easily in response to winds and currents. It is also vulnerable to ocean surface waves that may form in the open water, resulting from strong winds and frequent storms, and propagate into the ice field (Squire et al., 1995; Squire, 2007; Kohout et al., 2014). The ice tends to attenuate the incoming waves if any (Wadhams et al., 1988; Meylan et al., 2014), while the waves tend to bend ice repeatedly. The ice breaks if the stresses induced by bending exceed its flexural strength or if the repeated bending leads to fatigue failure (e.g., Langhorne et al., 1998).

Winds and currents may cause ice floes to break in the MIZ as well as in the IPI in the absence of waves if the wind- or current-induced ice internal stress exceeds the tensile or compressive strength of the ice. For example, during the ONR MIZ 2014 field campaign (http://www.apl.washington.edu/project/project.php?id=miz), floe breakup was observed during a period of low levels of wave activity (Arntsen et al., 2015). The floe breakup occurred in late July as melt accelerated and networks of melt ponds formed weak zones where modest wind forcing was sufficient for floes to break apart (Arntsen et al., 2015). This observation suggests the important role of wind forcing when wave forcing is not dominant, particularly in summer when the ice cover is weakened by thermodynamic forcing. The role of wind and current forcing in ice fragmentation is also reflected in the observations of ice floe fields in the IPI where there is likely little wave activity (Rothrock and Thorndike, 1984; Stern et al., 2014). Once an ice floe is broken, it becomes floes of smaller sizes, a process of floe size redistribution via mechanical forcing. The combination of wind-, current-, and wave-induced ice fragmentation makes the MIZ as well as much of the IPI full of ice floes with varying sizes as well as thicknesses (Rothrock and Thorndike, 1984; Wadhams, 1986; Stern et al., 2014).

Thus, to better describe the state of sea ice in a given area with floes of varying sizes and thicknesses, both a floe size distribution (FSD) and an ice thickness distribution (ITD) are needed. An ITD is able to describe the fractions of open water/leads and various ice thicknesses in that area (Thorndike et al., 1975; also see Hibler, 1980), while a FSD is able to represent the ice floes of various sizes ranging from meters to kilometers (Rothrock and Thorndike, 1984; Holt and Martin, 2001; Herman, 2010). The ITD is often described as the area fraction over a range of thicknesses (Thorndike et al., 1975; Hibler, 1980). The FSD may also be described as the area fraction over a range of floe sizes or caliper diameters — a measure of floe size (Rothrock and Thorndike, 1984). However, it is more often described equivalently as the number of floes per unit area over a floe size range (i.e., floe number distribution or FND) by many observational studies analyzing satellite images and aerial photographs of ice floes (Rothrock and Thorndike, 1984; Holt and Martin, 2001; Toyota et al., 2006; Steer et al., 2008). These studies indicate that the number of floes per unit area with caliper diameters not smaller than l, or the cumulative floe number distribution (CFND), can be described by a power law function N(l) ∝ l–α, where N(l) is the CFND and l is the caliper diameter of a floe. The power law function is a straight line in a log-log plot and –α is the slope of the line. Thus, the power law obeying CFND is characterized by a single exponent α over all floe sizes.

The Thorndike et al. (1975) ITD theory has been incorporated increasingly in operational forecast and climate models. However, no large-scale sea ice models have been able to simulate explicitly the evolution of FSD, not to mention simulating FSD and ITD jointly. The necessity of a better description of the state of sea ice in the MIZ has motivated research to develop coupled FSD and ITD theories for explicit representation of FSD and ITD in the MIZ (Zhang et al., 2015; Horvat and Tziperman, 2015). To our knowledge, the study by Zhang et al. (2015) is the first to develop a FSD theory coupled to the ITD theory of Thorndike et al. (1975) in order to simulate explicitly the evolution of both FSD and ITD simultaneously. The FSD theory includes a FSD function and a FSD conservation equation in parallel with the ITD equation of Thorndike et al. (1975). The FSD equation describes changes in FSD due to ice advection, thermodynamic growth, and lateral melting. It also includes changes in FSD because of mechanical redistribution of floe size due largely to ocean surface wave-induced ice fragmentation and wind- and current-induced ice ridging and lead opening calculated in the ITD equation. The FSD theory, tested in a simplified ITD and FSD sea ice model with idealized numerical experiments, is able to simulate FSD that follows a power law as observed by satellites and airborne surveys (Rothrock and Thorndike, 1984; Holt and Martin, 2001; Toyota et al., 2006; Steer et al., 2008; Perovich and Jones, 2014). The simulated values of the exponent of the power law are also generally in the range of the observations (Zhang et al., 2015).

However, the results of Zhang et al. (2015) are based on a zero-dimensional ITD and FSD sea ice model in idealized numerical experiments without ice dynamics and thermodynamics. Thus, it is necessary to expand the Zhang et al. (2015) study and incorporate the FSD theory into a two-dimensional dynamic–thermodynamic sea ice model. Here, we implemented the FSD theory into the Pan-arctic Ice–Ocean Modeling and Assimilation System (PIOMAS) to study the evolution of FSD in the MIZ and IPI of the Arctic Ocean. The sea ice component of PIOMAS employs the multicategory ice thickness and enthalpy distribution (TED) model (Zhang and Rothrock, 2001, 2003), a dynamic–thermodynamic model that also explicitly simulates sea ice ridging. Thus, PIOMAS is useful for simulating FSD and ITD jointly. This study is, to our knowledge, unique in its attempt to explicitly simulate the evolution of FSD in a dynamic-thermodynamic sea ice model on a pan-Arctic scale. It is our hope that this study will facilitate further model development and improvement in representing ice breakup processes. After briefly describing PIOMAS and the TED sea ice model, we present the implementation and parameterization of the FSD equation in the dynamic-thermodynamic TED model framework, followed by model configuration, forcing, initialization and simulations. An examination of model results, including model validation, leads to conclusions about the seasonal evolution of sea ice in both the MIZ and IPI.

PIOMAS is a large-scale coupled sea ice–ocean model. The PIOMAS ocean model is based on the Parallel Ocean Program (POP; Smith et al., 1992) that has been modified to impose open boundary conditions (Zhang and Steele, 2007). The PIOMAS TED sea ice model adopts the Thorndike et al. (1975) ITD theory (also see Hibler, 1980). In the ITD theory, the ice mass conservation is described by an ITD equation,

$∂gh∂t=−∇·(ugh)−∂(fhgh)∂h+Ψ+FL,$
(1)

where gh is the ice thickness distribution function, t is time, u is the ice velocity vector, fh is ice growth or melt rate, h is ice thickness, ψ is a mechanical thickness redistribution function for ridging, and FL is a source term for lateral melting (Hibler, 1980).

The ITD theory is augmented by an ice enthalpy distribution theory to conserve the thermal energy of ice (Zhang and Rothrock, 2001, 2003). As a result, the TED sea ice model can be used to integrate over multiple sub-grid categories each for ice thickness and ice enthalpy (e.g., Zhang et al., 2012). The TED sea ice model integration also includes multiple categories of snow depth following Flato and Hibler (1995; also see Zhang and Rothrock, 2003). In addition, the TED sea ice model employs a teardrop viscous-plastic ice rheology that determines the relationship between ice internal stress and ice deformation (Zhang and Rothrock, 2005), a mechanical redistribution function that determines ice ridging (Hibler, 1980; Rothrock, 1975; Thorndike et al., 1975), and a dynamics model to solve the ice momentum equation (Zhang and Hibler, 1997).

### a. Brief review of the FSD equation

Here, we develop a thickness, floe size, and enthalpy distribution (TFED) sea ice model by expanding the TED model to incorporate the Zhang et al. (2015) FSD equation in the PIOMAS framework in order to explicitly simulate the evolution of FSD. The FSD equation is given by

$∂gl∂t=−∇·(ugl)−∂(flgl)∂l+Φ,$
(2)

where gl(l) is the FSD function, l is caliper diameter, fl is the rate of change in floe size (caliper diameter l), and Φ is the mechanical floe size redistribution function. With gl(l) being the FSD function, gl(l)dl represents the area fraction of a floe size category centered at size l. The FSD function can be converted to the commonly used CFND function by $N(l)=∫0∞gl(l′)/(0.66l′2)dl′$, where 0.66l2 is the approximate area of a floe with caliper diameter l following Rothrock and Thorndike (1984). The second term in (2) describes the change in FSD due to ice advection. The third term describes the change in FSD due to thermodynamic growth or decay represented by freezing or lateral melting. The fourth term, the mechanical floe size redistribution function, describes the change in FSD due to open water creation or lead opening, ridging, and fragmentation.

The mechanical floe size redistribution function is separated into three terms, Φ = Φ0 + Φr + Φf, representing the mechanical changes in FSD due to open water creation (Φ0), ice ridging (Φr), and ice fragmentation (Φf). These terms are given in Zhang et al. (2015). Here we describe the ice fragmentation term Φf again for context. In the MIZ, ice fragmentation is largely induced by winds, currents, and ocean surface waves. Winds and currents may also force ice to deform and break in the IPI. The ice fragmentation term (Φf) is given by

$Φf=−Q(l)gl(l)+∫0∞β(l′,l)Q(l′)gl(l′)dl′,$
(3)

where Q(l) is the redistribution probability function and β(l1,l2) is a redistributor of FSD. The redistributor of FSD β(l1,l2) specifies how ice is transferred from one floe size category to another by breaking (Zhang et al., 2015). The redistribution probability function Q specifies whether ice breakup takes place and, if so, which floe size categories are to participate in the breakup processes and in the mechanical floe size redistribution. It is given by

$Q(l)=max[(1−∫0∞gl(l′)dl′/cb),0],$
(4)

where the constant cb is the floe size redistribution participation factor that specifies an area fraction of ice to participate in breaking. Thus the participation factor cb plays a prominent role in determining the redistribution probability function Q and therefore the magnitude of ice breakup and ensuing mechanical floe size redistribution (Zhang et al., 2015).

### b. Parameterization of floe size redistribution participation factor

Because winds and currents may cause ice to deform and break, the value of cb depends on winds which drive sea ice and upper ocean circulation. In the MIZ during periods of high levels of wave activity, the value of cb may also depend on wave conditions (wave energy, frequency, and direction). It would be ideal for wave conditions to be generated explicitly by wave models that are also able to capture wave–ice interactions. However, PIOMAS currently does not have a wave model component. Therefore it is necessary to parameterize cb by including the effect of waves in addition to those of winds and currents. Note that wave conditions depend on wind speed and fetch (the distance of open water over which winds are blowing) (e.g., Squire et al., 1995; Thomson and Rogers, 2014). They also depend on sea ice conditions because wave propagation and attenuation under ice are affected, through wave–ice interactions, by ITD and FSD which also control the flexural strength and hence the bending failure of sea ice in the MIZ (Wadhams et al., 1988; Squire et al., 1995, 2009; Kohout and Meylan, 2008; Dumont et al., 2011; Meylan et al., 2014). The ITD and FSD also control the tensile and compressive strength and hence the deforming failure of sea ice under wind and current forcing (Arntsen et al., 2015). In other words, in the MIZ the value of cb is a function of wind speed, fetch, ITD, and FSD. Such a function may be expressed by

$cb=kUa/max(hm,hc)exp[−a(1−fo)−b(1−lm/lmax)]Δt,$
(5)

where k, a, and b are fixed positive empirical constants (dimensionless), Ua is surface wind speed (m s−1), hm$(∫0∞gh(h)hdh)$ is mean ice thickness (m), hc is a cutoff ice thickness set to be 0.2 m (a typical ice thickness in the MIZ in winter; see Figures 5a–5b below), fo is the open water area fraction (averaged over each ocean grid cell and its eight surrounding ocean cells in the model to mimic fetch), lm$(=∫0∞gl(l)ldl)$ is mean ice floe size, lmax is the center of the largest floe size category at the high end of the FSD in the model (see table 1 in Zhang et al., 2015), and Δt is the model integration time step interval (s).

Table 1.
Numerical parameters used in the control and sensitivity simulations

aControl (CNTL) and sensitivity (SEN)

bConstants a and b in (5) and welding threshold wd

cParameter not applied, as FSD was excluded for comparison with CNTL

The formulation (5) makes the participation factor cb in a given area strongly dependent on local wind speed, open water fraction (related to fetch), ITD, and FSD. The value of cb increases with increasing wind speed, increasing open water fraction (i.e., our parameterization for fetch), and decreasing ice thickness. However, the value of cb decreases with decreasing floe size, reflecting the observations that larger floes are easier to break because they are subject to larger flexure-induced stresses or strains, while smaller floes are likely to ride with waves with little bending (e.g., Meylan and Squire, 1994; Squire et al., 1995). The value of cb is capped to the value of ice concentration in the event of a storm with high wind speeds that make the value of cb greater than ice concentration. This capping means that all ice floes are broken during such a strong storm.

Note that the expression of (5) as a function of wind speed, open water fraction, and sea ice thickness and floe size make it applicable also to the IPI where there are often low levels of wave activity. This applicability is because the ice fields in the central Arctic away from the MIZ, which appear more as a continuum, may also break because of wind- or current-induced ice deformation, resulting in leads or fractures interspersed (e.g., Wadhams, 1981; Hibler, 2001; Kwok et al., 2008). Here, the fractured ice fields are also described by the FSD equation, with a unified parameterization of cb.

### c. Change in FSD due to thermodynamics

The third term in (2) describes the change in FSD due to thermodynamic growth or decay represented by freezing or lateral melting. To calculate the thermodynamic term, it is necessary to determine the rate of change in floe size, fl. In the ice melting season, fl changes because of lateral melting. During the melting season, the ocean mixed layer temperature may increase with a gain of heat (per unit area), Qmix, because of ocean heat flux that enters the mixed layer from the deep ocean and surface heat flux that enters the mixed layer through leads (Perovich et al., 2008; Steele et al., 2010), which is described by

(6)

where Tmix is the mixed layer temperature, Tf is the freezing point, Cw is the volumetric heat capacity of water, and dmix is the mixed layer depth. This heat Qmix is available for both lateral melting and bottom melting such that

$Qbot=χQmix,$
(7a)

and

$Qlat=Qmix−Qbot,$
(7b)

where Qbot and Qlat are heat allocations for bottom melting and lateral melting, respectively, and χ is the fraction of ocean heat used for bottom melting, with the rest used for lateral melting. The value of the bottom heating fraction χ falls between 0 and 1.

The partition of the heat Qmix between lateral melting and bottom melting depends on the ratio of lateral to bottom floe areas (Steele, 1992), as well as the thermal stratification within leads and the upper ocean. As the mean floe size in a given area decreases, the bottom heating fraction χ decreases, while the lateral heating fraction 1 – χ increases. In late spring and especially summer when lateral and bottom melting occurs, the Arctic Ocean is generally in a regime of high stratification and low mixing because of various factors such as upper ocean warming, ice melting, and river runoff in the absence of a storm (e.g., Peralta-Ferriz and Woodgate, 2015). This condition suggests that a majority of the heat gain, Qmix, in the mixed layer is likely to come from surface heat flux that enters the mixed layer through leads (Steele et al., 2010). As a result, water temperatures at floe sidewalls (leads) are often higher than those at floe bottoms. This difference means that more heat should be allocated for lateral melting in summer conditions with generally smaller floes. Based on these considerations, we make the following parameterization for the bottom heating fraction:

$χ=max[(1−gh(0)dh)lm/lmax,0.2]$
(8)

where gh(0)dh is area fraction of open water, and the quantity (1 – gh(0)dh) is ice concentration (in area fraction).

Eq. (8) specifies that the heat allocation for bottom melting (lateral melting) decreases (increases) with decreasing mean floe size, reflecting the changes in the ratio of lateral to bottom floe areas. It also reflects, implicitly, the fact that water temperature in leads is likely higher than that at the ice bottom in summertime when ice floes are generally smaller than in other seasons, and in the MIZ, where ice floes are generally smaller than in the IPI. However, the bottom heating fraction χ is not allowed to drop below 0.2. Thus, when the simulated mean floe size becomes small, χ would be close to the value of 0.3 reported by Steele (1992) with floes of 30 m diameter. The bottom heating fraction χ would approach the value of ice concentration (1 – gh(0)dh) under the conditions of little ice breakup, when the mean floe size approaches the value of the maximum floe size allowed in the model, lmax (= 2502 m; Zhang et al., 2015), which generally occurs in the IPI. Ice concentration in the IPI often falls between 0.8 and 1.0., which means that χ would be close to the value of 0.8 reported by Steele (1992) with floes of 3000 m diameter.

Once the bottom heating fraction χ is parameterized, the heat allocations for bottom melting Qbot and lateral melting Qlat are determined. The term Qbot is then used to adjust the ice growth or melt rate fh in (1) such that$fh′=fh′dQbot/dt,$ where $fh′$is the adjusted ice growth or melt rate. The heat allocation for lateral melting, Qlat, is used to determine the lateral melting source term FL in (1) following Hibler (1980). It is also used to determine the rate of change in floe size fl in (2). To this end, we note that the area s and perimeter p of a floe with a caliper diameter l are given by

$s=λl2,$
(9)

and

$p=µl,$
(10)

respectively, where empirical constants λ = 0.66 and µ = 3.17 (Rothrock and Thorndike, 1984). The rate of change in floe size can be written as

$fl=dl/dt=µwlat/(2λ),$
(11)

where wlat is the sidewall melt rate of the floe such that ds/dt = wlatp (Steele, 1992). For a given FSD category l and ITD category h, there are $gldlλl2$ floes per unit area with ice thickness h. By taking into account (10) and (11), the heat needed to induce a sidewall melt rate of wlat for these floes per unit area is

$wlatpgldlλl2hQI=2flgldllhQI,$

where QI is the volumetric heat of fusion of ice. Given the rate of change in Qlat available for lateral melting of all floes of all thicknesses,

$2flQI∫0∞∫0∞gldllhghdh=−dQlat/dt.$
(12)

In (12), integrations are made over all floe sizes and all ice thicknesses. Based on the assumption of Zhang et al. (2015) that floes of all sizes have the same ITD in a given area and recalling that mean ice thickness is defined by $hm=∫0∞gh(h)hdh$ we can obtain the rate of change in floe size due to melting as

$fl=dQlat/dt2hmQI∫0∞(gl/l)dl.$
(13)

Interestingly, the derived fl does not depend on the empirical constants λ and µ.

In the ice growth season, ice forms not only from existing floes but also in the open water or leads/cracks between the floes. The ice formation that occurs between floes ‘welds’ floes together. Therefore, the rate of change in floe size fl during the growth season may become very large, depending on the rate of ice freezing. However, it is difficult to determine fl because of our lack of knowledge about the welding processes. Consequently, we avoid formulating fl here directly in order to circumvent the uncertainty in determining the rate of welding with floes of varying sizes. Alternatively, we allow ice floes of all sizes in a given area to be welded into the largest floe size category defined in the model, if the ice growth rate in that area exceeds a specified threshold for welding symbolized as wd (section 4). This allowance means that once the ice growth rate is above the welding threshold, we specify that the area fraction of the largest floe size category in the model (centered at lmax) takes on the value of the area fraction of ice concentration, while the area fractions of all the other floe size categories become zero. This simplified treatment is equivalent to making a large fl in a given area with ice growth exceeding the threshold.

Note that thermodynamically driven changes in FSD may also occur in conjunction with the formation of melt ponds (Arntsen et al., 2015). Networks of melt ponds are essentially weak zones that are easier to break apart under wind and current forcing. However, PIOMAS is unable to explicitly simulate the evolution of melt ponds because the model currently does not incorporate melt pond physics and parameterization. Consequently, the melt pond-induced changes in FSD are not explicitly simulated in the model. However, the first-order effects of melt ponds are captured in the model by the seasonal albedo evolution which in turn affects sea ice thickness and hence FSD via (5).

PIOMAS is configured to cover the seas north of 49°N, including the Arctic, Bering, Barents, and GIN (Greenland-Iceland-Norwegian) Seas and Baffin Bay (see Figure 1 in Zhang et al., 2014). Lateral open boundary conditions for PIOMAS along the 49°N southern boundaries are provided by a global ice–ocean model (Zhang, 2005). The modification to allow open boundaries in the POP ocean model (Zhang and Steele, 2007) makes it possible to nest PIOMAS into the global model. The model has 12 categories each for ice thickness, ice enthalpy, and snow depth, representing variable ice thicknesses up to 28 m (see Figure 1 in Zhang et al., 2000). The model also has 12 floe size categories partitioned following a Gaussian distribution to obtain a floe size mesh that varies smoothly in the l space, representing variable ice floe sizes up to 2828 m, with the largest floe size category centered at lmax = 2502 m (see table 1 in Zhang et al., 2015).

Figure 1.
Comparison of simulated cumulative floe numbers with TerraSAR-X data.

Cumulative floe number distribution (CFND) derived using TerraSAR-X images at six locations (a–f) in the Beaufort Sea on different days in 2014 (black lines), with corresponding model results (red lines). Image identification names, locations, and dates are marked.

Figure 1.
Comparison of simulated cumulative floe numbers with TerraSAR-X data.

Cumulative floe number distribution (CFND) derived using TerraSAR-X images at six locations (a–f) in the Beaufort Sea on different days in 2014 (black lines), with corresponding model results (red lines). Image identification names, locations, and dates are marked.

PIOMAS is driven by daily National Centers for Environmental Prediction/National Center for Atmospheric Research (NCEP/NCAR) reanalysis atmospheric forcing including 10-m surface winds, 2-m surface air temperature, and downwelling longwave radiation and cloud fraction (Kalnay et al., 1996). Surface air temperature and cloud fraction are used to calculate downwelling shortwave radiation following Parkinson and Washington (1979). Model spin-up consists of an integration of 30 years using 1948 reanalysis forcing repeatedly, initialized with a constant 2 m ice thickness in the areas of freezing surface air temperature, ocean temperature and salinity climatology (Levitus, 1982), and zero ice and ocean velocity. After this spin-up the model proceeds to simulate the period 1948–2010 without incorporating FSD. After that, the model continues to simulate the period 2011–2014, with FSD incorporated. During model calibration over the period 2011–2014, the empirical constants k, a, and b in (5) and the welding threshold wd are found by varying their values in a series of calibration runs and comparing the results with FSDs derived from Moderate Resolution Imaging Spectroradiometer (MODIS) images and power-law exponents (see section 5a). A best match with observations was obtained for the following values (with the range of values used in calibration provided in parentheses): k = 0.225/86400 (0.1/86400–0.5/86400), where 86400 is the seconds per day, a = 1.5 (1.0–5.0), b = 3.5 (1.0–5.0), and wd = 1.2×10−6 m s−1 (1.0–5.0 ×10−6 m s−1).

Note that PIOMAS is able to assimilate satellite sea ice concentration and sea surface temperature. However, no data assimilation is performed in this study so as to assess the simulated behavior of freely evolving FSD without constraint by observations. The focus is on the seasonal evolution of FSD in the Arctic Ocean, particularly in the MIZ, in 2014 when FSD observations derived from TerraSAR-X images are available for model calibration and validation. Some results for 2013 are used for a comparison with observations derived from MODIS images.

### a. Comparisons with observations

Before examining the seasonal evolution of Arctic sea ice FSD in the multicategory TFED sea ice model, model performance is assessed against MODIS and TerraSAR-X observations to verify that the model results are able to represent realistic natural conditions. MODIS on satellites Terra and Aqua is a multispectral sensor with a nominal resolution of 500 m. TerraSAR-X is an all-weather, day and night microwave radar sensor and its stripmap mode used here has a nominal ground resolution of 3 m. The MODIS data are used for model calibration and validation and the TerraSAR-X data only for model validation.

Model-simulated CFND is compared with available observations of CFND derived from the TerraSAR-X images at six locations in the Beaufort Sea in July and August 2014 (Figure 1). TerraSAR-X is a satellite of the German Aerospace Center (DLR, http://www.dlr.de) carrying an X-band synthetic aperture radar (SAR) instrument capable of all-weather day/night imaging of Earth’s surface in several operating modes and resolutions. The TerraSAR-X observations of CFND are derived using a combination of kernel graph cuts (Salah et al., 2011) and watershed and rule-based post-processing (Hwang and Ren, 2014), and are grouped into the model’s 12 floe size categories for comparison. The model results at these six locations indicate that the Zhang et al. (2015) FSD theory allows PIOMAS to obtain a CFND that is close to a power law, across space and time, as generally observed by satellites and aerial surveys (Rothrock and Thorndike, 1984; Holt and Martin, 2001; Toyota et al., 2006; Steer et al., 2008; Perovich and Jones, 2014).

The model results differ from the TerraSAR-X observations at the low end of the floe size range (< 200 m) because of the flattening of the observation curves (Figure 1). This flattening is due to the resolution limitation of the TerraSAR-X images, which results in a deviation from a power law for small floes. This deviation is common with satellite-derived CFNDs because of the resolution limitation (Holt and Martin, 2001). However, the model does not have this limitation because it includes floes as small as 0.1 m in caliper diameter. As a result, the simulated CFND follows a power law closely at the low end of the floe size range, and therefore differs from the satellite-derived CFND in the flattened zone.

Away from the flattened zone, the TerraSAR-X observed CFND displays a power law behavior (a straight line in the log-log plots), and the simulated CFND agrees reasonably well with the observations. The simulated floe numbers per square km are generally close to the observations in those floe size categories with mean caliper diameters greater than 150 m (Figure 1). On average over the six locations, the cumulative number of floes per square km with caliper diameter greater than 150 m is 2.8 from model results and 3.6 from TerraSAR-X observations (Table 2). While generally following a power law, the TerraSAR-X-derived CFND demonstrates a subtle tendency to deviate from a power law by showing a descent at the high end of the floe size range slightly steeper than a straight line in some of the log-log plots (Figure 1). Such a tendency is more obvious with the simulated CFND. Given that the model has a limited floe size range with the mean floe size in the largest floe size category being lmax = 2502 m, the model’s tendency to deviate from a power law at the high end in summertime is because of decreasing number, and ultimately disappearance, of floes of large sizes (up to the largest floe size category) as ice continues to break and melt. The fall-off from a power law at the high end is reported in other studies and often described by an upper truncated power law (Pickering et al., 1995; Burroughs and Tebbens, 2001; Lu et al., 2008). Here the upper truncated power law is replicated in PIOMAS.

Table 2.
Average cumulative number of floes derived from TerraSAR-X observations and corresponding model results

aControl (CNTL) and sensitivity (SEN) runs (see Table 1)

bAverage over the six locations with TerraSAR-X images with caliper diameter > 150 m

cPercentage model deviation from TerraSAR-X observations defined as (model-observation)/observation x 100

dAlso see Figure 1

As mentioned above, MODIS images were used for model calibration and validation. MODIS is a 36-band instrument on NASA’s Terra and Aqua satellites; we used the visible band images and extracted manually the cloud-free regions. We developed a database of FSD derived from more than 200 MODIS images collected throughout, mainly, the western Arctic during April through October in 2013 and 2014 (Stern et al., 2014). The MODIS-derived FSD is over a floe size range different from the model, with little overlap, because the MODIS images have relatively coarse resolution (pixel size 250 m). Thus it is not possible to compare PIOMAS-simulated CFND with MODIS-derived CFND directly, as is done using the TerraSAR-X data (Figure 1). Instead, we take advantage of the fact that the PIOMAS-simulated or MODIS-derived CFNDs generally obey a power law (also see Stern et al., 2014). The significance of a power law is its scale invariance such that the CFND may be characterized by a single exponent α over all floe sizes. This feature means that the power-law exponent may be a useful quantity for a model–observation comparison.

A comparison between model-simulated power-law exponents and those derived from MODIS images is given through a scatter plot (Figure 2a) and a map of point-by-point differences between model results and observations at the locations of observations (Figure 2b). There are 233 data points derived from the same number of MODIS images collected largely in the western Arctic during the period April–October in 2013 and 2014 (Stern et al., 2014; also see Figure 2b for locations). The MODIS-derived power-law exponents as well as model results generally fall within the range of 0.5–2.5 (Table 3). The comparison shows an overall mean model bias of –0.07 (or –5% against an observed mean value of 1.29), with a model–observation correlation of R = 0.57 (Figure 2a). This comparison suggests that the model captures about 32% of the variance of the observations. The model has a tendency to underestimate observed high exponents (cases with a higher proportion of small floes) and overestimate observed low exponents (cases with a lower proportion of small floes). This effect is because of the model’s tendency to create smoother fields than reality, as is also the case with the simulated ice thickness (see Schweiger et al., 2011, and Figure 3 below).

Figure 2.
Comparison of simulated power-law exponents of the cumulative floe number distributions with MODIS data.

(a) Model-simulated power-law exponents of the cumulative floe number distributions (CFND) and those derived using MODIS images collected from April through October in 2013 and 2014. The blue line indicates equality and the red line represents the best fit to the observations. The number of total observation points, observed and modeled mean values, model bias (mean model–observation difference), and model–observation correlation (R) are listed. (b) Point-by-point differences in power-law exponents between model results and observations at the locations of the observations; red (blue) color indicates that the modeled exponent is higher (lower) than the observed with more (fewer) floes of smaller size.

Figure 2.
Comparison of simulated power-law exponents of the cumulative floe number distributions with MODIS data.

(a) Model-simulated power-law exponents of the cumulative floe number distributions (CFND) and those derived using MODIS images collected from April through October in 2013 and 2014. The blue line indicates equality and the red line represents the best fit to the observations. The number of total observation points, observed and modeled mean values, model bias (mean model–observation difference), and model–observation correlation (R) are listed. (b) Point-by-point differences in power-law exponents between model results and observations at the locations of the observations; red (blue) color indicates that the modeled exponent is higher (lower) than the observed with more (fewer) floes of smaller size.

Table 3.
Mean power-law exponent derived from MODIS images and corresponding model results, with model–observation deviation and correlation

aControl (CNTL) and sensitivity (SEN) runs (see Table 1)

bDerived from more than 200 MODIS images (see Figure 2)

cPercentage model deviation from MODIS-derived mean value, defined as (model—observation)/observation x 100

Figure 3.
Comparison of simulated sea ice thickness with IceBridge data.

(a) Model-simulated sea ice thickness (m) and thickness data from the NASA IceBridge program (Kurtz et al., 2013) collected during March and April 2014 and available from the Unified Sea Ice Thickness Climate Data Record (Lindsay, 2010, 2013). The blue line indicates equality and the red line represents the best fit to the observations. The number of total observation points, observed and modeled mean values, model bias (mean model–observation difference), and model–observation correlation (R) are listed. (b) Point-by-point differences in ice thickness (m) between model results and observations at the locations of the observations; red color indicates that the modeled ice is thicker than the observed.

Figure 3.
Comparison of simulated sea ice thickness with IceBridge data.

(a) Model-simulated sea ice thickness (m) and thickness data from the NASA IceBridge program (Kurtz et al., 2013) collected during March and April 2014 and available from the Unified Sea Ice Thickness Climate Data Record (Lindsay, 2010, 2013). The blue line indicates equality and the red line represents the best fit to the observations. The number of total observation points, observed and modeled mean values, model bias (mean model–observation difference), and model–observation correlation (R) are listed. (b) Point-by-point differences in ice thickness (m) between model results and observations at the locations of the observations; red color indicates that the modeled ice is thicker than the observed.

PIOMAS is further evaluated using available sea ice (mean) thickness observations from the NASA IceBridge program (Kurtz et al., 2013) collected during March and April 2014 (Figure 3). The IceBridge ice thickness data are obtained from the Unified Sea Ice Thickness Climate Data Record in which sea ice thickness measurements from various sources are grouped into approximate 50-km aggregates in a unified, easy-to-use format (Lindsay, 2010, 2013) (http://psc.apl.uw.edu/sea_ice_cdr). The comparison between the available IceBridge ice thickness observations in 2014 (287 data points or aggregates in total) and the corresponding model results shows a low mean model bias of 0.05 m (or 2% against an observed mean value of 2.94 m), with a high correlation of R = 0.80 (Figure 3a). Thus PIOMAS captures about 64% of the variance of the observations taken along the IceBridge flight paths (Figure 3b). Note, however, that the accuracy of the IceBridge ice thickness data may be affected by the incorrect identification of side lobes of the impulse response of the snow radar, which impacts snow depth measurements used to estimate ice freeboard and hence ice thickness (Kwok and Haas, 2015). The induced uncertainty in the IceBridge ice thickness data may affect the model–observation comparison.

### b. Seasonal evolution of ice concentration and mean ice thickness

Before examining PIOMAS simulated seasonal evolution of FSD, it is useful, for the purpose of analysis, to examine the simulated seasonal evolution of ice concentration (1 – gh(0)dh) and (mean) ice thickness $(hm=∫0∞gh(h)hdh)$, two key measures of ITD. The seasonal evolution of the Arctic sea ice cover is obviously characterized by ice retreat from the subarctic seas to the Arctic Ocean in late spring through summer and advance again in fall (Figure 4). Thus the location of the ice edge (defined as 0.15 ice concentration) or the MIZ changes in space and time. For example, a large portion of the MIZ is located in the Bering, Barents, and Greenland Seas during much of the winter (January–March) and spring (April–June) (see the 0.15 and 0.80 ice concentration contours in Figure 5). In late spring and summer of 2014, much of the MIZ is well within the Arctic Ocean before the ice cover advances again in fall (October–December). The model generally overestimates ice extent (defined as the area with ice concentration greater than 0.15), in comparison with satellite observations of ice extent (Figure 4). However it underestimates ice coverage in some locations in summer (Figures 4d–4e). These discrepancies highlight the challenge to pinpoint the ice edge or MIZ locations accurately, even though the model captures much of the variance of the 2014 ice thickness observations along the IceBridge flight paths, away from the MIZ (Figure 3).

Figure 4.
Annual cycle of simulated sea ice concentration.

Simulated monthly averaged ice concentration for January (a), March (b), May (c), July (d), September (e), and November (f) of 2014. The white line indicates the satellite observed sea ice edge, defined as 0.15 ice concentration; the black line, the model-simulated ice edge. Satellite ice concentration data are from http://nsidc.org/data/nsidc-0081.html.

Figure 4.
Annual cycle of simulated sea ice concentration.

Simulated monthly averaged ice concentration for January (a), March (b), May (c), July (d), September (e), and November (f) of 2014. The white line indicates the satellite observed sea ice edge, defined as 0.15 ice concentration; the black line, the model-simulated ice edge. Satellite ice concentration data are from http://nsidc.org/data/nsidc-0081.html.

Figure 5.
Annual cycle of simulated ice thickness.

Simulated monthly averaged mean ice thickness $(hm=∫0∞gh(h)hdh)$ for January (a), March (b), May (c), July (d), September (e), and November (f) of 2014. The white line indicates the simulated ice edge; the black line, the simulated 0.8 ice concentration contour, representing the interior boundary of the marginal ice zone (MIZ).

Figure 5.
Annual cycle of simulated ice thickness.

Simulated monthly averaged mean ice thickness $(hm=∫0∞gh(h)hdh)$ for January (a), March (b), May (c), July (d), September (e), and November (f) of 2014. The white line indicates the simulated ice edge; the black line, the simulated 0.8 ice concentration contour, representing the interior boundary of the marginal ice zone (MIZ).

Figure 5 shows the simulated seasonal changes in ice thickness in ice-covered areas including the MIZ. In winter through mid-spring, ice is relatively thin in the MIZ located in the subarctic seas including the Bering, Barents, and Greenland seas (Figures 5a–5c). In summer, much of the MIZ moves to the Arctic Ocean, and the simulated ice in the MIZ is often relatively thick because most of the thin, first-year ice has melted and the ice that has survived is largely multi-year ice (Figures 5c–5d). In fall, the simulated ice cover expands rapidly, leading to relatively thin ice in a large area, both in and often north of the MIZ (Figure 5f). These patterns are reflected in Figure 6a, which shows a prominent peak in the average ice thickness over the MIZ in late spring and summer. This feature is in contrast to the average ice thickness over all ice-covered areas, which peaks earlier in April–May and decreases in summer, reflecting the annual cycle of ice growth and melt (e.g., Maykut and Untersteiner, 1971).

Figure 6.
Seasonal evolution of simulated ice thickness, floe size and power-law exponent, and observed wind speed.

Simulated 2014 seasonal evolution of mean ice thickness (a), mean floe size (b), and power-law exponent (c), along with NCEP/NCAR reanalysis surface wind speed (d), averaged over all ice-covered areas (with ice concentration greater than 0.15, solid line) and over the MIZ (with ice concentration between 0.15 and 0.8, dotted line) of the whole model domain covering the seas north of 49°N.

Figure 6.
Seasonal evolution of simulated ice thickness, floe size and power-law exponent, and observed wind speed.

Simulated 2014 seasonal evolution of mean ice thickness (a), mean floe size (b), and power-law exponent (c), along with NCEP/NCAR reanalysis surface wind speed (d), averaged over all ice-covered areas (with ice concentration greater than 0.15, solid line) and over the MIZ (with ice concentration between 0.15 and 0.8, dotted line) of the whole model domain covering the seas north of 49°N.

### c. Seasonal evolution of FSD

Figure 7 shows the simulated 2014 seasonal evolution of CFND at six different latitudes along the 150°W transect in the Beaufort Sea. These six locations are mostly in the seasonal ice zone: they are in the IPI in winter but mostly in open water or the MIZ in summer (July–September). At these locations, there are many fewer floes in winter than in summer because many of the floes are ‘welded’ together into bigger floes that enter the largest floe size category in freezing conditions. This effect is why the mean floe size averaged over all ice-covered areas is much larger in winter than in summer (Figure 6b). As a result, the simulated CFNDs are relatively flat with low slopes or exponents (Figure 7), and the simulated power-law exponent averaged over all ice-covered areas is much lower in winter than in summer (Figure 6c). Because of an increase in the number of big floes in the largest floe size category and a lack of increase in the number of smaller floes, the simulated CFNDs in winter with compact sea ice often deviate from a power law at the high end of the floe size range by showing a less steep descent than a straight line in a log-log space (Figure 7).

Figure 7.
Simulated cumulative floe number distribution along 150°W.

Simulated 2014 monthly mean CFND at six different latitudes along 150°W in the Beaufort Sea for January (black), March (blue), May (green), July, (green-yellow), September (orange), and November (pink). There is no orange curve in (a)–(c) because ice does not exist at those locations in September 2014.

Figure 7.
Simulated cumulative floe number distribution along 150°W.

Simulated 2014 monthly mean CFND at six different latitudes along 150°W in the Beaufort Sea for January (black), March (blue), May (green), July, (green-yellow), September (orange), and November (pink). There is no orange curve in (a)–(c) because ice does not exist at those locations in September 2014.

In May there is no deviation from a power law, even at the high end of the floe size range; all the CFNDs at the six different locations display a nearly straight line in a log-log space over all the floe size categories, indicating fully developed fields of ice floes across all sizes (Figure 7). As big floes are broken into smaller ones, the number of big floes decreases and the number of small floes increases more substantially. In July, enough of the big floes are broken such that the simulated CFNDs at the six locations all show a steeper descent than a straight line at the high end of the floe size range in a log-log space — an upper truncated power law behavior. By September, no CFNDs exist in the southern part of the 150°W transect because ice is completely melted away there (Figures 7a–7c). In the northern part of the transect where ice remains, floes in some of the largest size categories are completely depleted, which further increases the number of small floes (Figures 7d–7f). This effect is why the mean floe size averaged over all ice-covered areas is lower in summer than in winter (Figure 6b), while the power-law exponent is higher in summer than in winter (Figure 6c). In November, floes are again welded together into larger ones, and the CFNDs resemble those in March (Figure 7). Consequently, the simulated mean floe size (power-law exponent) averaged over all ice-covered areas increases (decreases) again (Figures 6b–6c).

While Figure 7 shows the seasonal evolution of the simulated CFND at several individual locations in the seasonal ice zone that is generally occupied by compact ice in winter and open water or the MIZ in summer, Figure 8 shows the seasonal evolution of CFND averaged over the entire MIZ, which moves back and forth between the subarctic seas and the Arctic Ocean in an annual cycle (Figure 4). Unlike those in the winter IPI of the Beaufort Sea, the simulated CFNDs over the winter MIZ show a large number of small floes (Figure 8a). This large number is because of enhanced ice breakup in the ice edge region in the subarctic seas induced by stronger winds (Figure 6d), and hence likely stronger waves, acting on thinner, weaker first-year ice (Figure 6a). As a result, the simulated mean floe size in the winter MIZ is not necessarily larger than that in the summer MIZ (Figure 6b), and the simulated power-law exponent is even greater than during some of the spring and summer months (Figure 6c).

Figure 8.
Simulated cumulative floe number distribution in the MIZ.

Simulated 2014 monthly mean CFND averaged over the MIZ for (a) January (red) and March (blue), (b) May (red) and July (blue), and (c) September (red) and November (blue).

Figure 8.
Simulated cumulative floe number distribution in the MIZ.

Simulated 2014 monthly mean CFND averaged over the MIZ for (a) January (red) and March (blue), (b) May (red) and July (blue), and (c) September (red) and November (blue).

In May through July (Figure 8b), the simulated CFNDs in the MIZ resemble an upper truncated power law, with some of the largest floes absent, broken into smaller ones. However, the number of small floes (< 200 m) in these months (Figure 8b) is actually lower than in the winter months (Figure 8a). This difference is because floes of small sizes or first-year ice are easily melted away in spring and summer in the MIZ. It may also be attributed to reduced ice fragmentation in the MIZ, linked to thicker and hence stronger (Figure 6a) ice after thinner, first-year ice is preferentially melted and to weaker winds during spring and summer (Figure 6d). Thus the number of medium-sized floes (> 200 m and < 1500 m) is generally higher during spring and summer (Figure 8).

In September, the biggest floes in the MIZ continue to break into smaller ones and the floes in the largest floe size category are depleted completely (Figure 8c). The number of small floes increases because of intensifying winds (Figure 6d) and small floes are less likely to be melted away at the end of summer. As a result, the mean floe size in the MIZ drops to a minimum (Figure 6b). In November, the number of larger floes increases again because of welding in freezing conditions (Figure 8c), leading to an increase in mean floe size and a decrease in the power-law exponent (Figure 6).

Figure 9 shows the simulated spatiotemporal evolution of mean floe size in 2014. The spatiotemporal evolution meets the expectation that smaller mean floe sizes are simulated in summer than in other seasons or in the MIZ than in the IPI all year long, with relatively high power-law exponents (Figure 10). In mid-fall through winter, large mean floe sizes are created in much of the IPI (Figures 9a–9b and 9f), with power-law exponents often below 1.0 (Figures 10a–10b and 10f). In mid-spring through summer, ice fragmentation occurs not only in the MIZ, but also in the IPI, leading to a decreased mean floe size (Figures 9c–9e) and increased power-law exponent (Figures 10c–10e) in the IPI. In the MIZ, however, there is no clear annual cycle in the variation of mean floe size (Figure 6b). This effect is because the MIZ changes in space and time (back and forth between the Arctic Ocean and the subarctic seas), and the simulated mean floe size (or power-law exponent) along the MIZ strip varies significantly (Figures 9 and 10), depending on local wind forcing and sea ice conditions.

Figure 9.
Annual cycle of simulated mean floe size.

Simulated monthly averaged mean floe size $(lm=∫0∞gl(l)ldl)$ for 2014. The white line indicates the simulated ice edge; the black line, the simulated 0.8 ice concentration contour, representing the interior boundary of the MIZ.

Figure 9.
Annual cycle of simulated mean floe size.

Simulated monthly averaged mean floe size $(lm=∫0∞gl(l)ldl)$ for 2014. The white line indicates the simulated ice edge; the black line, the simulated 0.8 ice concentration contour, representing the interior boundary of the MIZ.

Figure 10.
Annual cycle of simulated power-law exponent.

Simulated monthly averaged power-law exponent for 2014. The white line indicates the simulated ice edge; the black line, the simulated 0.8 ice concentration contour, representing the interior boundary of the MIZ.

Figure 10.
Annual cycle of simulated power-law exponent.

Simulated monthly averaged power-law exponent for 2014. The white line indicates the simulated ice edge; the black line, the simulated 0.8 ice concentration contour, representing the interior boundary of the MIZ.

### d. Model sensitivity

The results discussed up to this point are from a model run that is considered a ‘control’ simulation (denoted hereafter as CNTL). In addition, three sensitivity simulations were conducted in parallel to the CNTL run over the period 2011–2014. These sensitivity runs aimed at examining the effect of incorporating FSD, and at model sensitivity to the parameterizations of the participation factor cb and the ice growth rate threshold for welding wd (Table 1). The first sensitivity run (SEN1) does not incorporate FSD, which is used to assess how model results differ with and without simulating FSD (CNTL vs. SEN1). Like the CNTL run, the other two sensitivity runs (SEN2 and SEN3) incorporate FSD. Model sensitivity to cb described in (5) is represented by two different sets of constants a and b (CNTL vs. SEN2) (Table 1). Note that the value of cb is also dependent on the constant k in (5). However, given that the effect of changing k is more or less similar to that of changing a and b, the value of k is fixed at 0.225/86400 for all the simulations incorporating FSD, limiting the number of sensitivity simulations. In addition, model sensitivity to the welding threshold is also represented by two different values (CNTL vs. SEN3) (Table 1).

Incorporating FSD (CNTL) tends to reduce ice thickness all year round when compared to the run without (SEN1) (Figure 11). The reduction in ice thickness is generally greater in the MIZ than in the IPI. This result is because the MIZ has more small floes than the IPI, which enhances lateral melting, increases open water, and thus enhances ice–albedo feedback. The reduction is more noticeable in September (Figure 11e) when the mean floe size in the MIZ drops to a minimum (Figure 6b). Because of the generally low ice thickness in the MIZ (Figure 5), the percentage of ice thickness reduction there can be as high as 50% (Figure 12). The reduction in ice thickness in the MIZ throughout the year tends to retard ice expansion in fall and winter and accelerate ice retreat in late spring and summer.

Figure 11.
Annual cycle of simulated ice thickness difference.

Simulated 2014 monthly mean ice thickness difference between the CNTL and SEN1 model runs (see Table 1). The black line indicates the simulated 0.8 ice concentration contour for the CNTL run.

Figure 11.
Annual cycle of simulated ice thickness difference.

Simulated 2014 monthly mean ice thickness difference between the CNTL and SEN1 model runs (see Table 1). The black line indicates the simulated 0.8 ice concentration contour for the CNTL run.

Figure 12.
Annual cycle of percentage of simulated ice thickness difference.

Simulated 2014 monthly mean ice thickness difference in percentage between the CNTL and SEN1 model runs (see Table 1), defined as (CNTL – SEN1)/SEN1×100. The black line indicates the simulated 0.8 ice concentration contour for the CNTL run.

Figure 12.
Annual cycle of percentage of simulated ice thickness difference.

Simulated 2014 monthly mean ice thickness difference in percentage between the CNTL and SEN1 model runs (see Table 1), defined as (CNTL – SEN1)/SEN1×100. The black line indicates the simulated 0.8 ice concentration contour for the CNTL run.

Note that the model generally tends to overestimate ice coverage compared to satellite observations (Figure 4). This tendency is also reflected in the seasonal evolution of the whole ice extent (Figure 13a). The key point here is that incorporation of FSD has the tendency to improve the simulation of ice extent and ice edge locations, especially during summer and early fall (Figure 13b). The percentage improvement of the CNTL simulated ice extent over the SEN1 run is small during January–July and November–December, but is more noticeable during August–October, with a 10% improvement in September. However, in some locations where the model underestimates ice coverage, the incorporation of FSD is likely to increase the underestimation.

Figure 13.
Seasonal evolution of ice extent.

Satellite-observed 2014 seasonal evolution of ice extent (area with ice concentration greater than 0.15) and simulated extent in model (CNTL and SEN1) runs (a) and the percentage deviation of simulated from observed ice extent, defined as (model – observation)/observation×100), and between the SEN1 and CNTL runs (b).

Figure 13.
Seasonal evolution of ice extent.

Satellite-observed 2014 seasonal evolution of ice extent (area with ice concentration greater than 0.15) and simulated extent in model (CNTL and SEN1) runs (a) and the percentage deviation of simulated from observed ice extent, defined as (model – observation)/observation×100), and between the SEN1 and CNTL runs (b).

The SEN2 run has larger a and b and therefore smaller cb than the CNTL run under the same external forcing (Table 1). This means that the magnitude of ice breakup in SEN2 is generally smaller than CNTL. As a result, the SEN2 simulated CFNDs averaged over the MIZ show more large floes and fewer small floes than the CNTL results (Figure 14). The percentage deviations of SEN2 from the TerraSAR-X-derived CFNDs and MODIS-derived power-law exponents are noticeably larger than those of CNTL (Tables 2 and 3), indicating that the magnitude of ice breakup in SEN2 is less adequate than that in CNTL. This model result also suggests that the simulated CFND is sensitive to changes in cb. The SEN3 run, on the other hand, has a higher welding threshold than the CNTL run, indicating a slower pace in welding floes together in freezing conditions. As a result, the SEN3 simulated CFNDs averaged over the MIZ show fewer large floes and more small floes than the CNTL results in all seasons except summer (Figure 14). In summer, the effect of the welding threshold diminishes, as expected, and the SEN3 simulated CFNDs differ little from the CNTL results (Figures 14d–14e). The percentage deviations of SEN3 from the TerraSAR-X-derived CFNDs and MODIS-derived power-law exponents are somewhat smaller than those of CNTL, but its correlation with the MODIS power-law exponents is lower than the CNTL (Tables 2 and 3).

Figure 14.
Cumulative floe number distribution from sensitivity runs.

Simulated 2014 monthly mean CFND averaged over the MIZ for January (a), March (b), May (c), July (d), September (e) and November (f). Red, blue, and green lines indicate results from the CNTL, SEN2, and SEN3 model runs, respectively.

Figure 14.
Cumulative floe number distribution from sensitivity runs.

Simulated 2014 monthly mean CFND averaged over the MIZ for January (a), March (b), May (c), July (d), September (e) and November (f). Red, blue, and green lines indicate results from the CNTL, SEN2, and SEN3 model runs, respectively.

Changes in the participation factor and welding threshold impact the simulation of ice thickness, but the impact is relatively small (Figure 15). Against CNTL, SEN2 tends to slightly increase ice thickness, while SEN3 tends to slightly decrease ice thickness, all year long (Figure 15 is an example for September). The magnitude of the ice thickness difference is small when compared to that of ice thickness difference between CNTL and SEN1, even in the MIZ (Figure 11). This result suggests that adjusting the participation factor and welding threshold in a way that does not cause substantial CFND deviations from satellite observations does not substantially change the simulation of ITD.

Figure 15.
September ice thickness difference.

Simulated September 2014 mean ice thickness difference between the CNTL and SEN2 model runs (a) and the CNTL and SEN3 runs (b). The black line indicates the simulated 0.8 ice concentration contour for the CNTL run.

Figure 15.
September ice thickness difference.

Simulated September 2014 mean ice thickness difference between the CNTL and SEN2 model runs (a) and the CNTL and SEN3 runs (b). The black line indicates the simulated 0.8 ice concentration contour for the CNTL run.

A sea ice model of thickness, floe size and enthalpy distribution (TFED) has been implemented into PIOMAS to explicitly simulate, simultaneously, the evolution of Arctic floe size distribution (FSD) and ice thickness distribution (ITD). This effort was accomplished by coupling the FSD theory of Zhang et al. (2015) with the ITD theory of Thorndike et al. (1975) augmented by the ice enthalpy distribution theory of Zhang and Rothrock (2001) in the PIOMAS framework. Before model calibration and validation using remotely sensed observations, the implementation of the TFED model in PIOMAS involved three additional model developments essential to close the FSD equation and to couple it to the ITD equation. One was the parameterization of floe size redistribution participation factor based on the observations that wave conditions depend on wind speed and fetch as well as, through wave–ice interactions, ice conditions such as ITD and FSD that also control the flexural strength of sea ice in the MIZ. These observations helped to establish an empirical relationship linking the magnitude of ice fragmentation to winds, which drive sea ice and upper ocean circulation, open water fraction (or fetch in the MIZ), FSD, and ITD. Another was the parameterization in partitioning the heat in the ocean mixed layer for bottom and lateral melting based on FSD, and the determination of the rate of change in floe size based on the heat allocation for lateral melting. The third was to introduce a mechanism to allow ice floes to aggregate into larger ones under freezing conditions.

The results of these modeling developments indicate that simulated cumulative floe number distributions (CFNDs) generally follow a power law that holds across space and time, as observed by satellites and aerial surveys. In particular, they generally agree with the CFND observations over medium and large floe sizes (> 200 m) derived from TerraSAR-X images in the Beaufort Sea in summer 2014. Over small floe sizes (< 200 m), the TerraSAR-X-derived CFNDs are flattened due to the resolution limitation of the satellite images. However, because the model floe size categories are partitioned to include floes as small as 0.1 m, the model is able to simulate the dynamically and thermodynamically important small floes, following a power law closely in the low end of the floe size range. The simulated power-law exponents are also correlated (R = 0.57) with those derived using numerous MODIS images from, mainly, the western Arctic, with a low mean bias of –2%. Note that a power-law exponent comparison is a higher level comparison than a CFND comparison because the exponent is the absolute value of the slope of a power-law obeying CFND and therefore more sensitive to model or data uncertainty. Nevertheless, the model manages to capture 32% of the variance of the MODIS-derived exponents. The reason the simulated CFND generally matches the observations is because the parametrization reasonably represents the role of wind forcing in ice breakup during periods of either high or low levels of wave activity. In addition, the simulated (mean) ice thicknesses are correlated (R = 0.80) with the NASA IceBridge observations in 2014, with a mean bias of 2%.

In the marginal ice zone (MIZ), the simulated ice is, surprisingly, thinner in fall and winter than in summer. This model result is attributed to rapid expansion (melting) of first-year ice in fall and winter (summer). The simulated CFNDs show a large number of small floes in winter because of stronger winds and, implicitly, waves acting on thin, weak first-year ice in the ice edge region, causing large-scale ice fragmentation even though there is a tendency for floes to aggregate in freezing conditions. In mid-spring and summer, the ice is thicker in the MIZ because of the disappearance of thin, first-year ice and the simulated CFNDs resemble an upper truncated power law, with the largest floes mostly broken into smaller ones, a process of mechanical floe size redistribution. However, the number of small floes is lower than in winter because floes of small sizes or first-year ice are melted away easily.

In the ice pack interior (IPI), there are fewer floes in late fall and winter than in summer because many of the floes are welded together into larger floes in freezing conditions, leading to relatively flat CFNDs with low power-law exponents. The increase in the number of large floes and stable number of smaller floes result in a CFND that deviates from a power law at the high end of the floe size range by showing a less steep descent than a straight line in a log-log space, which is the opposite of an upper truncated power law that has a steeper descent, which is often seen in the simulated CFNDs in summer.

The simulated mean floe size fields show smaller mean floe sizes in the MIZ than in the IPI all year long, as expected. In mid-fall through winter, much of the IPI has large mean floe sizes because of floe aggregation. In mid-spring through summer, ice breakup also occurs in the IPI, although at a generally lower magnitude than in the MIZ, leading to a decrease in mean floe size. As a result, the mean floe size averaged over all ice-covered areas shows a clear annual cycle. However, there is no obvious annual cycle with the simulated mean floe size averaged over the MIZ, although the averaged mean floe size falls to a minimum in September. The mean floe size fields, together with ice concentration and thickness fields, provide a better description of the state of sea ice, which may be useful for the planning and management of economic, transportation, and subsistence activities in the Arctic.

The simulated CFND is found to be sensitive all year to changes in the floe size redistribution participation factor, but not sensitive in summer to changes in the welding threshold. The simulated ice thickness is found to be insensitive to moderate changes in the participation factor and welding threshold. Increasing the participation factor or the welding threshold tends to generally increase the number of small floes and therefore reduce ice thickness because of enhanced lateral melting, and vice versa. However, the changes in ice thickness are relatively small. When compared to a sensitivity simulation in which FSD is not incorporated, the control simulation that incorporates FSD leads to reduced ice thickness, mainly in the MIZ. The reduction in ice thickness in the MIZ can be as high as 50% in summer, which generally improves the model simulation of ice extent and yields an earlier ice retreat.

Although this study is the only attempt, to our knowledge, to explicitly simulate the evolution of FSD in a dynamic thermodynamic sea ice model on a pan-Arctic scale, much needs to be done to further improve model physics and parameterization. For example, it is necessary to explicitly simulate the evolution of melt ponds as they are found to play a role in changes in FSD (Arntsen et al., 2015). Additional improvement is possible as the ability to explicitly simulate FSD as a prognostic variable opens a door for incorporating additional model physics in the future. For example, model physics may be improved by parameterizing the influences of FSD on the mechanical properties of the ice and therefore its response to winds and ocean waves and currents (e.g., Shen et al., 1987; Feltham, 2005). FSD, together with ITD, may also be useful for improving the surface exchange of heat, mass, and momentum. For example, floe edges and pressure ridges in the Arctic all obstruct the flow of air or water past the ice and therefore are a source of form drag (Tsamados et al., 2014), which may be parameterized based on the spatiotemporal evolution of FSD and ITD. The empirical parameterization of the floe size redistribution participation factor should be considered to be an interim solution until the effects of wind, current, and wave forcing on ice fragmentation can be modeled explicitly. For example, the empirical parameterization provides a general framework for more direct incorporation of the effects of waves, wave–ice interactions, and wave-induced ice fragmentation, which are likely to be simulated explicitly in the next generation sea ice models coupled to wave models.

Finally, we want to point out that the properties of the FSD for smaller floes (< 200 m) are not well observed because of the resolution limitation of many satellite images (e.g., the MODIS and TerraSAR-X data analyzed here). There are published results that suggest a ‘regime’ change at some small floe sizes where the power law exponent may take different values (Toyota et al., 2006), as the small floes may only become smaller primarily by melt rather than by flexure. To shed more light on changes in FSD for small floes, high-resolution observations and image analyses are sorely needed, in addition to modeling efforts. High-resolution observations of FSD for small floes will also be of considerable value for model development and validation.

Model results are available by contacting JZ (zhang@apl.washington.edu).

© 2016 Zhang et al. This is an open-access article distributed under the terms of the Creative Commons Attribution License, which permits unrestricted use, distribution, and reproduction in any medium, provided the original author and source are credited.

We gratefully acknowledge the support of the Office of Naval Research and the NASA Cryosphere Program. We thank two anonymous reviewers for their constructive comments.

Arntsen
AE
,
Song
AJ
,
Perovich
DK
,
Richter-Menge
JA
.
2015
.
Observations of the summer breakup of an Arctic sea ice cover
.
Geophys Res Lett
42
: doi: .
Burroughs
SM
,
Tebbens
SF
.
2001
.
Upper-truncated power laws in nature systems
.
Pure Appl Geophys
158
:
741
757
.
Comiso
JC
.
2012
.
Large decadal decline of the Arctic multiyear ice cover
.
J Clim
25
:
1176
1193
.
Dumont
D
,
Kohout
A
,
Bertino
L
.
2011
.
A wave-based model for the marginal ice zone including a floe breaking parameterization
.
J Geophys Res
116
:
C04001
. doi: .
Feltham
DL
.
2005
.
Granular flow in the marginal ice zone
.
Philos T Roy Soc A
363
:
1677
1700
.
Flato
GM
,
Hibler
WD III
.
1995
.
Ridging and strength in modeling the thickness distribution of Arctic sea ice
.
J Geophys Res
100
:
18611
18626
.
Herman
A
.
2010
.
Sea-ice floe-size distribution in the context of spontaneous scaling emergence in stochastic systems
.
Physical Rev E
81
: 066123.
Hibler
WD III
.
1980
.
Modeling a variable thickness sea ice cover
.
Mon Wea Rev
1
:
1943
1973
.
Hibler
WD III
.
2001
.
Sea ice fracturing on the larger scale
.
Eng Fract Mech
68
:
2013
2043
.
Holt
B
,
Martin
S
.
2001
.
The effect of a storm on the 1992 summer sea ice cover of the Beaufort, Chukchi, and East Siberian seas
.
J Geophys Res
106
(
C1
):
1017
1032
.
Horvat
C
,
Tziperman
E
.
2015
.
A prognostic model of the sea ice floe size and thickness distribution
.
The Cryosphere Discuss
9
:
2955
2997
. doi: .
Hwang
B
,
Ren
J
.
2014
.
Monitoring spring-to-summer sea ice floe breakup using high resolution SAR and autonomous observing platforms
.
AGU Fall Meeting, San Francisco, 15–19 December
.
Kalnay
E
,
Kanamitsu
M
,
Kistler
R
,
Collins
W
,
Deaven
D
, et al
1996
.
The NCEP/NCAR 40-year reanalysis project
.
B Am Meteorol Soc
77
:
437
471
.
Kohout
AL
,
Meylan
MH
.
2008
.
An elastic plate model for wave attenuation and ice floe breaking in the marginal ice zone
.
J Geophys Res
113
: C09016. doi: .
Kohout
AL
,
Williams
MJM
,
Dean
S
,
Meylan
MH
.
2014
.
Storm-induced sea ice breakup and the implications for ice extent
.
Nature
509
:
604
607
.
Kurtz
NT
,
Farrell
SL
,
Studinger
M
,
Galin
N
,
Harbeck
JP
, et al
2013
.
Sea ice thickness, freeboard, and snow depth products from Operation IceBridge airborne data
.
The Cryosphere
7
:
1035
1056
. doi: .
Kwok
R
,
Cunningham
GF
.
2015
.
Variability of Arctic sea ice thickness and volume from CryoSat-2
.
Philos T Roy Soc A
373
: 20140157. doi: .
Kwok
R
,
Haas
C
.
2015
.
Effects of radar sidelobes on snow depth retrievals from Operation IceBridge
.
J Glaciol
61
(
227
): doi: .
Kwok
R
,
Hunke
EC
,
Maslowski
W
,
Menemenlis
D
,
Zhang
J
.
2008
.
Variability of sea ice simulations assessed with RGPS kinematics
.
J Geophys Res
113
: C11012, doi: .
Langhorne
PJ
,
Squire
VA
,
Fox
C
,
TG
.
1998
.
Break-up of sea ice by ocean waves
.
Cold Reg Sci Tech
27
:
438
442
.
Levitus
S
.
1982
.
Climatological atlas of the world ocean
.
Prof. Pap. 13: 173 pp
.
Silver Spring, MD
:
NOAA
.
Lindsay
RW
.
2010
.
Unified Sea Ice Thickness Climate Data Record
.
University of Washington
:
Polar Science Center, Applied Physics Laboratory
.
digital media
. http://psc.apl.washington.edu/sea_ice_cdr,.
Lindsay
RW
.
2013
.
Unified Sea Ice Thickness Climate Data Record Collection Spanning 1947–2012
.
:
National Snow and Ice Data Center
. http://dx.doi.org/10.7265/N5D50JXV.
Lindsay
RW
,
Schweiger
A
.
2015
.
Arctic sea ice loss determined using subsurface, aircraft, and satellite observations
.
Cryosphere
9
:
269
283
. doi: .
Lu
P
,
Li
ZJ
,
Zhang
ZH
,
Dong
XL
.
2008
.
Aerial observations of floe size distribution in the marginal ice zone of summer Prydz Bay
.
J Geophys Res
113
: C02011. doi: .
Maykut
GA
,
Untersteiner
N
.
1971
.
Some results from a time-dependent thermodynamic model of sea ice
.
J Geophys Res
76
:
1550
1575
.
Meier
WN
,
Hovelsrud
GK
,
van Oort
BEH
,
Key
JR
,
Kovacs
KM
, et al
2014
.
Arctic sea ice in transformation: A review of recent observed changes and impacts on biology and human activity
.
Rev Geophys
51
: doi: .
Meylan
MH
,
Bennetts
LG
,
Kohout
AL
.
2014
.
In situ measurements and analysis of ocean waves in the Antarctic marginal ice zone
.
Geophys Res Lett
41
:
5046
5051
.
Meylan
MH
,
Squire
VA
.
1994
.
The response of ice floes to ocean waves
.
J Geophys Res
99
:
891
900
.
Parkinson
CL
,
Washington
WM
.
1979
.
A large scale numerical model of sea ice
.
J Geophys Res
84
:
311
337
.
Peralta-Ferriz
C
,
Woodgate
RA
.
2015
.
Seasonal and interannual variability of pan-Arctic surface mixed layer properties from 1979 to 2012 from hydrographic data, and the dominance of stratification for multiyear mixed layer depth shoaling
.
Prog Oceanogr
134
:
19
53
.
Perovich
DK
,
Jones
KF
.
2014
.
The seasonal evolution of sea ice floe size distribution
.
J Geophys Res Oceans
119
:
8767
8777
. doi: .
Perovich
DK
,
Richter-Menge
JA
,
Jones
KF
,
Light
B
.
2008
.
Sunlight, water, and ice: Extreme Arctic sea ice melt during the summer of 2007
.
Geophys Res Lett
35
: L11501. doi: .
Pickering
G
,
Bull
JM
,
Sanderson
DJ
.
1995
.
Sampling power-law distribution
.
Tectonophysics
248
:
1
20
.
Rothrock
DA
.
1975
.
Energetics of plastic-deformation of pack ice by ridging
.
J Geophys Res
80
(
33
):
4514
4519
.
Rothrock
DA
,
Thorndike
AS
.
1984
.
Measuring the sea ice floe size distribution
.
J Geophys Res
89
:
6477
6486
.
Salah
MB
,
Mitiche
A
,
Ayed
IB
.
2011
.
Multiregion image segmentation by parametric kernel graph cuts
.
IEEE Image Processing
20
(
2
):
545
557
. doi: .
Schweiger
A
,
Lindsay
RW
,
Zhang
J
,
Steele
M
,
Stern
H
, et al
2011
.
Uncertainty in modeled Arctic sea ice volume
.
J Geophys Res
116
: C00D06. doi: .
Shen
HH
,
Hibler
WD III
,
Leppäranta
M
.
1987
.
The role of floe collisions in sea ice rheology
.
J Geophys Res
92
:
7085
7096
.
Smith
RD
,
Dukowicz
JK
,
Malone
RC
.
1992
.
Parallel Ocean General-Circulation Modeling
.
Physica D
60
(
1–4
):
38
61
.
Squire
VA
.
2007
.
Of ocean waves and sea-ice revisited
.
Cold Reg Sci Technol
49
(
2
):
110
133
.
Squire
VA
,
Dugan
JP
,
P
,
Rottier
PJ
,
Liu
AK
.
1995
.
Of ocean waves and sea-ice
.
Ann Rev Fluid Mech
27
:
115
168
.
Squire
VA
,
Vaughan
GL
,
Bennetts
LG
.
2009
.
Ocean surface wave evolvement in the Arctic Basin
.
Geophys Res Lett
36
: L22502. doi: .
Steele
M
.
1992
.
Sea ice melting and floe geometry in a simple ice-ocean model
.
J Geophys Res
97
:
17729
17738
.
Steele
M
,
Zhang
J
,
Ermold
W
.
2010
.
Mechanisms of summertime upper Arctic Ocean warming and the effect on sea ice melt
.
J Geophys Res
115
: C11004. doi: .
Steer
A
,
Worby
A
,
Heil
P
.
2008
.
Observed changes in sea-ice floe size distribution during early summer in the western Weddell Sea
.
Deep-Sea Res Pt II
55
:
933
942
.
Stern
H
,
Schweiger
A
,
Stark
M
,
Zhang
J
,
Steele
M
, et al
2014
.
The floe size distribution in the marginal ice zone of the Beaufort and Chukchi Seas
.
C11A-0349 - December
.
San Francisco
:
AGU
.
Strong
C
,
Rigor
IG
.
2013
.
Arctic marginal ice zone trending wider in summer and narrower in winter
.
Geophys Res Lett
40
:
4864
4868
. doi: .
Thomson
J
,
Rogers
WE
.
2014
.
Swell and sea in the emerging Arctic Ocean
.
Geophys Res Lett
41
:
3136
3140
. doi: .
Thorndike
AS
,
Rothrock
DA
,
Maykut
GA
,
Colony
R
.
1975
.
The thickness distribution of sea ice
.
J Geophys Res
80
:
4501
4513
.
Toyota
T
,
Takatsuji
S
,
Nakayama
M
.
2006
.
Characteristics of sea ice floe size distribution in the seasonal ice zone
.
Geophys Res Lett
33
(
2
): L02616.
M
,
Feltham
D
,
Schroeder
D
,
Flocco
D
,
Farrell
S
, et al
2014
.
Impact of variable atmospheric and oceanic form drag on simulations of Arctic sea ice
.
J Phys Oceanogr
44
:
1329
1353
. doi: .
P
.
1981
.
Sea-ice topography of the Arctic Ocean in the region 70°W to 25°E
.
Philos T Roy Soc
,
302
:
45
85
.
P
.
1986
. The seasonal ice zone, in
Untersteiner
N
, ed.,
The Geophysics of Sea Ice
.
New York
:
Plenum
: pp.
825
991
.
P
,
Squire
VA
,
Goodman
DJ
,
Cowan
AM
,
Moore
SC
.
1988
.
The attenuation rates of ocean waves in the marginal ice zone
.
J Geophys Res
93
:
6799
6818
.
Zhang
J
.
2005
.
Warming of the arctic ice-ocean system is faster than the global average since the 1960s
.
Geophys Res Lett
32
: L19602. doi: .
Zhang
J
,
Ashjian
C
,
Campbell
R
,
Hill
V
,
Spitz
YH
, et al
2014
.
The great 2012 Arctic Ocean summer cyclone enhanced biological productivity on the shelves
.
J Geophys Res
119
:
1
16
. doi: .
Zhang
J
,
Hibler
WD III
.
1997
.
On an efficient numerical method for modeling sea ice dynamics
.
J Geophys Res
102
:
8691
8702
.
Zhang
J
,
Lindsay
RW
,
Schweiger
A
,
Rigor
I
.
2012
.
Recent changes in the dynamic properties of declining Arctic sea ice: A model study
.
Geophys Res Lett
39
: L20503. doi: .
Zhang
J
,
Rothrock
DA
.
2001
.
A thickness and enthalpy distribution sea-ice model
.
J Phys Oceanogr
31
:
2986
3001
.
Zhang
J
,
Rothrock
DA
.
2003
.
Modeling global sea ice with a thickness and enthalpy distribution model in generalized curvilinear coordinates
.
Mon Wea Rev
131
(
5
):
681
697
.
Zhang
J
,
Rothrock
DA
.
2005
.
The effect of sea-ice rheology in numerical investigations of climate
.
J Geophys Res
110
:
C08014
. doi: .
Zhang
J
,
Rothrock
DA
,
Steele
M
.
2000
.
Recent changes in Arctic Sea ice: The interplay between ice dynamics and thermodynamics
.
J Clim
13
:
3099
3114
.
Zhang
J
,
Schweiger
A
,
Steele
M
,
Stern
H
.
2015
.
Sea ice floe size distribution in the marginal ice zone: Theory and numerical experiments
.
J Geophys Res Oceans
120
: doi: .
Zhang
J
,
Steele
M
.
2007
.
The effect of vertical mixing on the Atlantic layer circulation in the Arctic Ocean
.
J Geophys Res
112
:
C04S04
. doi: .

Funding was provided by the Office of Naval Research (grants N00014-12-1-0112, N00014-12-1-0359, and N00014-12-1-0448) as part of the Marginal Ice Zone, Department Research Initiative, and the NASA Cryosphere Program (NNX15AG68G).

## Competing Interests

The authors have no competing interests to declare.

This is an open-access article distributed under the terms of the Creative Commons Attribution License, which permits unrestricted use, distribution, and reproduction in any medium, provided the original author and source are credited.