Advances in Modeling Interactions Between Sea Ice and Ocean Surface Waves

Recent field programs have highlighted the importance of the composite nature of the sea ice mosaic to the climate system. Accordingly, we previously developed a process‐based prognostic model that captures key characteristics of the sea ice floe size distribution and its evolution subject to melting, freezing, new ice formation, welding, and fracture by ocean surface waves. Here we build upon this earlier work, demonstrating a new coupling between the sea ice model and ocean surface waves and a new physically based parameterization for new ice formation in open water. The experiments presented here are the first to include two‐way interactions between prognostically evolving waves and sea ice on a global domain. The simulated area‐average floe perimeter has a similar magnitude to existing observations in the Arctic and exhibits plausible spatial variability. During the melt season, wave fracture is the dominant FSD process driving changes in floe perimeter per unit sea ice area—the quantity that determines the concentration change due to lateral melt—highlighting the importance of wave‐ice interactions for marginal ice zone thermodynamics. We additionally interpret the results to target spatial scales and processes for which floe size observations can most effectively improve model fidelity.


Introduction
Sea ice is a critical component of the coupled polar climate system. An aggregate material, sea ice is composed of floes that range in horizontal scale from meters to many kilometers. These floes are particularly important in the marginal ice zone, a spatially and temporally variable region separating the open ocean from the interior pack ice, where ocean surface waves propagate into the sea ice.
Sea ice floes modify the energy spectrum of ocean surface waves as they travel through ice, dissipating and scattering wave energy as they interact with the sea ice cover Meylan & Squire, 1994;Montiel et al., 2016;Squire, 2007). Simultaneously, ocean surface waves cause sea ice to fracture (Mellor, 1986;Toyota et al., 2006) and prevent the formation of extensive sheets of flat ice that would otherwise occur in quiescent conditions (Roach, Smith, & Dean, 2018;Shen, 2004;Shen et al., 2001). The fragmentation of sea ice changes its mechanical properties (Feltham, 2005;Shen et al., 1987), melting at the edge of sea ice floes (Horvat & Tziperman, 2017;Steele, 1992), the ocean mixed-layer eddy field (Horvat et al., 2016;Horvat & Tziperman, 2018), and air-sea heat and momentum transfer (Asplin et al., 2014;Birnbaum & Lüpkes, 2002;Tsamados et al., 2014), which feed back strongly to the polar climate system Asplin et al., 2014;Roach, Horvat, et al., 2018). Because fragmented sea ice is easier to navigate through, sea ice fracture statistics are also of interest to polar shipping, maritime activities, and native communities (Meier et al., 2014).

