Sea Ice Formation in a Coupled Climate Model Including Grease Ice

Sea ice formation processes occur on subgrid scales, and the detailed physics describing the processes are therefore not generally represented in climate models. One likely consequence of this is the premature closing of areas of open water in model simulations, which may result in a misrepresentation of heat and gas exchange between the ocean and atmosphere. This work demonstrates the implementation of a more realistic model of sea ice formation, introducing grease ice as a wind and oceanic stress‐dependent intermediary state between water and new sea ice. We use the fully coupled land‐atmosphere‐ocean‐sea ice model, HadGEM3‐GC3.1 and perform a three‐member ensemble with the new grease ice scheme from 1964 to 2013. Comparing our sea ice results with the existing ensemble without grease ice formation shows an increase in sea ice thickness and volume in the Arctic. In the Antarctic, including grease ice processes results in large local changes to both simulated sea ice concentration and thickness, but no change to the total area or volume.


Introduction
Large-scale climate models struggle to accurately calculate Arctic sea ice volume (Shu et al., 2015) and thickness (Langehaug et al., 2013;Stroeve et al., 2014) and to capture trends in Antarctic sea ice extent (Turner et al., 2013). Various processes represented in the models have been investigated to explain this, including natural variability (Zunz et al., 2013), winds (Holland & Kwok, 2012), and melting ice shelves (Pauling et al., 2017). The latest generation of climate models include a more detailed representation of sea ice processes , and here we build on those advances by implementing a more sophisticated representation of subgrid-scale sea ice formation processes in a fully coupled climate model that includes feedbacks between changes to the atmosphere, ocean, and sea ice. We calculate historical climate simulations using our modified model and compare the representation of the simulated sea ice to that simulated from the same model without our modifications.
An important mode of sea ice formation results from ocean supercooling. Where the temperature of ocean water is lower than its salinity-dependent freezing temperature, supercooling occurs and frazil crystals may form. These small buoyant ice crystals rise to the surface and freeze to the underside of existing sea ice, thickening it. This mechanism for sea ice thickening depends on the supercooling of water at the ocean-sea ice interface and is different to congelation growth, which is the downward growth of ice crystals into the underlying water (Maksym, 2012). While most climate models broadly agree on the total sea ice mass budget, the partitioning of growth between these two processes varies between models and is likely to change as the climate warms, particularly in the Arctic where the increased open water area will make frazil growth increasingly important (Keen et al., 2020). It is important that both processes are represented appropriately in models so that the simulated sea ice continues to be realistic as the partitioning evolves. In completely calm conditions, frazil crystals can form a continuous, thin, flexible layer of sea ice over the open water surface, known as nilas and observed, for example, by Smedsrud and Skogseth (2006) and Winsor and Björk (2000). If there is wind-or ocean-driven turbulence, the frazil crystals mix with surface waters to form a slushy mix referred to as grease ice.
Supercooling may arise in response to river drainage or the mixing of ocean masses with different salinities (Foldvik & Kvinge, 1974;Martin & Kauffman, 1981). More commonly, it follows from a buoyant freshwater flux at depth, for example, beneath Antarctic ice shelves (Lewis & Perkin, 1986, 1983, or in response to extreme atmospheric surface cooling, for example, in leads and polynyas, where frazil created at the surface is mixed downward by wind-generated turbulence (Morales Maqueda et al., 2004). Polynyas are holes in the sea ice (or areas where ice does not form), created primarily by either strong offshore winds or by the creation of "hot spots" driven by warm waters upwelling or, in the Arctic, by solar heating (Morales Maqueda et al., 2004), while leads are fractures in sea ice caused by internal stresses. It is also possible for snow to be blown from the sea ice surface into areas of open water, creating a slush that is distinguishable from a frazil-formed layer of grease ice only through isotope analysis in a laboratory (Smedsrud & Skogseth, 2006;Weeks, 2010).
At present, most climate models remove supercooling from the surface of the ocean by transforming the energy deficit to a volume of new sea ice using the latent heat of freezing for ice (e.g., The Los Alamos National Laboraory sea ice model, CICE, Hunke et al., 2015, which constitutes the sea ice component of several climate models). This means that no grease ice is created and new sea ice forms instantly in response to surface supercooling, regardless of whether conditions are calm or turbulent. In reality, if grease ice forms, then it may persist for several days before atmospheric cooling causes the water fraction to solidify to create sea ice (Smedsrud & Skogseth, 2006). Exposure to a warm atmosphere may cause the solid fraction to melt and the grease ice may reduce or disappear altogether without ever forming new sea ice. This could mean that sea ice in climate models forms too fast, and areas of open water may close more quickly in the models than is realistic. Where grease ice forms close to sea ice in the real world, the grease ice may be "herded" against the sea ice edge by atmospheric and oceanic stress, leading to an uneven grease ice thickness distribution, and sometimes leaving part of the water area free from grease ice (Martin & Kauffman, 1981;Skogseth et al., 2009;Smedsrud & Skogseth, 2006;Smedsrud, 2011). This herding effect cannot be represented without grease ice being represented in the model, and omitting it could contribute further to the premature freezing over of leads and polynyas, and result in new sea ice created in the model being too thin, with impacts for sea ice and for atmospheric and oceanic processes that are coupled to it. These processes include heat and gas exchange between the polar ocean and atmosphere, which generally cools the upper ocean and warms the lower atmosphere (Morales Maqueda et al., 2004), and is inhibited by sea ice cover (Stephens & Keeling, 2000). Polynyas and leads therefore directly affect ocean and atmosphere heat and carbon dioxide cycles, as well as impacting the planetary radiation budget by reducing the surface albedo as dark holes in the relatively reflective ice (Brandt et al., 2005). In the Antarctic, coastal polynyas are a major source of sea ice production as extreme atmospheric cooling and strong offshore winds drive supercooling and the formation of grease ice, which is driven away from the coast by the strong winds, solidifying into new sea ice (which is also transported by the wind), and exposing the polynya surface water to further cooling (Morales Maqueda et al., 2004). Appropriate representation of polynyas in coupled climate models is therefore important for a realistic representation of sea ice formation and for realistic realizations of the atmospheric and oceanic processes that are coupled to it.
The importance of representing biological processes in climate simulations is increasingly recognized (Duarte et al., 2017;A. A. Sellar et al., 2019), and moden climate models, known as Earth System Models (ESMs), include these, for example, A. A. Sellar et al. (2019). Open water areas within sea ice are important for plankton (Arrigo et al., 1999) and macrofauna (Stirling, 1997), making the appropriate representation of leads and polynyas increasingly important. Sea ice one of the largest known biomes on Earth (Dieckmann & Hellmer, 2010), and high-latitude biological processes impact on global marine ecosystems and on global climate through interactions between biota incorporated in the sea ice and the atmosphere and ocean (Melnikov et al., 2002;Vancoppenolle et al., 2013). Realistic representation of these processes in ESMs requires grease ice to be represented because, to be incorporated into sea ice, biota are scavenged from the ambient ocean water and first incorporated into grease ice (Gradinger & Ikävalko, 1998), which then freezes to become sea ice. A coupled climate model includes wider atmospheric and oceanic processes that are likely to largely determine the volume of sea ice produced in the model, and biases in these are likely to dominate over any biases in the detailed sea ice formation calculations. For example, overestimation of the Arctic sea ice mass budget has been attributed to mainly atmospheric processes (Keen et al., 2020). However, a more physically realistic representation of sea ice formation makes the location and rate of sea ice growth more realistic. At present, no model that accounts for feedbacks between sea ice and atmosphere and ocean includes grease ice. Few field observations are available for grease ice because of the logistical difficulties of reaching and working in areas where it forms. This paucity of data is partly why grease ice processes are generally not represented in large-scale global climate models. Another reason is the computational expense of including subgrid-scale processes in a relatively coarse global model. Despite these challenges, a parameterization has been proposed to represent grease ice processes within leads in large-scale models (Smedsrud, 2011). The method has been demonstrated for partially ice-covered cells only in a sea ice model (Wilchinsky et al., 2015), and in a coupled sea ice-ocean model (Smedsrud & Martin, 2015). These previous works suggest that including grease ice processes in climate simulations results in thicker sea ice and an altered spatial distribution. However, because the models were not fully coupled, the atmospheric forcings experienced by the grease ice were unaffected by any changes to the sea ice. Therefore, the atmospheric forcings may not have been consistent with the simulated sea ice, and feedbacks between different components of the climate system were not considered. Here we extend those works to include a representation of grease ice processes in the coupled land-atmosphere-ocean-sea ice model, HadGEM3-GC3.1. In contrast to previous works, atmospheric and oceanic feedbacks that result from changes to the sea ice, driven by the inclusion of grease ice processes, are accounted for. Further differences from previous studies are that grease ice is considered for grid cells that are both ice-free and partially ice-covered, and either grease ice or nilas may form depending on whether conditions are calm or turbulent. In addition to making sea ice formation more realistic in the model calculations, the changes implemented here facilitate future ESM development, outside the scope of this work, to account for biological processes in sea ice environments. We assess the effect of implementing the grease ice scheme on modeled sea ice concentration and thickness, using an ensemble of historical simulations from 1964 to 2013, and use data derived from observations for the latter part of the same period as reference where possible. We note that the inclusion of new physical processes in a climate model generally necessitates retuning of the model, since previous tuning is likely to have accounted somewhat for the missing physics (Hourdin et al., 2017;Flato et al., 2013;Notz, 2015). Therefore, despite making model calculations more physically realistic, including a new process rarely results in closer agreement between observations and the simulated climate without further development effort to address the tuning. Nonetheless, it is preferable for climate models to include a physically accurate representation of processes where possible because tuning is unlikely to allow for the full range of physically possible responses to changes in other parts of the coupled climate system.

Model Description
The new scheme is implemented in the coupled atmosphere-land-ocean-sea ice model, HadGEM3-GC3.1, the physical core of the UK and New Zealand Earth System Models Williams et al., 2017). For the ocean component, GO6, based on NEMO3.6 (Madec & NEMO team, 2016), see Storkey et al. (2018), and for the sea ice component, GSI8.1, based on CICE5.1 (Hunke et al., 2015), see Ridley et al. (2018). The atmosphere component is provided by the Unified Model, using the GA7.1 configuration, and the land component is the JULES model, configured as GL7.1 (Walters et al., 2019). The ORCA1 grid (nominally 1°resolution) was used for the sea ice and ocean, with 75 vertical ocean layers, and the atmosphere model was run at 1.875°by 1.25°resolution, with 85 vertical levels. Simulations were implemented on the global domain.
In the model, sea ice is assigned to one of five thickness categories, and may move to a different category as it thins or thickens. Sea ice belonging to different categories can coexist in the same grid cell, and the sum of the concentration for the different categories gives the total ice concentration for the cell (ice concentration is the fraction of the grid cell covered by sea ice). In the standard scheme, if the ocean surface temperature is below its salinity-dependent freezing temperature, new sea ice forms. The amount of supercooling is transformed to an equivalent ice volume using the latent heat of freezing for ice. If there is no open water in the cell, then the new ice volume freezes to the existing sea ice, proportioned between the different ice categories according to their relative concentrations in the cell. If there is open water in the cell, then the new ice volume forms a layer of new sea ice of uniform thickness over the open water portion of the cell, with a minimum thickness of 5 cm ("collection depth"), and a maximum thickness of 60 cm (the minimum thickness requirement means that it may only partially fill the open water part of the cell). If the volume of new ice is greater than can be accommodated in the open water part of the cell, then the water fraction is covered with 60 cm thick new sea ice, and the remaining new ice is distributed between the categories of existing sea ice.
An appropriate collection depth is important for appropriate partitioning of sea ice growth between frazil and congelation processes, which is likely to become more important as the climate warms and the partitioning changes (Keen et al., 2020). Congelation growth is generally stronger for thinner ice, and so the collection depth, which determines new sea ice thickness, largely determines the transition from frazil to congelation growth. In general, a thicker collection depth results in a greater proportion of frazil growth (Keen et al., 2020). An earlier study using a forced model found a similar increase in frazil-driven sea ice growth when the collection depth was raised from 5 to 30 cm and when a grease ice scheme was implemented (Wilchinsky et al., 2015), but it is not clear that the same behavior would occur in a coupled model. That study also showed that grease ice thickness, when coupled to wind and ocean properties and to sea ice thickness, is more variable than a fixed collection depth, making the resulting sea ice thickness distribution more realistic than that calculated using a fixed collection depth (Wilchinsky et al., 2015). Including grease ice with a thickness coupled to ocean, wind and sea ice thickness makes the parameterization more physically complete and future proofs the representation of sea ice formation in the model, ensuring it responds appropriately to changes elsewhere in the climate system. Representing grease ice in the model has the additional advantage of facilitating future ESM development to include biological processes, as described in section 1.

