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 are 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 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, Wave-Induced 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.

1 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 (FSD) as a core aspect of their model physics (Horvat and Tziperman, 2015; Zhang et al., 2015; Roach et al., 2018b; 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. The integral of f(r) over all positive r is equal to the sea ice concentration, c.

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 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, impacting up to half of the sea ice cover depending on the season and hemisphere (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 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., 2012, 2014). Floe fracture by waves is the dominant process driving changes in floe perimeter in summer months when the most lateral melting occurs (Roach et al., 2019).

In a climate model, coupled wave–ice feedbacks are related to two sub-grid-scale distributions: the FSD, f(r), and the ocean surface wave spectrum, S(λ), where λ is the wavelength and $\int S\left(\mathit{\lambda }\right)\mathrm{d}\mathit{\lambda }=E$ (units m²) is the wave energy per square meter, with $E=\mathrm{4}\sqrt{H_{\mathrm{s}}}$ for H_s the significant wave height (Michel, 1968). 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., 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 sea ice fracture. This method (hereafter called SP-WIFF) is broadly analogous to a “superparameterization” (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., 2018b; Hunke et al., 2019) and has been coupled to an ocean surface wave model (Roach et al., 2019). Running SP-WIFF incurs high computational costs – in the case of its incorporation in CICE, it increases computation times by an order of magnitude (see Sect. 4). As it is stochastic SP-WIFF is not bitwise reproducible: two identically initialized and forced simulations using SP-WIFF will not produce identical output up to the level of machine precision, which is often a necessary feature for climate model development.

To provide a computationally inexpensive, flexible, and bitwise reproducible sub-grid-scale parameterization of wave-induced ice fracture, we 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), 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. The model is trained using 5.1 million input and output vectors taken from coupled CICE–WAVEWATCH III simulations (see Sect. 3). Implemented in free-running coupled simulations in a year not corresponding to the training dataset, 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 %.

2 SP-WIFF: a superparameterized wave fracture scheme

The superparameterization of sea ice fracture by ocean waves was derived in Horvat and Tziperman (2015) (Sect. 2.3); its statistical properties were explored in Horvat and Tziperman (2017) (Sect. 3.4), and it was implemented with offline wave–ice interactions in a sea ice model in Roach et al. (2018b) (Sect. 2.4) and 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 in which ocean surface wave energetics are described by a discrete unidirectional wave spectrum:

$$\begin{array}{}\text{(1)}& S\left({\mathit{\lambda }}_i\right)\mathrm{d}{\mathit{\lambda }}_i=\underset{\mathrm{0}}{\overset{\mathrm{2}\mathit{\pi }}{\int }}S({\mathit{\lambda }}_i,\mathrm{\Theta })\mathrm{d}{\mathit{\lambda }}_i\mathrm{d}\mathrm{\Theta }.\end{array}$$

Considering only those floes with horizontal size between r and r+dr, a fraction of the domain $\frac{f\left(r\right)}{\mathit{\tau }}\mathrm{\Omega }(r,t)\mathrm{d}r\mathrm{d}t$ (unitless) is broken by ocean surface waves over a period dt. The parameter τ is a prescribed timescale over which the floe fracture takes place. In the Roach et al. (2018a) implementation, τ is determined via an adapting time-stepping algorithm to satisfy the CFL criteria for fast-propagating waves and small-area FSD categories, respectively (see Horvat and Tziperman, 2017, Appendix A, for further details). The fractured area has its own floe size distribution – the fraction of Ω(r,t)dr dt that now belongs to floes with size between s and s+ds is F(r,s) ds (unitless), with $\underset{\mathrm{0}}{\overset{\mathrm{\infty }}{\int }}F(r,s)\mathrm{d}s=\mathrm{1}$. Generically, the time rate of change of area of floes of size r due to fracture by ocean surface waves is

$$\begin{array}{}\text{(2)}& \begin{array}{rl}{\displaystyle \frac{\partial f(r,t)}{\partial t}}& ={\displaystyle \frac{\mathrm{1}}{\mathit{\tau }}}\left[-f(r,t)\mathrm{\Omega }(r,t)\right.\\ & \left.+\underset{r}{\overset{\mathrm{\infty }}{\int }}\mathrm{\Omega }(s,t)f(s,t)F(s,r,t)\mathrm{d}s\right].\end{array}\end{array}$$