10.1029/2019MS001836
To first order, ocean surface waves affect sea ice by modifying the horizontal scale of individual floes (their "size")-a state variable that is not yet included in most sea ice models used for climate study (e.g., CICE5; Hunke et al. (2015))-rather than sea ice concentration or thickness. The attenuation of waves by sea ice depends strongly on floe size (Meylan & Squire, 1994;Montiel et al., 2016). Thus, information about the floe size distribution (FSD; Rothrock and Thorndike (1984)), the fraction of the ocean surface covered by floes at each size, is required to understand the coupled evolution of sea ice and ocean waves in polar regions.
The FSD varies widely across space and time as floes are influenced by many different physical processes, with floe sizes spanning orders of magnitude even within the area of a typical climate model grid cell (e.g., Steer et al. (2008)). A number of diagnostic parameterizations of the FSD have been developed in recent years (Bateson et al., 2019;Bennetts et al., 2017;Boutin et al., 2019;Williams et al., 2013a;Williams et al., 2013b;Zhang et al., 2015;Zhang et al., 2016). In contrast, building on the theory developed by Horvat and Tziperman (2015), Roach, Horvat, et al. (2018) presented a sea ice model with a fully prognostic subgrid-scale FSD parameterization. In this representation, floes change size as a result of lateral melting, lateral growth, fracture by surface waves, and, in freezing conditions, forming on the surface of the ocean and welding with one another. The model simulations were initialized without sea ice, and thus, this was the first global model to simulate an emergent, freely varying FSD arising from different physical processes. Waves were included as an external forcing outside the ice and propagated into the ice using a simple exponential attenuation scheme.
We build upon this earlier work to show results from a new coupled ocean wave and sea ice model that includes a prognostic FSD and emergent wave-ice feedbacks. Motivated by the observations presented by Roach, Smith, and Dean (2018) and the underestimation of average floe sizes in wave-free regions in Roach, Horvat, et al. (2018), we implement a new parameterization for ice formation that determines the size of new sea ice floes based on the properties of the ocean surface wave field. The ocean surface wave model and FSD-resolving sea ice model are coupled every 24 hr through a new scheme for floe size-dependent attenuation (Meylan et al., 2019), such that the wave field simulated in the wave model is functionally dependent on the evolving FSD, and in turn, floe fracture and new ice formation in the sea ice model are dependent on the wave field.
The majority of previous FSD studies using both models and observations Toyota et al., 2006;Williams et al., 2013a;Zhang et al., 2015) have approximated the FSD as a power law distribution, the so-called "power law hypothesis" (Clauset et al., 2009). This simplified approach assumes the true multiscale FSD is a straight line in log-log space and, thereby, condenses all degrees of freedom to a reduced number of free parameters, such as the power law slope and minimum floe size. However, the applicability of the power law to FSD observations has been challenged repeatedly (Herman, 2010;Horvat et al., 2019;, and assumption of such scaling behavior can significantly bias model results (Horvat & Tziperman, 2017). In contrast, the fully emergent model described here provides for scale-variant behavior, and we use metrics that do not require assumptions as to the distribution's shape. In this model, the FSD primarily impacts the polar climate system through its influence on lateral melt, which modifies sea ice area, thus participating in ice-albedo feedback. As described in section 3, lateral melt is determined by the perimeter of floes exposed to the ocean. Therefore, to quantify the evolution of the FSD and its interaction with waves, we use the lateral perimeter of floes per unit area of sea ice, P i .
In section 4, we investigate the spatial and seasonal variability of P i and explore the sensitivity of our results to different model versions. We further describe the additional observations required for comprehensive model validation in section 5 and conclude in section 6.

Experiment Configurations
To understand the evolution and sensitivity of the FSD scheme to the model changes we implement here, we perform a number of simulations where the underlying Roach, Horvat, et al. (2018) model is altered or new physics are added. The FSD scheme is implemented in CICE5 (Hunke et al., 2015) and is publicly available (see Acknowledgments). In the standard version of CICE5, without an FSD (NOFSD), the state variable is the ice thickness distribution, whereas in the Roach, Horvat, et al. (2018) FSD scheme, the state variable in CICE instead is the joint floe size and thickness distribution (FSTD). Integrating the FSTD over all ice  Roach, Horvat, et al. (2018).
thicknesses gives the FSD, which is analyzed in this study. Similarly, integrating the FSTD over all floe sizes gives the ice thickness distribution.
The FSD model determines floe sizes prognostically by the interaction of five key physical processes: new ice formation in open water, welding of floes in freezing conditions, lateral growth and melt, and fracture of floes by ocean surface waves. The FSD is included as a tracer variable in CICE, such that floes of all sizes are advected according to the standard sea ice momentum equation and floes of all sizes are affected equally by mechanical redistribution. The welding together of floes when the sea surface temperature is at or below the freezing point occurs according to the probability that pairs of floes placed randomly on the domain overlap, at a rate determined by field observations reported by Roach, Smith, and Dean (2018). Lateral melt and growth, respectively, reduce and increase both floe sizes and sea ice concentration. The lateral melt rate is consistent with the standard version of CICE, without a FSD (Hunke et al., 2015;Maykut & Perovich, 1987). The lateral growth rate is computed from the fraction of freezing potential determined to be near floe edges (Horvat & Tziperman, 2015). The sizes of floes fractured by ocean surface waves are computed from a "superparametrization." Realizations of sea surface height are discretized in space at high resolution (1 m), and the sizes of fractured floes are set to distances between extrema of the sea surface height field when the strain across triplets of these extrema exceeds a critical strain. For further details, see Roach, Horvat, et al. (2018). For all experiments here and in Roach, Horvat, et al. (2018), the FSD is discretized into 12 floe size categories, with category widths that follow a Gaussian distribution ranging from a lower limit of 0.5 m to an upper limit of 950 m.
A concise list of the numerical experiments and the naming convention we adopt are given in Table 1. Below, we outline the major changes to the Roach, Horvat, et al. (2018) model: • Changes to new ice production. In the Roach, Horvat, et al. (2018) model, all new sea ice grown initially forms as pancakes and is assigned to the smallest floe size category. This leads to artificially small floes in regions where waves are insignificant. We implement a new method to assign the floe size of new ice production as a function of the surface wave field, with additional detail given in section 2.2. We denote the original ice-production scheme as FSDv1 and the new wave-dependent scheme as FSDv2. • Changes to ocean surface waves. The majority of our experiments has ocean surface waves prescribed in the sea ice component, as in Roach, Horvat, et al. (2018). A new experiment has sea ice coupled to the ocean surface wave model Wavewatch III v5.16 (WAVEWATCH III Development Group, 2016), which simulates wave attenuation in both sea ice covered and open ocean regions. Wavewatch and CICE are coupled each 24 hr, with Wavewatch passing CICE the integrated ocean surface wave spectrum (which drives the FSD fracture response and wave-dependent new ice growth) and CICE passing Wavewatch grid cell-averaged sea ice concentration, sea ice thickness, and mean floe size (which drive the attenuation of ocean surface waves according to the empirical formulation of Meylan et al. (2019)). The coupling frequency is a compromise determined by the sensitivity of respective states and the computational expense of our initial coupling method. The attenuation is the sum of two parts. The first is an empirical fit to a multiple-floe scattering theory by Meylan et al. (2019) that accounts for floe size, thickness, and sea ice concentration. This theory finds that the scattering resembles the linear Boltzmann equation derived by Masson and Leblond (1989) with an error correction that accounts for multiple scattering. Attenuation from this theory becomes negligible when wavelengths are much greater than floe sizes; therefore, we include an additional term based on observations to attenuate waves like the inverse of the wave period squared. We denote the simulation with the coupled wave approach with the suffix "WAVE."Simulations without a suffix have a

10.1029/2019MS001836
prescribed ocean surface wave spectrum derived from a Wavewatch integration conducted independently of the sea ice model and subsequently propagated into the sea ice using the scheme described in Roach, Horvat, et al. (2018). This scheme damps the prescribed wave spectrum according to a simple parameterization of exponential attenuation in sea ice based on scattering (Kohout & Meylan, 2008) with distance into the sea ice along lines of constant longitude. • Change to the ocean component. One experiment (FSDv1-NEMO) has CICE5 coupled to a full-depth "dynamical" ocean model (Nucleus for European Modeling of the Ocean [NEMO]) based on the configuration of Rae et al. (2015) and was previously presented in Roach, Horvat, et al. (2018). We included it here to aid comparison with new experiments. All other experiments are run with a slab ocean model (SOM), as described in Bitz et al. (2012). The SOM has a single layer at a specified depth, h, that is constant in time and an ocean heat transport convergence field, Q flx , that is annually periodic and interpolated from a monthly mean climatology. Both h and Q flx are diagnosed from a long control run of the Community Climate System Model Version 4 (CCSM4). The NEMO-CICE run uses a 1 • tripolar grid, and the SOM-CICE runs use a displaced-pole 1 • grid (gx1v6). The surface wave component has no direct influence on the ocean component.
Sea ice, ocean, and surface waves in all experiments are forced by the atmospheric reanalysis JRA55 (Japan Meteorological Agency, 2013), for both the hindcast and coupled experiments.

New Ice Production
In the first version of our model (FSDv1), new ice that forms over open water is initially assigned to the smallest floe size category to represent the process of pancake ice formation. In reality, pancake floes typically form only in the presence of surface waves and/or winds; in the absence of such forcing, new ice forms as large thin sheets called nilas (Weeks & Ackley, 1986). In the second version of the model (FSDv2), we use the ocean surface wave field to differentiate between new ice formation of both types, a wave-dependent new ice growth.
When ocean surface waves are present during the formation of new sea ice, the curvature of the ocean surface exerts both a tensile (stretching) stress and a differential (bending) stress on newly formed floes as they freeze together (Shen et al., 2001). Because these stresses scale with the size of the object across which they are applied, a maximum floe size can be determined as the maximum lateral size of a floe that does not fracture through either tensile or bending failure mode. In both failure modes, the maximum floe size decreases with wave amplitude and increases with wavelength (Shen et al., 2001). Laboratory experiments using monochromatic waves (Shen, 2004) and field observations (Roach, Smith, & Dean, 2018) suggest that floe size during freezing conditions is primarily controlled by the tensile failure mode.
Following Shen et al. (2001), tensile failure limits floes to a maximum diameter D max , where is the wavelength, W A is the wave amplitude, g is gravitational acceleration, and i is ice density. The tensile stress mode parameter, C 2 (units: kg·m −1 ·s −2 ), is determined by local conditions. It is expected to be a function of temperature and salinity, but the functional relationship is not yet known, so here we estimate a constant C 2 = 0.167 kg·m −1 ·s −2 from Roach, Smith, and Dean (2018).
New ice is assigned the floe-size category containing D max calculated from equation (1). As equation (1) is derived assuming a monochromatic wave field, we assume that the wave fracture process responds to the most energetic wave frequency. From the wave spectrum as a function of frequency, E(f), we define the peak frequency f p as the value of f that maximized E(f). Then using the deep-water surface gravity wave dispersion relation, which gives the "peak wavelength," p .
For consistency with Roach, Smith, and Dean (2018), we assume the wave amplitude is W A = H m0 ∕2, where the spectral height parameter, H m0 , is given by for the wave energy spectrum as a function of frequency f , E(f) (Michel, 1999). The wave parameters used to compute D max are calculated in each grid cell when freezing occurs.

FSD Model Evaluation Metrics
A size of a floe, r, is related to its area, A(r), via A(r) = r 2 . If floes were perfect circles, the geometric factor = . As outlined in Rothrock and Thorndike (1984), a fit that represents the general noncircularity of floes is = 2.64 ≈ 0.84 , which we select here.
The areal FSD, f(r), is then defined for a given region of sea ice such that f(r)dr (dimensionless) equals the area of floes per unit ocean surface area for floe sizes between r and r+dr. The fractional sea ice concentration in a grid cell, c (dimensionless), is therefore obtained via an integral of f(r) over all floe sizes where we use the notation ∫  dr to denote an integral over all floe sizes and ⟨r i ⟩ to denote the i th moment of f(r).
From f(r), we also define the number FSD, F N (r), where F N (r)dr (units: m −2 ) is the number of floes per unit ocean surface area with size between r and r + dr. The areal FSD and number FSD are related, with Assuming the perimeter of an individual floe is equal to 2 r, we may define the perimeter of floes per unit area of ocean, P o , (units: m −1 ) as where we also define the notation ⟨r i ⟩ N as the ith moment of the number FSD.
The total perimeter of floes per unit area of ocean determines the rate of lateral melting and freezing of existing sea ice floes (Steele, 1992). As discussed previously, in the Roach, Horvat, et al. (2018) model, this is the only mechanism by which changes in the FSD directly impact sea ice concentration. Let floes change their size at a rate G r (units: ms −1 ), which is constant as a function of size and where the subscript r indicates that the rate changes floes in the radial rather than the vertical. Then the evolution of the FSD subject to thermodynamic melting,  T , is where the first term describes how floes move between size classes and the second term describes how changing the size of floes also increases or decreases the area of sea ice (Horvat & Tziperman, 2015).
Integrating equation (6) over all floe sizes, the first term vanishes, such that This reveals the total change of sea ice area due to thermodynamic melting and freezing along floe boundaries to be linearly related to the perimeter of floes per unit area of ocean, P o .
For comparison, the concentration change due to lateral melt with constant rate G r , for a region of floes with constant diameter L, is given by Steele (1992) as

10.1029/2019MS001836
The standard version of CICE5 without an FSD uses this equation with L = 300 m (Hunke et al., 2015), giving c t ≈ [0.0159 m −1 ]G r c. In this study we will principally consider the perimeter of floes per unit area of sea ice, P i ≡ P o ∕c, to emphasize how changes in the FSD arise in a way that is not affected by changes to sea ice concentration through other processes. Then, c t = P i G r c, and lateral melt schemes with and without an FSD (equations (9) and (10)) can be directly compared. P i is proportional to the ratio of FSD moments The standard CICE5 model therefore assumes a fixed value P i = 0.0159 m −1 . For a given area of sea ice, larger values of the perimeter per unit sea ice area implies the sea ice cover is composed of a greater number of smaller floes and is hence more fragmented.
Previously, Roach, Horvat, et al. (2018) defined the sea ice representative radius, r, to quantify how different processes that act on the FSD contribute to its evolution. Because r is a positive moment of the FSD, larger sizes attain a higher weight in the calculation of r. On the other hand, P i is a negative moment of the FSD, weighted more heavily by smaller sizes, so P i is more relevant for thermodynamic melting and freezing of floes. Because the FSD model implemented here has been developed in part to study how regions with small floes, such as the marginal ice zone, affect sea ice thermodynamics, floe perimeter per unit ice area will prove a more appropriate metric for the purposes of this study. Representative radius may be a more relevant choice in other situations, such as in the computation of parameterizations for wave propagation into the ice pack (Meylan et al., 2019;Williams et al., 2013a) and large-scale fracture processes. Note that both representative radius and floe perimeter can be integrated over surface areas to obtain meaningful hemispheric quantities, unlike other shape characteristics of the FSD such as the assumed power law exponent (Horvat & Tziperman, 2017).

The Role of Model Physics on Sea Ice Perimeter
We begin by examining spatial and seasonal variability in floe perimeter per unit ice area, P i , in both hemispheres, in model experiments that include different physics. Recall that all model experiments are listed in Table 1 and are outlined in section 2.  Figures 1k and 1l) assumes that all floes have a constant diameter of 300 m, resulting in a fixed floe perimeter per unit area of sea ice of 15.9 km −1 , (section 3). Experiments including an FSD have considerably higher values of spatially averaged P i than this fixed value. All experiments that include an FSD have similar modes of seasonal and spatial variability in floe perimeter per unit ice area: larger in winter (more fragmented ice) and smaller in summer (less fragmented ice) (Figures 1k and 1l). As we shall see below, this winter-summer difference occurs principally because smaller floes are frequently produced when waves encounter sea ice during winter growth. Annual variability in P i is typically smaller than intermodel variability, reflecting the large variation arising from modifications to model physics. Below, we summarize the changes resulting from each of the model physics variations outlined in section 2.   Roach, Horvat, et al. (2018) coupled to the NEMO ocean model, (c) same as (b) but using a slab ocean model, (d) same as (c) but now including the effects of wave-dependent ice growth, and (e) same as (d) but including the fully coupled Wavewatch III model. The seasonal cycle of Arctic P i , integrated over the full Northern Hemisphere, is shown in (k). (f)-(j) repeat (a)-(e) but for the Antarctic in September, with the seasonal cycle in Antarctic P i shown in (l). On (a)-(j), contours show 15% and 80% sea ice concentration, and perimeter is shown for sea ice concentrations greater than 15%.