Grease Ice Scheme
The new scheme is outlined in Figure 1. The surface ocean freeze-melt potential is converted to an ice volume as in the standard model. If the cell is ice covered, then the new ice volume freezes to, and thickens, the existing sea ice as in the standard scheme. If there is any open water in the cell, then the magnitude of the combined wind and ocean stress is calculated. If the net stress is zero, then no grease ice forms, and the new ice volume constitutes new sea ice, which forms as nilas. If the stress magnitude is greater than zero, and there is open water present, then we implement the grease ice scheme. Under the scheme, the new ice volume is not immediately considered to be new sea ice. Instead, some of it constitutes a volume of frazil ice, which makes up the solid fraction of a layer of grease ice in the open water part of the cell, comprising 25% frazil and 75% sea water. This follows the convention set by previous model studies (Heorton et al., 2017;Smedsrud & Martin, 2015;Wilchinsky et al., 2015). If there is grease ice in the cell persisting from the previous time step, then this is added to the new grease ice volume. Note that, for the purposes of this work, "grease ice" is distinct from "sea ice" and does not contribute to the values for sea ice concentration or volume (unless/until it freezes to become new sea ice, at which point it is no longer considered to be grease ice). We continue to refer to the open water fraction of a grid cell as open water, regardless of whether the water contains grease ice or not, that is, the sea ice concentration and water concentration sum to unity. Note also that although the volume of grease ice is preserved between time steps, the concentration and thickness of the grease are recalculated each time step.
In most cases, supercooling giving rise to frazil formation is driven by atmospheric cooling, and so it may be reasonable to assume that frazil produced in a partially ice-covered cell is concentrated in the open water fraction of the cell. However, it is also possible for supercooling to result from mixing of waters with different salinities, for example, where ocean masses meet, or where rivers or ice shelf melt provide freshwater fluxes (Martin & Kauffman, 1981). In these cases, it is unrealistic for all the frazil created in a cell to be concentrated in the open water, and in some cases this assumption could lead to problems. For example, if sea ice concentration is high, then forcing all frazil to be concentrated in a relatively small area of open water could lead to the formation of an unrealistically thick grease ice layer. Therefore, in partially ice-covered cells, not all of the frazil produced from surface supercooling is used to create grease ice. Instead, grease ice is created from only a proportion of the total frazil that is equal to the cell open water concentration. The remainder of the frazil thickens the existing sea ice (for ice-free cells, all frazil produced in the cell becomes part of the grease ice). Ideally, the origin of the supercooling (and hence of the frazil) would be determined from other model parameters and used to determine whether the frazil should be concentrated in the open water or not. For example, a full mixed layer model could be used to create frazil crystals in the water column, as demonstrated in Wilchinsky et al. (2015); however, that mixed layer model was not compatible with a coupled ocean model. The scheme presented here represents an improvement over the standard configuration (where no grease ice forms at all) but may underestimate the volume of grease ice in many cases.

