Impact of sea ice floe size distribution on seasonal fragmentation and melt of Arctic sea ice

. Recent years have seen a rapid reduction in the summer Arctic sea ice extent. To both understand this trend and project the future evolution of the summer Arctic sea ice, a better understanding of the physical processes that drive the seasonal loss of sea ice is required. The marginal ice zone, here deﬁned as regions with between 15 % and 80 % sea ice cover, is the region separating pack ice from the open ocean. Accurate modelling of this region is important to understand the dominant mechanisms involved in seasonal sea ice loss. Evolution of the marginal ice zone is determined by complex interactions between the atmosphere, sea ice, ocean, and ocean surface waves. Therefore, this region presents a signif-icant modelling challenge. Sea ice ﬂoes span a range of sizes but sea ice models within climate models assume they adopt a constant size. Floe size inﬂuences the lateral melt rate of sea ice and momentum transfer between atmosphere, sea ice, and ocean, all important processes within the marginal ice zone. In this study, the ﬂoe size distribution is represented as a power law deﬁned by an upper ﬂoe size cut-off, lower ﬂoe size cut-off, and power-law exponent. This distribution is also deﬁned by a new tracer that varies in response to lateral melting, wave-induced break-up, freezing conditions, and advection. This distribution is implemented within a sea ice model coupled to a prognostic ocean mixed-layer model. We present results to show that the use of a power-law ﬂoe size distribution has a spatially and temporally dependent impact on the sea ice, in particular increasing the role of the marginal ice zone in seasonal sea ice loss. This feature is important in correcting existing biases within sea ice models. In addition, we show a much stronger model sensitivity to ﬂoe size distribution parameters than other parameters used to calculate lateral melt, justifying the focus on ﬂoe size distribution in model development. We also ﬁnd that the attenuation rate of waves propagating under the sea ice cover modulates the impact of wave break-up on the ﬂoe size distribution. It is ﬁnally concluded that the model approach presented here is a ﬂexible tool for assessing the importance of a ﬂoe size distribution in the evolution of sea ice and is a useful stepping stone for future development of ﬂoe size modelling.

Abstract. Recent years have seen a rapid reduction in the summer Arctic sea ice extent. To both understand this trend and project the future evolution of the summer Arctic sea ice, a better understanding of the physical processes that drive the seasonal loss of sea ice is required. The marginal ice zone, here defined as regions with between 15 % and 80 % sea ice cover, is the region separating pack ice from the open ocean. Accurate modelling of this region is important to understand the dominant mechanisms involved in seasonal sea ice loss. Evolution of the marginal ice zone is determined by complex interactions between the atmosphere, sea ice, ocean, and ocean surface waves. Therefore, this region presents a significant modelling challenge. Sea ice floes span a range of sizes but sea ice models within climate models assume they adopt a constant size. Floe size influences the lateral melt rate of sea ice and momentum transfer between atmosphere, sea ice, and ocean, all important processes within the marginal ice zone. In this study, the floe size distribution is represented as a power law defined by an upper floe size cut-off, lower floe size cut-off, and power-law exponent. This distribution is also defined by a new tracer that varies in response to lateral melting, wave-induced break-up, freezing conditions, and advection. This distribution is implemented within a sea ice model coupled to a prognostic ocean mixed-layer model. We present results to show that the use of a power-law floe size distribution has a spatially and temporally dependent impact on the sea ice, in particular increasing the role of the marginal ice zone in seasonal sea ice loss. This feature is important in correcting existing biases within sea ice models. In addition, we show a much stronger model sensitivity to floe size distribution parameters than other parameters used to calculate lateral melt, justifying the focus on floe size distribution in model development. We also find that the attenuation rate of waves propagating under the sea ice cover modulates the impact of wave break-up on the floe size distribution. It is finally concluded that the model approach presented here is a flexible tool for assessing the importance of a floe size distribution in the evolution of sea ice and is a useful stepping stone for future development of floe size modelling.
The Arctic is currently in a state of transition Stroeve and Notz, 2018). Multi-year sea ice fraction has decreased by more than 50 % with an increasing proportion of the ice cover now seasonal first-year ice (Kwok, 2018;Maslanik et al., 2007). First-year ice does not have the same surface roughness or the same mechanical or thermophysical (salinity, conductivity, permeability) properties as ice that has developed over multiple years. In particular, first-year ice is thinner and weaker  and hence more vulnerable to fracture in response to external stress . Similarly, the region of the Arctic identified as the marginal ice zone (MIZ), generally defined as the region where ocean waves are able to significantly influence the dynamics of the sea ice (Strong et al., 2017), is projected to increase in extent . An alternative definition of the MIZ, and the one that will be used in the present study, is the region where the concentration of the sea ice extends between 15 % and 80 %. This definition of the MIZ is often more practical for modelling and observational studies where sea ice concentration data are more readily available than information about wave behaviour in sea ice.
Modelling the MIZ is a significant challenge due to its complexity; it is a region in which there is strong coupling between the sea ice, ocean, and atmosphere (Lee et al., 2012;McPhee et al., 1987). The sea ice cover in this region is significantly broken up and fragmented by the waves that define the MIZ (Liu et al., 1992). Wave intensity and storm frequency are projected to increase, which will strengthen wave-sea ice interactions (Casas-Prat et al., 2018;Day and Hodges, 2018). This continues a trend already observed over the past few decades (Stopa et al., 2016). Such interactions are even more prominent around Antarctica due to the dominance of seasonal sea ice in the region (Parkinson and Cavalieri, 2012) and large and increasing wave fetch (Young et al., 2011).
Floe size is a key parameter in describing the evolution of the MIZ (Rothrock and Thorndike, 1984). As sea ice floes become smaller, the available perimeter per unit area of sea ice cover increases, enhancing the lateral melt rate (Steele, 1992). Increased lateral ice melt increases the area of exposed ocean, allowing the input of more heat into the ocean mixed layer from solar insolation. Warming of the upper mixed layer also re-stratifies the ocean. These two processes increase heat available for ice melt through basal and lateral ice melting mechanisms. The former is a well-known mechanism, the albedo feedback (Curry et al., 1995). As the MIZ expands, the lateral ice melting is expected to become an increasingly significant driver of seasonal ice loss.
Currently climate models either assume a fixed and constant characteristic floe size across the Arctic cover, for all types of sea ice (Hunke et al., 2015), or they ignore floe size entirely. This approach does not allow for regional or temporal variations in floe size. Multiple sea ice processes depend on floe size. Lateral melt rate is a function of floe size; the melt rate is proportional to the perimeter per unit area of sea ice. A recent study has found that the basal melt rate may also be influenced by floe size (Horvat and Tziperman, 2018). Floe size can also impact the propagation of waves under the sea ice (Boutin et al., 2018;Meylan and Squire, 1994;Squire, 2007). The assumption of a fixed floe size also prevents sea ice models from accurately representing the impact of processes on the sea ice evolution that act via the perturbation of floe size such as lateral melting and wave-induced fragmentation of floes. Whilst these assumptions are significant, the use of a variable floe size within models will need to be justified against the increased computational cost. The most suitable modelling approach will be context dependent; for example, high-resolution regional sea ice models would be expected to require a higher complexity of floe size treatment than large-scale climate models.
There have been several observational studies aiming to characterise the floe size distribution (FSD) using techniques including satellite imagery and in situ studies (Stern et al., 2018a). FSD data are generally fitted to a power law (Rothrock and Thorndike, 1984). Values have been reported for the magnitude of the exponent of this power law ranging from 1.5 to over 3.5 between different datasets (Stern et al., 2018a). Comparing these observations is complicated by the fact that some studies report a value for the probability distribution of floe size and some for the cumulative floe size distribution. It has been recently pointed out that if a distribution adopts a power law for a probability distribution, it will have a tailing off for larger floes when plotted as a cumulative distribution (Stern et al., 2018a). Furthermore, a recent study (Stern et al., 2018b) found evidence to suggest that the exponent of the power-law FSD evolves throughout the year and is not fixed. This same study was also able to use two satellite datasets with different resolutions but operating over the same region to show floes from as small as 10 m and as large as 30 000 m follow power laws. Other studies find different values for these limits, for example Toyota et al. (2016) showed a power law extending to 1 m (using data collected in situ from a ship), whereas Hwang et al. (2017) found a tailing off from the power law around 300-400 m. As each study operates over a different spatial extent, with a different resolution and different algorithms used to extract the FSD, it is not trivial to identify whether the cut-offs in each scenario are physical or a product of limited resolution or spatial extent. Alternative approaches to a single power law have been proposed including the use of two power laws over different size ranges, with smaller floes found to have a smaller exponent (Steer et al., 2008). The Pareto distribution has also been discussed (Herman, 2010); it is analogous to a power law but with a non-constant exponent. To fully understand and characterise the FSD across the Arctic sea ice, good spatial and temporal coverage is required. Novel techniques, particularly those using autonomous platforms and robotic instruments, are enabling increased high-resolution data capture of sea ice and ocean conditions that can be used The Cryosphere, 14, 403-428, 2020 www.the-cryosphere.net/14/403/2020/ alongside time series of up to 1 m resolution FSD data obtained through remote sensing to better understand the factors driving FSD evolution (Thomson and Lee, 2017). These data could be applied within an approach analogous to that of Perovich and Jones (2014), who used aerial photography alongside simple parameterisations for lateral melting and floe fragmentation by waves, assuming the floe size cumulative distribution adopts a power law, to explore whether these processes could result in the observed changes to the FSD. There are also efforts to characterise the floe size distribution resulting from individual processes, such as laboratory analogues to the wave break-up of ice . Future Arctic expeditions including "Multidisciplinary drifting Observatory for the Study of Arctic Climate" (MOSAiC; Dethloff et al., 2016), planned to last 1 year within the central Arctic, should contribute to the existing FSD datasets.
Modelling studies have used contrasting approaches to represent floes as a distribution. A very simple approach is the use of a semi-empirical relationship between floe size and sea ice concentration (Lüpkes et al., 2012;Tsamados et al., 2015). Although this approach involves a simple amendment to the code and has a negligible computational cost, it is unable to respond to fragmentation processes. It will not capture the desired feedbacks during events such as storms that are expected to produce significant fragmentation of the sea ice cover. Furthermore, the parameters used within the relationship were constrained by a set of observations from a specific region and season and might not be applicable across the whole sea ice extent and full seasonal cycle.
Extending beyond using this simple dependency of floe size on sea ice concentration, Zhang et al. (2015) introduced a thickness, floe size, and enthalpy distribution. This model aims to represent the impacts on floe size of advection, thermodynamic growth, lateral melting, ice ridging, and ice fragmentation. However, the impacts of wind, current, and wave forcing are represented by an empirically parameterised floe size distribution factor.  focus on the incorporation of a physically realistic wave-induced break-up model (Williams et al., 2013a, b).  assume that the FSD follows a split power law, with a change in exponent at some critical diameter. The wave component of this model assumes steady-state conditions over a time step and uses a Bretschneider spectrum defined by a significant wave height and a peak period for computational efficiency and propagates it in the mean wave direction. The propagation directions are calculated from averages of the wave directions entering the neighbouring cells and weighted according to the respective wave energy. The model implementation also assumes floe sizes to be assigned to a minimum representative diameter if ice is too thin and compliant to be broken by waves. A recent study by Boutin et al. (2019) also considers the interactions between floe size and waves within the MIZ. This study includes a fully coupled ocean surfacewave model and is unique in considering the impact of mo-mentum transfer to the sea ice from the waves via the wave radiative stress.
There has also been a significant drive to develop a physically derived prognostic floe size-thickness distribution Tziperman, 2015, 2017;Roach et al., 2018a). A recent approach by Roach et al. (2018a) includes the representation of five processes: new ice formation, welding of floes, lateral growth, lateral melt, and fracture by ocean surface waves. This model has the advantage that it does not involve any assumptions about the form of the distribution. Provided the model incorporates good physical representations of the processes which impact floe size, the model should respond accurately to localised extremes in behaviour (such as the large waves associated with storms) or future changes (e.g. changing wind speeds). It is also possible to model floe evolution at the floe by floe scale, for example Herman (2018) uses a discrete-element model to investigate the wave-induced behaviour of floes.
For this study, a single power law will be applied to describe the FSD within a stand-alone sea ice model coupled to a prognostic mixed-layer model, hereafter referred to as the WIPoFSD model (Waves-in-Ice module and Power law Floe Size Distribution model). The distribution is defined by three parameters: d min , lower floe size cut-off for the distribution; d max , upper floe size cut-off; and α, the power-law exponent. α, d min , and d max are set to fixed values. We also introduce a new floe size tracer, l var , which evolves between fixed limits in response to four key processes: wave-induced break-up, lateral melting, advection, and a restoring mechanism in freezing conditions. The WIPoFSD model has been selected as it is able to respond to processes that influence floe size without the computational expense of a full prognostic FSD model. The model allows an assessment of how a power-law distribution of floes will impact the sea ice cover and by what mechanisms these changes occur. Furthermore, it provides a simple framework to explore the model sensitivity to the three parameters used to define the WIPoFSD. A series of additional experiments are also possible within this framework including imposing a variable exponent, changing the parameters that define the impact of waves on sea ice, and comparing the model sensitivity of the floe size parameters to other parameters that influence the lateral melt rate. A stand-alone sea ice model has been selected over a coupled approach to limit model complexity so that the physical impacts and feedbacks of imposing the WIPoFSD model can be more easily identified and to permit more sensitivity studies. The WIPoFSD model is coupled to a prognostic mixed layer so that mixed-layer feedbacks can also be considered.
In this study we present results to understand the thermodynamic response of the sea ice to a power-law-derived FSD and the individual impacts of wave-floe size and lateral melting-floe size interactions. Our focus will be on the impact of this FSD on the seasonal sea ice retreat and variability rather than on longer-term changes and trends.  (Hunke et al., 2015). This is a dynamic and thermodynamic sea ice model designed for inclusion within a climate model. CICE includes a large choice of different physical parameterisations; see Hunke et al. (2015) for details. Section 2.1 outlines the features pertinent to this study. Our local version also includes some state-of-the-art parameterisations not included within the general CICE distribution, also described in Sect. 2.1. The WIPoFSD model that we have implemented into standalone CICE is adapted from an implementation developed at the National Oceanography Centre of the UK within a coupled sea ice-ocean framework, called the NEMO-CICE-Waves-in-Ice (WIM) model (Hosekova et al., 2015;NERSC, 2016). This approach was originally developed to understand the impact of waves on the MIZ and the upper ocean via the thermodynamic and dynamic response with applications for the operational forecasting of the MIZ and largescale coupled sea ice-ocean global modelling, where assuming a power law is particularly practical. The model includes the wave attenuation and floe break-up model based on the Waves-in-Ice Model from the Nansen Environmental and Remote Sensing Center (NERSC) Norway (Williams et al., 2013a, b). An overview of this scheme is given in Sect. 2.2. Floe size is assumed to follow a single power law within the WIPoFSD model. Three new global parameters and one tracer are required to define this power law. The global parameters are d min , lower floe size cut-off for the distribution; d max , upper floe size cut-off; and α, the power-law exponent. The introduced variable FSD tracer, l var , is a function of several processes that change floe sizes: lateral melting, wave break-up of sea ice, advection, and freeze-up. We also introduce a new floe size metric l eff to characterise the FSD, the effective floe size. Section 2.3 outlines how the imposed FSD is defined and describes amendments made to model thermodynamics to account for the change in floe size treatment. This section also provides a definition of l eff . Further details about the treatment of floe size and how l var evolves are given in Sect. 2.4.

