A prognostic model of the sea-ice floe size and thickness distribution

. Sea ice exhibits considerable seasonal and longer-term variations in extent, concentration, thickness and age, and is characterized by a complex and continuously changing distribution of ﬂoe sizes and thicknesses, particularly in the marginal ice zone (MIZ). Models of sea ice used in current climate models keep track of its concentration and of the distribution of ice thicknesses, but do not account for the ﬂoe size distribution and its potential effects on air-sea exchange and 5 sea-ice evolution. Accurately capturing sea-ice variability in climate models may require a better understanding and representation of the distribution of ﬂoe sizes and thicknesses. We develop and demonstrate a model for the evolution of the joint sea-ice ﬂoe size and thickness distribution that depends on atmospheric and oceanic forcing ﬁelds. The model accounts for effects due to multiple processes that are active in the MIZ and seasonal ice zones: freezing and melting along the lateral 10 side and base of ﬂoes, mechanical interactions due to ﬂoe collisions (ridging and rafting) and sea-ice fracture due to wave propagation in the MIZ. The model is then examined and demonstrated in a series of idealized test cases.

few decades, Arctic sea ice has become thinner, less extensive, and more seasonal (Cavalieri and Parkinson, 2012). Regions that were once covered by ice year-round now are ice-free in the summer (Stroeve et al., 2012), and the Arctic marginal ice zone, defined as either the region of the ocean over which waves lead to the fracture of ice (e.g. Williams et al., 2013b), or as the area of ice with 25 concentration between 15% and 80%, which has been widening during the summer season (Strong and Rigor, 2013). High-latitude storms are capable of breaking thinning pack ice into smaller floes, changing ocean circulation and air-sea exchange (Asplin et al., 2012;Zhang et al., 2013;Kohout et al., 2015), with evidence suggesting that these storms will become more prevalent in the future (Vavrus et al., 2012). 30 Sea-ice cover is heterogeneous, composed of a distribution of floes of different areas and thicknesses. Floes can vary dramatically in size, ranging from newly-formed frazil crystals millimeters in size to pack ice in the Canadian Arctic with floes up to ten meters thick in places and hundreds of kilometers wide. The most dramatic intra-annual variability in sea ice cover is found in the MIZ, and in seasonal ice zones, regions which range from being ice-covered to ice-free over the year. As 35 summer sea-ice cover becomes thinner and more fractured, these regions will become larger, and the distribution of these floes and their size, shape, and properties may change. Events that generate surface waves, such as a fortuitously observed Arctic cyclone in 2011, the so-called "Great Arctic Cyclone" of 2012, and an energetic wave event observed in the Barents sea, can lead to the fracturing of floes (Asplin et al., 2012;Zhang et al., 2013;Collins et al., 2015). The fractured sea-ice 40 cover has increased floe perimeter, which may lead to enhanced melting and a more rapid reduction in sea-ice area compared to an unfractured sea-ice cover. Steele (1992) indeed demonstrated an increasing sensitivity of the ice cover to lateral melting with decreasing floe size, finding that below 30 m lateral melting was critically important. Smaller floe sizes may additionally lead to changes in the mechanical response of the sea-ice cover to forcing from the ocean and atmosphere, as floe 45 size is a parameter in collisional models of ice rheology (Shen et al., 1986(Shen et al., , 1987Feltham, 2005Feltham, , 2008. As sea ice attenuates wave energy, the diminished ice fraction may lead to further surface wave propagation into the ice field, enhancing fracturing farther from the sea-ice edge, and leading to further sea-ice area loss in a positive feedback loop (Asplin et al., 2014). Floe sizes can also affect the surface drag coefficient and therefore air-sea fluxes (Birnbaum and Lüpkes, 2002). Along floe 50 edges, ocean eddies may be generated due to the gradient in surface heat and stress boundary conditions between ice edge and open water (Niebauer, 1982;Johannessen et al., 1987). These eddies may more rapidly mix air-sea heat flux absorbed by open water to underneath sea-ice floes when floe sizes are comparable to the eddy length scale, but not when floe sizes are much larger. This in turn may have consequences for ice melt rates and ocean circulation (Horvat and Tziperman, 2014). 55 Given that it is not computationally practical to simulate all individual floes, properties of the ice cover can instead be described using statistical distributions. This approach was pioneered by Thorndike et al. (1975), who developed a framework for simulating the thickness distribution (ITD), g(h), defined such that g(h)dh is the fractional area of the sea surface covered by ice with thickness between h and h + dh. The Thorndike model evolves the prognostic equation where u is the horizontal ice velocity, G h is the rate of change of ice thickness due to melting and freezing (thermodynamics), and ψ, the "redistribution function", describes the creation of ice of thickness h by mechanical combination of ice of different thicknesses. Measurements of ice thickness are made possible by a variety of remote sensing techniques such as submarine sonar, fixed 65 moorings, helicopter borne electromagnetic induction, and satellite measurements (Bourke and Garrett, 1987;Yu and Rothrock, 1996;Renner and Gerland, 2014), which may be used to test model skill. Variants of the Thorndike model have been implemented in several general circulation models (GCMs, Bitz, 2008;Hunke et al., 2013), and have been used to understand sea ice behavior and predictability (Bitz et al., 2001;Chevallier and Salas-Mélia, 2012).

70
Modern approaches to modeling sea ice in GCMs, such as the community ice model (Hunke et al., 2013), generally approximate ice cover as a non-Newtonian fluid with a vertically layered thermodynamics, and simple thickness distribution (Thorndike et al., 1975;Semtner, 1976;Hibler, 1979). This approximation may not suffice, because it does not account for the distribution of floe sizes and therefore for the above mentioned related effects. 75 We aim to describe the sub-grid scale variability of the sea-ice cover by extending the ice thickness distribution to a joint distribution that includes both ice thickness and floe size. Rothrock and Thorndike (1984) were among the first to describe the distribution of lateral floe sizes, defining the floe size distribution (FSD) n(r) dr as the fractional area of the sea surface covered by floes with lateral size between r and r + dr. The size of a floe with area a is represented by its effective ra-80 dius, r = a/π, which represents floes as cylinders of radius r. Modeling of the lateral floe size distribution is hampered by the difficulty of measurement, as floe sizes vary over many orders of magnitude. Even with sufficient imagery, algorithms that identify and measure floes must overcome many obstacles, such as submerged floes, melt ponds, and clouds. In spite of these challenges, many observations of the floe size distribution have been made, often using helicopter or ship-board cam-85 eras, notably in the Alaskan and Russian Arctic (Holt and Martin, 2001), Sea of Okhotsk (Toyota and Enomoto, 2002;Toyota et al., 2006), Prydz Bay (Lu et al., 2008), and Weddell Sea (Herman, 2010;Toyota et al., 2011). These studies have focused on deriving and fitting scaling relationships measured distributions, leading to power-law (Toyota et al., 2006), Pareto (Herman, 2010), or joined power-law (Toyota et al., 2011) distributions of floe sizes. The temporal evolution of the floe size 90 distribution has been examined in a small number of observational studies (Holt and Martin, 2001;Steer et al., 2008;Perovich and Jones, 2014), that analyzed the change in the floe size distribution over several weeks or seasonally, but these observations, particularly in the marginal ice zone, are limited.
Herman (2010) modeled the FSD as a generalized Lotka-Volterra system, which admits as a solu-95 tion a Pareto distribution of floe sizes, and suggested that this distribution might fit observed FSDs. Toyota et al. (2011) showed that observed FSDs in the Weddell Sea may be fit by a power law and, that such a scaling relationship may be obtained by assuming that ice fracture is a self-similar process, following a renormalization group method. Zhang et al. (2015) developed a model for the floe size distribution evolution, assuming that all floes of different sizes have the same ITD. The present 100 paper, however, develops a model for the joint floe size and thickness distribution, allowing for different ice thickness distribution for each horizontal size class. The Zhang et al. (2015) paper shares many of our goals and we refer to it below, further elaborating on additional differences between the two studies in the treatment of thermodynamics, mechanical interactions and wave fracturing. Other modeling studies involving the temporal evolution of the floe size distribution have mainly focused 105 on understanding ocean wave propagation and attenuation in the marginal ice zone (Dumont et al., 2011;Williams et al., 2013a, b). These studies developed models of ocean wave propagation, attenuation and associated ice breakage, and modeled the FSD using the renormalization group method of Toyota et al. (2011).
The purpose of the present paper is to develop and demonstrate a framework for modeling the joint 110 distribution of floe sizes and thicknesses (referred to below as the FSTD) f (r, h), with f (r, h) dr dh being the fraction of the ocean surface area covered by floes of thickness between h and h + dh and lateral size between r and r + dr (a list of variable names and descriptions are provided in Table 1).
The ice thickness distribution g(h) and floe size distribution n(r) are obtained by integrating over The prognostic equation for the joint floe size and thickness distribution has the form, where r = (r, h), and ∇ = ( ∂ ∂x , ∂ ∂y ) is the two-dimensional Laplacian. The two dimensional spatial 120 domain may be thought of as corresponding to a single grid cell of a climate model, on the order of tens of km on a side. The term ∇ · (f (r) u) describes advection of the floe size distribution by the flow of ice. L T is the time rate of change of the floe size distribution due to thermodynamic effects. L M is the time rate of change due to mechanical interaction (rafting and ridging of floes).
L W is the time rate of change due to floes being fractured by surface ocean waves. We parameterize 125 each of the above processes, forced by grid-scale atmospheric and oceanic forcing fields. The major contributions of this paper are, first, that it presents the first treatment of the joint floe size and thickness distribution. In addition, each of the terms in equation (2) as developed below contains a novel formulation of the corresponding process that is physically based and less heuristic than used in previous studies.

130
The paper proceeds as follows: we first develop explicit representations for the different processes affecting the joint floe size and thickness distribution in response to atmospheric and oceanic forcing in section 2. The model response to individual forcing fields, in the form of air-sea heat fluxes, ice flow that leads to floe collisions, and surface waves, is analyzed in section 3. We conclude in section 4. size and thickness will change, but the total number of floes will not. Suppose that the only source or sink of ice volume is due to freezing and melting of existing floes, which causes them to change their size at a rate we denote as G r and thickness at a rate G h , and we define G ≡ (G r , G h ). Let N be the number distribution, such that N (r) dh dr is the number of floes in the range (h, h + dh), (r, r + dr) (a list of the variables used to describe FSTD thermodynamics is provided in Table 2).

145
The cumulative number distribution is defined as with ∂ 2 ∂r∂h (C) = N (r) = f (r)/πr 2 , and it obeys the conservation equation, since floes with a finite size and thickness r = (r, h) are, by assumption, neither created nor destroyed 150 by thermodynamic growth and melting. Expanding the right hand side and rearranging in the limit as dt → 0 leads to the time rate of change of the cumulative number distribution, where ∇ r = ( ∂ ∂r , ∂ ∂h ) is the vector of partial derivatives in (size, thickness) space. Changes to the cumulative number distribution are due to the transfer of ice to larger or smaller sizes by thermo-155 dynamic growth and melting. We next make the assumption that thickness changes due to melting and freezing do not depend on the floe radius, and that horizontal size changes do not depend on the thickness, i.e., ∂ ∂h (G r ) = ∂ ∂r (G h ) = 0. The time evolution of the floe size distribution solely due to freezing and melting of existing floes is derived by taking derivatives with respect to both thickness and size of (3), Without loss of generality, consider the interpretation of this equation for the case of freezing in which existing floes get thicker and larger. This implies that some of the area f (r) now moves to larger ice classes, represented by the first term in (4). Note that the integral over all size classes and 165 thickness of the first term vanishes, and therefore it does not describe ice area growth. The total ice area added or removed that belongs to floes of size r, N (r)d/dt(πr 2 ), equal to N (r)2πr G r , which is equal to the second term in (4). Zhang et al. (2015) include the effects of melting and freezing on the FSD, in a way that depends on the lateral growth rate (our G r ), but without evaluating this rate in terms of thermodynamic 170 forcing. Their formulation seems to lack the second term on the rhs of (4). The formulation presented here is for the joint FSTD, and therefore depends on both G r and G h . We further evaluate these rates below in terms of air-sea fluxes.
In addition to melting and freezing of existing floes we must also consider the rate of growth of pancake ice,Ȧ p , due to the flocculation of frazil crystals in patches of open water away from existing 175 floes. Pancakes are assumed to be created by freezing at the smallest size and thickness accounted for in the model, with an effective radius r p and thickness h min . The full expression for the rate of change of the floe size and thickness distribution due to thermodynamics, L T , is therefore, The floe size and thickness change rate vector G = (G r , G h ) is determined using the balance of is shown as the blue and white regions in Fig. 1, (see also Parkinson and Washington, 1979). The total lead area, A lead , is approximated as, where φ is the open water fraction, and the above integration is over the entire ranges of effective radius and thickness represented in the model.
The lead region heat flux, Q lead , is further partitioned into a part that leads to basal freezing or parameterization with a quadratic dependence on open water fraction Q l,l ∝ A 2 lead (Parkinson and Washington, 1979), and diffusive and molecular-sublayer parameterizations based on the tempera-205 ture of the surface waters (Steele, 1992;McPhee, 1992). While these parameterizations have been tested in some detail (Harvey, 1990;Steele, 1992), sensitivity analyses in previous studies have fixed (either explicitly or implicitly) the floe size distribution, and the impact of this assumption on the results is unclear. We choose to simply assume that the lead heat flux is mixed uniformly over the exposed surface of a floe, partitioned according to the ratio of ice basal and lateral surface areas, 210 where it contributes to ice growth or melt. The total fractional lateral surface area (that is, the area of the vertical edges of ice floes, per unit ocean area) is where N is the number distribution introduced above, 2πrh is the lateral area of one floe, and 2h/r represents an average over all ice floes, weighted by the floe size and thickness distribution. The 215 above result depends on the model including an explicit joint FSTD, without which this estimate for the lateral area would not be possible to obtain. The total basal ice surface area per unit ocean area is the ice concentration, c. The partitioning of heat flux from the lead region between the ice base and ice edges is therefore, The rate of change of ice thickness can be found using a model of ice thermodynamics, given the above derived open-water air-sea flux contribution Q l,b to the heat budget at the ice base. For example, ignoring ice heat capacity, ice thickness changes due to melting and freezing are related to the net heat flux into the ice from the surface above, Q surf (defined negative upward), and from below (where negative flux means ocean cooling), The rate of change of the lateral floe size is calculated from the corresponding contribution of the air-sea heat flux from the lead region Q l,l , The above equations can now be used to express the thermodynamic floe growth rate vector, G = 230 (G r , G h ).

Mechanical interactions
Wind and ocean currents can drive individual floe collisions, and therefore merge them together.
When one floe overrides another while remaining intact, the interaction is referred to as rafting. If the ice at the point of contact disintegrates into a rubble pile, forming a 'sail' and a 'keel', and the two  Table 3). The 240 integral of f (r) over all floe sizes and thicknesses, including open water, is equal to one. Therefore, ignoring thermodynamic and wave effects, we integrate (2) over a range of floe sizes that includes a vanishingly small interval of sizes around r = (r, h) = 0, The integral of f (r) over all floe sizes and thicknesses, but excluding open water (r = 0), is equal to the ice concentration, c. Integrating (2) as before but now excluding r = 0, The above definition of operator D M /Dt implies that D M (1)/Dt = ∇·u. The subscript M indicates that this operator represents concentration changes due to mechanical interactions only. D M c/Dt is equal to the total sea-ice area which is eliminated due to the collisions of floes per unit time.
Subtracting (8) from (9), This result implies that L M (r) has a δ(r) component due to open water creation in floe collisions, or the integral on the infinitesimally small range near zero size would have vanished. Note that the function δ(r) is the two-dimensional delta function: (9) suggests that there should be another term in L M (r) that, when integrated over all sizes leads to 260 D M c/Dt. This suggests the following form, where L c (r) is yet unspecified except that its integral over all sizes is one, and it is non-singular at The factor L c (r) quantifies the relative fraction of the total concentration lost due to collisions at each floe size. The terms in (10) Defining the deviatoric strain tensor,˙ ij =˙ ij − δ ij ∇ · u/2, equal to the divergence-free part of˙ ij , two relevant invariants may be written as E = ( I , II ) = (∇ · u, 2| −˙ ˙ ˙ | 1/2 ). The first invariant is the flow divergence and the second is calculated from the determinant of the deviatoric strain rate tensor, and is equal to the maximal shear strain rate. Given these definitions, we parameterize the rate of ice area loss due to collisions as, which allows us to write the mechanical interaction term in the FSTD equation as, This formulation is exactly equivalent to that of Thorndike et al. (1975), see appendix for details. In the case of ice flow characterized by pure divergence, E = (∇ · u, 0) and ∇ · u > 0, the mechanical The rearrangement of floe area in response to a unit amount of open water formation, L c (r), is represented using a collision kernel K(r 1 , r 2 ; r). Let K(r 1 , r 2 ; r) dr 1 dr 2 dr be equal to the number 305 of collisions per unit time between floes in the range (r 1 , r 1 + dr 1 ) and floes in the range (r 2 , r 2 + dr 2 ), that form floes in the range (r, r + dr), per unit area of open water formation. In general, the floe number distribution subject to mechanical combination of floes evolves according to

310
where the notation r dr is taken to mean an integral over all floe sizes and thicknesses resolved by the model. The factor of 1/2 prevents double-counting: since K is symmetric with respect to its first two arguments, each interaction pair (r 1 , r 2 ) is counted twice in the integral in (15). This represents the rate of change in the number of floes of size r 3 due to mechanical interactions. In reality, some floe collisions may lead to a rebound and erosion of floe edges rather than to a merging of the floes, 315 yet we do not account for such a process. The first term on the right-hand side of (15) represents the increase in floe number at size r due to collisions between floes of other sizes, and the second term represents the loss in floe number at size r due to combination of floes of size r with other floes.
Equation (15) is a generalization of the Smoluchowski coagulation equation that has been previously used to model the sea-ice thickness distribution (Godlovitch et al., 2011). If we multiply equation 320 (15) by the area of a floe of size r, we obtain the rate of change of the fractional area covered by floes of size r due to mechanical interactions, which is nothing but the definition of L M (r), We already concluded above that away from r = 0 we have L M (r) = L c (r). Therefore the above eqn gives, where ∂N/∂t is taken from (15). We represent the kernel K(r 1 , r 2 , r) as the product of two factors.
The first is the probability of collision via ridging or rafting of two floes of size r 1 and r 2 , termed P coll (r 1 , r 2 ) where the subscript "coll" is either "ridge" or "raft", and the probabilities are to be defined more specifically shortly.

330
The second factor is a delta function, δ(r − R(r 1 , r 2 )), that limits the pairs of collision partners to only those that form a floe of size r = R(r 1 , r 2 ), specified below, and whose area is smaller than the area of the two colliding floes combined. Noting again that the number distribution and area distribution are related through N (r) = πr 2 f (r), we combine (17) and (15) to find, The coefficient L * c is a normalization constant ensuring that the integral over L c (r) is one (11). In the discretized version of equation (18), two floe classes of discrete size r d 1 and r d 2 which combine to form floes of discrete size r d do not necessarily satisfy π(r d 1 ) 2 h d 1 + π(r d 2 ) 2 h d 2 = π(r d ) 2 h d . Ice volume conservation that is independent of the discretization is achieved by determining the newly 340 formed area of the new floes, in each time step, using the constraint that volume must be conserved, where ∆f (r) is the area change at size r in a single timestep due to the mechanical interaction considered here. Thus the total volume lost by floes at size r d 1 and r d 2 (lhs) is equal to the corresponding volume gained at size r d 3 (rhs).

Probability of collision
We choose the functions P coll (r 1 , r 2 ) to be proportional to the probability that two floes of size r 1 and r 2 will overlap if placed randomly in the domain, and they are calculated in a similar manner for both mechanical processes (rafting or ridging). We consider such an overlap as an indication that mechanical interaction has occurred. The area of each floe that may be deformed due to mechanical 350 interactions is restricted to a small region near the edge of the floe, represented in our model by a narrow annulus, which we term a "contact zone", of width δ cz = δ ridge or δ cz = δ raf t at the floe edge, which depends on the floe size and the interaction type; we also term the interiors of floes "cores" (Fig. 1). The area of a single floe of size s is therefore broken down as, πs 2 = A core (s) + A cz (s) = π(s − δ cz ) 2 + π(2δ cz s − δ 2 cz ).

355
The above defined probability of collision between floes of size r 1 and r 2 is proportional to the product of contact zone areas divided by the open ocean area, A, not including the core areas, The above probability that two floes will collide is based on geometric constraints. However, the rate of collisions depends also on the ice strain rate tensor˙ ˙ ˙ as explained above, and this tensor depends 360 on external forcings such as the strength of the prevailing winds and currents (Shen et al., 1987;Herman, 2011Herman, , 2013Bennetts and Williams, 2015), but the determination of that relationship is not a focus of the FSTD model presented here.
Data of the morphology and width distribution of ridges and rafts as a function of the size of the combining ice floes are scarce, though there are indications that rafts can be substantially larger than 365 ridges (Hopkins et al., 1999). We crudely define the width of the contact zone in ridging to be 5 meters, or the size of the smaller of the two combining floes, whichever is smaller, δ ridge (r 1 , r 2 ) = min(5 m, r 1 , r 2 ).
For rafting, we assume a larger portion of the smaller floe may be uplifted, up to 10 meters, δ raf t (r 1 , r 2 ) = min(10 m, r 1 , r 2 ).

370
Both choices lead to larger ridges and rafts as the size of the interacting floes increases. Given observations of these processes one can refine the above choices, to which our model is not overly sensitive. Finally, we assume that ridging occurs for floes thicker than 0.3 m, and rafting occurs when both floes are thinner than 0.3 m, consistent with the study of Parmerter (1975), with a smooth transition between the two regimes implemented by a coefficient γ(h) which tends to one for thick-375 nesses that are prone to rafting and to zero for ridging,

New floe size 380
The ice area lost in an interaction is different for rafting and ridging. In rafting, the entire contact zone is replaced by ice whose thickness is the sum of that of the original floes. In ridging, the contact zone is increased in thickness by a factor of 5, compressing its area by a factor of 1/5 (Parmerter and Coon, 1972). Given that our model assumes each floe has a uniform thickness, we treat floes formed by ridging or rafting to be of uniform thickness, chosen to conserve volume. This choice eliminates 385 the need for keeping track of sea-ice morphology. Observations (Collins et al., 2015;Kohout et al., 2015) have indicated that floes may break up along ridges, in which case equation (18) may be used to provide information about the ridge density. This is a potential future extension of the present work.
Assuming without loss of generality that r 1 ≤ r 2 , the area of the newly formed floes is therefore 390 given by the sum of the areas minus the area lost to either ridging or rafting. We then divide this area by π and take the square root to find the size of the newly formed floes. The thickness of the formed floe is calculated from volume conservation. We therefore have, where V (r) = V ([r, h]) = hπr 2 is the volume of an ice floe.

Swell fracture
Sea surface height variations due to surface ocean waves strain and possibly break sea-ice floes into 400 smaller floes of varying sizes. Since this process does not create or destroy sea-ice area, the response of the FSTD to fracture of sea ice by waves obeys the conservation law, where L W (r) is the time rate of change of floes of size and thickness r = (r, h) due to fracture of ice by surface waves in 2, and the integral is over all sizes and thicknesses (a list of the variables used 405 to describe the response of the FSTD to ice fracture by waves is provided in Table 4). Suppose that an area of floes Ω(r, t) dr with sizes between r and r + dr is fractured per unit time. Let new floes resulting from this process have the floe size distribution F (r, s) ds, equal to the fraction of Ω(r, t) that becomes floes with size between s and s + ds. The rate of change of area of floes of size r due to fracture by ocean surface waves is then, The first term is the loss of fractional area of size r that is fractured per unit time, and the second is the increase in the area occupied by floes of size r due to the fracture of floes of larger sizes. (2008) Fig. 6) to a quadratic function of the period and mean 420 thickness (Fig. 2). Kohout and Meylan (2008) only report an attenuation coefficient for wave periods longer than 6 seconds and thicknesses less than 3 meters (red box in Fig. 2), so we extrapolate to shorter periods and higher thicknesses using this fit when necessary. We convert the attenuation coefficients from a function of wave period to a function of wavelength using the deep-water surface gravity wave dispersion relation λ = gT 2 /2π. We consider for simplicity a one dimensional domain and assume floes flex with the sea surface height field η(x), experiencing a strain = h 2 ∂ 2 η ∂x 2 (Dumont et al., 2011, p. 4). If the maximum strain, which occurs at the trough and crest of a wave, exceeds an empirically defined value crit , 445 the floe will break. For a monochromatic swell wave of wavelength λ, this leads to floes of size λ/2. For a discretization into N λ spectral lines with spacing ∆λ, spectral amplitudes are defined as