Dynamical Versus SOM
the Antarctic). The SOM uses a single-layer ocean depth and ocean heat transport convergence prescribed from a fully coupled run of CCSM4, developed independently of NEMO, and therefore, structural differences between these two oceanic configurations result in differing sea surface temperatures. FSD evolution is coupled to oceanic freezing/melting potential, which depends on sea surface temperatures, due to its role in lateral melt and floe freezing processes. The consistently lower hemispheric perimeter in FSDv1-NEMO relative to FSDv1 may reflect cooler sea surface temperatures and a higher freezing potential with corresponding increased floe welding. Qualitatively, FSD results from the full ocean and slab ocean experiments are similar, particularly in the Arctic, which guided our decision to proceed using only SOM experiments.

Changes to New Ice Production
The addition of the wave-determined new ice growth scheme (FSDv2) results in a large change in floe perimeter per unit ice area relative to the original model (FSDv1) year-round (Figures 1k and 1l). The year-round spatial average P i decreases fourfold in the Arctic and by approximately 50% in the Antarctic. Because floes are permitted to initially form as sheet ice in the absence of waves, FSDv2 is less fragmented in the center of the ice pack compared to experiments that assign all new ice into the smallest size category (Figures 1c, 1d, 1h, 1i, 1k, and 1l), particularly in the Arctic, where the simulated waves do not often reach the ice interior. This introduces much greater spatial heterogeneity in simulated floe sizes. In the Arctic, changes to new ice production result in a well-defined transition between wave-affected regions, often defined as the marginal ice zone, and the interior sea ice. In the Antarctic, where the ice cover is more seasonal and the modeled influence of waves is greater, many floes still form as pancakes initially. However, the floe perimeter per unit ice area decreases in areas of multiyear sea ice, such as the Weddell Sea.