Description of standard model physics
Within the CICE v5.1.2 model we use the incremental remapping advection scheme (Lipscomb and Hunke, 2004), an ice thickness redistribution scheme (Lipscomb et al., 2007), along with five ice thickness categories (Hunke et al., 2015). The default elastic-viscous-plastic (EVP) rheology is used (Hunke and Dukowicz, 2002) along with an ice strength formulation (Rothrock, 1975). The frictional energy dissipation parameter is set to 12. A topologically based melt pond scheme is used (Flocco et al., 2012) in conjunction with a delta-Eddington radiation scheme (Briegleb and Light, 2007). The atmospheric and oceanic neutral drag coefficients are assumed constant in time and space. An ocean heat flux formulation is used at the ice-ocean interface (Maykut and McPhee, 1995).
The rate of thermodynamic ice loss is calculated as follows: where A refers to the sea ice concentration, H to the ice thickness, L to the floe diameter (300 m in the default set up), and α shape a geometrical parameter to represent the deviation of floes from having a circular profile (0.66 in the default setup). The terms w top , w bas , and w lat refer to the melt rate at the floe upper surface (top melt), base (basal melt), and sides (lateral melt). The lateral melt rate is calculated as follows: Here m 1 = 1.6×10 −6 m s −1 K −m 2 and m 2 = 1.36 (Perovich, 1983). T is the elevation of the surface water temperature above freezing. The basal and top melt rates are not explicitly calculated, but instead expressed as changes in height derived from a consideration of fluxes over the top and bottom floe surfaces (Hunke et al., 2015). Both lateral and basal melting are reliant on there being sufficient heat flux from the ocean to the sea ice to produce the predicted melting. The model calculates a melting potential term, F frzmlt , for the upper ocean layer. If F frzmlt < 0 in a grid cell where sea ice is present, lateral and basal melting will occur. F frzmlt is proportional to the difference between the sea surface temperature and sea ice freezing temperature (up to a maximum limit of 1000 W m −2 ). If the total heat flux required to produce the calculated basal and lateral melt exceeds the value permitted by the melting potential, then both values will be reduced proportionally such that the total heat flux required equals F frzmlt . Note that H stays constant with respect to lateral melt; so discarding the w top and w bas terms in Eq. (1) we have an expression for the rate of sea ice concentration loss In these simulations, the default CICE fixed slab ocean mixed layer (ML) is not used, and instead a prognostic mixed-layer model is used wherein the temperature, salinity, and depth of the layer are all able to evolve with time . These variables evolve based on surface fluxes and entrainment-detrainment at the base of the ML. The ML entrainment rate is calculated based on the mechanical energy input by wind forcing and surface buoyancy fluxes and profiles of water properties beneath the mixed layer (Kraus and Turner, 1967). This implementation also includes a minimum ML depth, set to 10 m. The prognostic mixed-layer model used here cannot capture the full extent of ocean variability; however it is sufficient to represent sea ice-mixed-layer feedbacks via the mixed-layer properties. Tsamados et al. (2015) have previously compared the performance of the prognostic ML model used here to observations (Peralta-Ferriz and Woodgate, 2015). The mixed layer was found to be generally realistic, though it shows a bias towards too shallow mixedlayer depths through the melting season. A number of amendments are made to CICE version 5.1.2 based on recent work by . The maximum meltwater added to melt ponds is reduced from 100 % to 50 %. This produces a more realistic distribution of melt ponds (Rösel et al., 2012). Snow erosion, to account for a redistribution of snow based on wind fields, snow density, and surface topography, is parameterised based on Lecomte et al. (2015) with the additional assumptions described by . The "bubbly" conductivity formulation of Pringle et al. (2007) is also included, which results in larger thermal conductivities for cooler ice.

