The sea ice component of GC5: coupling SI³ to HadGEM3 using conductive fluxes

We present an overview of the UK's Global Sea Ice model configuration version 9 (GSI9), the sea ice component of the latest Met Office Global Coupled model, GC5. The GC5 configuration will, amongst other uses, form the physical basis for the HadGEM3 (Hadley Centre Global Environment Model version 3) climate model and UKESM2 (UK Earth System Model version 2) Earth system model that will provide the Met Office Hadley Centre/UK model contributions to CMIP7 (Coupled Model Intercomparison Project Phase 7). Although UK ocean model configurations have been developed for many years around the NEMO (Nucleus for European Modelling of the Ocean) ocean modelling framework, the GSI9 configuration is the first UK sea ice model configuration to use the new native NEMO sea ice model, SI³ (Sea Ice modelling Integrated Initiative). This replaces the CICE (Community Ice CodE) model used in previous configuration versions. In this paper we document the physical and technical options used within the GSI9 sea ice configuration. We provide details of the implementation of SI³ into the Met Office coupled model and the adaptations required to work with our “conductivity coupling” approach and provide a thorough description of the GC5 coupling methodology. A brief evaluation of sea ice simulated by the GC5 model is included, with results compared to observational references and a previous Global Coupled model version (GC3.1) used for CMIP6, to demonstrate the scientific credibility of the results.

1 Introduction