Coupled Versus Prescribed Ocean Surface Waves to the Sea Ice Model
Compared to the changes arising from the new ice growth scheme, the differences in spatial variability in floe perimeter per unit ice area arising from the addition of wave coupling (FSDv2-WAVE) are smaller (Figures 1d, 1e, 1i, and 1j). The WAVE experiment uses a comprehensive scheme for attenuation within sea ice (Meylan et al., 2019), with increased physical fidelity relative to the uncoupled runs that use a simple  wave scheme that propagated waves from outside the sea ice into the sea ice by following meridians. Relative to uncoupled runs, the coupled WAVE experiment is more fragmented in the central Arctic due to increased wave propagation and less fragmented in the central Antarctic as a result of reduced wave propagation. The WAVE experiment determines the attenuation of waves in sea ice as a function of floe size, thickness and concentration, and permits wave generation within the sea ice zone, resulting in significant regional differences in perimeter compared to the uncoupled runs. These results are physically plausible, given the limitations of the initial attenuation scheme, which are particularly apparent for the land-locked central Arctic. The WAVE experiment is the most advanced and physically faithful simulation in the SOM configuration, and thus, we discuss it further below. Figure 3b shows how different floe-size-modifying processes contribute to the seasonal cycle in floe perimeter per unit ice area, shown in Figure 3a, for the FSDv2-WAVE experiment in the Arctic. This is the area-weighted average monthly tendency in floe perimeter per unit ice area from the five floe-size-modifying processes. Recall that floe perimeter per unit ice area is weighted more heavily by smaller floe sizes, while representative radius is weighted more heavily toward larger floe sizes. The seasonal cycle is mostly driven by new ice growth, except in summer, when the dominant processes driving changes in floe perimeter per unit ice area are wave fracture and lateral melt. Therefore, in the model, wave-controlled new ice growth controls floe perimeter before lateral melt begins, and fracture of floes by waves becomes more important. This highlights the role of wave-ice interactions in determining basin-wide sea ice floe statistics. Figure 3c shows the seasonal cycle in floe representative radius, where changes naturally have the opposite sign to floe perimeter (low perimeter indicates that sea ice is less fragmented, suggesting larger floes). The seasonal cycle in representative radius is driven by any process that impacts sea ice area, not just floe-size-modifying processes, and is driven principally by standard sea ice area changes. Representative radius peaks in September, when smaller floes have melted away, leaving only large floes. Of the floe-size-modifying processes, welding, wave fracture, and new ice growth are the most important for determining representative radius. Welding stands out more for representative radius compared to the tendencies for floe perimeter, because it has a greater effect on the large floes. Conversely, lateral melt has a minimal impact on representative radius, as it principally affects the smaller floes. Results from the Antarctic are similar but show a sharper increase (decrease) in perimeter per unit ice area (representative radius) at the end of the summer driven by enhanced formation of new small floes and enhanced wave fracture relative to the Arctic.