Waves-in-ice module
The full details of this module are described in Williams et al. (2013a, b), to which the reader is referred for details; here we provide an overview of the elements pertinent to our study alongside developments unique to the WIPoFSD model. The waves-in-ice module described here reproduces wave conditions near the sea ice edge within the MIZ. Local wind direction determines the direction of wave propagation with adjustments made for attenuation imposed by the sea ice cover. This is a compromise dictated by availability of forcing data, lack of observational studies, and the coarse resolution of the CICE model.
The module operates using its own internal time step defined by where c is the Courant-Friedrichs-Lewy (CFL) condition, here set to 0.7, x min is the size of the smallest grid cell, and c g,max is the highest available group velocity. This is necessary due to the high wave speeds observed in the Arctic. Over each module time step, the wave field is advected, attenuation of waves is calculated, and any ice-breaking events are identified. Note also the forcing fields within each module time step are interpolated between the prior reading and the subsequent reading to ensure smooth variations in the field (note this only applies if the grid cell remains ice-free over this period). We construct the wave energy spectra using H s , the significant wave height (m), and T p the peak wave period (s). These parameters are obtained from the ERA-Interim reanalysis dataset (Dee et al., 2011). The forcings are updated at 6 h intervals, but only for locations where the sea ice is at less than 1 % coverage, i.e. grid cells where there will be negligible wave-ice interactions. The ocean surface wave spectra, S (m 2 s −1 ), are then constructed using the two-parameter Bretschneider formula, Here ω is the frequency (rad s −1 ). H s and T p are used rather than the full wave energy spectra for consistency with Williams et al. (2013a, b). Once the wave field S is defined, it needs to be advected into the ice-covered regions. In the first instance this involves defining the directional space of advection. A principal direction is defined as that of the boundary surface stress component of the ocean. This is generally close to the atmospheric wind direction; however, sea ice also contributes to the boundary surface stress. The waves are advected in five directions spaced equally around the principal direction, with the total angular size of the surface wave spread equal to 90 • . The energy is distributed amongst the bins according to 2/π(cos θ ) 2 , where θ is the deviation from the principal wave direction. The wave energy spectra are then discretised into 25 individual frequencies from a minimum wave period of 2.5 s and a maximum of 23 s. The wave energy spectra are then advected in each defined direction using an upwind advection scheme with each individual spectrum advected separately using its group velocity c g (ω). This advection process is necessary because the wave forcing, derived from the ERA reanalysis data, does not cover areas with a sea ice cover. Furthermore, due to differences between the modelled sea ice edge and observations, there can exist ice-free regions within the model for which no wave forcing data are available.
The decision to use the ocean surface stress to define the primary direction of wave propagation rather than the Stokes drift direction was made because the Stokes drift direction data were not available within the sea ice field at the time of model development. The use of ocean surface stress will be sensible for wind-driven seas, but not for swell-driven seas where the Stokes drift is a more appropriate choice. Stopa et al. (2016) discuss wave climate in the Arctic bewww.the-cryosphere.net/14/403/2020/ The Cryosphere, 14, 403-428, 2020 tween 1992 and 2014 and they find that regions exposed to the North Atlantic wave climate will be strongly influenced by swells generated within the North Atlantic Ocean. Semienclosed and isolated seas, e.g. Laptev and Kara seas, are more event driven and have an equal mix of wind-driven and swell-driven waves. The results presented in this study should therefore be considered in the context that the direction of wave propagation is a significant approximation. Furthermore, we are only able to represent the impacts of waves generated externally to the sea ice cover within this set-up. The choice of surface wave spread is also non-trivial. Wadhams et al. (2002) showed that a wave propagating into the MIZ could experience significant wave spreading until it was essentially isotropic. However, a distinction was found between wind seas where the isotropic state could be achieved within a few kilometres and swell seas where spreading occurs much more slowly, if at all. Wave spreading has been shown to be dependent on the wavelength. Montiel et al. (2016) found that shorter wavelengths experienced spreading and longer wavelengths did not with a transition between these two regimes defined by the maximum floe size. This is consistent with the observed behaviour of wavedriven regimens and swell-driven regimes. Using a fixed surface wave spread across a limited number of categories is a significant simplification of the rather complex spreading behaviour of waves; however it represents a balance between short wave periods that quickly achieve an isotropic state and longer wave periods that propagate much further into the MIZ before they experience significant spreading. After advection, the attenuation of waves over each wave time step is calculated. This will be calculated for each individual wave energy spectrum: where S at is the wave spectrum after attenuation (m 2 s −1 ), α dim is the dimensional attenuation coefficient (m −1 ), t wav is the module time step (s), and other variables are as previously defined. α dim can also be described as the rate of exponential attenuation per metre. It is here modelled as a sum of the linear wave scattering at floe edges in addition to a viscosity term. It is also updated discontinuously when the wave energy is large enough to cause ice breakage. α dim effectively becomes a function of mean floe size, sea ice concentration, ice thickness, and wave period (see Williams et al., 2013a, for further details). After attenuation, the wave energy spectra within each grid cell are reconstructed as a discretised function of ω by summing the advected spectra from each of the five incident directions. The final spectra, S(ω), can then be advected using the process described above for subsequent time steps. If we assume that the sea surface elevation follows a Gaussian distribution, i.e. non-linear affects that can cause asymmetry are neglected, we can calculate the following properties of interest from the wave energy spectra: the mean square surface elevation of the ocean, η 2 ; the mean square surface elevation of the sea ice, η 2 ice ; the mean square strain for the sea ice (modelled as a thin elastic plate), ε 2 ; and the representative wave period, T W . Each of these metrics requires the computation of integrals over frequency, here approximated using Simpson's rule (see Williams et al., 2013a, for further details). H s , a model output, can be calculated as 4 η 2 (Laing et al., 1998).
The floe fragmentation scheme used is identical to Williams et al. (2013a), which should be referred to for a detailed description of the scheme. An overview of this scheme is presented here. Ice-breaking events occur when the probability that the breaking strain amplitude, E s , exceeds the breaking strain, ε c , becomes larger than a critical probability, P crit : We assume that the spectrum is narrow enough to be considered monochromatic. In this case P crit = e −1 and the criterion reduces to E s > ε c √ 2. E s is defined as 2 ε 2 , i.e. twice the standard deviation in strain. ε c is calculated as σ c /Y * , where σ c is the flexural strength and Y * the effective Young's modulus for the sea ice. σ c and Y * are calculated using empirically derived expressions, where both are dependent on the brine volume fraction.
T W is used to calculate the representative wavelength, λ W , required to update the FSD after a wave fragmentation event (see Sect. 2.4 for details on how the FSD is changed). λ W is calculated as 2π/k W , where k W = k ice (2π/T W ). Here k ice (ω) is the positive real root of the dispersion equation for a section of ice-covered ocean.

Floe size distribution model
We employ a number-weighted FSD, N (x), where x is the floe diameter. N (x) is fitted to a power law as shown in Fig. 1. It is described by the following equation: Here N has units of reciprocal metres, d min and l var have units of metres, and α is unitless. l var is the variable FSD tracer, also in metres. l var evolves independently in each grid cell as a function of physical processes between the upper and lower floe size cut-offs of the distribution, d max and d min respectively. d max also has units of metres. d min , d max , and α can all be defined independently for each grid cell; however in this study they will be fixed across the sea ice cover within an individual simulation. l var can be considered to represent the history of a given area of sea ice in terms of physical processes that affect the FSD. The model is initiated with l var set to d max in all grid cells where sea ice is present. The floe size number distribution factor, C, is determined such that the total area of individual floes, N α shape x 2 , sum to equal the total sea ice area, Al 2 The Cryosphere, 14, 403-428, 2020 www.the-cryosphere.net/14/403/2020/ Within individual grid cells the variable FSD tracer, l var , varies between these two limits. l var evolves through lateral melting, wave break-up events, freezing, and advection. Not shown is how changes in l var will also impact the floe size number distribution factor, C. Panel (b) shows how d min , d max , and α can all be varied to produce floe size number distributions. All axes within both panels are logarithmic to base 10.
(where l 2 is the total grid cell area): It should be noted this treatment of N means that in this model the sea ice cover consists only of floes between the limits of d min and d max in diameter. There are no floes with sizes outside these limits. It is useful here to define an additional floe size parameter, l eff , the effective floe size. l eff is defined as the floe size of a distribution of identical floes that would produce the same lateral melt rate in a given instant to a distribution of nonuniform floes, when under the same conditions with the same total sea ice cover. Equation (3), used to calculate the lateral melt rate, can be adapted for use within the WIPoFSD model: The lateral melt rate of a given area of sea ice is proportional to the total perimeter of that sea ice. It is therefore also useful to introduce a second parameter called perimeter density, ρ P , which is the length of the ice edge per unit area of sea ice cover. l eff is hence the constant floe size which produces the same ρ P as an FSD. First, Eqs. (8) and (9) can be used to give an expression for the total sea ice area, Al 2 : The total ice edge length, P fsd , within a grid cell, can also be expressed in terms of the WIPoFSD parameters: We can then divide the second expression by the first to give ρ fsd P , which is P fsd divided by the total ice area in the grid cell, Al 2 : Whilst perimeter density has not been a standard parameter to report from observations, it can be easily calculated from available FSD data. A similar value has been reported by Perovich (2002), though this was reported per unit area of domain size (i.e. ocean plus sea ice area). We can then also define ρ con P , the perimeter density for a distribution of floes of constant size, using an analogous approach: L corresponds to the constant floe size; hence for the 300 m case we would get a perimeter density of 0.0159 m −1 . Setting the perimeter density expressions for both a constant floe size and power-law FSD to be equal, and noting that this defines L = l eff , we obtain Note that Eqs. (13) and (15) are not valid where α = 2 or 3. For these cases, α is taken to be 2.001 and 3.001 to maintain code simplicity with only a negligible cost to accuracy.