Kohout and Meylan
Let the width of the domain to which the FSTD model is applied be D (e.g., the width of a GCM grid cell which borders on open water). A realization of the sea surface height η(x) is generated according to, where x ranges from 0 to D, the random phases φ i are drawn from a uniform distribution between 0 and 2π, and α(λ i ) is the attenuation coefficient for waves of wavelength λ i .
If the strain is calculated locally from η(x) the critical strain is reached almost everywhere for a realistically-generated wave field (see supplement, Fig. S10). Instead, a floe is assumed to fracture 455 when it is strained between three successive local extrema of η, where points are defined to be extrema if they are a local maximum or minimum over a distance of 10 m on both sides, based on the observations of Toyota et al. (2011) who find this to be the order of the smallest floe size affected by wave fracture. For a triplet of successive extrema (max, min, max; or min, max, min) of η, (x * i−1 , x * i , x * i ), the strain felt by the floe at x * i is calculated by a finite difference approximation 460 (see supplement, section S3). When the magnitude of this strain exceeds the critical strain, crit = 3 × 10 −5 , the floe will break. This determines a set of points at which a floe of thickness h will fracture, X * i (h). From this set of points we define the size of the fractured floe as X * i+1 − X * i . We form a histogram R(r, h) of the number of occurrences of each fracture of size r, which is normalized so that rR(r, h s ) dr = D. In this way, R(r, h s )dr is equal to the number of fractures 465 with size between r and r + dr and thickness h s when waves affect a fully ice-covered domain of length D. We assume that a floe of size s will fracture only when X * i+1 − X * i = r < s, and that the number of fractures of size r is either proportional to R(r) (for r < s), or zero (for r >= s). The total length of fractures of size r is thus proportional to rR(r), or zero, for r > s. The floe size distribution formed by the fracture of a floe of size s, F (s, r) is therefore equal to the total length of floes of size 470 r that are formed by this fracturing of a floe of size s, normalized such that ∞ 0 F (s, r)dr = 1, i.e., The upper limit of the normalization integral in the denominator is truncated to s because the integrand vanishes for larger values of r as explained above. The delta function δ(h − h s ) represents the fact that fracture does not change ice thickness, i.e., any floes formed from the fracture of ice with 475 thickness h s will also have thickness h s .
The function Ω(r, t)dr is the fractional area that belongs to floes of size between r and r + dr that is fractured per unit time. It is set equal to the the area fraction covered by floes of size r, f(r), multiplied by the fraction of the domain reached by waves of group velocity c g per unit time, c g /D, multiplied by the probability that floes of size r will fracture by waves. To calculate this probability, 480 we note that r R(r ) is the total length of the domain covered by waves that can break floes into size r . Integrating this over r from zero to a size r we find the total width of the domain covered by waves that can produce floes smaller than r, which is the same as the length of the domain covered with waves that can break floes of size r into smaller sizes. Normalizing by the domain width D, we find the final factor in the expression for Ω, The group velocity is taken to be that of the mean zero-crossing wavelength, c g = λzg 8π . Observations of wave propagation in ice (Collins et al., 2015) have suggested that the propagation speed of fracture in ice may be slower than the group velocity of surface waves. With more data, the above choice for c g may be re-evaluated.