Interactions With Ocean Surface Waves
The FSDv2-WAVE experiment presented here includes two-way interactions between prognostically evolving waves and sea ice on a global domain. Figure 4 shows significant wave heights simulated by Wavewatch in the FSDv2-WAVE simulation. The 2000-2014 average spatial fields (Figures 4a, 4b, 4d, and 4e) show waves penetrating into the ice cover, some to large distances but entirely contained within the 80% sea ice concentration contour, sometimes used to define the marginal ice zone (e.g., Strong et al., 2017). Probability distributions of all monthly mean significant wave heights in sea ice of concentrations greater than 15% reveal larger wave heights in the Antarctic than the Arctic, reducing as a function of concentration.

Impact on Lateral Melt
The addition or modification of model physics in our range of experiments results in strongly differing floe perimeters. Figure 5 shows the annual-mean hemispheric-average floe perimeter per unit ice area for all experiments listed in Table 1, with the corresponding total lateral melt that occurs for each pole throughout the year.
The total amount of lateral melt differs substantially between the experiment without an FSD, where all floes are assumed to be 300 m in diameter, and the various versions of the FSD model, increasing near-linearly with sea ice floe perimeter per unit ice area (Figures 5a and 5b). The highest floe perimeter case, FSDv1, has a sixfold increase in lateral melt relative to the simulation without an FSD. In all cases, some compensating decrease in basal melt occurs (Figures 5c and 5d). The FSDv1-NEMO experiment lies outside the apparent relationship in the Arctic for lateral melt and in the Antarctic for basal melt. Recall that the change in concentration due to lateral melt depends on the product of the lateral melt rate and floe perimeter. The NEMO-CICE simulation has different sea surface temperatures, and therefore lateral melt rates, so we would expect a different gradient for a linear relationship than the SOM experiments.  Figure 6 shows the difference in 2000-2014 mean sea ice volume per unit area between each slab ocean configuration including a FSD and the NOFSD simulation, in the months where lateral melt is greatest. There are major reductions in sea ice volume per unit area in all three FSD experiments in the Southern Hemisphere and toward the ice edge in the Northern Hemisphere, corresponding to areas of increased floe perimeter in Figure 1. The assumption of all new floe growth as pancakes in FSDv1 leads to the highest perimeter among all of our experiments, and we suspect that FSDv1 overestimates the reductions in sea ice. The changes to the sea ice volume in FSDv2 and FSDv2-NEMO seem more plausible. FSDv2 shows overall increased ice volume in the central Arctic, while FSDv2-WAVE shows some patches of increased ice volume. These spatial differences in the areas of increased ice volume result from different propagation of waves in the sea ice zone between the simulations with waves included as forcing versus the wave-coupled simulation. Addition of the FSD also results in reductions in sea ice concentration of up to around 30%, with a broadly similar pattern to sea ice volume, but with some differences in the central Arctic (not shown).
A full mechanistic investigation of these impacts is beyond the scope of the present study. We stress that these results are obtained using prescribed atmospheric forcing-which constrains the location of the sea ice edge-and a slab ocean-which neglects any negative or positive feedbacks by ocean dynamics. The climate impact of changes to floe perimeter should be investigated using a fully coupled modeling system in future work. We include Figure 6 as an indication of sea ice impacts purely from the sea ice physics changes in the absence of ocean and atmospheric coupled interactions. The sensitivity shown in Figure 5 and potential impacts shown in Figure 6 highlight the requirement for observations to constrain simulated floe perimeter, the assumed lateral melt rate, and the partitioning between basal and lateral melt, as discussed further in section 5. which requires identifying and analyzing many individual floes. Perimeter is weighted toward smaller floe sizes, as noted above, which are better resolved in the model than in observations. Perovich (2002) made observational estimates of floe perimeter per unit ice area from aerial photographs captured during the 1998 Surface Heat Budget of the Arctic Ocean (SHEBA) field campaign, on 10 dates between June and September. These estimates represent areas covering a 50 km square centered on the SHEBA site, a Lagrangian drifter, and resolve floes to 10 m-larger than the smallest floes in our model but at finer scales than many other observational estimates (e.g., Stern, Schweiger, Stark, et al., 2018).