Processes that impact l var
In our model there are four ways in which the floe size distribution can be perturbed: lateral melt, break-up of floes by ocean waves, advection of floes, and restoration due to freezing. Changes in l var impact the entire FSD via the floe number distribution factor, C, which is also a function of l var , as defined in Eq. (9). Note that C is also a function of sea ice concentration and therefore, for processes such as lateral melting, changes in both l var and A will contribute to changes in the floe number distribution. It should be noted here that the WIPoFSD model is not intended to represent the impact of physical processes on the details of the floe size distribution; it is indeed not possible to do so in a framework where a power law is imposed. Instead the impact of the different processes considered here is represented via parameterisations, here expressed in terms of the model variable, l var . As lateral melt involves the loss of ice volume from the sides of floes, it can be expected to reduce floe size. To represent this in the model, we set the reduction in l 2 var from lateral melting to be proportional to the reduction in A, the sea ice concentration, from lateral melting: If we then express A final in terms of A initial and A lm , the reduction in sea ice concentration from lateral melting, we obtain The act of reducing l var alone acts to redistribute sea ice area attributed to floes larger than l var to floes smaller than l var . However, the change in A also independently acts to reduce C, as described above. The combined effect is to decrease the number of floes across the whole distribution. Previous studies, such as that by Horvat and Tziperman (2017), have shown that lateral melting causes stronger deviation from the power law for smaller floes than larger floes. However lateral melting also results in floes smaller than d min that will contribute to an even higher lateral melt relative to the floe size.
Hence the behaviour of this lateral melt scheme compensates between these two expected changes to the distribution. Section 2.4 outlines the conditions necessary to trigger the break-up of floes by waves. If these conditions are fulfilled, l var is updated according to the following expression: where λ W is the representative wavelength, as defined in Sect. 2.2. Here l var can be considered a fragmentation length scale, defining the transition from a regime where floes are broken up by waves to a regime where the number of floes is increasing due to this break-up of larger floes.
There are three processes thought to be the main drivers of floe formation and growth during freezing conditions: lateral growth, welding of floes, and formation of new floes (Roach et al., 2018a). The focus of this study is on the seasonal melt and fragmentation of sea ice rather than the winter evolution; hence a simple floe growth restoration scheme is used. During conditions when the model identifies frazil ice growth, l var is restored to its maximum value according to the following expression: where T rel is a relaxation time which relates to how quickly the ice floes would be expected to grow to cover the entire grid cell area. It is set to 10 d as standard. In grid cells that transition from being ice-free to having a sea ice cover, l var is initiated with its minimum value, i.e. d min . The behaviour of the full floe number distribution depends not only on l var but also on A, the sea ice concentration. During periods of freezing when the sea ice concentration increases significantly, both C and l var will increase in value, leading to increases in the number density across all sizes of floes. This is consistent with a scenario where lots of new floes are being formed. During periods of freezing where the sea ice concentration does not increase significantly (e.g. where the sea ice area fraction is already close to 1), then l var will increase and C will decrease. This represents a shift in the distribution from smaller floes to larger floes. It corresponds physically to a scenario where floe welding is the dominant process driving changes in the FSD. l var is transported using the horizontal remapping scheme with a conservative transport equation, the standard within CICE for ice area tracers (Hunke et al., 2015). An amendment to the usual scheme involves calculating a weighted average of the l var over ice thickness categories after advection and the subsequent mechanical redistribution. This is necessary as the tracer is not defined independently for each thickness category unlike other tracer fields. It is useful here to comment on the choice of advection scheme. Firstly, properties that scale to the root of the sea ice area, such as the floe diameter, cannot be advected as an area tracer. Secondly, it has been shown that normalised or mean properties relating to the FSD also do not advect as an area-conserved property (Horvat and Tziperman, 2017). Here, l var is a parameter assigned to areas of sea ice to represent the prior history of that sea ice area in terms of processes that can affect the FSD. l var is not a parameter attributed to individual floes and it is calculated independently to the FSD and is not a diagnostic property calculated from the distribution. Hence, it is appropriate to treat l var as an ice area tracer.
It is worth commenting here on the limitations of the modelling approach to floe size used in this study. The use of a power-law distribution with a fixed exponent to describe the FSD is a valuable simplification to explore the impact of floe size on the Arctic sea ice. The tracer l var is an internal model The Cryosphere, 14, 403-428, 2020 www.the-cryosphere.net/14/403/2020/ tool used to enable parameterisations of how individual processes impact the FSD within this constrained framework. The parameterisations described in this section are necessarily approximations of how these processes might impact the FSD and should not be considered exact physical descriptions.

Methodology
Our modified version of CICE is run over a pan-Arctic domain with a 1 • tripolar (129 × 104) grid. The surface forcing is derived from the 6-hourly NCEP-2 reanalysis fields (Kanamitsu et al., 2002). The mixed-layer properties are restored over a timescale of 5 d to a monthly climatology reanalysis at 10 m depth taken from the MyOcean global ocean physical reanalysis product (MYO reanalysis;Ferry et al., 2011). This restoring is needed to effectively represent advection within the mixed layer. The deep ocean post detrainment retains the mixed-layer properties; however it is restored over a timescale of 90 d to the winter climatology (herein meaning the mean of 1 January conditions from 1993 to 2010) from the MYO reanalysis. All simulations are spun up between 1 January 1990 and 31 December 2004 using the standard set-up described in Sect. 2.1 with a constant floe size of 300 m (without the WIPoFSD model included). Simulations are initiated on 1 January 2005 using the output of the spin-up and evaluated for 12 years until 31 December 2016. Results are all taken from the period 2007-2016 to allow 2 years for the model to adjust to the addition of the WIPoFSD model. A reference run is also evaluated over this period using the standard setup and a 300 m constant floe size. Figure 2 shows this model simulates the climatological monthly sea ice extent realistically for this period. All further simulations are evaluated over the same time period using the same initial model state, however with the WIPoFSD model imposed. Some simulations have additional modifications made to the model as described.

Results
Results are presented for the pan-Arctic domain with a focus on the melting season. All plots compare the mean behaviour over 10 years from 2007 to 2016 against the reference simulation, referred to as ref, which uses a constant floe size of 300 m. The results for 2005 and 2006 are discarded to allow 2 years for the model to adjust to the imposed FSD. In this study we are trying to understand the impact of the FSD and associated processes on the seasonal sea ice loss. The years 2007-2016 have been selected as the baseline for these simulations as they will capture the current climatology of the Arctic, including the record September minimum sea ice extent observed in 2012.  (Comiso, 1999, marked Observations, green ribbon, large dashes). The ribbon shows, in each case, the region spanned by the mean value plus or minus 2 times the standard deviation for each simulation. This gives a measure of the interannual variability over the 10-year period. Results show the new model performs either comparably to or better than the previous optimum set-up throughout the year. In addition, the mean CICE-ML sea ice extent falls within the interannual variability of the observations between June and December, i.e. most of the melting season, suggesting this reference state is suitable for studies focusing on this period.