Grease Concentration and Thickness
For ice-free cells, the grease ice volume, V g , is distributed evenly over the cell, giving grease ice concentration C g = 1, and grease thickness H g = V g /C g (grease ice concentration is the fraction of the grid cell area covered by grease ice). Since the surface area covered by a grid cell varies at high latitudes in the ORCA grid, calculations in the sea ice model are carried out with respect to concentration (which is unitless) rather than area. Volume is calculated in the model as the product of concentration and depth, and therefore has units of meters (grid cell area is generally used in postprocessing and analysis to convert this to cubic meters).

Journal of Advances in Modeling Earth Systems
For partially ice-covered cells, all open water is assumed to represent leads in the sea ice. A lead-sea ice element is conceptualized as extending the full width of the cell, with length Y, made of sea ice length L i and lead length L l , see Figure 2a, where H i and H g are the sea ice and grease ice thicknesses, respectively. The grease ice has span L g , which may not equal L l if conditions are conducive to herding as described below. Sea ice in each thickness category that is present in the cell makes up the lead walls for a fraction of the lead length proportional to that category's relative concentration, see Figure 2b. We set the length of the lead-sea ice element, Y = 5 km, sea ice length, L i = YC i and lead length, L l = Y − L i , following Wilchinsky et al. (2015) (C i is ice concentration). This means that for a cell with C i = 0.9, L l = 500 m.

No Herding
If there is insufficient open water to create leads of at least 10 m length, that is, L l < 10, then herding does not occur (Heorton et al., 2017;Smedsrud & Skogseth, 2006) and the grease ice is spread in a layer of uniform thickness over the lead surface ( Figure 2a). If the grease layer is thicker than the sea ice for any part of the lead, then the grease ice thickness is reduced to match the sea ice for that lead section, H g = H i . The solid fraction of the grease ice that is thereby removed from the lead thickens the sea ice in this category, and the water fraction drains to the ocean.
Where C i is low, some grease ice may form at large distances from ice floes and is unlikely to all be herded against ice edges, or to all overflow onto ice floes (note that there is no distinction in the model between underflowing and overflowing). For cells with C i < 0.05, any grease ice therefore forms a uniformly thick layer in the open water part of the cell. The grease ice thickness is compared to the sea ice thickness for the different categories of existing ice, and if H g > H i , then some grease ice overflows onto the sea ice of that category. For each existing ice category, the volume of overflowed grease ice is V o g ¼C i ðH g − H i Þ, and the volume of grease ice in the water is updated: g . This treatment means that grease ice remaining in the open water may be thicker than the existing sea ice, but it is considered more realistic than piling grease ice formed over a large open water area onto small ice floes, potentially thickening them by an unrealistically high amount.