Toward Model-Observation Comparisons
Because the SHEBA data represent just a handful of point measurements, these perimeter observations cannot be used to evaluate the modeled seasonal cycle. Still, they may be useful to compare to model output. All SHEBA locations have ice concentrations greater than 80%, indicating pack-ice-type conditions. The mean perimeter per unit ice area from the observations is approximately 20 km −1 . This value is of the same order of magnitude of floe perimeter simulated by FSDv2 and FSDv2-WAVE in the interior ice pack (Figures 1d  and 1e). A full model validation would require observations and model output over comparable temporal periods, ice types, and spatial scales. Horvat et al. (2019) provides estimates of FSDs over a larger spatial and temporal scale, covering the entire Northern Hemisphere and all nonsummer months from 2010-2018. Such a spatial and temporal scale is sufficient to evaluate a seasonal cycle and regional change. These estimates, however, are made using the CRYOSAT-2 altimeter and therefore have minimum resolved floe sizes at the subkilometer scale, where their measurement tool experiences significant resolution uncertainty. As the model simulations here focus on floes smaller than 1 km and floe perimeter is most sensitive to these small floes, such a comparison is best made in terms of representative radius and for experiments using extended floe size categories, as in Horvat et al. (2019).

Discussion
Development of models of the sea ice FSD has thus far been heavily limited by the paucity of observations suitable for model validation. Sea ice floe sizes can be estimated from images captured by cameras carried by aerial vehicles, visible-band satellites (e.g., Landsat, MODIS, and MEDEA), and synthetic aperture radar (e.g., TerraSAR-X and RADARSAT-2) using image processing tools to separate individual floes. The spatial and temporal coverage of existing observations is limited , and the range of floe sizes resolved varies between studies (Stern, Schweiger, Stark, et al., 2018). Naturally, observational uncertainty is introduced by (and not limited to) (1) lack of distinct floe edges, (2) floes at sizes below the image resolution, (3) floes extending beyond the image frame, and (4) visual interference such as cloud cover in visible-band imagery. Some studies use similar methods (Stern, Schweiger, Stark, et al., 2018), but a comprehensive intercomparison of image processing methods for floes has not yet been carried out.
Nearly all studies fit power laws to either the cumulative or noncumulative floe number density distribution (Stern, Schweiger, Stark, et al., 2018). In most cases, the value of the exponent is obtained from a least squares regression on a log-log plot, which may lead to erroneous overconfidence in the applicability of the power law hypothesis, and depends on the binning of the data (Clauset et al., 2009;Stern, Schweiger, Stark, et al., 2018). Although power law exponents are often discussed, the data do not universally show straight lines in log-log space (Herman, 2010;Horvat et al., 2019;Steer et al., 2008;Wang et al., 2016); and very few have been assessed using statistical tests of the power law hypothesis (Stern, Schweiger, Stark, et al., 2018). In our experience of FSD model development, power law exponents have not proven useful for understanding the physical processes that affect sea ice floes and, with limited understanding of the several heuristic parameters that set their scaling, are not easily compared to model results. To advance modeling capabilities, we suggest a series of observations to improve the capacity of floe-sensitive sea ice models: 1. Estimates of floe perimeter per unit area in different seasons and over a wide spatial domain. Floe perimeter is the key variable that determines the sea ice concentration change that occurs due to lateral melt and, if models are designed to examine the thermodynamic impacts of fragmented ice, must thus be validated directly. By definition, floe perimeter per unit area is weighted toward smaller floe sizes (section 3), and so high-resolution observations are required to compare with the small floes resolved in the model. We show steps toward an example comparison in section 4.4, but much more data are required for model validation.
2. Estimates of significant wave heights in sea ice during different seasons and at the basin scale. We have shown above that sea ice floe perimeter differs significantly depending on the description of ocean surface wave attenuation. A data set of calibrated wave heights from remote sensing would permit calibration of the wave attenuation scheme and a comparison of the methods currently available in Wavewatch III (WAVEWATCH III Development Group, 2016). In situ calibration of attenuation has been an active area of research (e.g., Meylan et al., 2014), but only satellite products can achieve the coverage required for model development at this early stage. The possibility of obtaining probability distributions of significant wave heights in ice from synthetic aperture radar observations has recently been successfully explored (Stopa et al., 2018), but achieving such measurements at the basin scale is still not feasible. 3. In situ observations of basal and lateral melt rates. Our results depend on the standard sea ice model assumption of a lateral melt rate, w lat = m 1 (T w − T ) m 2 , a function of the difference between the sea surface temperature T w and the freezing point T f , with coefficients m 1 and m 2 the best fit to data from a single field site quoted by Maykut and Perovich (1987). Furthermore, we, along with Bennetts et al. Besides observational validation of simulated floe perimeter, determining the long-range climate impacts of changes to the partitioning of sea ice melt requires fully coupled, centennial-scale simulations. Investigation of the sea ice response to fragmentation, for example, during storms, on seasonal and subseasonal time scales would require a coupled wave model that exchanges information at higher temporal resolution than we have implemented here. We hope that such modifications to the modeling configuration are made in future work to determine the relevance of enhanced lateral melt for the polar climate system.