General impact of an imposed distribution
The WIPoFSD model introduces new parameters that can be constrained through observations. Stern et al. (2018b) were recently able to show a region of floe sizes could be described by power laws over a size range from 10 to 30 000 m. This is the largest range of floe sizes that a power law has produced a good fit to; hence these are set as the standard values for d min and d max in this study. A collated analysis of observations (Stern et al., 2018a) shows that α can adopt values generally ranging from 1.6 to 3.5 (when the FSD is reported as a probability distribution). A standard exponent value of α = 2.5 is adopted as an intermediate value over this range, noting in addition that this value is consistent with the ranges reported by Stern et al. (2018b). The simulation using these standard FSD parameters, α = 2.5, d min = 10 m, and d max = 30 000 m, will be referred to as stan-fsd (see Table 2). Figure 3 displays the percentage difference in sea ice extent and volume for stan-fsd compared to ref. In addition, it shows the spread of twice the standard deviation of these simwww.the-cryosphere.net/14/403/2020/ The Cryosphere, 14, 403-428, 2020 ulations as a measure of the interannual variability. The impact on the pan-Arctic scale is small, with sea ice extent and volume reductions of up to 1.2 %. The difference in sea ice area reaches a maximum in August whereas the difference in sea ice volume peaks in September. The differences in both extent and volume evolve over an annual cycle, with minimum differences of −0.1 % and −0.2 % observed respectively between December and January for ice area and April and May for volume. The annual cycles correspond with periods of melting and freeze-up and are a product of the nature of the imposed FSD. The interannual variability shows that the impact of the WIPoFSD model with standard parameters varies significantly depending on the year. In some years the difference between the stan-fsd and ref set-ups can be negligible, and in other years it can be up to 2 %. Lateral melt rates are a function of floe size but freeze-up rates are not, and hence model differences only increase during periods of melting and not during periods of freeze-up. The difference in sea ice extent decreases rapidly during the freeze-up conditions; this is a consequence of the fact this lateral freeze-up behaviour is predominantly driven by ocean surface properties, which are strongly coupled to atmospheric conditions in areas of low sea ice extent. In comparison, whilst atmospheric conditions initiate the vertical sea ice growth, this atmosphere-ocean coupling is rapidly lost due to insulation of the warmer ocean from the cooler atmosphere once sea ice extends across the horizontal plane. Hence a residual difference in sea ice thickness and therefore volume propagates throughout the winter season. The difference in sea ice extent shows an additional trough in June. This feature is something also seen consistently within the data for individual years and can most likely be attributed to particular weather patterns that occur during the spring season. Figure 4 shows the absolute difference in the mean cumulative annual melt components between the two simulations. The plot shows lateral, basal, top, and total melt (as defined in Sect. 2.1). A large increase can be seen in the lateral melt; however the change in total melt is negligible. This is because the lateral melt increase is largely compensated by a reduction in basal melt. The top melt also shows a negligible change. Figure 4 also shows the change in basal melt in stan-fsd only accounting for the loss of basal surface area available for melting. To explain how this is calculated, imagine for a given time step the sea ice fraction for that grid cell in the stan-fsd simulation is 0.81 and in the ref simulation it is 0.90. If this physical reduction is the only factor causing changes to the total basal melt, then the basal melt rate per unit grid cell area would also reduce by the same factor of 10 % from ref The Cryosphere, 14, 403-428, 2020 www.the-cryosphere.net/14/403/2020/ to stan-fsd. The reduction in the total basal melt volume can then be calculated for this grid cell accounting only for the reduction in sea ice fraction as the product of 0.1, the basal melt rate per unit grid cell area, and the area of the grid cell. This process can be repeated over every grid cell to obtain the total reduction in basal melt volume accounting only for reduction in sea ice concentration. The agreement (to within 1 standard deviation) between this synthetic reduction in basal melt and the actual reduction in basal melt suggests that the loss of ice area by lateral melt is sufficient to explain most of the basal melt compensation effect. Figure 5 shows the spatial distribution for the predicted reduction in basal melt from stan-fsd to ref, the actual reduction in basal melt, and the difference between the actual reduction and predicted reduction in basal melt. These map plots are presented as monthly averages for March, June, and September averaged over  shows that the predicted basal melt can capture the regional distribution of the changes in basal melt from ref to stan-fsd, not just the area-integrated quantity. Figure 6 explores the spatial distribution in the changes in ice extent and volume for 3 months over the melting season, March, June, and September. Data are shown only for regions where the sea ice cover exceeds 5 % of the total grid cell. These results show the differences increase in magnitude through the melting season. Although the pan-Arctic differences in extent and volume are marginal, Fig. 6 shows distinct regional variations in sea ice area and thickness metrics. Reductions in the sea ice concentration and thickness are seen both within and beyond the MIZ with reductions of up to 0.1 and 50 cm observed respectively in September. Within the pack ice, increases in the sea ice concentration of up to 0.05 and ice thickness of up to 10 cm can be seen. In September the biggest increases in thickness are directed along the North American coast, particularly within the Beaufort Sea. To understand the non-uniform spatial impacts of the FSD, it is useful to look at the behaviour of l eff . Regions with an l eff greater than 300 m will experience less lateral melt than the equivalent location in ref (all other things being equal), whereas locations with an l eff below 300 m will experience more lateral melt. The distribution of l eff is shown in Fig. 6 where in general we see a transition from larger floes to smaller floes moving from the pack ice into the MIZ, with the transition to an l eff of a size less than 300 m observed within the MIZ. Most of the sea ice area must therefore experience less lateral melting compared to ref. This result shows that the increase in lateral melt observed in Fig. 4 is localised to regions where the sea ice concentration is around 50 % or below.

Exploration of the parameter space
It has been previously discussed that the floe size parameters used within the WIPoFSD model are poorly constrained by observations. In this section experiments are performed using different permutations of these parameters to assess For the first study the α is changed from 2.5 to 3.5, previously identified as the most extreme value within a reasonable observed range for the power-law exponent. This simulation will be referred to as (A). Figure 7 is analogous to Fig. 4, comparing the component and total melt evolution for an FSD with an α = 3.5 compared to one with an α = 2.5 (with d min and d max set to standard values). The plot shows an increase in the cumulative lateral melt, as seen before for stan-fsd compared to ref. Now, however, the basal melt is www.the-cryosphere.net/14/403/2020/ The Cryosphere, 14, 403-428, 2020 less effective at compensating the lateral melt, resulting in a significant increase in the total melt. There is also now a non-negligible reduction in the top melt, with the interannual variability showing the increase in total melt and reduction in top melt are consistently produced for each year of the simulations. The difference in cumulative total melt reaches a maximum in August and subsequently decreases slightly. This suggests that increasing the magnitude of α results in an earlier melting season and a correspondingly reduced melt in the late season. The predicted change in basal melt based on the reduced sea ice area is again plotted and is able to account for 90 % of the actual reduction in basal melt. This is in contrast to Fig. 4, where the predicted reduction in basal melt was too high compared to the simulated reduction. The interannual variability shows that this underprediction of the reduction in basal melt is consistent throughout individual years. This implies the presence of additional mechanisms such as albedo and other mixed-layer feedbacks causing nonnegligible changes in the basal melt rate; however reduction in the sea ice concentration remains the leading-order impact. Figure 8 shows difference map plots between the two simulations. The ice area and thickness are reduced across the sea ice cover with reductions of over 5 % and 0.5 m respectively seen in particular locations during September. However, even in March, after the freeze-up period, reductions of 0.1 m or more in sea ice thickness can be seen within the ice pack. The response of sea ice can once again be understood through the behaviour of the l eff . l eff is below 30 m across the entire ice cover throughout all 3 months studied, leading to increased lateral melt rates across the sea ice.
A further 17 sensitivity studies using different permutations of the parameters have been completed. These are formed by varying the three key defining parameters of the FSD shown in Fig. 1 in order to span the range of values reported in observational studies: for α values of 2, 2.5, 3, and 3.5 to span the general range of values reported in observations (Stern et al., 2018a); for d min values of 1, 20, and 50 m are selected. These have been selected to reflect the different behaviours reported in studies, with some showing powerlaw behaviour extending to 1 m (Toyota et al., 2006) and others showing a tailing off at an order of tens of metres (Stern et al., 2018b). A further limitation for d min is the smallest floe size where individual floes can be distinguished i.e. the transition from a floe regime to a brash ice regime. For the upper cut-off, d max , values of 1000, 10 000, 30 000, and 50 000 m are selected, again to represent the distributions reported in different studies. The largest value is taken as 50 km for d max as this serves as an upper limit to what can be resolved within an individual grid cell on a CICE 1 • grid. In addition, this model does not account for processes that are expected to be The Cryosphere, 14, 403-428, 2020 www.the-cryosphere.net/14/403/2020/ Figure 8. As Fig. 6 except now the difference between (A) compared to stan-fsd is given, i.e. the impact of changing α from 2.5 to 3.5 with the other FSD parameters held at standard values. l eff is reported for the simulation with the higher magnitude α. In general, the plots show a reduction in the sea ice concentration and ice thickness across the sea ice cover. This corresponds to the behaviour of the l eff , with the l eff 30 m or below across the sea ice cover.
important for the evolution of floes at the kilometre scale and above, such as wind stresses and melt ponds (Arntsen et al., 2015;Wilchinsky et al., 2010). A total of 14 of the 17 permutations for these sensitivity studies are generated by selecting all the different α-d min permutations (except the two already investigated). Each of these simulations has d max = 30 000 m. The further three simulations vary d max with the α and d min fixed to 2.5 and 10 m respectively. Figure 9 shows the change in mean September sea ice extent and volume relative to ref plotted against mean annual l eff , averaged over the sea ice extent. The impacts range from a small increase in extent and volume to large reductions of −22 % and −55 % respectively, even within the parameter space defined by observations. Furthermore, there is almost a one-to-one mapping between mean l eff and extent and volume reduction. This suggests l eff is a useful diagnostic tool to predict the impact of a given set of floe size parameters. The system varies most in response to the changes in the α, but it is also particularly sensitive to d min . The colour of the marker indicates the value of the α, the shape indicates the value of d min , and the three experiments using standard parameters but different d max (1000, 10 000, and 50 000 m) are indicated by a crossed red square. The parameters are selected to be representative of a parameter space for the WIPoFSD that has been constrained by observations. Model response ranges from small increases in the sea ice extent and volume to reductions of over 20 % and 50 % respectively. The mean l eff is shown to be a good predictor of the response of the sea ice extent and volume.