Herding
For cells with C i > 0.05 and L l > 10, grease ice may be subject to herding, that is, may be piled up against (and overflow onto) the sea ice by atmospheric and oceanic stress, forming a wedge shape ( Figure 2a) rather than being distributed in a layer of even thickness (Heorton et al., 2017;Smedsrud & Martin, 2015;Smedsrud, 2011;Wilchinsky et al., 2015). We follow Wilchinsky et al. (2015) and project the stress onto the leads, somewhat arbitrarily assuming all leads to be orientated at 30°to the stress direction (HadGEM3-GC3.1 contains

Journal of Advances in Modeling Earth Systems
no information on subgrid-scale lead orientation). Using the projected stress to implement the model proposed by Smedsrud (2011), we calculate the concentration and thickness of the herded grease ice in the lead, and the volume of any grease ice that overflows onto the sea ice, V o g for each lead part, that is, for each part of the lead that has walls corresponding to a specific ice thickness category (Figure 2b).
Assuming the thick end of the grease ice wedge has thickness H i , Equation 1 determines the maximum possible span of grease ice that the lead can accommodate, L max g , from the stress, τ, and the granular resistance of the grease ice, k r (Figure 3). Note that this is the maximum span available for the grease ice to occupy, and the actual span (calculated later) may be smaller if there is insuffcient grease ice to fill this span. The granular resistance, k r , can be thought of as the resistance of the grease ice to a solid wall moving through it, with units Nm −3 . If the wall exerts force (per unit length of wall), F, over a grease ice depth, H g , then the wall will move with constant speed, that is, the resistive force from the grease ice will match F, if F¼k r H 2 g (Smedsrud, 2011). If the lead is not wide enough to accommodate L max g (i.e., L max g > L l ), then the wedge is truncated ( Figure 3b). In this case, we set L max g ¼L l , and calculate the thickness of the thin end of the wedge, H min g , from Equation 2. If the wedge is not truncated (Figure 3a), then H min g ¼0.
Having defined the wedge shape corresponding to the largest grease ice span allowed by τ, H i , k r and L l , we follow Wilchinsky et al. (2015) and use Equation 3 to calculate the corresponding grease ice volume, V max g . This is the maximum grease ice volume that can be accommodated in the lead without overflowing.
If V max g is greater than the actual volume of grease ice, V g , then all the grease ice can be accommodated in the lead and there is no overflowing. Continuing to follow Wilchinsky et al. (2015), the actual span of the grease ice, L g , is then given by Equation 4. . Cross section of the wedge shape occupied by herded grease ice in a lead with the maximum possible grease ice span, L max g : (a) when the maximum span can be accommodated in the lead; (b) when the wedge shape is truncated so that the maximum grease ice span can be contained within the lead.

Journal of Advances in Modeling Earth Systems
If V max g is less than V g , then the excess grease ice volume overflows onto the sea ice, V o g ¼V g − V max g , and the volume of grease ice remaining in the lead is updated: V g ¼V max g . The solid part of the overflowed grease ice thickens the existing ice, and the water part drains to the ocean. Equations 1-4 are carried out separately for each lead part (lead parts are defined by different values of H i , e.g., for the parts of the lead associated with different ice thickness categories in Figure 2b). The grease ice concentration, C n g , and thickness, H n g , are now determined for each lead part, n, from Equations 5 and 6, where V g is now the updated grease ice volume. The total grease ice concentration is C g ¼∑ n C n g , where C n i is the concentration of sea ice in the thickness category corresponding to lead part n, so C i ¼∑ n C n i .

New Sea Ice
Once grease ice concentration and thickness have been calculated, the atmosphere to ocean heat flux, Q a → o , determines whether any, or all, of the grease ice melts back into the ocean, freezes to become new sea ice (depending on the sign of Q a → o ), or persists as grease ice to the next model time step.
The latent heat of freezing is used to calculate the volume of water that can be frozen by Q a → o (or the volume of ice that can be melted). This is converted to the equivalent grease ice volume, accounting for the fact that only the water fraction of the grease ice can freeze, and only the solid fraction can melt.
The concentration and thickness of any new sea ice is then determined from the concentration and thickness of grease ice, and by the magnitude of Q a→o , following Wilchinsky et al. (2015). The latent heat associated with this freeze (or melt) is added to (or subtracted from) the ice to ocean heat flux that is returned from the sea ice model to the ocean model. For cells with low ice concentration (C i < 0.05), any freezing or melting occurs from the surface downward, that is, the surface of the grease ice layer freezes (or melts) first, and the volume of grease ice to be frozen (or melted) determines the depth to which freezing (or melting) occurs. The concentration of the new sea ice is then the grease ice concentration, C g , and the thickness of the new sea ice is the depth to which the grease ice froze. Conversely, for cells with C i > 0.05, we assume that the grease ice occupies leads in the sea ice, and freezing and melting occur laterally at the lead walls, that is, the full depth of the grease ice layer freezes (or melts), and the volume of grease ice to be frozen (or melted) determines the concentration of grease ice that freezes (or melts). In this case, the thickness of the new sea ice is the thickness of the grease ice, and the concentration is the concentration of the grease ice that froze. The grease ice volume is reduced by the volume of grease ice that has melted or frozen into new sea ice. If not all of the grease ice has frozen or melted, then the remainder persists to the next model time step where the solid fraction is added to any new ice volume created from the surface ocean freeze-melt potential (Figure 1).

Transport
When sea ice in a cell is transported in the dynamics part of the sea ice model, any grease ice that remains in the cell after the freeze-melt steps above is transported with it as a passive tracer (grease ice only exists in the sea ice component of the model, not the ocean component). To avoid grease ice in ice-free cells remaining static, cells containing grease ice are required to have C i > 0.00005 in at least one sea ice category (in both the standard HadGEM3-GC3.1 model, and in our modified version that includes grease ice, sea ice is melted from categories with concentration lower than this prior to the dynamics calculations to avoid numerical instabilities). For cells where the sea ice concentration in every thickness category is too low, any sea ice volume in the thinnest sea ice category (which may be up to 60 cm thick), is "spread out" over the cell in an attempt to achieve C i > 0.00005. If this does not result in a layer of sea ice that is at least 20 cm thick, then the solid fraction of some of the grease ice is considered to be new sea ice and removed from the persisting 10.1029/2020MS002103

Journal of Advances in Modeling Earth Systems
grease ice volume (note that in standard HadGEM3-GC3.1 model simulations, the solid fraction of the whole grease ice volume would be considered new sea ice).