Conclusions
In this study, we discuss advances in modeling wave-ice interactions using a sea ice model and a wave model, which are both widely used individually for applications from seasonal forecasting to climate studies. We make use of the floe perimeter per unit of sea ice area as a metric for model assessment due to its role in determining the sea ice concentration change arising from lateral melt. We show results from a series of model experiments that increase in physical fidelity through addition of a scheme that applies a tensile stress limitation to the size of new ice floes and coupling to an ocean surface wave model. The difference in floe perimeter per unit ice area arising from incorporation of different physical processes is substantial at this stage of model development. When validation products are available, future work should more robustly test impact of varied floe size category resolution.
The final experiment (FSDv2-WAVE) includes the new wave-dependent ice growth scheme and wave model coupling, making it our most physically faithful and state-of-the-art modeling system. Its simulated seasonal cycle in floe perimeter per unit ice area has a similar magnitude to existing observations in the interior pack ice and suggests characteristics that could be tested with further observations on an extended spatial and temporal domain. Examining how different FSD processes contribute to floe perimeter, we find that new ice growth is most significant in nonsummer months and wave fracture is the most significant during summer months, closely followed by lateral melt. Thus, during the months when lateral melt occurs, wave fracture is the dominant FSD process determining sea ice floe perimeter per unit ice area in our model, making uncertainty in the representation of sea ice fracture by ocean surface waves particularly significant for simulated lateral melt. We also demonstrate properties of the simulated wave field in the marginal ice zone, which are in line with intuition.
The total amount of lateral melt in our simulations depends strongly on the simulated floe perimeter, which in turn depends strongly on the properties of ocean surface waves as they travel through sea ice. Although individual model processes have been constrained by observations (Meylan et al., 2019;Roach, Smith, & Dean, 2018), new, large data sets from remote sensing for sea ice floe perimeter and wave heights in ice are required to make further progress on model calibration and evaluate which of our configurations is most realistic. Our recent advances in modeling of the sea ice FSD prioritize specific observational requirements and will enable new opportunities for seasonal forecasting and numerical investigations of subgrid-scale sea ice processes in the future. NA16NWS4620043, Years 2017-2021, with the National Oceanic and Atmospheric Administration (NOAA) and U.S. Department of Commerce (DOC). S. D. was funded via Marsden Contract VUW-1408. Michael Meylan provided the floe-size-dependent wave attenuation scheme, described in Meylan et al. (2019). We are grateful to Donald Perovich for sharing observational data from Perovich and Jones (2014). The authors also wish to acknowledge the contribution of New Zealand eScience Infrastructure (NeSI) high-performance computing facilities to the results of this research. We are grateful to the two anonymous reviewers for suggestions that strengthened the manuscript. Model output used in this manuscript is publicly available via Zenodo (https:// doi.org/10.5281/zenodo.3463580). Model code is publicly available online (https://github.com/lettie-roach).