490
The effects of the fracture of ice by waves on the FSD is represented by Zhang et al. (2015) based on an expression similar to (19), assuming that only floes with horizontal size larger than a specified threshold break, that a fractured floe is equally likely to form any smaller size within a specified range, and that all floes in a given size class have the same ITD. In the representation in the present paper of the effects of ice fracture by waves on the joint FSTD, the wave spectrum plays a central 495 role in determining the resulting floe sizes, as well as the propagation distance over which ocean waves are attenuated by the ice field. Information about the specific thickness of individual floe sizes informs the strain rate failure criterion and therefore determines which floes will be fractured.

Model results
To demonstrate and understand the model's response to a variety of forcing scenarios, we first ex-500 amine its response over a single time step in three runs with idealized forcing fields. Each of these scenarios applies one of the following forcing fields: a net surface cooling Q = −100 W m −2 which induces ice growth, a rate of ice flow convergence of ∇ · u = −5 × 10 −9 s −1 which induces floe collisions, and a surface gravity wave field of a single wavelength λ = 56 m and amplitude of 1 m, leading to ice fracture. The model is initialized with a size and thickness distribution composed of 505 two Gaussian peaks (Fig. 3a). The first (referred to as size I below) has a mean size of 90 m and a mean thickness of 0.25 m. Ice at this size and thickness is susceptible to fracture by surface waves and rafting. The second peak (size II) has a mean size of 15 m and a mean thickness 1.5 m. Ice at this size and thickness tends to ridge rather than raft, and is not susceptible to fracture given our specified wave field. This second point is important, as it demonstrates a possible scenario in which knowledge 510 of the ITD and FSD, separately, would not be sufficient to evolve the FSTD, as some floes, independent of their thickness, will not fracture. The initial sea-ice concentration is 75%. The domain width is D = 10 km, and the width of the lead region is set to be r lw = r min = 0.5 m, the smallest floe size resolved in this model. The critical strain amplitude for flexural failure, crit , is set to 3 × 10 −5 in line with other studies (Kohout and Meylan, 2008;Dumont et al., 2011). Williams et al. (2013a) 515 formulated a more complex expression for the critical failure limit, and this was found to have a significant effect on wave fracturing (Williams et al., 2013b). We examine the model sensitivity to some of the main parameters used in these model simulations in the supplement (Sec. S1).
When two floes of size r and s combine due to rafting or riding interactions, they form a new floe with effective radius r > max(r, s). For an arbitrary floe size discretization into size bins, this new 520 size may not lie within a bin representing a size larger than those of the two interacting floes. As a result, interacting floes may accumulate at a single bin size rather than move into bins representing larger sizes. The minimum bin resolution necessary to avoid this problem is set by the interaction of two floes that are the same size r, with r smaller than the ridge width δ ridge . When two such small floes interact via ridging in our model, one of them becomes 5 times thicker and its area is reduced 525 by a factor of 5. They therefore form a floe of size 6/5r. We select a variable discretization, with r n+1 = 6/5r n , with 64 floe sizes between 0.5 and 156 meters. There are 14 thickness categories, 13 of which are equally spaced between 0.1 m to 2.5 m. To conserve volume when thick floes combine or grow due to freezing, the 14th thickness category incorporates all thicknesses greater than 2.5 m. We examine the numerical convergence of the model in the supplement (Sec. S2) finding 530 that increasing this resolution does not significantly alter the numerical results.
The difference between the model state after a single one-hour time step and the model initial conditions is shown in Figs. 3b-d. Cooling leads to growth in both thickness and size (Fig. 3b) with the impact of lateral growth being less visible than the change in thickness. The shift in thickness is seen by the negative tendency (blue shading) for thicknesses smaller than the maximum of the 535 initial distribution, and positive tendency at sizes larger than the initial maximum (red shading).
These tendencies correspond to the shifting of floes from thinner to thicker floes due to the freezing.
The shift in horizontal size is less apparent in the figure, due to the separation of scales between size and thickness: lateral growth rates are comparable to vertical growth rates (1 cm/day), but given that there is more than an order of magnitude difference between the floe size and thickness, the 540 size change corresponds to a smaller relative change than the thicknesses change. The size response would be more apparent for smaller initial floe sizes not included in this idealized model experiment.
Mechanical interactions (Fig. 3c) lead to growth at three distinct clusters of size and thickness. Swell fracture (Fig. 3d) leads to the fracturing of many of the floes of size I, shown as a negative tendency at the eliminated size class. Floes of size II are not affected because they are smaller than twice the wavelength of the specified surface gravity wave field. Since the specified wave field is monochromatic, the area of floes of size I that are broken is shown as a positive tendency at a floe 555 size equal to half of the wavelength of the surface gravity wave, λ/2 = 28 m. Ice thickness does not change when the ice is fractured.
Next, two one-month simulations are performed using the same initial distribution to show the behavior of the model forced by two different fixed strain rate scenarios (Fig. 4). The first (Fig. 4a,b) simulates convergence of fixed magnitude ( I = −10 −7 , II = 0) s −1 , and the second (Fig. 4c,d) 560 simulates shear of fixed magnitude ( I = 0, II = 10 −7 ) s −1 . When there is no convergence, the rate of open water formation due to collisions (13) is 0.5×10 −7 s −1 , equal to the magnitude of the strain rate tensor divided by two, When there is no shear, and only convergence, the amount of open water formation due to collisions 565 is 10 −7 s −1 , equal to the magnitude of the strain rate tensor, In both scenarios the norm of the strain rate tensor is the same, ||E|| = 10 −7 s −1 . In the case of only shear (Fig. 4c,d), ice concentration is diminished by a factor of roughly 18%, corresponding to a 22% increase in mean ice thickness, and with no change in ice volume. In contrast, in the case of 570 convergence only (Fig. 4a,b), ice concentration is diminished by 36%, with a corresponding 56% increase in mean ice thickness, again with no change in ice volume. Thus shear motions lead to collisions and the combinations of floes with one another, but at a reduced rate when compared to convergence of ice flow, for the same strain rate tensor norm. In the case of shear only, the two initial peaks in the FSTD are smeared out over a range of floe sizes and thicknesses (Fig. 4b), with 575 the variety of floe sizes and thicknesses increasing in number over time. Since there is twice as much open water formation in the case of convergence only, and therefore an increased number of mechanical interactions, the distribution of floe sizes and thickness is smeared more rapidly, and over a larger range (Fig. 4c).
Fig . 5 shows the response of the joint floe size and thickness distribution to a single-week ex-580 periment that simulates a seven-day period of ice fracture by surface waves, using a wave spectrum that leads to ice breaking into a broader range of floe sizes. The experiment uses the Bretschneider (Michel, 1968, p. 24) surface wave spectrum as function of period T , S(T ) dT , where H s = 2 m is the significant wave height (the mean wave height of the 1/3 highest surface 585 waves), and T z = 6 s is the mean time interval between zero-crossings of the observed wave record.
We use the surface gravity wave dispersion relation λ = gT 2 /2π to write S(T ) dT as a wavelength spectrum S(λ) dλ. The wavelength bins are spaced to correspond uniquely to floe size bins, and there is a one-to-one relationship between a wave's wavelength and the floe size of new floes formed through fracture of existing floes by that wave. The peak wavelength of the wave spectrum is at 590 T ≈ 6.75 s, corresponding to λ ≈ 70 m. As before, the domain width D is set to 10 kilometers. Large floes (size I) are rapidly fractured, with the fractional area corresponding to these floes decreasing, and the distribution shifts towards smaller sizes (Fig. 5a, gray lines). After one week, the fractional area belonging to floes in the range from 75-125 m decreases from 37% to 0%, with mean floe size decreasing by 67% (Fig. 5b, blue line). As a consequence, the total lateral surface area rises as floes 595 are broken and their lateral sides are exposed, increasing by 63% over the week (Fig. 5b, blue line).