Sensitivity runs to explore specific model components and additional relevant parameters
A series of sensitivity studies have been performed to explore the behaviour of the WIPoFSD model and understand how it interacts with other model components. Table 1 defines the important parameters considered in this section and Table 2 provides a summary of the sensitivity experiments performed. The first two entries in Table 2 α Power-law exponent within the WIPoFSD model. Standard value fixed at 2.5.
l eff The effective floe size is defined as the floe size of a distribution of identical floes that would produce the same lateral melt rate in a given instant to a distribution of non-uniform floes, when under the same conditions with the same total ice cover. See Eq. (15).
The dimensional attenuation coefficient, as used in Eq. (6).

P crit
The critical probability that must be exceeded for wave-breaking events to occur, as used in Eq. (7).
T rel The floe restoration rate, as used in Eq. (19). Set to 10 as default.
α shape Floe shape parameter to account for the deviation of floes from a perfect circle. Standard value of 0.66 (Rothrock and Thorndike, 1984).
m 1 Melt rate parameter, as used in Eq.
m 2 Melt rate parameter, as used in Eq.
scribed above and a constant floe size of 300 m respectively. Studies (A)-(C) are a selection of the simulations described in Sect. 4.2 to allow a comparison between model sensitivity to the parameters that define the FSD and model sensitivity to other relevant parameters and components within the WIPoFSD model. In the following section a bracketed letter will follow descriptions of sensitivity studies, which correspond to the letter assigned in Table 2. Table 3 reports key metrics for the sensitivity studies described in Table 2, plus a selection of the different sensitivity studies described in Sect. 4.2. For each experiment the September sea ice extent and volume size are reported for both the full sea ice extent and MIZ only (taken as a mean between 2007 and 2016), with the MIZ defined here as regions with between 15 % and 80 % sea ice cover. In addition, the mean cumulative lateral, basal, top, and total melts until September are reported in each case, and the September mean l eff and mean sea ice perimeter per square metre of ocean area are both reported averaged over the MIZ. For each value reported (except for the l eff ) the difference from stanfsd is also stated. Cells highlighted in italic and bold font deviate by 1 and 2 standard deviation(s) respectively from the stan-fsd mean value (the standard deviation is calculated from the set of 10 annual values for each metric).

Imposing a variable exponent on the floe size distribution
The shape of the FSD between its limiting values is defined by α. Recent evidence suggests this may not be constant in time or space (Stern et al., 2018b). We have investigated the impact of this behaviour through the use of two alternative modelling approaches. The first approach imposes a sinusoidal annual cycle on α (D): Here d refers to the current day of the year (for example 45 would refer to 14 February) and d ann is the total number of days in the year (here taken to be 365). This curve was selected as a reasonable fit to the observations of Stern et al. (2018b), though it should be noted that these observations were taken from the Beaufort and Chukchi seas so should not be assumed to be representative of the entire Arctic Ocean. The second sensitivity experiment assumes that α is a function of sea ice concentration, A (E). This is derived from the observation that α increases in magnitude as the melting season advances and in locations of lower sea ice concentration: The Cryosphere, 14, 403-428, 2020 www.the-cryosphere.net/14/403/2020/  (I) Weak ice P crit is reduced by a factor of 10.
(J) Weaker wave attenuation α dim is reduced by a factor of 10.
(K) Reduced floe growth rates T rel is increased from 10 to 365.
(L) Fewer circular floes α shape is reduced from 0.66 to 0.44.
(M) Perfectly circular floes α shape is increased from 0.66 to 0.79. This is the approximate value of this parameter for a perfect circle.
(N) Reduced lateral melt rate The parameters m 1 and m 2 are reduced by 10 % each to 1.44 × 10 −6 m s −1 K −m 2 and 1.22 respectively.

(O) Increased lateral melt rate
The parameters m 1 and m 2 are increased by 10 % each to 1.76× 10 −6 m s −1 K −m 2 and 1.48 respectively.
(P) Shallow mixed layer The minimum mixed-layer depth is reduced from 10 to 7 m.
(Q) Deep mixed layer The minimum mixed-layer depth is increased from 10 to 20 m.
The limits were selected to try and capture the variability of the exponent seen within observations. The results in Table 3 show imposing the time-varying α (D) has a very small impact on the sea ice cover, whereas the spatial-varying α (E) causes a moderate reduction in September ice extent and volume of about 3 % and 5 % respectively. It is worth noting that the mean l eff over the MIZ does not correlate well with the size of the response of the system in these cases compared to simulations with a fixed α, with l eff being much higher than expected given the size of the sea ice extent and volume reduction. The value of the sea ice perimeter averaged over the MIZ is more consistent with the observed changes in sea ice extent and volume, particularly for experiment (E). This shows that it is useful to have mulwww.the-cryosphere.net/14/403/2020/ The Cryosphere, 14, 403-428, 2020 Table 3. A summary of the metrics for each of the sensitivity studies described in Table 2. Metrics are reported for sea ice extent, MIZ extent, total sea ice volume, MIZ volume, mean l eff within the MIZ, mean sea ice perimeter per square metre of ocean area within the MIZ, and cumulative melt top, basal, lateral, and total melt. All metrics are reported for September, except the cumulative melt, which is reported for all months up to and including September and given as an average between 2007 and 2016. Means for l eff and ice perimeter are taken as averages over the MIZ with each grid cell equally weighted. The values within the parentheses give the change from stan-fsd. Cells highlighted in italic and bold font deviate by 1 and 2 standard deviation(s) respectively from the stan-fsd mean value (the standard deviation is calculated from the set of 10 annual values for each metric). www.the-cryosphere.net/14/403/2020/ tiple approaches to collapsing the FSD into a representative value. Whilst map plots of l eff can be very useful for understanding the regional impacts of an FSD, as in Fig. 6, the mean value can be misleading. Figures 10 and 11 show how α and the resultant l eff respectively evolve in these two simulations averaged over both the overall sea ice cover and the MIZ. The region spanned by twice the standard deviation of individual years within the simulation is also shown. Whilst l eff in both regions behaves in corresponding ways for the simulation with a time-varying α (D), experiment (E) shows the mean α and hence l eff within the MIZ are small and approximately constant throughout the year, despite the overall sea ice pack showing strong seasonal variability for these quantities. During the peak melting period between May and August the mean l eff is lower for experiment (D) within the pack ice and experiment (E) within the MIZ. Given the much stronger changes seen for experiment (E) compared to experiment (D) relative to stan-fsd, this supports previous findings that the impact of the WIPoFSD model is primarily dependent on the behaviour of the FSD within the MIZ. (D) shows the strongest interannual variation in l eff between March and May, whereas for (E) it is strongest in the peak melting season between July and August. Figure 11 also includes the annual evolution of l eff for the stan-fsd simulation. Unlike (D) and (E), stan-fsd shows no strong annual oscillation in the l eff across the overall pack ice.