For over a decade, the Met Office has been developing global ocean model configurations based on the NEMO (Nucleus for European Modelling of the Ocean) ocean modelling framework (Madec et al., 2022). Standard UK configurations (see Guiavarc'h et al., 2024; Storkey et al., 2018; Ridley et al., 2018) are developed in collaboration with partners at the National Oceanography Centre (NOC), British Antarctic Survey (BAS), and the Centre for Polar Observation and Modelling (CPOM) as part of the UK's Joint Marine Modelling Programme (JMMP). These global configurations, amongst other purposes, provide the ocean and sea ice components of the Met Office Global Coupled (GC) model, which is used for modelling across a range of timescales, from short-range forecasting to centennial climate projections, as part of the Met Office seamless forecasting approach (Brown et al., 2012).

This paper is focussed on the latest version of the Met Office Global Coupled configuration, GC5 (Xavier et al., 2024), which will form the physical basis for Met Office Hadley Centre/UK model contributions to CMIP7 (Coupled Model Intercomparison Project Phase 7) with the HadGEM3 (Hadley Centre Global Environment Model version 3) physical climate model and UKESM2 (UK Earth System Model version 2). The GC5 coupled model is comprised of a combined Global Atmosphere and Land component (GAL9; Willett et al., 2024), coupled using the OASIS3-MCT coupler (Ocean Atmosphere Sea Ice Soil Model Coupling Toolkit; Craig et al., 2017) to a combined Global Ocean and Sea Ice component (GOSI9; Guiavarc'h et al., 2024). The atmosphere and land components of GC5 run on a staggered latitude–longitude grid using the Met Office Unified Model (MetUM, hereafter simply UM) and JULES (Joint UK Land Environment Simulator; Best et al., 2011) modelling systems respectively. The ocean and sea ice components run on a tripolar grid and are built around the NEMO ocean modelling framework. The sea ice component of GC5 is the new NEMO sea ice model, SI³ (Sea Ice modelling Integrated Initiative; Vancoppenolle et al., 2023). This is a change from previous model configurations, which used the CICE (Community Ice CodE) sea ice model (Hunke et al., 2015). The modelling systems used for the atmosphere, land, and ocean components of GC5, however, are still the UM, JULES, and NEMO respectively, consistent with previous GC model versions, e.g. the GC3.1 configuration used for CMIP6 (Williams et al., 2017).

This paper provides an in-depth description of the “GSI9” (Global Sea Ice version 9) sea ice model configuration, which forms the sea ice component of the GC5 model. A brief evaluation of the configuration is provided using sea ice model output simulated by the fully coupled GC5 system. Documentation and wider performance of the GC5 model meanwhile are discussed in the GC5 paper (Xavier et al., 2024), whilst a more detailed analysis of the ocean in the context of forced ocean–sea ice experiments can be found in the GOSI9 system description paper of Guiavarc'h et al. (2024).

As well as documenting the GSI9 sea ice model configuration, which is built around the SI³ model, this paper also describes the steps undertaken to incorporate SI³ into the Met Office coupled model – including adaption of SI³ to work with the “conductivity coupling” scheme used by the Met Office (West et al., 2016; Ridley et al., 2018). The coupling in the original implementation of HadGEM3 is thoroughly documented in Hewitt et al. (2011). However, with subsequent changes to the Global Coupled model made over more than a decade, some aspects require updating. Moreover, the change from CICE to SI³ means that many of the detailed schematics and processes described in Hewitt et al. (2011) are now out of date. We therefore provide a complete documentation of the coupling within GC5, with a particular focus on the sea ice exchanges.

This paper is organised as follows: Sect. 2 describes the GSI9 sea ice configuration used within GC5; Sect. 3 provides details on the integration of SI³ within the Met Office coupled model, including modification of SI³ to work with the Met Office conductivity coupling, and presents a detailed overview of sea ice coupling within GC5; Sect. 4 provides a brief evaluation of GC5 sea ice output, with comparison to the GC3.1 model used within CMIP6; Sect. 5 ends the paper with some discussion and future plans.

2 GSI9 sea ice configuration

The GSI9 sea ice model configuration is based on the native NEMO sea ice model, SI³ (Vancoppenolle et al., 2023), which was developed from the LIM3 model of Rousset et al. (2015) with some functionality merged from CICE. SI³ was first made available at NEMO version 4.0; it is fully embedded in the code and invoked from within the surface boundary code (SBC) module. The version of SI³ used for GSI9 is based on the NEMO version 4.0.4 release as described in Guiavarc'h et al. (2024). NEMO is the ocean modelling framework of choice for the UK and has formed the basis of global ocean model components for over a decade. Met Office is part owner of NEMO, and therefore SI³, as one of two UK NEMO consortium members (alongside NOC). Use of SI³, therefore, offers considerable efficiencies related to management overheads, technical development of the code, and integration into Met Office systems. For example, an advantage of using the sea ice model native to NEMO is that the coupling is simplified; interpolation of velocity points required between NEMO (Arakawa C-grid) and CICE (Arakawa B-grid) at previous configurations (Hewitt et al., 2011) is no longer necessary.

2.1 Model structure

Aside from the change in sea ice model, the sea ice physics used within GSI9 remains similar to the previous, CICE-based, GSI8.1 configuration documented in Ridley et al. (2018). Like CICE, SI³ is a dynamic–thermodynamic continuum sea ice model that includes an ice thickness distribution (ITD; see Thorndike et al., 1975), conservation of horizontal momentum, an elastic–viscous–plastic (EVP) rheology, and energy-conserving halo-thermodynamics (Vancoppenolle et al., 2023). SI³ is run on the same grid as the NEMO ocean model component and on every ocean time step; the sea ice “levitates” above the modelled ocean surface, rather than being embedded within it.

For the GSI9 configuration, five thickness categories are used to model the subgrid-scale ITD, and an additional ice-free category represents open water. The bounds of the thickness categories are determined using a function of domain-mean ice thickness, which is specified as 2.0 m (namelist variable rn_himean). This sets the maximum thickness category bounds to 0.00, 0.45, 1.13, 2.14, 3.67, and 99.0 m. As with previous model versions, the ice–atmosphere exchange is undertaken separately for each ice thickness category using a conductivity coupling scheme in which surface exchanges are calculated externally within JULES. More details on the model coupling can be found in Sect. 3 below.

The sea ice namelist values used in the GSI9 configuration are provided in Appendix A, with SI³ options in Table A1 and JULES options in Table A2. Departures from the SI³ default options are highlighted in the tables, and descriptions of these changes are included. Full details of the GSI9 sea ice configuration are provided in the following subsections, covering dynamics, thermodynamics, and radiation components.

2.2 Dynamics

Horizontal sea ice velocities are calculated by solving the momentum equation of Hibler (1979), which includes terms for internal ice stress, wind and ice–ocean stresses, sea surface tilt, and Coriolis effects. The ice properties of the different thickness categories are advected following the horizontal velocity field, using the second-order scheme of Prather (1986).

After the advection has been performed, mechanical deformation and lead opening convert thinner ice to either thicker ice or open water, by redistributing the global ice state variables into the different ice thickness categories. Following Thorndike et al. (1975), the redistribution function is separated into three components: (i) dynamical inputs (opening and net closing rates); (ii) the participation function, which describes the amount of ice with a given thickness participating in the mechanical deformation; and (iii) the transfer function, which determines to where in thickness space the ice is transferred as a result of deformation. The opening and net closing rates are determined following Flato and Hibler (1995), using a formulation that includes energy dissipation by shear and convergence and a deformation term which relates this to the rheology. The EVP rheology of Hunke and Dukowicz (2002) is used, which employs an elastic wave modification to improve the computational efficiency of a viscous–plastic (VP) rheology. Although the “adaptive” version of the EVP rheology (aEVP; Kimmritz et al., 2016) is the default in SI³, this is not used in the GSI9 configuration owing to issues arising from interaction with the ice shelf basal melt parameterisation in the ocean model (Guiavarc'h et al., 2024; Storkey et al., 2018). Super-cooled water at the edge of the Ronne–Filchner ice shelf led to continuous sea ice growth in the Weddell Sea using aEVP, where a much higher proportion of stationary ice is simulated than for EVP, associated with the improved convergence of the aEVP rheology. These issues will be addressed in future configurations. Ice strength is parameterised following Hibler (1979), which represents a departure from previous CICE configurations where the strength scheme of Rothrock (1975) was utilised.

The participation function for the mechanical redistribution is that of Lipscomb et al. (2007), which favours the closing of open water and deformation of thin ice over the deformation of thicker ice. Participation is independent of whether the deformation process is rafting or ridging. The transfer function considers rafting and ridging separately: rafting doubles the ice thickness, and newly ridged ice is linearly redistributed to the new thickness categories using a function based on Hibler (1980). Mechanical redistribution in SI³ is formulated to conserve ice area, with ice volume remaining constant under rafting. Under ridging, mass from the ocean is added to the sea ice since observations show that newly formed pressure ridges are porous (Leppäranta et al., 1995; Høyland, 2002). Additionally, during deformation, a fraction of snow falls into the ocean, set here to 50 % (via namelist parameters rn_fsnwrdg and rn_fsnwrft).

Ice–ocean momentum exchange is calculated following a standard approach, using a simple bulk formula with a constant exchange coefficient and rotation angle and using the ocean current velocity provided by NEMO. The neutral drag coefficient has been increased for GSI9 from the SI³ default of 5.0 × 10⁻³ to 1.0 × 10⁻², consistent with the previous GSI configurations. This was partly motivated by improvements to the sea ice shown by Roy et al. (2015) when using increased drag. In order to increase the ocean time step, semi-implicit ice–ocean drag was implemented within the ocean component of GC5 (Guiavarc'h et al., 2024). Wind stress is provided as an external forcing, calculated in JULES as part of the surface exchange scheme. A parameterisation of ice–atmosphere form drag based on Lüpkes et al. (2012) is used within JULES, following Renfrew et al. (2019), including stability dependence and floe size as part of the form drag parameterisation (Lock et al., 2022). This scheme is different from that used in previous GSI configurations and results in a net reduction of atmosphere–ice drag, as detailed in Renfrew et al. (2019).

2.3 Thermodynamics

Thermodynamic growth and melt of the sea ice are modelled after the dynamics are applied, using a multilayer scheme based on Bitz and Lipscomb (1999). For each thickness category, as in previous GSI configurations, SI³ models a snow–ice column consisting of four vertical layers of ice, plus an optional snow layer above. The SI³ default is for two ice layers, but Vancoppenolle et al. (2023) suggest that between two and five is suitable. The thermodynamics scheme has been modified following West et al. (2016) to utilise ice surface temperatures and conductive heat fluxes into the ice provided by the surface exchange scheme in JULES. Full details are given in Sect. 3. Ice–ocean sensible heat flux is calculated following Maykut and McPhee (1995) as a function of local turbulent friction velocity and temperature difference between the ice and ocean, with a transfer coefficient specified as 5.7 × 10⁻³. Additionally, the thermal conductivity of snow has been increased from the default 0.31 to 0.50 $\mathrm{W}{\mathrm{m}}^{-\mathrm{1}}{\mathrm{K}}^{-\mathrm{1}}$ to generate thicker ice in winter and improve summertime sea ice area and extent.

The thermodynamics also includes a lateral melting scheme that reduces the ice concentration if sea surface temperature (SST) is above freezing. The method imposes a lateral melt rate as a function of ice concentration and SST, following Bitz et al. (2001) and, unlike in previous GSI configurations, includes a parameterisation for floe size distribution. Further details are provided in Vancoppenolle et al. (2023). The default behaviour in SI³ of using heat in leads for basal melting of the sea ice before heating the ocean has been turned off in GSI9 to allow ocean warming in the marginal ice zone and at the ice edge. This in turn resulted in a large increase in lateral melting, which necessitated tuning of that scheme. The beta exponent and the minimum floe diameter, used in the lateral melting scheme to describe the relationship between ice concentration and floe diameter (see Eqs. 26 and 27 of Lüpkes et al., 2012), were respectively increased from the default 1.0 to the maximum recommended value of 1.2 and from the default 8 m to the maximum recommended size of 10 m. This had the effect of reducing the lateral melting and increasing the basal melting.

Once the thermodynamic melt and growth rates have been calculated, ice properties are exchanged between neighbouring thickness categories. As described by Rousset et al. (2015), this is analogous to a transport in thickness space, where the velocity is equal to the net ice growth rate and is achieved using the semi-Lagrangian linear remapping scheme of Lipscomb (2001).

Unlike in previous GSI configurations, the vertically averaged bulk ice salinity in SI³ evolves in time, for each thickness category, as a function of salt uptake during ice growth, gravity brine drainage, and flushing (following Vancoppenolle et al., 2009). Salinity is assumed to have a linear vertical distribution with a profile shape dependent on the evolving bulk salinity. Salinity is used for both freshwater exchange and in the calculation of all sea ice thermodynamic properties including specific heat, thermal conductivity, enthalpy, and freezing/melting temperatures. Full details of the salinity scheme are given by Vancoppenolle et al. (2023).

2.4 Radiation

Under the conductivity coupling approach used within GC configurations (see Sect. 3 below), the radiation is calculated externally to the sea ice model as part of the surface exchanges in JULES. The CCSM3 (Community Climate System Model) scheme from CICE (Hunke et al., 2015) is used within JULES for computing albedo and radiative fluxes over sea ice. The scheme remains similar to that used in the previous model version described by Ridley et al. (2018), with the impact of melt ponds on albedo calculated in JULES using ponds modelled by the topographic melt pond formulation of Flocco et al. (2010, 2012) within SI³. A further modification has been made here to allow the penetration of visible light into the sea ice. More details on these aspects are provided in Sect. 3 below as part of the coupling documentation.

3 Coupling SI³ within GC5

As SI³ is called from within the NEMO ocean model, the atmosphere–ice coupling is inherited from the atmosphere–ocean coupling framework, OASIS3-MCT (hereafter referred to as simply OASIS). The ocean (NEMO-SI³) and atmosphere (UM-JULES) models run in parallel for each coupling period, which contains several (typically between 2 and 4) ocean and atmosphere time steps. At each coupling instance, all ice variables required by the atmosphere are passed first to the ocean model, where they are sent to the OASIS coupler along with other ocean model variables to be remapped to the atmosphere grid. UM-JULES reads these variables, which function as the bottom boundary condition for the ensuing coupling period. UM-JULES then outputs to OASIS its own variables which are required by NEMO-SI³, and these are also remapped to the ocean grid by OASIS. Upon being read by NEMO, variables required by the ice are passed to SI³, functioning as the top boundary condition for the ensuing coupling period (concurrent with that of the atmosphere). This architecture is consistent with the framework documented for the initial implementation of HadGEM3 by Hewitt et al. (2011) and used for previous GC model configurations.

Figure 1 Schematic overview of the coupling between UM-JULES and NEMO-SI³ used in GC5. The general coupling approach is illustrated in the form of a time series across the bottom of the schematic; the upper 80 % of the figure shows the detail for a single coupling instance, including a full list of variables passed in both directions with the key arithmetic operations performed on them. Here $a_t^{\text{OCN}}$ denotes the ice area fraction field sent by NEMO-SI³ at the coupling instant shown, $a_t^{\text{ATM}}$ the re-gridded ice area fraction field received by UM-JULES after being passed through OASIS, and $a_{t-\mathrm{1}}^{\text{ATM}}$ the ice area fraction field from the previous coupling instant, used by UM-JULES for the preceding coupling period. Horizontal dashed red lines are used to denote the coupling interface between NEMO-SI³ and UM-JULES via OASIS. Solid coloured arrows denote the passage of information between the various components.

Table 1 List of variables passed from NEMO-SI³ to UM-JULES each coupling cycle. Note that not all field indices are documented here; field indices 12–16 are intentionally omitted as they correspond to variables that were used previously in the coupling but are not used in GC5.

Table 2 As Table 1 but for variables passed from UM-JULES to NEMO-SI³. Note that field indices 57–58 are intentionally omitted as they correspond to variables that were used previously in the coupling but are not used in GC5.

Table 3 Typical coupling frequencies with corresponding model time steps and nominal resolution used by the atmosphere–land and ocean–sea ice components for the low- and medium-resolution configurations of the GC5 coupled model.

The coupling framework is demonstrated schematically in Fig. 1. A complete list of variables passed is shown in Table 1 (ocean to atmosphere) and Table 2 (atmosphere to ocean). The various re-gridding methods used for the different variables are indicated. Energy or mass flux variables are passed using conservative remapping. Where possible, second-order conservative remapping is used for increased accuracy. However, most fields contain sharp, irregular horizontal gradients and are bounded above or below by physical constraints (e.g. sea ice fraction must be between 0 and 1), properties that make second-order conservative remapping undesirable owing to the potential for overshoot. Therefore first-order conservative remapping is used in most cases. Dynamical variables such as wind speed, for which conservation is not required, are passed using simple bilinear remapping. Typical coupling frequencies, model time steps, and resolutions used for the GC5 model configuration can be found in Table 3.

3.1 Conductivity coupling in GC5

SI³ is coupled to UM-JULES using a conductivity coupling framework, in which surface variables are calculated within the atmosphere model surface exchange scheme, and the thermodynamic interface between the two models is placed below the surface within the top layer of the ice–snow column. This contrasts to standard bulk formulae coupling in which this interface is placed above the ice–snow surface, and the sea ice thermodynamics solves for the surface variables (i.e. temperature and energy fluxes) as well. In previous GC versions, CICE was coupled to UM-JULES using the same method, based on the implicit coupling framework described in Best et al. (2004).

The conductivity coupling framework is used instead of the traditional framework to enable surface processes to respond more quickly to changes in the atmospheric boundary layer (West et al., 2016). By placing the surface exchange in the atmosphere model, the surface temperature and surface flux can respond instantly to changes in near-surface atmosphere conditions (and vice versa), whereas in the standard framework there would be a delay equal to the coupling period length before either could respond (see West et al., 2016). The sea ice temperatures still experience a delayed response, but this is considered a lesser problem as these are already subject to a damped, slowed response in reality. A secondary reason for using the conductivity coupling framework is to maintain consistency with all land surface types for which UM-JULES also calculates surface variables and exchanges (Best et al., 2011). The coupling methodology is described briefly below and is discussed in more detail in West et al. (2016) and Ridley et al. (2018).

In the conductivity coupling framework, the surface exchange calculations over sea ice are carried out in the atmosphere model using the implicit scheme of Best et al. (2004) and not within the sea ice models SI³ or CICE. The following four energy fluxes are output from the surface exchange and sent through the OASIS coupler to NEMO-SI³ to be used as the forcing for the sea ice thermodynamics:

1. top conductive flux, the flux of conduction from the surface of the snow–ice column to the middle of the top thermodynamically active layer

2. penetrating solar flux, the flux of penetrating solar radiation from the surface of the snow–ice column to the middle of the top thermodynamically active layer

3. top melting flux

4. net sublimation flux.

The sea ice thermodynamics scheme then solves for new temperatures in the ice and snow layers using this forcing but does not solve for surface variables. At the end of each coupling period, the temperature and effective conductivity (conductivity divided by half the layer thickness) of the top thermodynamically active layer are passed through OASIS for the atmosphere to use as a bottom boundary condition for the surface exchange in the ensuing coupling period.

The conductivity coupling framework used is similar to that employed by the CMIP6 GC3.1 configuration (Ridley et al., 2018) and in the initial implementation of HadGEM3 described by Hewitt et al. (2011) (using a much more basic zero-layer thermodynamic configuration). However, its implementation in GC5 differs in a number of ways to the documentation provided in Hewitt et al. (2011) as described below:

1. All variables relating to sea ice are now passed separately for each thickness category. This enables the surface exchange scheme to make full use of the ice thickness distribution, allowing better simulation of rapid sea ice growth in areas of thin ice (for example, as documented by Holland et al., 2006).

2. The use of multilayer thermodynamics in GC5 entails the passing of two additional variables from ocean to atmosphere: the temperature and effective conductivity of the top layer of the snow–ice column. If snow thickness is zero, the top layer is taken to be the top ice layer; as snow thickness increases from 0 to a threshold hs_min, the top layer temperature and conductivity passed to the atmosphere change linearly from the values in the top ice layer to those in the snow. Note that this is an update from GC3.1 (Ridley et al. (2018), where quantities used changed abruptly from the top ice layer to a snow layer as snow thickness crossed the threshold hs_min.

3. GC5 uses semi-implicit coupling to pass the four atmosphere–ice energy fluxes described in Sect. 3.1. In this formulation, UM-JULES does not pass grid-box-mean fluxes to the ocean. Instead, it divides these by ice concentration upon receiving the new values from NEMO-SI³ to create “pseudo-local” fluxes. These fields are passed through OASIS to NEMO-SI³ where they are multiplied by the same ice concentration field. The resulting grid-box-mean fluxes are provided to the sea ice model for use over the ensuing coupling period. This formulation is necessary to both globally conserve energy and force the ice thermodynamics with energy fields proportional to the amount of ice in each grid cell. A full description of, and justification for, the semi-implicit coupling approach is given in Ridley et al. (2018).

4. The radiative melt-pond scheme used by GC5 entails the passing of additional variables in each direction. Surface temperature is passed from the atmosphere to the ocean to be used in the melt-pond scheme to determine growth and melt of pond refrozen lids, whilst melt-pond effective area fraction (the fraction of sea ice covered by radiatively active melt ponds) and melt-pond depth are passed from ocean to atmosphere to be used in the radiation scheme for calculating albedo.

5. Penetrating solar radiation is now modelled. Hence, in addition to the three atmosphere–ice fluxes previously included (top conductive flux, top melt flux, and net sublimation flux), a fourth energy flux, solar radiation penetrating the sea ice surface, is passed from the atmosphere to the ocean. The proportion of penetrating solar radiation is calculated as an extension of the Semtner (1976) scheme used in previous GC configurations (Ridley et al., 2018). Visible light that penetrates the sea ice and is not scattered back out is passed through the coupler to be used in the sea ice model.

Several other sea ice variables are passed to the atmosphere beyond those already discussed: ice and snow thickness, which are used in the albedo calculations; ice area, which is fundamental in quantifying the contribution of the sea ice surface exchange to the whole grid cell; and combined ice and ocean velocity, which is used both in dynamic boundary layer calculations and in calculating turbulent fluxes in the surface exchange.

3.2 Adaption of SI³ for use with conductivity coupling

As detailed above, conductivity coupling has been used with all previous GC versions, for which CICE was adapted to work with this method. To implement conductivity coupling in SI³ for GC5, two major modifications were required. Firstly, the NEMO coupling interface has been changed to allow the top conductive flux to be received from the OASIS coupler and used as the top boundary condition for the thermodynamic solver, rather than the surface exchange boundary conditions (downwelling radiative fluxes, air temperature, specific humidity). Likewise, the coupling has been modified to send the new ice–atmosphere variables (temperature and effective conductivity) from the topmost thermodynamically active layer through the OASIS coupler to be used in the surface exchange calculation. Secondly, the thermodynamic routine (icethd_zdf_bl99) was modified to solve for only internal snow and ice layer temperatures, leaving the surface temperature equation to be calculated elsewhere. This option is controlled in SI³ by a logical (ln_cndflx) within the surface boundary condition namelist. In addition, there is an option to “emulate” the conductivity coupling approach in cases where no external surface exchange scheme is available (e.g. when testing or running NEMO-SI³ in forced-atmosphere mode). This is controlled by an additional namelist logical (ln_cndemulate), where SI³ calculates the conductivity fluxes needed for the surface boundary condition from the usual bulk formulae input using its own surface exchange calculation.

3.3 Assessment of energy conservation in the GC5 coupling

To compare energy conservation across the coupler in GC5, global area averages of top conductive flux, top melt flux, and net sublimation flux sent to the OASIS coupler were compared to global area averages of the re-gridded flux fields received by the ocean. Over the course of a 1 d simulation, average errors were of the order of 2 × 10⁻³ W m⁻² for top conductive flux, 2 × 10⁻⁵ W m⁻² for net sublimation flux, and 5 × 10⁻⁶ W m⁻² for top melt flux, approximately 0.05 %, 0.005 %, and 0.05 % of the absolute flux fields respectively. These errors are similar in magnitude to those reported in Sect. 4 of Hewitt et al. (2011) for HadGEM3-AO, as well as for the GC3.1 CMIP6 configuration (not shown).

4 Model evaluation

In this section we present a brief evaluation of the sea ice simulated by the GC5 coupled model using the above-documented GSI9 sea ice model configuration. The intention is not to provide a thorough assessment of the sea ice performance in GC5 but rather a sanity check that the configuration documented here is performing sensibly. As per previous coupled model versions, GC5 has been developed with traceable science across a hierarchy of model resolutions, with ocean (atmosphere) resolution ranging from 1° (130 km) to $\mathrm{1}/\mathrm{12}$° (25 km) (e.g. Guiavarc'h et al., 2024; Storkey et al., 2018; Roberts et al., 2019). The sea ice model options used are identical across all the different model resolutions, and so we limit our attention here to the medium-resolution model, which uses eORCA025 (nominal $\mathrm{1}/\mathrm{4}$°) ocean–sea ice resolution and N216 (∼ 60 km in midlatitudes) atmosphere–land resolution.

The simulated sea ice evaluated in this section is the last 50 years from 100-year “present-day” control runs forced by greenhouse gases and aerosols from the year 2000 (see Williams et al., 2017). Simulated sea ice from GC5 is compared against the GC3.1 CMIP6 configuration (Williams et al., 2017; Ridley et al., 2018), along with reference datasets including observations of sea ice concentration from the Hadley Centre Sea Ice and Sea Surface Temperature dataset version 2 (HadISST.2.2.0.0; Titchner and Rayner, 2014), observations of sea surface temperature from the ESA CCI SST L4 dataset of Good et al. (2019), and sea ice thickness from the Pan Arctic Ice Ocean Modeling and Assimilation System (PIOMAS; Schweiger et al., 2011) reanalysis, which assimilates observations of sea ice concentration but not sea ice thickness. The present-day simulations employed here represent how the climate would evolve if emissions were fixed at the year-2000 level for 100 years, which of course can differ from observed conditions for the period 1990–2009. Thus, our comparisons with observations can be considered more of a benchmark than a direct assessment.

Figure 2 (a) The model seasonal cycle of sea ice area (10⁶ km²) in the Arctic for GC5 (red) and GC3.1 (blue). Area estimates from the HadISST.2 sea ice dataset are included in grey, with ± 2 standard error (SE) shading and error bars and ± 20 % indicated with chain lines. Panels (b)–(e) show simulated mean sea ice fraction with HadISST.2 0.15 contours added in orange from GC5 (b, c) and from GC3.1 (d, e) for September and March respectively. GC5 and GC3.1 data are the last 50 years from 100-year model simulations using year-2000 continuous forcing, whilst HadISST.2 data are from the period 1990–2009.

The Arctic seasonal cycle of sea ice area (Fig. 2a) demonstrates that GC5 and GC3.1 are both within 20 % of HadISST.2, a criterion that has been commonly used when evaluating climate models in the context of model selection (e.g. Massonnet et al., 2012), apart from in August. The sea ice area in GC5 is closer than GC3.1 to the HadISST.2 observations apart from late spring and summer when the two are very similar. In GC5 the ice concentration is greater than GC3.1 in the marginal seas in winter (Fig. 2c and e), particularly the Sea of Okhotsk, where it is now closer to HadISST.2. Meanwhile, summer concentration increases are seen in the central Arctic north of Siberia and in the Canadian Archipelago (Fig. 2b and d). When compared to HadISST.2 observations, the mean spatial pattern of ice concentration (Fig. 2c and e) shows excess marginal ice cover in the Greenland Sea and reduced cover in the Labrador Sea, for both models. These biases are related to the preference of the ocean model for convective overturning in the Labrador Sea rather than the northeast Atlantic, as described by Megann et al. (2014) and Guiavarc'h et al. (2024). The bias of an early Arctic summer minimum in August is present in both model configurations. Having an areal minimum for Arctic sea ice in August is not uncommon for models; Roach et al. (2020) showed that around a quarter of CMIP6 models have lower average sea ice area in August than in September. In September Arctic sea ice is still melting at the peripheries but, owing to the onset of the polar night at higher latitudes, starting to freeze up in the centre of the pack. The evolution of September sea ice area is therefore dependent on these two competing processes. Whether September area is higher or lower than August will depend on the timing of when the Arctic transitions from net melting to net growth. It is likely that the timing of this transition is out by a few days in the model. This competition between lower-latitude melting and high-latitude refreeze is not captured in the extent metric, which has a minimum in September for both these model configurations (e.g. see Rae et al., 2015). It is also worth noting here that uncertainty in passive-microwave satellite-derived observations of sea ice is considerably high in summer, where the presence of surface melt ponds can lead to an overestimation of up to 25 % in concentration (see Kern et al., 2020). Therefore, a non-negligible portion of the offset between model and observations could be related to errors in the observations.

Figure 3 March mean sea ice thickness (m) from the last 50 years of 100-year model simulations using year-2000 continuous forcing for (a) GC5 and (b) GC3.1 and (c) from the PIOMAS reanalysis for the period 1990–2009. Panel (d) shows mean sea ice thickness difference between GC5 and GC3.1 (GC5–GC3.1) with hatched areas identifying differences that are significant at the 95 % level, calculated using a Welch t test. Panels (e) and (f) show corresponding 50-year March mean sea ice velocity for GC3.1 and GC5 respectively with coloured shading depicting ice speed and velocity arrows overlain in black.

Although the spatial pattern of winter mean sea ice thickness in GC5 is similar to GC3.1 (Fig. 3a and b), the ice thickness in GC5 has increased across the whole Arctic, with largest increases north of the Canadian Archipelago (Fig. 3d), meaning that the spatial distribution of Arctic sea ice thickness is now more comparable to the PIOMAS reanalysis (Fig. 3c). This improvement in the central Arctic thickness pattern is associated with changes in the sea ice dynamics. The strength of the Beaufort Gyre is lower in GC5 than in GC3.1 (Fig. 3e and f), meaning a longer residence time and subsequent thickening of the multi-year sea ice piled up north of Greenland and the Canadian Archipelago. This reduction in Beaufort Gyre speed is very likely linked to the change in atmosphere–sea ice drag scheme, whereby a net reduction in drag (Renfrew et al., 2019) will have led to reduced sea ice velocity, along the lines of that described by Johns et al. (2021).

Figure 4 (a) The model seasonal cycle of sea ice area (10⁶ km²) in the Antarctic for GC5 (red) and GC3.1 (blue). Area estimates from the HadISST.2 sea ice dataset are included in grey, with ± 2 standard error (SE) shading and error bars and ± 20 % indicated by chain lines. Panels (b)–(e) show simulated mean sea ice fraction with HadISST.2 0.15 contours added in orange from GC5 (b, c) and from GC3.1 (d, e) for September and March respectively. GC5 and GC3.1 data are the last 50 years from 100-year model simulations using year-2000 continuous forcing, whilst HadISST.2 data are from the period 1990–2009. Panels (f) and (g) show the annual mean SST difference (K) from the ESA CCI SST L4 dataset of Good et al. (2019) for GC5 and GC3.1 respectively (model − observations).

The Antarctic sea ice area in GC5 has increased considerably from GC3.1, and the seasonal cycle is now comparable with HadISST.2 (Fig. 4a), albeit with a slight phase lag suggesting that ice growth is too slow in early winter (May–July) and ice melt is too slow in early summer (December). It has been suggested that the suppressed rate of growth in winter is associated with the temporal pattern of insolation (Goosse et al., 2023), with other climate models displaying the same issue (e.g. DuVivier et al., 2020). Antarctic sea ice concentration has increased, and extent has expanded, in GC5 compared with GC3.1 at all times of the year, as illustrated for September and March respectively in Fig. 4b and d and Fig. 4c and e. These are now much closer to the HadISST.2 dataset. This is because of a considerable improvement in GC5 of Southern Ocean surface temperatures (see Fig. 4f and g and Storkey et al., 2024), where previously a warm ocean bias led to low sea ice area (Ridley et al., 2018; Rae et al., 2015). Despite the general increase in Antarctic sea ice cover, as illustrated in Fig. 4b and d, there is a minor reduction in the Weddell Sea compared to GC3.1. This is associated with the emergence of multi-year sensible heat polynyas, a feature common to many climate models (Heuzé et al., 2015) including Met Office configurations (Megann et al., 2014; Ridley et al., 2022).

5 Conclusions

In this paper we have presented the UK's GSI9 sea ice model configuration, used within the latest version of the Met Office Global Coupled model configuration GC5, which will form the physical model basis for UK contributions to CMIP7 with HadGEM3 and UKESM. GC5 includes a change to the sea ice model component compared with earlier GC versions with the implementation of the new NEMO SI³ model in place of CICE. We have described how SI³ has been adapted to work with the conductivity coupling used in Met Office models and provided a thorough documentation of sea ice (and wider) coupling in GC5. A brief evaluation of the GC5 sea ice using continuous year-2000 climate forcing has been presented, which shows that the sea ice simulated by this configuration compares well with observational references. A comparison was also performed with the CMIP6 model version, GC3.1, which shows that the mean state and variability of the GC5 sea ice are generally improved compared to GC3. This is particularly so in the Antarctic where the sea ice is much improved throughout the year in response to the reduction of warm biases in the ocean, as described in Storkey et al. (2024).

Future development of the Global Sea Ice configurations will include exploring

• the ridging sea ice strength formulation of Rothrock (1975) and the exponential ITD transfer function of Lipscomb et al. (2007), which were used in previous GC configurations. Although these schemes have since been included in SI³ (under the EU IS-ENES3 project), they were not available in time to be used in the GSI9 configuration.

• the land–fast ice modelling scheme of Lemieux et al. (2016), which is already available as an option in SI³.

• alternative sea ice rheology schemes, including the adaptive elastic–viscous–plastic (aEVP; Kimmritz et al., 2016) and elastic–anisotropic–plastic (EAP; Tsamados et al., 2013) rheologies, which are already included in SI³ (through the EU IMMERSE and IS-ENES3 projects).

• improved representation of ice–ocean and ice–atmosphere drag using the form-drag scheme of Tsamados et al. (2014), which has been included in JULES through the EU-APPLICATE project and is currently being ported into SI³; another option for improving the ice–ocean drag is to adopt the methodology outlined in Roy et al. (2015).

• the inclusion of floe–size distribution and wave–ice interaction as discussed by Bateson et al. (2022).

• the adaptation of the radiation scheme to include penetrating shortwave radiation into the sea ice under melt ponds or simulate the freshwater impacts of melt ponds.

Appendix A: GSI9 model namelists

Table A1 SI³ namelist and hard-wired parameters of scientific significance. Variables that are changed from the SI³ defaults are highlighted with an asterisk with default values given below in brackets.

Table A2 JULES sea ice namelist updated from Ridley et al. (2018).