The sea ice model component of HadGEM3-GC3.1

. A new sea ice conﬁguration, GSI8.1, is implemented in the Met Ofﬁce global coupled conﬁguration HadGEM3-GC3.1 which will be used for all CMIP6 (Coupled Model Intercomparison Project Phase 6) simulations. The inclusion of multi-layer thermodynamics has required a semi-implicit coupling scheme between atmosphere and sea ice to ensure the stability of the solver. Here we describe the sea ice model component and show that the Arctic thickness and extent compare well with observationally based data.


Introduction
HadGEM3-GC3.1 is the global coupled model configuration to be submitted for physical model simulations to CMIP6 (Coupled Model Intercomparison Project Phase 6) by the Met Office Hadley Centre (Williams et al., 2017).It is comprised of global atmosphere, GA7.1, and land surface GL7.0 components (Walters et al., 2017), coupled to global ocean GO6 (Storkey et al., 2018) and global sea ice GSI8.1 components.This paper describes the global sea ice component (GSI8.1)which is embedded in the NEMO (Nucleus for European Modelling of the Ocean; Madec, 2008) ocean configuration (GO6) and uses a tripolar grid, while the atmosphere model (GA7.1), a configuration of the Met Office Unified Model (MetUM), and land surface (GL7.0), a configuration of JULES (Joint UK Land Environment Simulator; Best et al., 2011), use a staggered latitude-longitude grid.The communication between the two components is through the OASIS-MCT coupler (Ocean Atmosphere Sea Ice Soil Model Coupling Toolkit; Valcke, 2006).

Model description
Relative to its predecessor GC2, GC3.1 has new modal aerosol and multi-layer snow schemes, a newly introduced multi-layer sea ice scheme (described below), and a number of parameterization changes in all the model components, including a set relating cloud and radiation and revision to the numerics of convection (Williams et al., 2017).The GSI8.1 global sea ice configuration builds on the previous version used in HadGEM3-GC2.0,GSI6 (Rae et al., 2015), and is based on version 5.1.2 of the Los Alamos sea ice model CICE (Hunke et al., 2015).The CICE model determines the spatial and temporal evolution of the ice thickness distribution (ITD) due to advection, thermodynamic growth and melt, and mechanical redistribution/ridging.At each model grid point the sub-grid-scale ITD is modelled by dividing the ice pack into five thickness categories, with an additional icefree category for open-water areas.The initial implementation of CICE within the HadGEM3 coupled climate model is described in Hewitt et al. (2011).The key differences between GSI6 and GSI8.1 are the replacement of zero-layer thermodynamics with a multi-layer scheme, the addition of prognostic melt ponds, and the coupling to the atmosphere on ice thickness categories.For GSI8.1 the ice-atmosphere coupling is undertaken by category with all thermodynamic fluxes (conduction, surface melt, and sublimation), as well as snow depth and melt pond fraction and depth, being cal-Published by Copernicus Publications on behalf of the European Geosciences Union.culated separately for each ice thickness category (ITC).Appendix A contains details of namelist options and parameters used in GSI8.1 and Appendix B the C pre-processing keys used to build the model.