Other parameters affecting the floe size distribution
The two processes currently represented in the model that actively reduce l var are lateral melting and wave-induced fragmentation of floes. Two simulations are undertaken where either waves are no longer able to influence l var (F) or lateral melting is no longer allowed to influence l var (G). An additional three simulations are performed to focus on how waves may be influencing sea ice via reductions in l var : the incident significant wave height at the point of entering the sea ice cover is increased by a factor of 10 (H), the floe breaking strain is reduced by a factor of 10 (I), and the wave attenuation coefficients under the sea ice are reduced by a factor of 10 (J). The results in Table 3 show that the wave-l var interaction is more important than the lateral melt-l var interaction in driving the increase in lateral melt observed by imposing the standard FSD. Study (F), where waves no longer reduce l var , shows a 3 % increase in MIZ volume compared to stanfsd, whereas study (G), where l var does not change as a result of lateral melt, shows an increase in MIZ volume of less than 1 %. For the three simulations performed to explore the behaviour of the wave advection model, i.e. (H), (I), and (J), the strongest response is produced by reducing the wave attenuation rate of the model (J). The weakest response is produced by increasing the ice vulnerability to wave fracture (I). Figure 12 shows difference plots of sea ice concentration and l eff The mean α is taken as the equally weighted average across all grid cells where the sea ice concentration exceeds 5 % (total extent) or is between 15 % and 80 % (MIZ only). The imposed annual oscillation in α is identical for all grid cells for (D); hence the MIZ behaviour has not been plotted as it will be identical to the annual oscillation in α across the total sea ice extent. The ribbon shows, in each case, the region spanned by the mean value plus or minus 2 times the standard deviation for each simulation. Both set-ups show an annual oscillation in the value of α averaged over the total sea ice extent. For experiment (E), no obvious annual trend in the mean value of α can be seen when averaged over the MIZ, though the interannual variation is at a maximum during the peak melting season between July and September. between stan-fsd and (J), where the attenuation rate of waves under sea ice is reduced. The plots show a reduction in the sea ice concentration of around 1 % across the MIZ throughout the year for (J). This can be attributed to the reduction of l eff in the same region by magnitudes of greater than 100 m.
The floe restoration rate is the parameter, T rel , used in Eq. (19). As a standard it is set to 10 d; however this value is not well constrained. This effectively means that l var is restored rapidly during freezing conditions, and hence the FSD is effectively initiated in each melting season with no memory of the previous year. There is not enough evidence available to either validate or invalidate the assumption that the FSD retains no memory of the previous melting or freeze-up season. An experiment (K) has been performed where T rel is increased from 10 to 365 d to explore the impact of interseasonal memory retention within the FSD model. The results in Table 3 show that, whilst this change to the model did reduce the l eff and increase the perimeter density metrics by significant amounts, it did not produce a significant www.the-cryosphere.net/14/403/2020/ The Cryosphere, 14, 403-428, 2020 Figure 11. Annual variation in mean l eff averaged over 2007-2016 for two simulations with variable α. The plots show the evolution of l eff throughout the year for a simulation with a time-dependent α (D, red ribbon) or a sea-ice-concentration-dependent α (E, blue ribbon). Also shown is the behaviour of l eff for a simulation with a fixed α of 2.5 (stan-fsd, grey ribbon). Results are shown for the total sea ice area (solid) and MIZ only (dashed). The mean l eff is taken as the equally weighted average across all grid cells where the sea ice concentration exceeds 5 % (total extent) or is between 15 % and 80 % (MIZ only). The ribbon shows, in each case, the region spanned by the mean value plus or minus 2 times the standard deviation for each simulation. The results show that introducing a variable α produces much larger intra-annual variations in l eff across the overall sea ice extent than with a fixed α. change in either the melt components or sea ice extent and volume.
In Fig. 13 we show the evolution of simulations stan-fsd, (F), and (K) over 2015 averaged over selected grid cells. The year 2015 has been chosen as representative over the 2007-2016 period. There are two subplots: the first gives l eff averaged over grid cells with a sea ice concentration within the MIZ on 31 August 2015, selected as the approximate date of the 2015 minimum sea ice extent in simulations. This set of grid cells is chosen to capture grid cells that are marginal for at least some of the year without also becoming ice-free, which would create an artificial seasonal cycle in l eff . For the second subplot, the same set is further constrained to grid cells with between 15 % and 30 % sea ice concentration on 31 August 2015. Figure 6 shows that significant reductions in l eff are generally seen at the outer edge of the sea ice extent, so further restricting the maximum sea ice concentration in this way will capture this region. The significant reduction of l eff by up to 120 m between (K) and stan-fsd in August and September shows that the wave break-up of floes is a significant component of both the floe size reduction and the subsequent reduction in sea ice concentration seen in Fig. 6 for these locations. The difference between (K) and the maximum possible l eff , of just over 540 m, during the melting season primarily captures the impact of lateral melting on floe size as floe restoration will not be active during this period. We see a reduction of up to 50 m for the more marginal set of grid cells, so whilst not insignificant the impact is a factor of around 3-4 times lower than the wave fragmentation in these regions. This suggests that mechanical break-up of floes is a necessary precondition for the lateral melting feedback on floe size to become significant. This effect will not be as strong for other selections of FSD parameters, particular those where l eff is below 50 m even when l var = d max . For these simulations we expect the much larger increase in lat- The Cryosphere, 14, 403-428, 2020 www.the-cryosphere.net/14/403/2020/ eral melt, as seen in Fig. 7, to produce a stronger lateral melt impact on the FSD. For (K), where l var restoration rates during freezing conditions are reduced, l eff is significantly lower throughout the year including during the melting season. l eff varies between 360 and 480 m for the full MIZ grid cell selection, significantly reduced from the 450-540 m seen for the stan-fsd simulation. We also see a well-defined seasonal cycle, unlike with stan-fsd.

Lateral melt parameters
The first-order impact of introducing a variable floe size is on the lateral melt volume. Equation (1) shows the lateral melt volume is calculated from several parameters beyond just floe diameter, L, including lateral melt rate, w lat , and floe shape, α shape . α shape is currently fixed to a constant value, 0.66. There has been significantly less interest in characterising how the shape of floes varies and to characterise a floe shape distribution, particularly given available evidence suggesting floe size and shape may be uncorrelated parameters (Gherardi and Lagomarsino, 2015). Two sensitivity studies are performed: one with α shape reduced to 0.44 (L), corresponding to 3 : 1 rectangular floes or similar distortions from a perfect circle, and one with α shape increased to 0.79, corresponding to approximately circular floes (M). w lat is a function of two parameters, m 1 and m 2 (see Eq. 2). These parameters have been estimated from observations and hence are subject to uncertainty. Experiments are undertaken with either both m 1 and m 2 reduced by 10 % (N) or both increased by 10 % (O). A reduction in these parameters reduces the lateral melt rate and an increase increases it. Table 3 shows that all four of these sensitivity studies did not produce a large model response in terms of the overall sea ice extent and volume. Reducing the floe shape parameter (L) produced the strongest response in the lateral melt volume, and more generally the model metrics were more sensitive to α shape than the melt coefficients, m 1 and m 2 . The much stronger model sensitivity to the floe size parameters justifies the focus on floe size as the main uncertainty in lateral melt volume calculation.

Minimum mixed-layer depth
The minimum ocean mixed-layer depth is a constant within the prognostic mixed-layer model required to prevent the mixed-layer depths reaching unrealistically small values. As a standard it is set to 10 m. The depth of the mixed layer is important for the strength of mixed-layer feedbacks, with a deeper mixed layer acting to damp any feedbacks via mixedlayer properties. These feedbacks include the albedo feedback mechanism and the negative feedback of increased lateral and basal melts (meltwater perturbs the mixed-layer properties towards less favourable melting conditions). Sensitivity studies are performed with the minimum mixed-layer depth both reduced to 7 m (P) and increased to 20 m (Q).
The challenge with this set of experiments is that, unlike the other sensitivity studies presented here, it acts to influence the evolution of the sea ice both via changes in the lateral melt and via the basal melt and sea ice freeze-up rates, determined by ocean properties. Experiment (P) shows a small increase in the total sea ice extent and volume and (Q) a small decrease; however both result in larger increases in the MIZ extent and volume. In comparison to other sensitivity studies, the changes in the lateral and basal melt are www.the-cryosphere.net/14/403/2020/ The Cryosphere, 14, 403-428, 2020 small, suggesting that mixed-layer feedbacks do not have a significant role in the impacts of the FSD found in stan-fsd compared to ref. It should be noted, however, that the evidence presented here is not enough to rule out the existence of multiple compensating feedback processes.

