WIFF1.0: A hybrid machine-learning-based parameterization of Wave-Induced sea-ice Floe Fracture

Ocean surface waves play an important role in maintaining the marginal ice zone, a heterogenous region occupied by sea ice floes with variable horizontal sizes. The location, width, and evolution of the marginal ice zone is determined by the mutual interaction of ocean waves and floes, as waves propagate into the ice, bend it, and fracture it. In previous work, we developed a one-dimensional “superparameterized" scheme to simulate the interaction between the stochastic ocean surface wave field and sea ice. As this method is computationally expensive and not bitwise reproducible, here we use a pair of neural 5 networks to accelerate this parameterization, delivering an adaptable, computationally-inexpensive, reproducible approach for simulating stochastic wave-ice interactions. Implemented in the sea ice model CICE, this accelerated code reproduces global statistics resulting from the full wave fracture code without increasing computational overheads. The combined model, WaveInduced Floe Fracture (WIFF v1.0) is publicly available and may be incorporated into climate models that seek to represent the effect of waves fracturing sea ice. 10


Introduction
Sea ice is a composite of individual pieces, called floes. Floes are discrete solid pieces of sea ice with horizontal geometric length scales (sizes) ranging from meters to tens of kilometers. Yet in climate model simulations, sea ice is treated as a continuum field (Hibler, 1979;Golden et al., 2020). To represent the effect of changes to sea ice floe geometry in climate models, modeling centers have started to incorporate prognostic sub-grid-scale parameterizations of the sea ice floe size distribution 15 (FSD) as a core aspect of their model physics (Horvat and Tziperman, 2015;Zhang et al., 2015;Roach et al., 2018;Boutin et al., 2018;Bateson et al., 2020;Hunke et al., 2019). The FSD is a probability distribution f (r), where f (r)dr is the percentage of a grid area comprised of floes with a size between r and r + dr.
One major advantage of simulating a sub-grid-scale FSD is it permits the representation of coupled interactions between sea ice floes and ocean surface waves. Ocean surface waves propagate into, are attenuated by, and fracture sea ice in regions 20 known as the marginal ice zone (MIZ, Wadhams et al., 1988;Langhorne et al., 1998;Squire et al., 1995). Wave-affected sea-ice-covered regions are observed to be several million square kilometers in size in both hemispheres (Horvat et al., 2020).
When waves propagate into and fracture sea ice, there is no direct change to sea ice concentration or thickness. Instead, this fracture process alters sea ice floe sizes, increasing the sensitivity of the sea ice to external thermodynamic or dynamic forcing (Steele, 1992;Feltham et al., 2006;Horvat et al., 2016;Horvat and Tziperman, 2017) and altering the attenuation of ocean 25 wave energy (Perrie and Hu, 1996;Meylan et al., 2021). Thus there is a hypothesized coupling between sea ice loss and wave activity -waves propagate into the ice, fracturing it, causing it to melt, and enhancing wave propagation (Kohout et al., 2011;Asplin et al., 2012Asplin et al., , 2014. Floe fracture by waves is the dominant process driving changes in floe perimeter in summer months when the most lateral melting occurs . In a climate model, coupled wave-ice feedbacks are related to two sub-grid-scale distributions: the FSD, f (r), and the ocean 30 surface wave spectrum, S(λ), where λ is the wavelength and S(λ)dλ = E is the total wave energy. Heuristic parameterizations have been developed to relate bulk properties of the ocean surface wave field to the FSD -generally assuming that fractured floes follow a power-law distribution (Williams et al., 2013;Zhang et al., 2015;Bateson et al., 2020). Yet there is conflicting evidence about whether power-law FSDs are observed in nature (Herman, 2013;Stern et al., 2018a, b;Horvat et al., 2019), and whether power-law size distributions are generated by the process of wave-induced floe fracture (Horvat et al., 35 2016;Herman et al., 2021).
One challenge in relating S(λ) and f (r) is that the response of sea ice is determined by the two-dimensional ocean height field, a stochastic (i.e. one of many possible) representation of the ocean wave spectrum. An approach designed by Horvat and Tziperman (2015) was to generate a high-resolution, one-dimensional ocean surface wave field in every ice-covered grid cell, and explicitly resolve the strain experienced by the ice, the attenuation of wave energy, and the resulting statistics of 40 sea ice fracture. This method (hereafter called SP-WIFF) is broadly analogous to a "super-parameterization" (SP) approach used in coupled climate models for parameterization of cloud effects (Randall et al., 2003;Grabowski, 2004) or oceanic deep convection (Campin et al., 2011). SP-WIFF is included as a component of the sea ice model CICE (Roach et al., 2018;Hunke et al., 2019), and has been coupled to an ocean surface wave model . Running SP-WIFF incurs high computational costs -in the case of its incorporation in CICE, it overall computation times by an order of magnitude (see 45 Sec. 4). The stochastic generation of wave heights also means SP-WIFF is not bitwise reproducible, often a necessary feature for climate model development.
To provide a computationally inexpensive, flexible, and bit-wise reproducible sub-grid-scale parameterization of waveinduced ice fracture, we here present an accelerated parameterization of wave-induced floe fracture using a pair of neural networks for input classification and sea ice fracture (NN-WIFF). The full code, WIFF1.0 (Wave-Induced Floe Fracture), 50 is publicly available, and contains SP-WIFF and NN-WIFF, along with code for retraining or adapting NN-WIFF to prescribed error thresholds, variable input/output variables, and model configurations. Trained using input from coupled CICE-WAVEWATCH 3 simulations, NN-WIFF produces global sea ice variability that is not statistically significantly different from those produced by SP-WIFF, while reducing its computational overhead by more than 95%.

55
The super-parameterization of sea ice fracture by ocean waves was derived in Horvat and Tziperman (2015) (Sec. 2.3), its statistical properties explored in Horvat and Tziperman (2017) (Sec. 3.4), was implemented with offline wave-ice interactions in a sea ice model in Roach et al. (2018) (Sec. 2.4), and was introduced into a fully-coupled wave-ice model in Roach et al. (2019). We summarize it here.
Consider a region corresponding to a climate model grid cell where ocean surface wave energetics are described by a uni-60 directional wave spectrum S(λ i )dλ i (integrated over all angles). For floes with horizontal size between r and r + dr, a fraction of the domain Ω(r, t) dr dt (unitless) is fractured by ocean surface waves over a period dt. That fractured area has its own floe size distribution -the fraction of Ωdr dt that now belongs to floes with size between s and s + ds is F (r, s) ds (unitless). The rate of change of area of floes of size r due to fracture by ocean surface waves is then, (1)

65
The first term in Eq. 1 is the loss of floe area at size category r per unit time, and the second is the increase in floe area in size category r due to the fracture of floes of larger sizes. The limits of integration reflect the fact that any given floe can not fracture into a larger floe, and therefore F (s, r) = 0 for s < r. Evaluation of Ω and F occurs as follows within each ice thickness category: S1: The discrete unidirectional wave spectrum S(λ)dλ is converted to a 1-dimensional ice strain field η(x) of length 10 km,

70
(see (Horvat and Tziperman, 2015, eq. 20-21)) S2: A collection of fracture lengths {L} is found by finding the distance between successive strain extrema that would fracture a floe. These are binned into a probability distribution A(r), where A(r)dr is the percentage of fracture lengths with length between r and r + dr. Note we replace H(r) in Horvat and Tziperman (2015) with A(r) here to avoid confusion with sea ice thickness, H.

75
S3: S1-S2 are repeated (increasing the size of {L}) until convergence, defined as a mean absolute difference between successively updated values of A(r) below 5×10 −4 .
The histogram A(r)dr is then used to compute Ω and F , 80 SP-WIFF therefore involves a stochastic component (S1) and an expensive peak-finding and histogram-generation component (S2). Its execution time is uncertain because the number of steps until convergence (S3) is not known in advance.
Previously, to circumvent these challenges, SP-WIFF was implemented in CICE using a fixed phase in step S1 and only iterated once in step S3. We refer to this single-iteration implementation as SP-WIFF 1 . Despite this simplification, CICE runtimes are increased by approximately a factor of 4 (see Sec. 4) over benchmark estimates without wave fracture. We seek to improve 85 the performance of SP-WIFF by replacing S1-S3 with a neural network trained on SP-WIFF input and output data, which we call NN-WIFF.

Accelerating the parameterization of coupled ice-wave interactions
The NN-WIFF code replaces the full super-parameterization (SP-WIFF) with a parameterization that replicates its statistics but has reduced computational overhead. It consists of two full-connected feedforward neural networks trained using the 90 Python Keras deep learning library, one 100x100 input classification scheme for determining whether NN-WIFF will run, and a second with five hidden layers of 100 nodes each for generating fracture histograms. Network architectures were chosen using a metalearning approach varying loss functions, and numbers of nodes and layers (we provide code for altering network size in the WIFF1.0 release). Training is performed using the 'adam' optimization scheme using the resilient backpropagation method. Rectified linear unit activation functions are used for each of the hidden neurons, and we use a softmax activation for 95 final output layer in both cases as we seek a binary value (for the classifier) and a histogram that sums to one (for the fracture histogram network).
Because SP-WIFF can be run on any arbitrary input, we may generate a training dataset of unlimited size. In this study, we use data from a single year (2009) of a coupled simulation of waves and sea ice (as in the simulation FSD-WAVEv2 in Roach et al. (2019), but here we use a fully-converged wave fracture parameterization instead of a single application of S1-S2, 100 i.e. SP-WIFF instead of SP-WIFF 1 ). We take 6-hourly wave spectra, sea ice concentration, sea ice thickness, and fracture histograms in both hemispheres, a total of 180 million possible input vectors (most without sea ice). We prune training data by including only locations where the local sea ice thickness is less than 10 m, sea ice concentration is greater than 0.01, and the significant wave height H s = 4 √ E is greater than 0.1 m. These thresholds are also checked before executing SP-WIFF in the version of WIFF released here. This eliminates spurious wave fracture calls in areas of anomalous sea ice conditions or low 105 wave energies. We find a total of 17.9 million sea ice points, of which 5.1 million meet the criterion above.
Thus in total, the full training data set comprises N = 5.1 million input vectors of 25 spectral amplitudes, 1 ice thickness, and 1 ice concentration. We assume these include the potential phase space of inputs to the SP-WIFF code, and the potential range of wave energies and sea ice states. These input vectors are identified with a set of output vectors from the "true" SP-WIFF output -in this case 5.1 million 12-category floe probability distributions A i dr i , of which 2.6 million are non-zero (i.e.

110
SP-WIFF leads to floe fracture). We choose to use the 12 category FSD presented in Roach et al. (2019), however this can readily be altered for different uses. In the WIFF1.0 release, we include code for: (see code and data availability statement) -Defining a custom training data set using prescribed input spectra and thicknesses.
-Defining a custom output histogram size (i.e. different FSD categories).
-Designing and training the input classification scheme.

115
-Designing and training the histogram generation scheme.
We also provide data that was used to train networks and produce the results shown here (see Code and Data Availability).
We used coupled model output of SP-WIFF for the 12-category floe probability distributions in our training dataset. These can also be generated offline, and we provide offline code (SP_WIFF_Standalone.m) for generating output histograms for any arbitrary input that is bit-for-bit the same as the code in SP-WIFF.

Input Classification
The histogram A i dr i is a collection of numbers that sums to one, and is the target of our parameterization acceleration. This motivates the use of a softmax activation layer (as the output sums to one). Yet when SP-WIFF is executed, sea ice does not always fracture, and a true output histogram could be a vector of zeros. We therefore require a classification layer that determines when a histogram will be produced by NN-WIFF, and when it will not, to reduce the potential for false positives. 125 We found the bulk of potential false positives could be eliminated by introducing the wave energy cutoff of H s = 4 √ E = 0.1m above, and therefore the classification network is implemented after that cutoff is applied.
The training data for the classifier is the 5.1 million input vectors with corresponding binary values (fracture/no fracture).
This is randomly partitioned 70/30 into training ( 3.6 million) and validation ( 1.5 million) datasets, with an even proportion in each set (51%) leading to fracture. The input classifier then takes an input vector and returns a single value, 0 < χ < 1, which We found that is minimized for χ crit = 0.54 at = 12.5%. This error rate corresponds to a false positive rate of 10.6% and a false negative rate of 14.0%.
135 Fig. 1 shows error characteristics of the classifier evaluated on the validation dataset. In Fig. 1(a-c), we show the relationship between false positive rate (blue) and false negative rate (green) as a function of (Fig. 1a) sea ice concentration, (Fig. 1a) sea ice thickness, and (Fig. 1c) log 10 E. We also plot the run rate (dashed black line), or the fraction of times NN-WIFF is executed, as well as the overall accuracy (black line), equal to 1-.
Generally, the input classifier has high accuracy across the range of ice concentrations and thicknesses, with a flat distribution the largest portion of potential input states, and the 12.5% error in classification does not contribute to statistically significant difference in climate model output. Still, improvements in the input classification may be achievable if desired, for example for higher wave energies, and we include code to re-train the classifier in WIFF1.0.

Neural Network for Generating Fracture Histograms
The fracture histogram network takes a length-27 input vector and outputs a vectorÂ i dr i which sums to one and approximates 150 the histogram A i dr i . As a result, sea ice floes will fracture if this input wave and sea ice statistics pass the input classification.
To assess its performance, we use a custom loss function that minimizes the "representative size error", which is the absolute difference in the predicted histogram in each floe size category, weighted by the floe size in that category, and has units of meters. Available training data is randomly partitioned into a training dataset of 70% of the input vectors, 155 and a validation dataset of the remaining 30% (∼800,000 input vectors). Training is performed until the validation loss fails to improve over 20 subsequent steps. Indexing training data by j, the overall network performance is the average error for the training data, where N is the total number of inputs used to train the network. Aidri andÂidri in each floe size category.
In Figure 2 we display aggregated fit characteristics for the histogram neural net. Figures 2(a-d) display number histograms (red, left axis) of sea ice concentration (C), ice thickness (H), and log 10 E from the validation dataset, as well as of the "true" representative radius R, where from validation dataset input. To evaluate the standardized neural net error, we also compute the standardized size error (SSE) 165 as, which weights errors by the expected mean floe size. Across the entire validation dataset, median SSE is 3.9%.
There is little dependency of SSE on sea ice concentration (2a) or sea ice thickness (2b). Median SSE errors are below 5% across all ice concentrations and thicknesses considered here, with the peak of the interquartile range below 10%. While 170 median errors for both wave energy and R also lie below 10% for all input wave energies and true output floe sizes, there are higher SSE values for low wave energies and high R values. Large values of R can occur for a small number of potential floe fractures, which typically occurs for low wave energies that do not repeatedly fracture the sea ice. In general, the neural network (red line, Fig. 2e) reproduces the mean fracture histogram A i (black line, Fig. 2e) accurately at each floe size category.
Median error between A i dr i andÂ i dr i is less than 6.2% in all categories, and across all 12 categories averages 3.3%.

175
To understand what regions of the phase space give rise to error, in Fig. 3(a-c) we plot 2-dimensional heat maps of SSE, evaluated on the validation dataset, for each of the six combinations of ice thickness, concentration, R, and log wave energy.
Colors are shown only for those coordinates with at least 0.005% of the data (40 points), and the whitepoint of the colormap is set for an error of 10%. Blue values are those with SSE < 10%. Fig. 2 suggests that the input vectors with largest error are those with low input wave energy and high resulting R. In general, we find similar results when examining these two-dimensional maps. In general, we find low error (SSE<=5%) across the range of input sea ice thickness and concentration (Fig. 3a). However expanding the error map in R as a function of thickness (Fig. 3b), concentration (Fig. 3f), and wave energy (Fig. 3c) reveals that the largest part of training error comes from regions with low input sea ice concentration, small thickness, and low wave energy. This combination of inputs gives rise to high output R values, which have the highest SSE values. Although locations with low sea ice concentration are not used as 185 WIFF input, we do not presently exclude locations with low sea ice thickness -that these regions produce pathological output suggests they might be isolated or excluded from future versions of NN-WIFF.

Sea Ice Model Results
We next examine the relationship between versions of the WIFF code when implemented in a free-running sea ice model.
Here 2. (SP-WIFF 1 ) the SP-WIFF model but using a single iteration of steps S1-S2 with fixed phase.
The histogram output A i dr i from CICE is a daily average over 24 hour timesteps. When wave fracture does not occur, A i dr i values are set to zero. Thus though A i dr i is normalized to one if wave fracture occurs, if a grid cell does not always fracture 200 during one of the averaging timesteps, the CICE history output will not be normalized to one. We choose to include such points in this analysis, renormalizing histogram outputs to sum to one.
Over the course of this year, we find that in the NN-WIFF simulations, there are 686,000 gridpoint-days during which NN-WIFF is executed, for a total of 12.3 million NN-WIFF calls overall. This compares to 746,000 and 14.1 million for SP-WIFF and SP-WIFF-1. The discrepancy between function calls is due to the use of the input classifier, which introduces false 205 negatives as discussed in Sec. 3.1.
In Fig. 4(a,b) we repeat Fig. 2   We plot histograms of the errors between SP-WIFF and SP-WIFF 1 (red) and NN-WIFF (blue) in both hemispheres as Fig. 4(c-d) from the same sample used in producing Fig. 4(a,b). Dashed lines are SSE and solid lines are sum of absolute error (SAE) in each category. In both metrics, and in both hemispheres, NN-WIFF outperforms SP-WIFF 1 when it is executed. In the Arctic, median SSE is 8.9% (NN-WIFF) vs 10.5% (SP-WIFF 1 ), and median SAE is 5.8% (NN-WIFF) vs 5.9% (SP-WIFF 1 ).
In Fig. 5(top row) we show monthly-average differences in sea ice concentration between SP-WIFF and NN-WIFF in March and September at each hemisphere. There is little difference between the two, with a maximum global area difference of 47,000 km 2 in January (0.3% of total sea ice area). We repeat this analysis on sea ice volume per unit area in Fig. 5(middle row). We again find little difference throughout the year, with a global maximum volume difference of 33 km 3 in January (0.12% of total 220 sea ice volume). We also examine the log ratio of SP-WIFF to NN-WIFF mean floe size (representative radius) as the bottom row in Fig. 5. Some areas have larger floe sizes in SP-WIFF compared to NN-WIFF -those where the discrepancy in mean floe sizes is a multiple of two or more, i.e. | log 10 (R N N /R SP )| > log 10 (2) have a maximum global area of 1100 km 2 in April, or 4.8% of total sea ice area. These regions are visible as locations lying at the boundary between the MIZ and sea ice pack, i.e. locations where NN-WIFF does not trigger wave fracture but SP-WIFF does. We evaluate these differences between model 225 integrations using the CICE quality-control check (Roberts et al., 2018). Although there is not bit-for-bit reproducibility, we find that the different implementations of WIFF are not "climate changing" according to the two-stage paired thickness test.
This demonstrates that differences in WIFF implementation do not have an emergent effect on sea ice model state.
NN-WIFF substantially accelerates the parameterization of wave-induced floe fracture. Integrating these simulations using on the Cheyenne supercomputer with daily I/O on a 384x320 displaced pole grid, using 1 node and 24 cores, the standard SP-

230
WIFF executes 1 year in approximately 24 hours. The single-iteration version SP-WIFF-1 takes 3h30. The NN-WIFF model takes 1h10 minutes. For comparison, a standard CICE integration with the FSD but no waves takes 56 minutes, and a CICE integration with no FSD and no waves takes 45 minutes. (Note that these runtimes do not include integration of a wave model, as this was run offline and provided as a forcing to the sea ice model.) Thus while execution times will vary with computer architecture and user-specified compiler settings, NN-WIFF reduces the overhead associated with simulating wave-ice fracture 235 without significant added computational cost.

Conclusions
Here we have presented WIFF1.0, code for performing wave-induced floe fracture in sea ice models that are coupled to an ocean surface wave model. While the full "super-parameterized" version of WIFF (SP-WIFF) is capable of simulating waveaffected marginal ice zones with active wave-ice interactions, it results in high computational expense which renders it difficult 240 to use in long climate integrations. Instead, we trained a pair of neural networks (NN-WIFF) to accelerate this parameterization using a large set of true climate model input/output pairs, which results in a median "standardized size error" of 3.9% when evaluated on a set of 800,000 input states, and 8.4% in free-running climate model simulations.
Because of the use of a classification layer, NN-WIFF does introduce some false positives and negatives -with a classifier accuracy of 87.5%. Yet the "overall" accuracy of calls to NN-WIFF is 96.5%. Because we now employ an energy threshold 245 for determining whether NN-WIFF would be called, the classifier network is only called on a reduced subset of sea ice points.
The energy threshold eliminates potential calls to NN-WIFF for approximately 72% of total sea ice points, thus the likelihood that NN-WIFF is inappropriately called at any given sea ice point is 0.125 × 0.28 = 3.5%.
Compared to a single-iteration version of SP-WIFF that was previously employed for computational reasons, NN-WIFF produces improved error characteristics, notably reduced sum and "representative" absolute errors, and is significantly less The network has been trained on present-day (2009) sea ice conditions from a single model, and therefore may have less success for different climates where the phase space of wave spectra and ice thicknesses changes. We attempted to use a broad 255 range of input model vectors that represent a full seasonal cycle of sea ice and wave statistics in both hemispheres. The appeal of parameterization acceleration is that new versions of NN-WIFF can be trained by directly using other climate model output, and model input variables and thresholds can be adjusted to provide better fits to this model. We provide code in the WIFF1.0 repository for running SP-WIFF in a standalone fashion, collating training data from a climate model, and (re)-training the network, as well as producing the plots used in this paper.

260
Because of the ease of obtaining training data from climate model output, this parameterization acceleration approach has, and could continue to, find a use in many other aspects of climate model simulations (Rasp et al., 2018), for similarly expensive parameterizations, especially those that are not the solution to primitive equations Pal et al. (2019); Brenowitz et al. (2020).
For example, sea ice rheological parameterizations generally rely on complex implicit solvers run to a specified tolerance, but may, much like SP-WIFF, be able to be decomposed into a series of neural-net "black boxes".

265
Our primary goal in developing NN-WIFF is to provide a method for accelerating the existing SP-WIFF paramterization, which will permit the implementation of wave-ice fracture parameterizations in climate-scale sea ice models. In general, NN-WIFF allows for a cost-effective implementation of wave-ice fracture in sea ice models with no significant increase in runtime.
Because a neural network is a simple set of elementary functions, it can be employed straightforwardly in any climate model that desires to simulate the fracture of sea ice by ocean surface waves.

270
Code and data availability. Code for CICE6 and Icepack is available at https://github.com/lettie-roach/CICE/tree/mlwave and https://github. com/lettie-roach/Icepack/tree/mlwave, and will be given a public DOI upon paper acceptance and incorporation into the core CICE code.