Note that for ease of interpretation, notation in Eq. (2) differs from the analogous equation set in Horvat and Tziperman (2015) (their Eqs. 19–22). Equation (2) is one tendency term in the evolution of the FSD, which responds to multiple external forcings (e.g., thermodynamic growth and melting; advection) as described in Horvat and Tziperman (2015). External forcing terms are computed at each model time step, but the FSD f(r,t) is prognostically evolved using the above-referenced adaptive time-stepping scheme. Wave–ice coupling is performed at each sea ice model grid step, allowing for feedbacks between the FSD and, for example, floe-size-dependent wave attenuation schemes (Meylan et al., 2021).

The first term in Eq. (2) 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 cannot fracture into a larger floe, and therefore $F(s,r)=\mathrm{0}$ for s<r. As Eq. (2) is a generic tendency equation for fractured sea ice, SP-WIFF then refers to a parameterization of both Ω and F, which in each ice thickness category are evaluated as follows.

• S1.

  The discrete one-dimensional wave spectrum S(λ)dλ is converted to a one-dimensional ice strain field η(x) of (arbitrarily chosen) length 10 km (see Horvat and Tziperman, 2015, their Eqs. 19–22).

• 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 length-weighted percentage of fracture lengths with length between r and r+dr. Note we have replaced r⋅R(r) in Horvat and Tziperman (2015) with A(r) here to simplify notation.

• 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 $\mathrm{5}\times {\mathrm{10}}^{-\mathrm{4}}$.

Each wave field and corresponding fracture length collection (steps S1 and S2) is performed independently, i.e., on an unbroken floe of length 10 km. We next use the histogram A(r)dr to compute Ω and F. First,

$$\begin{array}{}\text{(3)}& \mathrm{\Omega }(r,t)=\underset{\mathrm{0}}{\overset{r}{\int }}A(s,t)\mathrm{d}s,\end{array}$$

and Ω(r,t) dr is equal to the length-weighted fraction of all fracture lengths smaller than r, which assuming a random horizontal distribution of floes is also equal to the fraction of floes of size r that will be broken. The resulting distribution of floe sizes is determined by the histogram itself:

$$\begin{array}{}\text{(4)}& \begin{array}{rl}& F(r,s,t)\mathrm{d}s=\\ & \left\{\begin{array}{ll}A(s,t)\mathrm{d}s/{\int }_{\mathrm{0}}^{r}A(l,t)\mathrm{d}l=A(s,t)\mathrm{d}s/\mathrm{\Omega }(r,t)& \text{if }r\ge s,\\ \mathrm{0}& \text{if }r<s.\end{array}\right.\end{array}\end{array}$$

The denominator in Eq. (4) ensures normalization of $F(r,s,t)$ over s, and therefore $F(r,s,t)\mathrm{d}s$ is the probability distribution of floe sizes from a broken floe of size r.

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. Note that this convergence is related to the stochastic histogram generation code, which converges to a steady-state distribution of fractures, and is not related to feedbacks from the waves to the ice. 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₁. Despite this simplification, CICE runtimes are increased by approximately a factor of 4 (see Sect. 4) over benchmark estimates without wave fracture. We seek to improve 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.

3 Accelerating the parameterization of coupled ice–wave interactions

The NN-WIFF code replaces the full superparameterization (SP-WIFF) with a parameterization that replicates its statistics but has reduced computational overhead. It consists of two full-connected feed-forward neural networks trained using the Python Keras deep learning library, one 100 node × 100 node 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. We introduced a classifier layer as SP-WIFF frequently returns unfractured sea ice in low-wave regimes, and we wish to train and run the histogram-generating network only when the sea ice will fracture. Recent observations of universal threshold behavior for sea ice breakup (Voermans et al., 2020) may provide the opportunity in future work to replace this classification layer with a simple physically based threshold for when the sea ice fractures.

Network architectures were chosen using a meta-learning approach of varying loss functions and numbers of nodes and layers (we provide code for altering network size in the WIFF1.0 release). Training is performed with the “adam” optimization scheme using the resilient backpropagation method. Rectified linear unit (RELU) activation functions are used for each of the hidden neurons, and we use a softmax activation for the final output layer in both cases as we seek a binary value (for the classifier) and a histogram that sums to 1 (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, i.e., SP-WIFF instead of SP-WIFF₁). We take 6-hourly wave spectra, sea ice concentration, sea ice thickness, and fracture histograms in both hemispheres, which results in a total of 180 million possible input vectors (most without sea ice). While the coupled wave–ice simulation includes the full distribution of sea ice floe sizes and thicknesses, only mean floe size, mean floe thickness, and sea ice concentration are passed to the wave module, as they are required to compute wave attenuation, and therefore we use these parameters to build NN-WIFF. 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_{\mathrm{s}}=\mathrm{4}\sqrt{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 wave energies. We find a total of 17.9 million sea ice points, 5.1 million of which meet the criterion above.

Thus, in total, the full training dataset 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, as well as 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_idr_i, 2.6 million of which are nonzero (i.e., 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 the following (see “Code and data availability” statement).

• Defining a custom training dataset using prescribed input spectra and thicknesses.

• Defining a custom output histogram size (i.e., different FSD categories).

• Designing and training the input classification scheme.

• Designing and training the histogram generation scheme.

We also provide data that were 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.

Figure 1 Performance of input classification for NN-WIFF. (a) Percentage of times the classifier runs (run rate, dashed black), calls the histogram network incorrectly (false positive rate, blue), or fails to call the histogram network (false negative rate, green) and overall accuracy (black) for the input classifier as a function of sea ice concentration. (b, c) Same as (a) for (b) sea ice thickness and (c) log ₁₀ wave energy. (d–f) Heat map of average classifier accuracy as a function of (d) ice concentration and log ₁₀ wave energy, (e) ice concentration and ice thickness, and (f) ice thickness and log ₁₀ wave energy. Regions of phase space with fewer than 75 input vectors (less than 0.005 % of input data) are uncolored or empty.

3.1 Input classification

The histogram A_idr_i is a collection of numbers that sums to 1 and is the target of our parameterization acceleration. This motivates the use of a softmax activation layer (as the output sums to 1). 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. We found the bulk of potential false positives could be eliminated by introducing the wave energy cutoff of $H_{\mathrm{s}}=\mathrm{4}\sqrt{E}=\mathrm{0.1}$ m above, and therefore the classification network is implemented after that cutoff is applied.

The training data for the classifier are the 5.1 million input vectors with corresponding binary values (fracture$/$no fracture). This is randomly partitioned $\mathrm{70}/\mathrm{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, $\mathrm{0}<\mathit{\chi }<\mathrm{1}$, which executes the NN-WIFF scheme if χ<χ_crit. The classification threshold χ_crit is determined by optimizing the error properties of the classifier when evaluated against the validation dataset. We optimize the total error rate, ϵ, where

$$\begin{array}{}\text{(5)}& \mathit{ϵ}={\displaystyle \frac{\text{false positives + false negatives}}{\text{total input vectors}}}.\end{array}$$

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 %.

Figure 1 shows error characteristics of the classifier evaluated on the validation dataset. In Fig. 1a–c, we show the relationship between the 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 ₁₀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-ϵ. Of interest for future model and scientific development is the peak in SP-WIFF calls for wave energies between 0.01 and 0.1 m² (significant wave heights around 70 cm).

Generally, the input classifier has high accuracy across the range of ice concentrations and thicknesses, with a flat distribution of false positives and false negatives. Lower accuracy occurs at the highest range of wave energies, at which both false positives and false negatives are higher than 20 %. We plot two-dimensional maps of the overall accuracy in Fig. 1d (log ₁₀E vs. ice concentration), Fig. 1e (ice thickness vs. ice concentration), and Fig. 1f (log ₁₀E vs. ice thickness). Some regions of phase space have classifier accuracies as low as 42 %. The classifier accuracy is high in the regions with the highest run rates – high sea ice concentrations, thicknesses, and wave energies between 0.01 and 0.1 m². As we explore below, these regions represent the largest portion of potential input states, and the 12.5 % error in classification does not contribute to a 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 retrain the classifier in WIFF1.0.

3.2 Neural network for generating fracture histograms

The fracture histogram network takes a length-27 input vector and outputs a vector ${\widehat{A}}_i\mathrm{d}{r}_i$, which sums to 1 and approximates the histogram A_i dr_i. As a result, sea ice floes will fracture if this input wave and the 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,

$$\begin{array}{}\text{(6)}& \mathrm{RSE}=\sum _i\left|{\widehat{A}}_i-A_i\right|r_i\mathrm{d}{r}_i,\end{array}$$

and has units of meters. We choose this error statistic over more standard metrics (root mean square error, for example) because mean floe size is an observable quantity (Horvat et al., 2020) and a positive moment of r, and it is therefore less sensitive to errors at small floe sizes or changes in floe size discretization. In the WIFF1.0 release we also provide code for training the neural network using other standard error metrics. Available training data are randomly partitioned into a training dataset of 70 % of the input vectors 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,

$$\begin{array}{}\text{(7)}& P={\displaystyle \frac{\mathrm{1}}{N}}\sum _{\mathrm{1}}^N{\mathrm{RSE}}_j\approx \mathrm{30}\text{m},\end{array}$$

where N is the total number of inputs used to train the network.

Figure 2 Performance of NN-WIFF on validation data. Histogram (a, left axis) of sea ice concentration, C, in the validation dataset. Median (a, right axis) standardized size error (SSE, Eq. 9) as a function of C (solid line). Dashed lines are the interquartile distance of SSE. (b–d) Same as (a), but for (b) sea ice thickness, (c) log ₁₀ wave energy, and (d) representative radius ($\stackrel{\mathrm{‾}}{R}$, Eq. 8) of the “true” output histogram. (e) Mean histogram output A_idr_i for (black) the validation dataset and (red) parameterization output ${\widehat{F}}_i\mathrm{d}{r}_i$. Dashed lines are the standard deviation of error between A_idr_i and ${\widehat{A}}_i\mathrm{d}{r}_i$ in each floe size category. Crosses indicate floe size bin centers.

Figure 3 Examination of training errors in NN-WIFF. (a) 2-D map of median SSE as a function of (x axis) sea ice concentration and (y axis) log ₁₀E from the validation dataset. (b–f) Same as (a) but for (b) H and log ₁₀E, (c) $\stackrel{\mathrm{‾}}{R}$ and log ₁₀E, (d) C and H, (e) H and log ₁₀E, and (f) $\stackrel{\mathrm{‾}}{R}$ and C. Regions are not colored for less than 0.005 % (∼ 40 points).

In Fig. 2 we display aggregated fit characteristics for the histogram neural net. Figure 2a–d display number histograms (red, left axis) of sea ice concentration (C), ice thickness (H), and log ₁₀E from the validation dataset, as well as of the “true” representative radius $\stackrel{\mathrm{‾}}{R}$, where

$$\begin{array}{}\text{(8)}& \stackrel{\mathrm{‾}}{R}=\sum _i{r}_i{A}_i\mathrm{d}{r}_i,\end{array}$$

from validation dataset input. To evaluate the standardized neural net error, we also compute the standardized size error (SSE) as

$$\begin{array}{}\text{(9)}& \text{SSE}={\displaystyle \frac{\text{RSE}}{\stackrel{\mathrm{‾}}{R}}},\end{array}$$

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 (Fig. 2a) or sea ice thickness (Fig. 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 median errors for both wave energy and $\stackrel{\mathrm{‾}}{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 $\stackrel{\mathrm{‾}}{R}$ values. Large values of $\stackrel{\mathrm{‾}}{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 ${\stackrel{\mathrm{‾}}{A}}_i$ (black line, Fig. 2e) accurately at each floe size category. Median error between ${\stackrel{\mathrm{‾}}{A}}_i\mathrm{d}{r}_i$ and $\stackrel{\mathrm{‾}}{{\widehat{A}}_i}\mathrm{d}{r}_i$ is less than 6.2 % in all categories, and across all 12 categories it averages 3.3 %.

To understand what regions of the phase space give rise to error, in Fig. 3a–c we plot two-dimensional heat maps of SSE, evaluated on the validation dataset, for each of the six combinations of ice thickness, concentration, $\stackrel{\mathrm{‾}}{R}$, and log wave energy. Colors are shown only for those coordinates with at least 0.005 % of the data (40 points), and the white point of the color map is set for an error of 10 %. Blue values are those with SSE < 10 %.

Figure 2 suggests that the input vectors with the largest error are those with low input wave energy and high resulting $\stackrel{\mathrm{‾}}{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 $\stackrel{\mathrm{‾}}{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 $\stackrel{\mathrm{‾}}{R}$ values, which have the highest SSE values. Although locations with low sea ice concentration are not used as 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.

4 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 we use output from the global sea ice model CICE6.1.4 (Hunke et al., 2019) in stand-alone mode forced by the JRA55 atmospheric reanalysis (JRA-55 and Japan Meteorological Agency, 2013) at a nominal 1^∘ spatial resolution. We include time-varying ocean surface wave spectra as sea ice model forcing, taken from output of the CICE–WAVEWATCH III coupled integration described in Roach et al. (2019). The model is run over the year 2005 to provide different boundary conditions and initialization from the 2009 data used as input to the network training. We examine results from three separate methods for obtaining the output histogram A_idr_i:

1. (SP-WIFF) using the full SP-WIFF model but the converged solution described in Sect. 2,

2. (SP-WIFF₁) the SP-WIFF model but using a single iteration of steps S1–S2 with fixed phase, and

3. (NN-WIFF) the trained neural network described above.

The histogram output A_idr_i from CICE is a daily average over 24 h time steps. When wave fracture does not occur, A_idr_i values are set to zero. Thus, though A_idr_i is normalized to 1 if wave fracture occurs, if a grid cell does not always fracture during one of the averaging time steps, the CICE history output will not be normalized to 1. We choose to include such points in this analysis, renormalizing histogram outputs to sum to 1.

Over the course of this year, we find that in the NN-WIFF simulations, there are 686 000 grid point 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 negatives as discussed in Sect. 3.1.

Figure 4 WIFF error characteristics in a climate model. (a) Mean Arctic histogram output A_idr_i in a free-running climate simulation using (black) SP-WIFF, (blue) NN-WIFF, and (red dashed) SP-WIFF₁. Points are included in the histogram if all three cause wave fracture. The interquartile range of SP-WIFF outputs is given as dashed black lines. (b) Same as (a) for the Antarctic. (c) Histogram of errors between SP-WIFF and (blue) NN-WIFF or (red) SP-WIFF-1 in the Arctic. Solid lines are the sum of absolute error (SAE). Dashed lines are the standardized size error (SSE). (d) Same as (c) for the Antarctic.

In Fig. 4a and b we repeat Fig. 2e but now examining (a) Arctic and (b) Antarctic averages of A_idr_i values for each of the three methods. To estimate differences in fracturing, Fig. 4 shows data only from locations and times at which WIFF leads to sea ice fracture in all three of the simulations and at which the daily average wave energy exceeds 0.1 m². Dashed lines show the interquartile range of SP-WIFF values. Median values are similar between all three implementations.

We plot histograms of the errors between SP-WIFF and SP-WIFF₁ (red) as well as NN-WIFF (blue) in both hemispheres as Fig. 4c–d from the same sample used to produce Fig. 4a and b. Dashed lines represent SSE and solid lines are the sum of absolute error (SAE) in each category. In both metrics and in both hemispheres, NN-WIFF performs as accurately as or more accurately than SP-WIFF₁ when executed and with significantly reduced computational cost. In the Arctic, median SSE is 8.9 % (NN-WIFF) vs. 10.5 % (SP-WIFF₁), and median SAE is 5.8 % (NN-WIFF) vs. 5.9 % (SP-WIFF₁). In the Antarctic, median SSEs are 8.1 % (NN-WIFF) vs. 13.4 % (SP-WIFF₁), and median SAEs are 5.52 % (NN-WIFF) vs. 9.02 % (SP-WIFF₁).

In Fig. 5 (top row) we show monthly average differences in sea ice concentration between SP-WIFF and NN-WIFF in March and September in each hemisphere. There is little difference between the two, with a maximum global area difference of 47 000 km² 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³ in January (0.12 % of total sea ice volume). We also examine the log ratio of SP-WIFF to NN-WIFF mean floe size (representative radius) in the bottom row in Fig. 5. Some areas have larger floe sizes in SP-WIFF compared to NN-WIFF – those with a discrepancy in mean floe sizes that is a multiple of 2 or more, i.e., $\left|{\log }_{\mathrm{10}}\right(R_{\mathrm{NN}}/R_{\mathrm{SP}}\left)\right|>{\log }_{\mathrm{10}}\left(\mathrm{2}\right)\approx \mathrm{0.3}$, have a maximum global area of 1100 km² 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 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, in spite of the impact of wave–ice interactions on sea ice state (Roach et al., 2018b).

NN-WIFF substantially accelerates the parameterization of wave-induced floe fracture. Integrating these simulations on the Cheyenne supercomputer with daily I/O on a 384×320 displaced pole grid using 1 node and 24 cores, the standard SP-WIFF executes 1 year in approximately 24 h. The single-iteration version SP-WIFF-1 takes 3 h 30 min. The NN-WIFF model takes 1 h 10 min. For comparison, a standard CICE integration with the FSD but no waves takes 56 min, and a CICE integration with no FSD and no waves takes 45 min. (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 significantly reduces the overhead associated with simulating wave–ice fracture.

Figure 5 Implementation of WIFF in a climate model. Difference (top row) in monthly sea ice concentration, SIC, between SP-WIFF and NN-WIFF in March and September in each hemisphere. The middle row is the same as the top row for sea ice volume, V. The bottom row is the same for the log ₁₀ ratio of the representative radii.

5 Conclusions

Here we have presented WIFF1.0, which is code for performing wave-induced floe fracture in sea ice models that are coupled to an ocean surface wave model. While the full “superparameterized” version of WIFF (SP-WIFF) is capable of simulating wave-affected marginal ice zones with active wave–ice interactions, it results in high computational expense, which renders it difficult 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 false positives and negatives – with a classifier accuracy of 87.5 %. With an energy threshold for determining whether NN-WIFF is to 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 $\mathrm{0.125}\times \mathrm{0.28}=\mathrm{3.5}$ %, and the “overall” accuracy of calls to NN-WIFF is 96.5 %. When compared to a single-iteration version of SP-WIFF that was previously employed for computational reasons, NN-WIFF produces equal or better error characteristics and is significantly less expensive. Global patterns of sea ice area and volume are statistically indistinguishable between NN-WIFF throughout a year of climate model simulations, with minor differences in representative floe sizes at the boundary between the MIZ and pack ice.

The network has been trained on present-day (2009) sea ice conditions from a single model and may therefore have less success for different climates in which the phase space of wave spectra and ice thicknesses changes. This may present challenges as model projections indicate that the future state of sea ice and wave climate will differ substantially from the present (for example Casas‐Prat and Wang, 2020; Roach et al., 2020; SIMIP Community, 2020), and the typical input sea ice and wave states to NN-WIFF may depart from the phase space examined here. In our training of NN-WIFF, we used a broad range of input model vectors, which sample both geographic and temporal variability in sea ice and wave conditions. While the input dataset spanned only a single year (2009), our assumption is that the future trajectory of ice and wave states lies contained within the geographic and temporal envelope of current model conditions. For example, while thinning Arctic sea ice may lead to larger percentages of the Arctic being affected by waves (Aksenov et al., 2017), interactions between high wave states and thin sea ice are presently simulated along the periphery of the Arctic Ocean and in the Southern Ocean and are therefore contained within the phase space of the training data examined here.

Still, it is possible that future inputs are not sampled as a part of the present-day conditions we examine here. Existing versions of NN-WIFF can be adaptively retrained using other updated climate model output to broaden the phase space of the training dataset, and we provide code for directly producing fracture histograms for any arbitrary wave and sea ice input vectors in WIFF1.0. NN-WIFF model input variables and thresholds may also be adjusted to increase the accuracy of the method, for example the inclusion of a simpler breaking threshold as described by Voermans et al. (2020). We intend in future versions of WIFF to implement a hybrid functionality that permits input vectors that lie outside the training phase space to be flagged and directly call SP-WIFF to permit better handling of outlier inputs. We note that the phase space of potential sea ice concentration, thickness, and wave spectral states simulated by the CICE simulations considered here may not reflect the real phase space of these variables nor that simulated by other climate model configurations. While we assume that CICE–WAVEWATCH III simulations largely mirror this particular set of potential input states, implementation of NN-WIFF in other modeling systems may require retraining to reflect a broader or narrower set of input conditions.

The model of Horvat and Tziperman (2015) was developed under assumptions about how sea ice breaks, namely that it is free-floating and responds instantaneously to curvature in the ocean surface height field. Recent developments in observing wave-induced breakup may point to the use of a single threshold for sea ice floe fracture (Voermans et al., 2020) that could replace the classifier trainer here and improve the representation of floe fracture in SP-WIFF. The aim of developing WIFF1.0 is simply replacing SP-WIFF with an approximating neural-network-based code. Future versions of WIFF could be developed in light of new observations (e.g., Voermans et al., 2019) and incorporate more detailed feedbacks between the FSD and wave models. The requirement of passing the full wave spectrum between wave and ice models may be too cumbersome for some current-generation coupling schemes, and it may be possible that fewer model parameters are needed to develop accurate fracture histograms (e.g., the abovementioned threshold). We provide code in the WIFF1.0 repository for running SP-WIFF in a stand-alone fashion, collating training data from a climate model, and (re)training the network, as well as producing the plots used in this paper.

Because of the ease of obtaining training data from climate model output, this parameterization acceleration approach has found 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). It may have applications in sea ice modeling as well. 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”. This methodology could also be used to accelerate the computation of wave–wave interaction and “source” terms in the WAVEWATCH III code.

Our primary goal in developing NN-WIFF is to provide a method for accelerating the existing SP-WIFF parameterization, 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 with the aim of simulating the fracture of sea ice by ocean surface waves.