Caveats to the Grease Ice Scheme
This first attempt to represent grease ice in a fully coupled climate model makes the model more physically representative. More comprehensive observations of grease ice properties would allow some necessary simplifications to be addressed. For example, grease ice in our model does not alter the surface roughness or radiative properties of the ocean surface, despite these being different for water-and grease ice-covered surfaces.
We assume a fixed solid fraction of 25% for grease ice, which is within the range of reported values observed in situ and in laboratory experiments (Martin & Kauffman, 1981;Maus & De La Rosa, 2012;Smedsrud & Skogseth, 2006;Winsor & Björk, 2000). In reality, however, the solid fraction is likely to increase as grease ice solidifies into sea ice, as observed in Smedsrud and Skogseth (2006) and described by Maus and De La Rosa (2012); however, more observations are needed to define or parameterize a globally realistic rate for the increase. In our implementation of the grease ice scheme, brine rejection is only associated with the formation of new sea ice, and there is no change to ocean salinity when grease ice forms or melts. In reality, the gradual release of brine as grease ice solidifies into sea ice can result in a more gradual salinification of ocean surface waters (Skogseth et al., 2009), which may have implications for local hydrography in some places (Smedsrud & Skogseth, 2006). It is also likely that the water content of grease ice is more saline than the ambient ocean water (Heorton et al., 2017;Smedsrud & Skogseth, 2006), and should therefore be associated with a lower freezing temperature. A range of salinities have been observed for grease ice, for example, Smedsrud and Skogseth (2006), making it difficult to define an appropriate deviation from the ambient salinity. We therefore neglect this and assume the salinity-dependent freezing temperature of water in grease ice to be that for the ambient ocean water. This means that the freezing of water within the grease ice may be associated with a slightly smaller energy change than is realistic. This first implementation of grease ice processes in a coupled model could constitute a framework for a more sophisticated representation of grease ice that allows the solid fraction to evolve, gradually releasing brine and adjusting the ambient ocean salinity. Further field and/or laboratory studies would be needed to constrain these processes and the development is beyond the scope of this work, but is facilitated by it.
Where grease ice forms in partially ice-covered cells and is subject to herding, a value for its granular resistance is required, k r in Equation 1. This is a function of the internal friction angle (Lambe & Whitman, 1979), and the grease ice bulk porosity (Dai et al., 2004), which are not well known. We set k r = 866 Nm −3 , following Wilchinsky et al. (2015). Sensitivity tests in that study showed that higher values of k r result in less herding, which may mean that leads freeze over faster and newly formed sea ice is thinner. A similar sensitivity was shown for the assumed lead orientation: smaller angles relative to the stress direction result in less herding since the drag stress perpendicular to the lead is reduced (Wilchinsky et al., 2015). However, the value for this orientation angle is necessarily arbitrary since the model contains no information on the orientation of subgrid-scale leads. The width of the lead-sea ice element, Y, in Figure 2 is set at 5 km, following Wilchinsky et al. (2015); however, Heorton et al. (2017) show that the degree to which grease ice is herded against the lead walls is sensitive to this, and Smedsrud and Martin (2015) suggest that the square root of the grid cell area may be more appropriate. At high latitudes, this would mean a different value for different cells, with particularly large differences at the latitudes where sea ice advances and retreats each year. Setting this value too low, or assuming an inappropriately low angle for the lead orientation has an effect equivalent to that which results from setting k r too high (Wilchinsky et al., 2015), that is, less herding may result than is realistic. It is hoped that future observation-based studies will allow k r to be assumed with greater confidence and provide evidence for a statistical distribution of lead orientations.

Data to Assess the Impact of the Grease Ice Scheme
Changes were made to the ocean and sea ice components of HadGEM3-GC3.1 to implement the new scheme for global coupled land-ocean-sea ice-atmosphere simulations using historical forcings from January 1964 to December 2013 inclusive. There was no discernible impact on computation time. Sea ice area and thickness have a high degree of natural variability that is represented by high internal variability in a fully coupled 10.1029/2020MS002103

Journal of Advances in Modeling Earth Systems
climate model (Notz, 2015). We therefore use a three-member ensemble of simulations with grease ice included (GREASE), for comparison against an equivalent ensemble using the standard sea ice formation scheme (STANDARD). The simulations in STANDARD formed part of the UK submission to CMIP6 (Eyring et al., 2016), and the GREASE simulations were branched from the STANDARD simulations in 1964. STANDARD simulations were branched from the preindustrial control simulation (Menary et al., 2018) at points in that simulation when the climate was in a different state (A. Sellar et al., 2018), and STANDARD and GREASE use the same historical external forcings, described in Eyring et al. (2016). Fully coupled models have high internal variability and three ensemble members is considered the minimum required for an acceptable representation of likely climate in the CMIP6 experiments (Eyring et al., 2016). The spread of climate realizations simulated by the ensemble is likely to be greater than the impact of the grease ice scheme. The ensembles for GREASE and STANDARD are therefore likely to overlap but impacts from the grease ice scheme should be evident in any persistent difference between the ensemble means, as discussed in Eyring et al. (2016) and Flato et al. (2013).
Using historical forcings allows impacts from the scheme to be assessed in the context of data derived from observations. We present the difference between the sea ice area in STANDARD and GREASE alongside the total sea ice area derived from satellite-borne observations using two different algorithms: bootstrapping (Comiso, 2017) and the NASA Team algorithm (Cavalieri et al., 1996), referred to as BOOTSTRAP and NASATEAM, respectively. To assess any impact of the grease ice scheme on sea ice thickness and volume, we present data from the Pan-Arctic Ice Ocean Modeling and Assimilation System (PIOMAS) for the Arctic, which combines satellite-derived sea ice concentration and sea surface temperatures with modeling (Schweiger et al., 2011), and for the Antarctic, we use data from the Global Ice-Ocean Modeling and Assimilation System (GIOMAS), which combines satellite-derived sea ice concentration with modeling (Zhang & Rothrock, 2003). PIOMAS data agree well with some sea ice thickness observations in the Arctic (Stroeve et al., 2014), but a comparison study of different thickness data sets derived from remote observations, including PIOMAS, found all derived sea ice thickness data to be associated with reasonably high uncertainty (Wang et al., 2016). In particular, PIOMAS may underestimate the thickness of thick ice and underestimate the thickness of thin ice (Wang et al., 2016). A paucity of observations means that GIOMAS data have not been validated in the Antarctic to our knowledge, although they have been shown to agree reasonably well with observations in the Arctic (Zhang & Rothrock, 2003). Nonetheless, in the absence of spatially comprehensive Antarctic observations, GIOMAS provides a useful data set, derived partially from observations, for comparison with model results and has been used as such in other studies, for example, Shu et al. (2015).
Implementing the grease ice scheme does not affect the timing or magnitude of the seasonal cycle in total sea ice area (Figures 4 and 5). In the Arctic, the maximum and minimum occur in March and August, respectively, in agreement with the observation-derived data. In the Antarctic, the maximum and minimum occur in September and February, respectively, making the maximum a month later in GREASE and STANDARD

Journal of Advances in Modeling Earth Systems
than in the observation-derived data. For the purposes of this work, "winter" hereafter refers to March and September in the Northern and Southern Hemispheres, respectively, and "summer" refers to August and February in the Northern and Southern Hemispheres, respectively. The geographical areas marked in Figure 6 are referred to in the discussion of the spatial distribution of effects from the grease ice scheme in sections 5 and 6.