Albedo scheme
The albedo scheme used in GSI8.1 is based on the scheme used in the CCSM3 model (see Hunke et al., 2015) and has separate albedos for visible (< 700 nm) and near-infrared (> 700 nm) wavelengths for both bare ice and snow.The scheme is described in Sect.3.6.2 of the CICE User's Manual (Hunke et al., 2015).The penetration of radiation into the ice, as described by Hunke et al. (2015), is not included here.For this reason, following Semtner (1976), a correction is applied to the surface albedo to account for scattering within the ice pack.
This configuration includes the impact of surface melt ponds on albedo as an addition to the CCSM3 albedo scheme.The melt pond area fraction, f p (n), and depth, h p (n), for ice in thickness category n, are calculated with the CICE topographic melt pond formulation (Flocco et al., 2010(Flocco et al., , 2012;;Hunke et al., 2015).Where the pond depth, h p (n), on ice of thickness category n is shallower than 4 mm, the ponds are assumed to have no impact on albedo, and the albedo, α pi (n), of such ponded ice is simply equal to that of bare ice, α i .Where the pond depth is greater than 20 cm, the underlying bare ice is assumed to have no impact, and the ponded ice albedo is assumed equal to that of the melt pond, α p .For ponds deeper than 4 mm but shallower than 20 cm, the underlying bare ice is assumed to have an impact on the total pond albedo, and the bare ice and melt pond albedos are combined linearly (Briegleb and Light, 2007): Because the impact of melt ponds on albedo has been included explicitly, the reduction in bare ice albedo with increasing temperature (Hunke et al., 2015), which was intended to account for melt pond formation, is not included.However, the reduction in snow albedo, α s (n), with increasing surface skin temperature, intended to take account of the lower albedo of melting snow, has been retained and takes the form where α c and α m are the albedos of cold and melting snow respectively, T m is the snow melting temperature (i.e., 0 • C), T (n) is the surface skin temperature of ice in thickness category n, and T c is the threshold temperature, below T m , at which surface melting starts to affect the snow albedo.The scheme calculates the total grid-box albedo, α(n), of ice in thickness category n for each of the two wavebands by combining the albedo, α pi (n), of the ponded fraction, calculated as described here, with the albedos of bare ice, α i , and snow, α s (n), weighted by the melt pond fraction, f p (n), and the snow fraction, f s (n): (3) The snow fraction, f s (n), representing surface inhomogeneity due to windblown snow, for category n is empirically parameterized via a calculation based on snow depth, h s (n): where h snowpatch is a length scale parameter (Hunke et al., 2015).Note that this is different from the parameterization used in the previous configuration, GSI6.0, described by Rae et al. (2015).

Thermodynamics
GSI8.1 is the first sea ice configuration of the Met Office model to use multi-layer thermodynamics.Previously, the sea ice model used the zero-layer formulation described in the appendix to Semtner (1976), in which surface temperature reacts instantaneously to surface forcing and conduction within the ice is uniform.In the new formulation, the sea ice has a heat capacity, which depends on the temperature and salinity, and hence conduction can vary in the vertical.The ice is divided into four vertical layers, each with its own temperature and prescribed salinity (a fixed salinity profile); an additional snow layer is permissible on top of the ice (Fig. 1).The thermodynamics scheme is very similar to that described by Bitz and Lipscomb (1999), present in CICE5.1.2, in which the diffusion equation with temperature-dependent coefficients is solved by the iteration of a tridiagonal matrix equation.However, it is modified as described by West et al. (2016), with surface exchange calculations, carried out, separately for each thickness category, in the Met Office surface exchange scheme, JULES.The use of JULES allows near surface temperature to evolve smoothly on the atmosphere time step (West et al., 2016), which is short compared with the coupling frequency.The modular structure of JULES allows a consistent treatment of surface exchange (vegetation canopies, snow, soils, and sea ice) throughout the model (Best et al., 2011).All parameters passed between JULES and CICE at each coupling step are shown in the table in Appendix C. The diffusion equation is forced from above by the conductive flux from the ice surface into the top-layer interior, which for each ITC is calculated by the surface exchange and passed to the ice model.The category top-layer temperature, thickness, and conductivity then become the bottom boundary conditions for the next iteration of the surface exchange.

Semi-implicit coupling
The OASIS-MCT coupler used within GC3.1 does not have the functionality to regrid variables with time-varying weights.Consequently, for the atmosphere-ice coupling the conductive flux F condtop , calculated by the surface exchange for a single-atmosphere (low-resolution) grid cell, is divided amongst the underlying ocean model (high-resolution) cells in proportion to grid cell area.This means that ocean cells with a low ice fraction receive too much energy, while cells with a high ice fraction receive too little.The problems resulting are twofold: an unphysical "inverse imprint" of ice fraction occurs in the spatial pattern of conductive flux (as shown in Fig. 2e), and in a large number of instances the CICE temperature solver is forced with exceptionally high local conductive fluxes -resulting in the iterative temperature solver failing to converge.
In order to render the coupling more physically realistic, and thereby increase the reliability of thermodynamic convergence, the coupling was made semi-implicit.The sea ice fraction is now passed by first-order conservative regridding to the atmosphere at a coupling instant, and this new sea ice fraction used within JULES to apportion F condtop to produce a "pseudo-local" conductive flux.This new flux is then passed to the ocean model in the normal way, where it is multiplied by the ice fraction field on the ocean grid to produce the grid-box mean field that is implemented over the ensuing time step (Fig. 2).The grid-box mean field has the favourable properties of both conserving energy and of restricting incoming atmospheric fluxes to be proportional to the ocean grid underlying the ice fraction which improves the convergence of the CICE temperature solver (see Fig. 2f).Coupled fields can be shown to be exactly equivalent to the physically desirable solution that would be produced if fluxes were divided amongst underlying ocean grid cells in proportion to ice area.The semi-implicit coupling was found to conserve energy to a similar order of magnitude to the previous, explicit coupling, with an average grid cell error of under 10 −4 W m −2 across the Arctic.
With the implementation of the semi-implicit coupling there remained two cases in which the CICE temperature solver, forced by the JULES conductive flux, would fail to converge.In the first case, convergence becomes very slow for thin (< 0.2 m), melting ice; the surface exchange scheme would occasionally calculate large conductive fluxes for which the solver failed to converge within the required 100 iterations.To deal with this problem, a maximum threshold of 1000 h I W m −3 (where h I is ice thickness in metres) is specified for the conductive flux; any surplus conductive flux above this value is repartitioned to the base of the ice, and added to the ice-to-ocean heat flux.The second issue is a consequence of the way in which the CICE thickness distribution interacts with the coupling method.At high latitudes it is common for a large number of ocean cells each to underlie partially a single atmospheric grid cell.In cold winter conditions, conducive to strong ice growth, the fraction of ice in the thinnest ITC, a 1 , can be very small.With cold atmospheric temperatures, surface flux and conduction through ice in this category are necessarily strongly upwards; in some cells, random effects, perhaps dynamical, would cause conduction to be stronger than in others and also lead to lower top-layer temperatures.However, stronger conduction also promotes ice growth, which reduces the fraction a 1 in the grid cell, as it gets promoted to a 2 , rendering its top-layer temperature less visible to the atmosphere.In a small number of cases this was found to cause runaway cooling, with toplayer temperature in isolated cells cooling to below −100 • C, forced by high negative conductive fluxes calculated by a surface exchange scheme that was seeing much higher grid-box mean ice temperatures.This problem was solved by linearly reducing conductive flux to zero as top-layer ice temperature fell from −60 to −100 • C, with the excess flux passed directly to the bottom of the ice (and therefore helping to grow more ice, the effect that would be expected to occur in reality).These two processes were found to direct an average of 0.4 and −0.2 W m −2 respectively to the ice base over the course of 1 year in the Arctic.

Dynamics
The standard elastic-viscous-plastic rheology (EVP) for ice dynamics in CICE is used here (Hunke et al., 2015).However, NEMO is on a C-grid and CICE on a B-grid.This is dealt with though simple interpolation from CICE to NEMO as described in Hewitt et al. (2011).The remapping is the transport/advection algorithm scheme and ridging schemes which are the default in CICE.

Model evaluation
An example of the sea ice evaluation provided here is the Arctic multiannual mean winter (December, January, and February -DJF) ice thickness as diagnosed by the model (Fig. 3), both in its CMIP6 configuration of GC3.1 and its previous stable version GC2 (Williams et al., 2015).The model present-day control is forced by greenhouse gases and aerosols from the year 2000 for 100 simulated years.The evaluation data is CryoSat-2 satellite thickness (Tilling et al., 2016) inferred from freeboard measurements from 2011 to 2015 along with the 1990-2010 mean thickness from Pan-  Arctic Ice Ocean Modelling and Assimilation System sea ice reanalysis (PIOMAS; Zhang and Rothrock, 2003;Schweiger et al., 2011), inferred through the assimilation of observed sea ice fraction.The CryoSat-2 thickness retrievals are included solely as a guide to the ice thickness distribution as the ice has thinned since 2000.Figure 3 also show the model sea ice extent, depicted by the 15 % ice concentration contour, compared with the 1990-2009 mean from the HadISST1.2(Hadley Centre Sea Ice and Sea Surface Temperature data set) sea ice analysis (Rayner et al., 2003).The PIOMAS and GC3.1 Arctic ice thicknesses are comparable in spatial pattern save for PIOMAS depicting a larger area of thick ice adjacent to north Greenland and the Canadian Archipelago and GC3.1 depicting thicker ice in the central Arctic.CryoSat-2 depicts thinner ice in the western Arctic.Both GC3.1 and CryoSat-2 show thicker ice along the east Greenland coast than PIOMAS.The DJF ice extent in GC2 is low compared with the PIOMAS analysis and CryoSat-2 data but consistent with the low ice thickness.The extent compares well with the HadISST analysis in GC3.1; however, the ice is overly extensive in the Greenland and Norwegian seas.This is likely because the deep Atlantic water, in the ORCA025 configuration, is predominately formed in the Labrador Sea and there is very little convection, contrary to observations (Pickart et al., 2003), in the Greenland and Irminger seas.As a consequence the waters off east Greenland are rather static and the surface waters cool resulting in excess sea ice.The winter sea ice extent simulated by GC3.1 is much closer to the HadISST observations than was the case for GC2 in the Bering-Chukchi, Barents, and Labrador seas.The Antarctic sea ice extent has improved considerably between GC2 and GC3.1 (Fig. 4), the difference being due to a substantial, although not complete, reduction in the Southern Ocean warm bias (see below).
Figure 5 shows the mean seasonal cycle of volume for the GC3.1 model compared with that from GC2.The veracity of the model seasonal cycle of sea ice volume informs us whether the annual energy budget to the ice is well balanced.Unfortunately there are few observational means to assess this, and so here we use the PIOMAS model as a reference in the Arctic and satellite estimates from ICESat (Ice, Cloud,and land Elevation Satellite) for the Antarctic (Kurtz and Markus, 2012).It can be seen that, in agreement with Fig. 3, the GC2 Arctic ice volume is low, and that it is in closer agreement with PIOMAS at GC3.1.Both models have near identical annual cycles suggesting that, under presentday forcing, the new sea ice physics has not substantially altered the seasonal energy balance.The estimates for the 2003-2008 Antarctic ice volume from ICESat are for a minimum of 3357 km 3 in summer to a maximum of 11 111 km 3 in winter (Fig. 5).The GC3.1 volume compares well with ICE-Sat in the summer (2735 km 3 ) but is a little higher in winter (17 087 km 3 ).The GC2 configuration had a warm bias in the Southern Ocean, principally caused by excess solar insulation due to low cloud reflectivity (Williams et al., 2017), which was melting the Antarctic ice.This bias has been considerably reduced in CG3.1 resulting in a substantial increase in the Antarctic sea ice volume.

Summary
The GSI8.1 sea ice configuration of the Met Office Hadley Centre CMIP6 coupled model HadGEM3-GC3.1 has a number of physical enhancements compared to the previous version GSI6, including the introduction of multilayer thermodynamics and an explicit representation of the radiative impact of melt ponds.A semi-implicit coupling scheme refines the transposition of atmospheric fluxes to the sea ice, improving the stability of the thermodynamic solver.The final GC3.1 namelist options and pre-processor keys (see Appendix A and Appendix B) produce ice thickness and extent that are in good agreement with analyses.

Figure 1 .
Figure 1.Schematic demonstrating the time evolution of ice temperature, following an increase in downward short wave surface flux, for (a) the GSI6 zero-layer ice thermodynamics scheme; (b) the GSI8 multilayer scheme.

Figure 2 .
Figure 2. Demonstrating the implementation and effect, of semi-implicit coupling within the CICE-JULES surface exchange scheme.Panel (a) shows an example atmospheric grid cell overlying four ice/ocean grid cells calculates a sea ice conductive flux representing 210 J energy; (b) the division of energy between cells with standard coupling, in proportion to grid cell area; (c) the division with semi-implicit coupling, in proportion to ice area.(d) shows a conductive flux field on the atmosphere grid; the resulting flux field on the sea ice grid is shown for (e) standard coupling and (f) semi-implicit coupling.

Figure 3 .
Figure 3. Mean winter (December, January, and February) Arctic sea ice thickness from the HadGEM3-GC3.1 (50-year mean from the year 2000 equilibrium simulation) (b, d), from the PIOMAS (1990-2009) model reanalysis (a), and inferred from the CryoSat-2 (2011-2015) sea ice freeboard measurements (c).The orange and black lines show the 15 % ice concentration contours for the model simulations and the HadISST1.2sea ice analysis respectively.

Figure 5 .
Figure 5.The HadGEM3 model annual cycle of sea ice volume in the Arctic (a) and Antarctic (b).Volume estimates from the PIOMAS model reanalysis are included for the Arctic and ICESat estimates of volume in the Antarctic (grey dashed lines).