Discussion
We present here a series of simulations and additional sensitivity studies completed with the newly developed WIPoFSD model to explore the impacts of a variable power-law-derived floe size distribution model on the Arctic sea ice. It is useful to consider the physical mechanisms that drive the simulation results. It was previously noted that the increase in lateral melt observed when imposing the WIPoFSD model was compensated by a loss in basal melt, resulting in a more moderate increase in the total melt. Within the model there are three possible mechanisms causing the limited basal melt. Firstly, the increase in lateral melt will correspond to a reduction in available ice area for basal melting. It is shown in Figs. 4 and 7 that this mechanism is able to explain most of the reduction in basal melt, but the difference remains large enough that further mechanisms need to be considered. The second mechanism concerns the melting potential of the ocean. If there is a large enough increase in the lateral melt to result in insufficient melting potential, both the lateral and basal melt will be reduced proportionally, as described in Sect. 2.1. A simulation (not presented) to explore this impact shows it has only a limited impact on the basal melt, and not enough to explain the observed compensation effect. The third mechanism concerns lateral melt feedback on the basal melt rate via the perturbation of mixed-layer properties. Higher freshwater release from the increase in lateral melt will lower the temperature and salinity of the ocean mixed layer, which will reduce the basal melt rate. However, the lateral melt increase also reduces the ice concentration, lowering the albedo of the ice-ocean system. This increases the absorption of shortwave solar radiation into the mixed layer, raising the temperature of the mixed layer; i.e. it has the opposite effect of the increased freshwater input. These two competing feedbacks explain the overprediction of basal melt in Fig. 4 but underprediction of basal melt in Fig. 7. The increase in total melt observed in Fig. 7 will likely correspond to a more efficient use of the available melt potential and the aforementioned albedo-feedback mechanism. The interaction between the mixed layer and FSD is further explored through the (P) and (Q) sensitivity studies where the minimum mixed-layer depth was reduced and increased respectively. These studies provide further evidence that mixed-layer feedbacks are not a leading-order effect of the FSD, given the very small perturbations of the melt component from the stan-fsd simulation. Larger changes are seen for the sea ice extent and volume metrics. However, the same mixed-layer feedbacks that change the melt rates can also independently influence the freeze-up rate of sea ice; hence it is not possible to directly attribute the changes produced by varying the minimum mixed-layer depths specifically to WIPoFSD-related feedbacks. It should also be noted that the prognostic mixed-layer model used here provides a limited representation of sea ice-ocean interactions and feedbacks.
The strength of these interactions may increase within a fully coupled sea ice-ocean model (Rynders, 2017).
The series of sensitivity studies to both the floe size parameters and other aspects of the WIPoFSD model are useful to understand the limitations of the model. An important result is the limited sensitivity of the model to the m 1 , m 2 , and α shape parameters, i.e. experiments (L)-(O), with significant perturbations of these parameters reducing the sea ice extent by around 1 % or less. These are additional constants needed to calculate the lateral melt rate beyond floe size. If a strong sensitivity was found to these parameters, it would suggest that these should be considered as alternative targets rather than the FSD for future model development. Instead, these experiments support the focus on floe size as the primary uncertainty in lateral melt calculation. Experiment (K) showed very little model response to increasing the floe freeze-up timescale, T rel , from 10 to 365 d. This result suggests that the use of more physically derived parameterisations of the floe growth during freezing conditions (e.g. Roach et al., 2018b) would not have a significant impact within the model framework presented here. However, Fig. 13 shows that the seasonal l eff evolution is dependent on the floe restoration rate, and there may be specific events, such as strong winter breakup events, where accurate modelling of floe growth is important to then understand the sea ice evolution during the subsequent melting season.
The sensitivity studies also give insight into the impact of waves on the sea ice cover. In particular, the two sensitivity studies that switch off the lateral melt-l var (G) and wave-l var (F) feedback mechanisms respectively showed that the latter had a stronger influence on both the evolution of l eff and the changes in sea ice area and extent when imposing standard parameters on WIPoFSD. This impact was enhanced through various perturbations to the wave model. The increase in significant wave height (H) and reduction in ice strength (I) are representative of future Arctic conditions when the sea ice is expected to be thinner  with storms of increasing strength and duration (Basu et al., 2018). The results presented here suggest that these changes will have only a limited impact on sea ice extent and volume via the floe size feedback mechanism. The strongest response in sea ice extent and volume was observed with a reduction in the attenuation rate (J). It is important to note that the attenuation rate is a function of floe size, with smaller floes driving stronger attenuation. This creates a feedback where the fragmentation caused by one wave changes the way subsequent waves propagate through the MIZ. It should be noted that the wave component of the WIPoFSD model is a simplified representation of waves propagating into sea ice and involves a The Cryosphere, 14, 403-428, 2020 www.the-cryosphere.net/14/403/2020/ number of approximations. In particular, the directional behaviour of the waves will be more suitable for wind waves than swell waves. As discussed in Sect. 2.2, swell waves have been observed to have longer wavelengths and reduced attenuation rates, suggesting they would interact differently with the FSD than wind waves. More generally, modelling the propagation and energy loss of waves as they travel under sea ice is a complex problem and an area of active research (Meylan et al., 2017), and there have been recent efforts to produce coupled wave-sea ice models (Boutin et al., 2019;Herman, 2017). However, any increase in complexity in modelling the waves will result in increased computational cost. Further observations about wave attenuation in sea ice are needed to judge the complexity of the model approach required to produce sufficient accuracy. As stated above, the model shows a strong sensitivity to the floe size parameters with some selections of the WIPoFSD parameters showing moderate increases in the sea ice extent and volume, and other selections driving reductions of these values by over 50 % in September. The limited observational data available to constrain the selected parameters is therefore a significant challenge of this modelling approach. Furthermore, a not insignificant model response of the order of 5 % relative to ref has been observed to sensitivity experiment (E) performed here to explore the impacts of the nonuniform α. Sensitivity experiments (D) and (E) were performed on the basis of evidence from Stern et al. (2018a, b) that α is not a fixed value and instead evolves spatially and temporally. Whilst it would be interesting to explore the impact of a variable d min , especially considering the strong sensitivity of the model response to this parameter, we do not have an analogous set of observations focusing on how d min may vary in space and time.
The WIPoFSD model used here assumes a power-law distribution with the exponent α, lower cut-off d min , and upper cut-off d max all fixed at constant values. Each grid cell has a locally defined variable FSD tracer, l var , which evolves in response to wave break-up events, lateral melt, and freezing conditions. The use of l var to represent variability within the FSD puts limits on the physical fidelity of the parameterisations of processes that change the FSD in our model. However, if l var is not used to represent variability in the distribution, then within a power-law framework over a fixed floe size range the only other component of the system that can change is the exponent. The exponent in such a set-up becomes an emergent parameter rather than one determined from observations. An important component of this study is to perform sensitivity studies of the sea ice mass balance to the range of exponents seen in observations. An investigation of the evolution of the floe size distribution itself, power law or otherwise, is better approached with a prognostic model of the proximate physical processes, such as in the manner of Roach et al. (2018a). Future improved understanding of the FSD may then allow the development of improved parameterisations of floe size and related processes that do not re-quire the assumptions made in this study regarding the shape and floe size range of the distribution. Upcoming studies, including MOSAiC, should provide further observational evidence to develop these parameterisations. The longer-term aim is the development of a floe size model for use in climate models that can reasonably capture the physical impacts of the FSD on the complete sea ice-ocean-atmosphere system without the full complexity of the prognostic floe sizethickness distribution model. The identification of l eff as a useful floe size parameter may provide a method to report useful FSD information over a larger spatial and temporal scale, as this value can be calculated from the ice perimeter length within a unit area and avoids the need to report a full distribution. This would allow an assessment of the regional, intra-annual, and inter-annual variability of the FSD and identify the FSD parameters and components that best reproduce these desired features. There have been recent efforts to develop techniques to obtain a representative floe size metric from satellite imagery over large spatial and temporal scales, though so far these techniques have only been demonstrated at low resolution (Horvat et al., 2019).
The reference simulation (ref) used in this study underpredicts summer sea ice concentration in the pack ice but overpredicts the concentration at the sea ice edge, consistent with other studies that use the CICE sea ice model ). An analysis of the historically forced simulations used within of the Coupled Model Intercomparison Project Phase 5 (CMIP5) found that coupled models consistently performed poorly in capturing the regional variation in sea ice concentration, showing this problem is not specific to CPOM CICE simulations (Ivanova et al., 2016). This suggests that models currently underestimate the role of the MIZ in driving the seasonal sea ice loss. The WIPoFSD model is shown here to have a non-uniform impact on the sea ice cover, with an enhancement in lateral melt and a corresponding reduction in sea ice concentration within the MIZ, as shown in Fig. 6. Whilst the changes are generally small, it shows that the use of an FSD model, either in the described form or otherwise, may be an important step towards improving the accuracy of sea ice models.

Conclusion
Climate model representations of sea ice currently assume that the size of floes that make up the sea ice is constant; however, observations show that floes adopt a distribution of sizes. A power law generally produces a good fit to observations of the floe size distribution (FSD), though the size range and exponent reported for this distribution can vary significantly between different studies. A power-law-derived FSD model including a waves-in-ice module (WIPoFSD) has been incorporated into the Los Alamos sea ice model coupled to a prognostic mixed-layer model, CICE-ML. In the WIPoFSD model, the FSD is defined by a lower floe size cut-off, upper www.the-cryosphere.net/14/403/2020/ The Cryosphere, 14, 403-428, 2020 A. W. Bateson et al.: Impact of FSD on fragmentation and melt of Arctic sea ice floe size cut-off, and exponent. A variable FSD tracer is also introduced, which varies in response to lateral melting, wave break-up events, and freezing conditions. The lower and upper floe size cut-offs and exponent are set to fixed values. A standard set of parameters for the WIPoFSD model is identified from observations and the results of a sea ice simulation using these parameters is compared to one with a constant floe size of 300 m. Inclusion of the WIPoFSD model within CICE-ML results in increased lateral melt compensated by reductions in basal melt, resulting in only moderate impacts on the total melt. The primary mechanism by which the increased lateral melt reduces the basal melt is shown to be the reduction in available ice area for basal melt. The impact is not spatially homogeneous, with losses in sea ice area and volume dominating in the marginal ice zone (MIZ). These impacts partially correct existing model biases in the standalone CICE-ML model, suggesting the inclusion of an FSD is an important step forward in ensuring that models can produce realistic simulations of the Arctic sea ice. A series of sensitivity experiments explore the limitations of the model. The model does show a strong response to a reduction in wave attenuation rate, suggesting this is an important component in understanding wave-sea ice interactions. Different selections of parameters for the FSD show a large impact on the modelled sea ice state, with some showing a moderate increase in mean September sea ice extent and volume, with others reducing these metrics by over 20 % and 50 % respectively. A newly defined parameter, effective floe size, is found to be a good predictor of model response for simulations where the lower floe size cut-off and power-law exponent are fixed. The impact of a non-uniform exponent was also explored based on observations that these parameters evolve for a given region of sea ice. Results suggest that this parameter could further enhance the differential behaviour seen between pack ice and the MIZ in response to the imposition of an FSD. These sensitivity studies also showed that the choice of WIPoFSD parameters is a source of much larger model uncertainty than other constants used within the lateral melt parameterisation, justifying the focus on developing an FSD model as a priority for improved accuracy of sea ice modelling.
Whilst the model presented here does make a major assumption that the floe size distribution adopts a power law, this is consistent with most observations. Furthermore, it has been shown that the model can easily be modified to adapt to additional findings such as the inclusion of a non-uniform exponent. This means the WIPoFSD model is a useful tool for assessing the importance of the FSD in the evolution of sea ice, particularly the seasonal retreat. Climate models require an important balance to be maintained between physical fidelity and computational expense. The simplicity of the WIPoFSD model makes it a useful stepping stone for the development of new parameterisations of floe size within climate models that can reasonably capture the physical impacts of the FSD without a large computational cost. Planned ob-servational studies such as MOSAiC should help in the development of these novel parameterisations.
Code and data availability. Model output used in this paper is publicly available via the University of Reading research data archive (https://doi.org/10. 17864/1947.223;Bateson et al., 2019). Please contact the corresponding author to discuss access to the code.
Author contributions. LH, with support from YA, developed the original version of the WIPoFSD model with a coupled CICE-NEMO framework. DS adapted the model into the CPOM CICE stand-alone set-up. AWB further developed the WIPoFSD model and completed the simulations and analysis under the supervision of DLF, DS, LH, JKR, and YA. DS provided additional technical support. AWB composed the paper with feedback and contributions from all authors.
Competing interests. The authors declare that they have no conflict of interest.