Impacts in the Arctic
In the Arctic, GREASE and STANDARD capture the sea ice minimum area well according to BOOTSTRAP, but overestimate the magnitude of the winter maximum area according to both BOOTSTRAP and NASATEAM (Figure 4a). In this work we restrict ourselves to a discussion of the impact of the grease ice scheme, using the observation-derived data sets for context, rather than discussing differences between the model and observation-derived data more widely.
The range of Arctic ice thicknesses simulated in GREASE is broader than in STANDARD and includes thicker ice in both summer and winter (Figures 7a and 7b). Herding of grease ice against the sea ice edge, and the lateral growth of new sea ice forming in leads, means that new sea ice in GREASE may be thicker than new sea ice forming in STANDARD. In GREASE, new sea ice forming in a partially ice-covered cell may be as thick as the existing sea ice, whereas in STANDARD, new sea ice has a uniform thickness of up to 60 cm, which is exceeded only once the grid cell has become completely ice covered. The PIOMAS thickness distributions do not have the bimodal shape of the STANDARD and GREASE distributions (Figures 7a and 7b). The two modes represent single-year and multiyear ice. The latter is broadened when the grease ice scheme is implemented because grease ice herded in leads against the edge of thick multiyear sea ice persists and consolidates into new sea ice with a thickness that may match the multiyear sea ice thickness. The single mode in PIOMAS may reflect an underestimation of thick ice thicknesses in PIOMAS in summer and winter, combined with an overestimation of thin ice thicknesses in winter, as suggested in Wang et al. (2016). There is a cold bias in the standard HadGEM3-GC3.1 historical simulations that leads to an overestimation of Arctic sea ice thickness , which is likely to mean that there is too much thick Arctic sea ice in STANDARD, and this bias increases in GREASE when the grease ice scheme is implemented.
The total Arctic sea ice volume is greater in GREASE than in STANDARD for most of the simulated period (Figures 8a and 8b). In conditions of nonnegligible oceanic and/or wind stress, the grease ice scheme introduces a delay to the formation of new sea ice, as grease ice is first created, and new sea ice may not form until time step(s) after the surface becomes supercooled (in contrast to STANDARD, where surface supercooling is transformed instantly to new sea ice). Also, new sea ice formed from frozen grease ice is thicker than the nilas that is formed in STANDARD, as discussed above, and is less likely to cover the open water fraction of a grid cell. Areas of open water therefore take longer to freeze over in GREASE than in STANDARD, leaving the ocean subject to increased atmospheric cooling and driving the production of an increased volume of sea ice in GREASE, relative to STANDARD. The total Arctic sea ice volume simulated in GREASE and STANDARD becomes more similar toward the end of the time series, reflecting the warming of the ocean and atmosphere in recent decades (a warmer ocean requires a greater degree of cooling in order to freeze, and the cooling provided by the atmosphere is reduced as the atmosphere warms).
There are some small local differences in sea ice concentration between GREASE and STANDARD ( Figure 9), with the grease ice scheme giving a slight decrease in some areas in winter and a slight increase in summer. The winter decrease occurs at locations where the STANDARD concentration is higher than the NASATEAM climatology, and so brings the model slightly closer to the observation-derived data (Figure 9d). The summer concentration increases in GREASE, however, occur mainly at the edge of the summer ice pack, where concentration in STANDARD is already higher than in the climatology (Figure 9c). Implementing the grease ice scheme therefore pushes the summer concentration in the model further from the climatology; however, Figure 4 shows that the NASATEAM algorithm, from which the climatology is derived, underestimates ice area relative to the bootstrap algorithm, as also shown in Comiso et al. (1997), and the climatology may therefore show too small an ice-covered area.
The effect on Arctic sea ice thickness is much greater than the effect on concentration, and ice simulated in GREASE is thicker than in STANDARD for most of the Arctic Ocean in both winter and summer (Figure 10), similar to results presented in Smedsrud and Martin (2015). An exception to this is the northern Barents Sea, where the winter ice is slightly thinner in GREASE than in STANDARD (Figure 10d). The thickening in GREASE enhances what is already a positive thickness bias in STANDARD, relative to PIOMAS. This could be attributable to the assumptions made in the grease ice scheme that determine the degree to which the grease ice is herded against lead edges, or may be attributable to partially compensating biases elsewhere in the model. Alternatively, Arctic sea ice thickness may be underestimated in the PIOMAS data, as suggested

Journal of Advances in Modeling Earth Systems
by the single mode in the thickness distribution, which does not differentiate between single-and multiyear ice (Figures 7a and 7b).
In both GREASE and STANDARD, new sea ice forms thermodynamically in the Arctic Ocean and along coastlines in the Arctic (Figure 11). It is transported to the edges of the Arctic Ocean to where it melts,

Journal of Advances in Modeling Earth Systems
this can be seen in Figure 12a, where positive areas show ice moving into a grid cell, increasing the concentration in the cell, while negative values indicate either ice divergence (i.e., ice leaving a grid cell) or ridging within a grid cell (reducing the total concentration but leaving the ice volume unchanged). The effect of the grease ice scheme on these processes in the Arctic is very small and is mostly confined to the edges of the sea ice pack. This is not surprising since Figures 11 and 12 reflect