Conclusions
We developed a model that simulates the evolution of the FSTD, using as input large-scale oceanic and atmospheric forcing fields, which may be useful as an extension to sea-ice models presently used in global climate models, in particular in regions with a continuously varying FSTD, such as 600 the marginal ice zone. We included representations of the impact of thermodynamics (melting and freezing), mechanical interactions of rafting and ridging due to floe collisions, and of floe fracture by ocean surface waves, all processes that are active in marginal or seasonal sea-ice zones. We demonstrated the effect of these processes using model runs forced by external forcing fields including air-sea heat flux, ice flows leading to mechanical interactions, and specified surface wave field, and 605 considered the effects of these forcing fields individually and when combined. We demonstrated the effects of mechanical interactions in the presence of both shearing and straining ice flows, separately accounting for ridging and rafting. We studied the effect of surface waves, first for idealized single-wavelength wave fields, and then accounting for a more realistic surface wave spectrum. We examined the response to melting and freezing both along existing floe bases and lateral edges, and 610 in open water, leading to pancake ice formation.
While the present paper focuses on the development of parameterizations needed to represent the FSTD dynamics and to testing the model with individual forcing fields, we hope to next study the consequences of realistic forcing fields on the FSTD and compare model output to the few available observations. Another important future direction is the model development and testing that will 615 allow for implementation of this model into sea-ice models used in GCMs, allowing for realistic ice thermodynamics, constitutive stress-strain relationship, wave model, and ice motions driven by ocean currents and winds. At the same time, an implementation into a GCM would require making the model more efficient by replacing the high resolution we could afford to use here in floe size and thickness by a simplified approach, possibly assuming a functional form of the FSTD and simulating 620 only its moments as is often done in atmospheric models of the particle size distribution.
The study of FSTD dynamics, and the development of a prognostic FSTD model, are made difficult by the scarcity of observations of the floe size distribution and its seasonal and long term evolution. Such observations are required to constrain uncertain parameters used in the model developed here, and help determine the dominant processes which need to be included in FSTD models 625 to be incorporated in global climate models.
Appendix A: Comparison of rate constants in Eq. 14 to those in Thorndike et al. (1975) Thorndike et al. (1975 employed the following parameterization of the function ψ (1), which represents the rate of change of area belonging to ice of thickness h due to mechanical interactions: where ∞ 0 w r (h) = −1, and the coefficients α 0 and α c are, where θ = arctan( II / I ). Using the trigonometric identity, with ||E|| ≡ 2 I + 2 II , ψ may be rewritten as, Identifying w r = − h L c (r) dh, and 1 2 (||E|| − I ) = D M c Dt , recovers the floe-size-integrated form of (14).   Figure 2. The natural logarithm of the attenuation coefficient α calculated by Kohout and Meylan (2008) (dash, inside the red box) and a quadratic fit to this attenuation coefficient that is used in section 2.3 (solid).

Initial Distribution
Floe Thickness (m)    Tz Zero-crossing period for wave record 3 Table 4. Variables used in the representation of the fracture of ice by surface waves in the FSTD model