Journal of Advances in Modeling Earth Systems
changes in sea ice concentration, which Figure 9 shows to be only slightly affected by the grease ice scheme in the Arctic.
The widespread thickness changes in response to the grease ice scheme are greater than may be expected from the increased cooling that follows the increased open water area ( Figure 9d shows only a slight increase in open water area). It is possible that leads remain open for longer in GREASE, leading to production of an increased ice volume as discussed above, but that when they freeze over the ice that fills them is thicker, meaning that the internal stresses are less likely to result in new leads (HadGEM3-GC3.1 uses the elastic-viscous-plastic ice rheology from Hunke (2001) to parameterize internal stresses resulting from sea ice deformation rather than explicitly representing ice fracture). The sea ice concentration may therefore remain largely unchanged, but leads in GREASE may persist for longer and occur less frequently. A smaller number of longer-lived leads may result in more ocean cooling than a greater number of short-lived leads if the longer opening time allows convection to develop in the exposed near-surface waters. As cooled (denser) surface water sinks, it drives an upwelling of warmer water, allowing the ocean to lose more heat to the atmosphere, resulting in greater ocean cooling than occurs if just the exposed surface water cools and freezes. The increased cooling drives increased frazil production, and therefore increased sea ice production. This effect can occur even for relatively shallow convection depths. Figure 13c shows some deepening of the winter mixed layer under the pack ice in the central Arctic Ocean in GREASE, relative to STANDARD, which Hatching marks areas not significant at the 95% confidence level following a Student's t test. Note that only sea ice is included here, not grease ice, and the ring-like features are numerical artifacts from the model advection scheme.
suggests increased convection and so supports this explanation (note that the strong deepening of the mixed layer in the Barents Sea is not statistically significant). There is a shallowing of the winter mixed layer south of Iceland in GREASE, indicating increased stratification driven by the higher volumes of meltwater exported out of the Arctic through the Denmark Strait, following the greater volume of sea ice, relative to STANDARD. There is also a deepening of the winter mixed layer on the northeastern side of the Greenland-Scotland ridge and a shallowing on the southwestern side in GREASE, relative to STANDARD (Figure 13c). Atmospheric cooling creates dense water that sinks on the northeastern side of the ridge, and then flows south, rising to cross the ridge before sinking below lighter, warmer water carried northward by the North Atlantic Current. Increases in convection on the north side of the ridge, and in stratification on the south side, may indicate increased atmospheric surface cooling on the north side of the ridge, but the mixed layer is already reasonably deep here in STANDARD, making the anomaly relatively small. It is difficult to attribute this effect to changes in model sea ice formation processes and further investigation is beyond the scope of this study, but may be worthwhile for future works.

Impacts in the Antarctic
Both GREASE and STANDARD underestimate the maximum Antarctic sea ice area, relative to BOOTSTRAP, but agree well with NASATEAM, although as noted earlier the simulated maximum Hatching marks areas not significant at the 95% confidence level following a Student's t test. Note that only sea ice is included here, not grease ice, and the ring-like features are numerical artifacts from the model advection scheme.
occurs around a month later in GREASE and STANDARD (Figure 4b). The trend in simulated maximum Antarctic sea ice area is largely unaffected by the implementation of the grease ice scheme, and both GREASE and STANDARD overestimate the rate of decline relative to BOOTSTRAP and NASATEAM (Figures 5c and 5d). The distributions of sea ice thicknesses in GREASE and STANDARD are similar for both summer and winter, but include thicker ice than GIOMAS in summer (Figures 7c and 7d). The decreasing trend in summer and winter sea ice volume is also unaffected by the implementation of the grease ice scheme, and disagrees with the slightly increasing trend in both of these fields in GIOMAS (Figures 8c and 8d).
It may be anticipated that the impact of the grease ice scheme would be greater in the Antarctic (Smedsrud & Martin, 2015), where frazil ice accounts for a greater proportion of sea ice production than in the Arctic (Maksym, 2012). However, the grease ice scheme has little impact on the total Antarctic sea ice area, volume or the overall distribution of sea ice thicknesses (Figures 4b, 5c, 5d, 7c, 7d, 8c, and 8d), but it does create large local differences in Antarctic sea ice concentration and thickness between GREASE and STANDARD, particularly in winter (Figures 14 and 15).
In summer, differences in sea ice concentration between GREASE and STANDARD are small (Figure 14c), but sea ice around the Antarctic coast is generally thicker in GREASE (Figure 15c), as open water at the coast remains open for longer, exposed to increased atmospheric cooling which drives increased sea ice Hatching marks areas not significant at the 95% confidence level following a Student's t test.
production. In winter, sea ice concentration in the Amundsen Sea is much lower in GREASE than in STANDARD and is similarly decreased (although more weakly) everywhere around the northern sea ice edge, except in the Western Pacific, where the sea ice concentration is much higher in GREASE than in STANDARD (Figure 14d). There are also changes to winter sea ice thickness, with a large increase in the 10.1029/2020MS002103

Journal of Advances in Modeling Earth Systems
Western Pacific and around the Antarctic Peninsula, a decrease in the Amundsen Sea, and some smaller areas of decrease, for example, at the location associated with the Weddell Sea polynya (Figure 15d). Differences in thickness are generally weaker than in the Arctic because Antarctic sea ice is thinner, and so the maximum thickness for new sea ice in GREASE is also thinner, making it closer to the maximum thickness allowed in STANDARD (note that the scale in Figure 10 is different from that in Figure 15). Hatching marks areas not significant at the 95% confidence level following a Student's t test. 10.1029/2020MS002103

Journal of Advances in Modeling Earth Systems
Sea ice concentration increases thermodynamically in Antarctic coastal polynyas in both GREASE and STANDARD (Figure 11b), and the ice is then transported offshore by sea ice divergence (Figure 12b). Coastal polynyas are areas of low concentration and thickness in Figures 14a, 14b and 15a, 15b. In GREASE, grease ice production has replaced at least some of the sea ice that forms instantly in the polynyas in STANDARD, meaning that these areas of open water are likely to remain exposed to atmospheric cooling for longer, increasing sea ice production. Grease ice is transported by the same wind and ocean stresses that drive sea ice divergence in Figure 12b, so there may be a decrease in sea ice formation at the coast, and an increase slightly north of the coast, where the transported grease ice freezes. This happens in the Western Pacific, where coastal sea ice formation is reduced in GREASE, relative to STANDARD, creating a negative anomaly in Figure 11d. Sea ice divergence at the Western Pacific coast then also reduces, since grease ice is transported instead of sea ice, creating a positive anomaly in Figure 12d. North of these coastal anomalies, thermodynamic sea ice production increases in GREASE where the transported grease ice freezes (Figure 11d). The increase in sea ice concentration and thickness in the Western Pacific in GREASE, relative to STANDARD (Figures 14d and 15d) therefore follows from the enhanced surface cooling at the coastal polynyas, despite the increase being displaced from the coast. The increased volume of sea ice forming in this area leads to an increase in the sea ice divergence that transports ice to the northern sea ice edge in the Western Pacific (Figure 12d), resulting in increased melt, which drives a shallowing the surface mixed layer in the Western Pacific in Figure 13d in GREASE (note that Figure 11 shows changes in ice concentration, so the melting of equal areas of thick and thin ice appears the same).
Similar processes explain the reduction in sea ice concentration and thickness in the Amundsen Sea in GREASE, relative to STANDARD (Figures 14d and 15d). There is a decrease in thermodynamic sea ice production toward the northern ice edge here in GREASE, relative to STANDARD, where the surface mixed layer depth is high in both STANDARD and GREASE (Figures 13b and 13d), indicating convection. In areas of convection, warmer water rises to the surface where it cools and sinks, driving an overturning and often maintaining a polynya within the sea ice cover. In STANDARD, the supercooled surface water is transformed to sea ice, but in GREASE, grease ice is produced instead. The negative anomaly just south of the Amundsen Sea northern ice edge in Figure 11d shows that at least of some of the grease ice does not freeze at the production location. In STANDARD, divergence transports sea ice from the convective region to the northern ice edge where it melts. In GREASE, some grease ice is transported to the ice edge where it melts, without ever having frozen to form sea ice. This production, transport and melt of grease ice, rather than sea ice, creates the negative-positive anomaly pairs close to the northern ice edge in the Amundsen Sea in Figures 11d and 12d. The former shows a decrease in the production and melt of sea ice in GREASE, since grease ice is produced and melts instead, and the latter shows a reduction in sea ice divergence to the northern ice edge in GREASE, since grease ice is transported instead.
The production of grease ice in place of at least some of the sea ice that forms in STANDARD means that open water areas freeze over less readily in GREASE, enhancing atmospheric surface cooling and driving increased convection. This results in a deepening of the mixed layer in the Amundsen Sea in GREASE, relative to STANDARD (Figure 13d). Ordinarily, increased surface supercooling is associated with increased sea ice production. However, the proximity of this area to the northern ice edge means that if grease ice is produced instead of sea ice, then at least some of it is transported to the northern ice edge where it melts without ever having formed sea ice. This reduces the sea ice concentration and thickness in the outer Amundsen Sea in GREASE, relative to STANDARD (Figures 14 and 15). This reduction is roughly equal in magnitude to the increase in the Western Pacific following the implementation of the grease ice scheme, and we therefore do not see the same increase in total sea ice volume in the Antarctic that we see in the Arctic.

Summary and Concluding Remarks
We have demonstrated a framework whereby grease ice formation and grease ice herding processes can be represented in sea ice formation calculations in a fully coupled global climate model. Whereas in the standard sea ice formation scheme, sea ice forms instantly in response to ocean surface supercooling, it may take several model time steps (depending on the ocean to atmosphere heat flux, this could be several days) for new sea ice to form when the grease ice scheme is implemented. This, and the nonuniform thickness distribution of grease ice (following herding by the wind against sea ice edges), which may freeze to form a nonuniform distribution of sea ice, means that areas of open water persist for longer when the grease ice scheme is implemented, prolonging the ocean's exposure to atmospheric cooling and driving increased frazil ice production. For the three-member ensemble used here, this increased frazil production drives an increase in Arctic sea ice volume. In the standard sea ice formation scheme, frazil ice is considered to be sea ice, whereas in the new scheme presented here, it forms grease ice, which may be transported from the supercooling location before freezing to form new sea ice.
In both hemispheres, implementing the grease ice scheme results in some local redistribution of sea ice. In general, new sea ice in areas of partial ice cover is thicker when the grease ice scheme is implemented, following herding of the grease ice against the sea ice edge, and the lateral growth of new sea ice, which closes leads laterally (rather than forming a cap across the upper surface of the lead). This means that new sea ice may be as thick as any existing sea ice in partially ice-covered grid cells. This thickening effect is greater in the Arctic, where sea ice is generally thicker, than in the Antarctic, although the grease ice scheme does drive a thickening of summer sea ice in Antarctic coastal areas.
In the Antarctic, changes in winter concentration, and to a lesser extent thickness, are associated with the production of grease ice, rather than sea ice, in polynya regions. The increased surface cooling when the grease ice scheme is implemented drives an increase in both sea ice concentration and thickness in the Western Pacific, as grease ice is transported away from the areas of supercooling at the coast and freezes into the ice pack, leaving the polynya surfaces exposed to further cooling and further frazil production. In the Amundsen Sea, grease ice forms in an area of convection relatively close to the northern ice edge, and is transported northward to where it melts without ever having frozen to form new sea ice. Sea ice concentration and thickness therefore decrease in the Amundsen Sea when the grease ice scheme is implemented. These two regions dominate the sea ice response to the grease ice scheme in the Antarctic and are of roughly equal magnitude, leaving little net change to total Antarctic sea ice volume.
We have shown that the implementation of a more detailed sea ice formation scheme results in some changes to the spatial distribution of sea ice, particularly in the Antarctic winter, but minimal change to the total area in either summer or winter in either hemisphere. Including grease ice drives an increase in volume and thickness for simulated Arctic sea ice, and causes local changes (thinning and thickening) to simulated Antarctic sea ice. Sea ice volume represents latent heat and so is important to ocean energy balance, but is difficult to estimate from observations because it requires reliable and spatially comprehensive sea ice thickness measurements. It is important that model calculations of sea ice thickness take into account all physical processes that impact on it, so that the sea ice volume response to current and future changes in the ocean and atmosphere can be reliably estimated, despite a lack of observations. Our results concur with earlier works (Smedsrud & Martin, 2015;Wilchinsky et al., 2015) in finding model sea ice thickness and volume to be sensitive to the representation of sea ice formation, emphasizing the need for grease ice processes to be included in models.
The grease ice scheme presented here makes the representation of sea ice formation more physically plausible. This implementation contains some necessary assumptions which previous works have shown are likely to impact the results, and which could be better constrained if more observations of grease ice properties were available from field and laboratory studies. More observations of sea ice thickness are also needed, particularly in the Antarctic, to constrain model development work and to assess model biases in simulated sea ice volume.
While our knowledge of the precise magnitude of the increase in Arctic sea ice thickness is limited by the small ensemble size, its scale demonstrates that grease ice formation is a relevant and important process for climate simulations. Although in this work the thickening effect does not lead to an improved (more realistic) thickness because HadGEM3-GC3.1 overestimates historical Arctic sea ice, this is connected to the cold bias in the historical simulations and should not be interpreted to mean that the processes leading to increased ice growth are unrealistic. The state of sea ice in a climate model depends on the interaction of all model components and therefore the implementation of a new process, such as grease ice formation, generally requires further tuning steps for the other model components. The changes seen here to result from the inclusion of grease ice processes in the model, including increased thermodynamic growth in areas where there is high ice divergence and/or thick partial ice cover, local effects such as those seen in the Amundsen Sea, and greater differences between seasonal and multiyear ice thicknesses, provide a more realistic description of sea ice in those areas, and this can be used to inform appropriate tuning for other processes in the model.

Data Availability Statement
Model data produced and analyzed in this work are publicly available under Creative Commons Attribution-NonCommercial 4.0 International license at this site (Mackie et al., 2020).