The contribution of melt ponds to enhanced Arctic sea-ice melt during the Last Interglacial

HadGEM3 is the first coupled climate model to simulate an ice-free Arctic during the Last Interglacial (LIG), 127 000 years ago. This simulation appears to yield accurate Arctic surface temperatures during the summer season. Here, we investigate the causes and impacts of this extreme simulated ice loss. We find that the summer ice melt was predominantly driven by thermodynamic processes: atmospheric and ocean circulation changes did not significantly contribute to the ice loss. We 5 demonstrate these thermodynamic processes were significantly impacted by melt ponds, which formed on average 8 days earlier during the LIG than during the pre-industrial control (PI) simulation. This relatively small difference significantly changed the LIG surface energy balance, and impacted the albedo feedback. Compared to the PI simulation: in mid-June, of the absorbed flux at the surface over ice-covered cells (sea-ice concentration>0.15), ponds accounted for 45-50%, open water 35-45%, and bare ice and snow 5-10%. We show that the simulated ice loss led to large Arctic sea surface salinity and temperature changes. 10 The sea surface temperature and salinity signals we identify here provide a means to verify, in marine observations, if and when an ice-free Arctic occurred during the LIG. Strong LIG correlations between spring melt pond and summer ice area indicate that, as Arctic ice continues to thin in future, the spring melt pond area will likely become an increasingly reliable predictor of the September sea-ice area. Finally, we note that models with explicitly modelled melt ponds seem to simulate particularly low LIG sea-ice area. These results show that models with explicit (as opposed to parameterised) melt ponds can simulate very 15 different sea-ice behaviour under forcings other than the present-day. This is of concern for future projections of sea-ice loss.


Introduction
Interglacials are periods of globally higher temperatures which occur between cold glacial periods (Sime et al., 2009;Otto-Bliesner et al., 2013;Fischer et al., 2018). Glacial-interglacial cycles are largely driven by changes in the Earth's orbit which affect incoming radiation. The Eemian period, here called the Last Interglacial or LIG, occurred 130,000-116,000 years ago. 20 At high latitudes orbital forcing led to summertime top-of-atmosphere short-wave (TOA SW) radiation 60-75 W m −2 greater for the LIG, compared to the preindustrial (PI) period. This drove differences between the LIG and PI surface energy balance.
Whilst the significance of these PI to LIG surface energy balance differences vary between climate models Otto-Bliesner et al., 2013), proxy records of the LIG indicate that mean Arctic summer land temperature was + 4-5°K higher than the PI (CAPE members, 2006). 25 Prior to 2020, most climate models simulated LIG temperatures which were too cool compared with LIG temperature data (Otto-Bliesner et al., 2013;IPCC, 2013). Recently, Guarino et al. (2020b) found that the loss of Arctic sea-ice in the summer likely drove these warm Arctic temperatures. As discussed below, Guarino et al. (2020b) suggested that melt ponds may have been a significant driver of ice loss, and that previous climate models, with a less comprehensive representation of melt ponds, may not have simulated the loss of enough Arctic sea-ice during the LIG. Kageyama et al. (2021) explored the ocean core based 30 proxy records of LIG Arctic sea-ice change. They found that sea-ice changes are more difficult to determine than temperature changes, with some conflicting interpretations of proxy data from the available records, and imprecision in dating materials from cores in the high Arctic. This makes it difficult to determine the mechanisms or distribution of sea-ice loss during the LIG from these preserved biological data.
In terms of understanding mechanisms that drove Arctic sea-ice change during the LIG, three main factors are known to 35 affect summer sea-ice behaviour: albedo feedbacks (Curry et al., 1995), cloud cover feedbacks (Kay and Gettelman, 2009), and ocean heat transport changes (Årthun et al., 2019;Auclair and Tremblay, 2018). As well as these, changes to sea-ice distribution caused by changes to wind patterns and ocean circulation may also affect summer sea-ice extent (Kageyama et al., 2021;Deser and Teng, 2008;Wang et al., 2009). Albedo feedbacks are strongly influenced by melt ponds, systems of pools that form from meltwater and begin to collect on the Arctic ice surface in spring. Pond-covered ice has a lower albedo, at 0.1-0.5, 40 than bare ice, at 0.6-0.65, or snow, at 0.84-0.87 (Perovich et al., 2002b;Eicken et al., 2004;Perovich and Tucker, 1997).
Pond-covered ice thus absorbs a higher fraction of incident solar radiation, and transmits a greater fraction of incident radiation to the ice and ocean below. This difference accelerates the melting of the ice beneath ponds, with melt-rates of pond-covered ice up to 2-3 times the melt-rate of bare ice (Fetterer and Untersteiner, 1998). Over the last decades, melt ponds have played a key role in reducing the surface albedo (Eicken et al., 2004;Maslanik et al., 2007;Perovich et al., 2007); throughout melt 45 season, nearly 60% of the summer sea-ice area may be covered by ponds (Fetterer and Untersteiner, 1998;Eicken et al., 2004).
In spite of their importance, melt ponds have only rather recently been explicitly included in CMIP models. In CMIP6 models, the most common approach is to implicitly parameterise melt ponds by reducing the ice/snow albedo when surface ice temperatures approach 0°C (e.g. Collins et al., 2006;Curry et al., 2001). This tuning has been relatively successful for reproducing realistic melt rates for the present day (Collins et al., 2006;Curry et al., 2001). However, pond formation is 50 affected by sea-ice processes throughout melt season (e.g. evolving topography and snow cover), which this tuning does not represent (Kwok et al., 2009). Therefore, in recent years, there has been increasing interest in incorporating more detailed melt pond models into GCMs. HadGEM3 includes one of the most comprehensive, to date, melt pond schemes in its sea-ice component CICE5.1 (Hunke et al. (2015), detailed in Section 2), a result of a series of developments and improvements (Taylor and Feltham, 2004;Lüthje et al., 2006;Skyllingstad et al., 2009;Scott and Feltham, 2010;Flocco et al., 2012).

55
The Coupled Model Intercomparison Project Phase 6 (CMIP6) Palaeoclimate Model Intercomparison Project Phase (PMIP4) or CMIP6-PMIP4 LIG experimental protocol prescribes differences between the LIG and PI in orbital parameters, as well as differences in trace greenhouse gas concentrations (Otto-Bliesner et al., 2017). This standardised climate modelling protocol enables the community to use models to explore these mechanisms using a multi-model approach.
16 models ran the CMIP6-PMIP4 LIG simulation. All 16 models showed a substantial reduction in LIG Arctic sea-ice 60 compared to the PI (Kageyama et al., 2021). They yielded a minimum Arctic sea-ice area which ranged between 0.2 -5.7 million km 2 . Whilst inter-model differences were variously attributed to differences in the albedo feedback, ocean circulation and heat transport, atmospheric circulation, and cloud cover, these aspects have not yet been fully analysed for all the models which ran the simulation (Guarino et al., 2020b;Kageyama et al., 2021;Otto-Bliesner et al., 2020).
Of these 16 models, Guarino et al. (2020b) showed that the model HadGEM3 gave a good match with proxy temperatures 65 (related to its complete simulated loss of summer sea-ice): the average LIG temperature anomaly in HadGEM3, for all locations with observations, was +4.9 ± 1.2 K compared with the observational mean of +4.5 ± 1.7 K. Guarino et al. (2020b) indicated that the HadGEM3-simulated summer LIG sea-ice loss was highly influenced by albedo changes, and suggested this was due to HadGEM3's detailed representation of sea-ice, and specifically melt pond, physics. This model also had a good match with all, except one, marine core sea-ice datapoint (Kageyama et al., 2021). This, alongside the complete summer sea-ice loss, makes 70 this model of interest for aiding our understanding of sea-ice loss mechanisms, particularly those related to melt ponds, both for the LIG and for the future (Guarino et al., 2020b).
Here, we analyse this first simulation of an ice-free Arctic during the LIG using HadGEM3 (Guarino et al., 2020b), with explicit melt pond dynamics (CICE 5.1): we examine in detail how melt ponds contributed to Arctic sea-ice loss during the LIG, particularly the enhanced LIG summer sea-ice loss. In this more in-depth follow-on study to Guarino et al. (2020b), we 75 analyse potential thermodynamic and dynamic contributors, and provide multi-model context that was not available to Guarino et al. (2020b). We investigate possible ice loss drivers, quantifying the impact of thermodynamic or dynamic processes for the enhanced loss. We then investigate thermodynamic processes in detail, and study what drove LIG surface albedo changes, and in particular the impact of melt ponds on the LIG-PI albedo difference. Finally, we study the predictability of summer sea-ice loss from the spring melt pond area, and compare results between the two simulated periods. This addresses key gaps in our 80 understanding of how melt ponds impact sea-ice behaviour in warm climates, and in particular, answers the question of how HadGEM3's melt pond scheme led to the simulated ice-free LIG Arctic.

The HadGEM3 model
All the simulations analysed in this study use the low resolution version of the latest UK physical climate model, HadGEM3-85 GC31-LL, hereafter HadGEM3 (Williams et al., 2018). HadGEM3 is a fully coupled climate model that uses the Unified Model (UM) (Walters et al., 2017) for the representation of the atmosphere, the Joint UK Land Environment Simulator (JULES) for the representation of land surface processes (Walters et al., 2017), and the NEMO3.6 (Madec et al., 2015) and the CICE5.1 (Hunke et al., 2015) models for the representation of the ocean and the sea-ice, respectively. a regular latitude-longitude grid for the atmosphere. For the ocean, an orthogonal curvilinear grid with a grid-spacing of approximately 1 • is used. Note that the grid-spacing for the ocean model decreases down to 0.33 • between 15 • N and 15 • S of the equator, as described by Kuhlbrodt et al. (2018). For the vertical discretization, the UM atmospheric model utilizes 85 pressure levels (terrain-following hybrid height coordinates) while the NEMO ocean model uses 75 depth levels (rescaledheight coordinates).
The standard elastic-viscous-plastic rheology (EVP) has been applied for ice dynamics with default CICE remapping advection algorithm and ridging schemes (Hunke et al., 2015). Ice thermodynamics are based on Bitz and Lipscomb (1999) with 4 ice and 1 snow layer. A semi-implicit coupling scheme between atmosphere and sea-ice has been introduced to ensure the stability of the solver (West et al., 2016). The evolution of the sea-ice is separately calculated for 5 ice thickness categories 100 within each grid-cell. For our study it is important to mention that the albedo calculation is based on the scheme used in the CCSM3 model (Hunke et al., 2015), but includes surface melt ponds by applying the explicit topographic melt pond model of Flocco et al. (2012) and Flocco et al. (2010). Meltwater, formed as a result of snow melt, ice melt and precipitation, runs downhill under the influence of gravity and collects on sea-ice starting at the lowest surface height applying the sub-grid scale ice thickness distribution. The evolution of pond fraction and depth as well as the formation of ice lids are calculated. The 105 albedo of ponds of depth >20cm is 0.27 (Ridley et al., 2018b), significantly smaller than the albedo of ice and snow, which varies through the year in the range 0.5-0.9 and is calculated as described in Hunke et al. (2015). The albedo of open water is 0.07±0.03 at the latitudes and over the months -above 70°N, over the months March-July -for which the open water albedo is considered in our analysis (Ridley et al., 2018b). In other sea-ice models without an explicit pond scheme, the ice albedo is reduced when the surface temperature approaches freezing temperature (see e.g. Hunke et al., 2015) to indirectly account for 110 the impact of melt ponds. This adjustment has been removed here not to double count for the impact of ponds on albedo.
For full details on model configuration, performance and improved physics compared to older model versions see Williams et al. (2018).

Simulations
The Pre-industrial (PI) simulation used in this study was prepared and run by the UK Met Office as part of the sixth coupled 115 model intercomparison project CMIP6 (Eyring et al., 2016). This simulation uses invariant solar, greenhouse gas (GHGs), ozone, tropospheric aerosol, volcanic and land-use forcing for the year 1850, see Menary et al. (2018) for details. The climate system took about 615 model-years of spin up to attain a steady state. These years are not used in our analysis. Of the subsequent 500 model-years of production run (Menary et al., 2018), the first 200 are used here in our analysis.
The LIG simulation analysed in this study was first presented by Guarino et al. (2020b); it constitutes the UK's PMIP4 LIG 120 contribution, as part of the wider CMIP6 project. The LIG experiment fully complies with the standard PMIP4 experimental protocol for Last Interglacial climate simulations, as described by Otto-Bliesner et al. (2017). In more detail, this simulation is a time-slice of the Earth's climate at 127,000 years ago (i.e. 127k). The Last Interglacial climate was forced using 127 k constant astronomical parameters based on Berger and Loutre (1991), and constant atmospheric trace GHG concentrations derived from ice core records (see Otto-Bliesner et al. (2017)  The LIG simulation was initialized from the end of the 615 years of PI spin-up. A further 350 model-years of LIG spin-up were required for the atmosphere and the (upper-) ocean to reach quasi-equilibrium. See Williams et al. (2020) for details on how the LIG spin-up was evaluated and what metrics were used to assess the atmospheric and oceanic equilibria. After having 130 attained quasi-equilibrium, the simulation was continued for further 200 years of production run. This length of simulation has been shown to be long enough to capture model internal variability (Guarino et al., 2020a). The first 35 years of LIG spin-up, and all 200 years of the LIG production run, are used in our analysis.

Analysis
For a given ice-covered cell (defined as a cell with sea-ice concentration> 0.15), the area of the cell at a subgrid scale of (i) 135 bare ice and snow (ii) melt ponds on ice and (iii) ocean that is exposed to the atmosphere and not covered by ice (

Enhanced sea-ice loss during the LIG
Here, we examine the HADGEM3-simulated LIG sea-ice and upper ocean, in order to identify factors that contributed to LIG summer sea-ice loss. We first compare LIG and PI sea-ice area, and sea surface temperatures and salinities. We then consider the LIG production run and spin-up to determine the primary drivers of LIG ice loss.
150 Figure 1 compares the annual cycle of the Arctic sea-ice area between the PI and LIG period from the HadGEM3 simulations.
The LIG winter ice area was slightly lower than the PI, with the smallest difference in area at 0.61±0.56 million km 2 in early April ( Fig. 1). Compared to the PI, an enhanced rate of LIG sea-ice melt from early May until late June led to a consistently ice-free LIG Arctic from early August until early October (Figs 1 and 2). The minimum LIG ice area in September was 0.09 ± 0.11 million km 2 , while the PI minimum ice area was 5.54 ± 0.98 million km 2 . Long-term mean ice thickness during the 155 LIG was also thinner than during the PI in all months, as seen in Fig. 3. The earlier sea-ice retreat resulted in a warmer and less salty ocean surface during July, August and September during the LIG (Figs. 4 and 5). The differences reached more than 5°C (LIG-PI SST anomaly) and around 1-2ppt for temperature and salinity respectively. The interaction of ocean surface conditions and sea-ice extent in the model simulations potentially allow sea surface temperature (SST) observations from proxy records to be used as a signature of ice conditions during the LIG. In 160 particular it may be informative for the marine core community to search for large summer SST differences from the PI, of > 4°C , at latitudes near the expected LIG ice edge (e.g. 70°N; east of Greenland) and smaller summer SST differences from the PI, of >2°C, in the central Arctic/Beaufort sea. This could help identify an sea-ice-free summer LIG Arctic in marine observations. However, we note that a partially ice-covered summer Arctic, with thin LIG sea-ice, might yield a similar signature to this but with lower magnitude temperature anomalies everywhere, as very thin summer ice cover has a small insulating effect (Serreze 165 and Barry, 2011). Thus it is additionally useful to consider sea surface salinity (SSS) changes (Fig. 5). These have a clearer signature than SST. While the salinity patterns were similar in May for both simulations, the stronger melt during the LIG caused the region with LIG winter ice cover to become significantly fresher than the PI, by around 1-2ppt in June and July.
Additionally, a difference of 0.5-1.5ppt was retained at least from March-September. Beyond helping identify the observational signature of a sea-ice-free Arctic, we aim to understand: what caused the large differences in summer sea-ice in this model? 170 Particularly, why did the spring melt increase so significantly in spite of the lower atmospheric CO 2 concentration in the LIG?
Processes with a significant impact on summer sea-ice include: thermodynamic processes such as ice-albedo feedbacks, ocean heat transport, and cloud cover feedbacks. Ice-albedo feedbacks may be significantly impacted by the presence of melt ponds, as discussed in Section 1. The preconditioning of winter sea-ice may also lead to a reduced sea-ice extent the following summer (Parkinson and Comiso, 2013;Williams et al., 2016). In addition, summer sea-ice extent may be affected by changes to sea-ice distribution caused by changes to wind patterns and ocean circulation (Kageyama et al., 2021;Deser and Teng, 2008;Wang et al., 2009) further discussed in Section 3.1.1. It is currently unknown which process was most significant for the enhanced LIG sea-ice loss.
We first look at the spin-up simulation for the LIG, with a particular focus on the first year of the spin-up (Figs. 6 and 7). Fig. 6 shows that the winter Arctic sea-ice retained a similar area to the PI control period over the first 15 years of the spin-up 180 period, beyond which it decreased only slightly (<1 million km 2 ) from the PI control. However, an ice-free summer state was reached after only 4 years (Guarino et al., 2020b); we find August sea-ice area halved during the first year of the spin-up run from nearly 6 million to around 3 million km 2 . Differences between PI and LIG developed rapidly during this first summer, within a few months of the January switch from PI to LIG forcings: this suggests that winter preconditioning did not play the dominant role in the enhanced melting of sea-ice during the LIG spin-up. Similarly, this rapid loss of sea-ice within the first 185 few spin-up years cannot be explained by change to ocean heat transport, which is often linked to reduced sea-ice (Holland et al., 2006;Steele et al., 2010): upper ocean heat transport takes decades to equilibrate in numerical models after a significant perturbation, and deeper ocean heat transport, centuries to millenia (Kantha and Clayson, 2000). Therefore, as can be seen from  Guarino et al. (2020b) period. In order to deduce the importance of two of the remaining processes, ice-albedo and cloud cover feedbacks, we consider the surface energy budget (Fig. 7). This is because the Arctic Ocean's surface heat balance is closely linked to its ice mass balance (Uttal et al., 2002;Untersteiner, 1961): at the surface, shortwave (SW) flux in particular is known to play a dominant role in summer sea-ice melt, and is linked to ice-albedo feedbacks (Maykut and Perovich, 1987); longwave flux is related to the longwave cloud radiative forcing, and cloud-cover feedbacks (Gupta et al., 1992;Niemelä et al., 2001). Summer positive TOA 195 radiation (Kageyama et al., 2021) led to a positive net absorbed SW flux anomaly, of up to 75 Wm −2 in June, that contributed nearly all of the net absorbed surface heat flux anomaly (Fig. 7a). Unlike changes related to preconditioning and ocean heat transport, this net SW flux anomaly was already present in the first year of LIG spin-up run, reaching 55 Wm −2 in June (Fig.   7b). This immediate response, coupled with the immediate halving of August sea-ice area, suggests the absorbed shortwave radiation played a key role in summer ice loss, and thus that ice-albedo feedbacks may also have been important. We note also 200 that the longwave anomaly accounted for < 5 Wm −2 of the total surface heat flux anomaly in Fig. 7b, so longwave forcings and feedbacks related to cloud cover were likely not dominant contributors to the enhanced LIG sea-ice loss, as confirmed by April-September is shown, computed as the time-average from monthly model output. Guarino et al. (2020b). Other differences between the LIG and PI surface heat budget contributed < 20 Wm −2 monthly to the surface energy budget (Fig. 7). This suggests that thermodynamic processes that led to the enhanced LIG summer sea-ice melt must predominantly have resulted from this surface SW anomaly.

205
Therefore, by elimination, the dominant contributors to the enhanced LIG summer sea-ice loss were (i) thermodynamic processes, driven by this surface SW anomaly and likely related to ice-albedo feedbacks, and/or (ii) changes to the summer ice distribution. In the next section, we investigate the relative importance of these processes.

Thermodynamic versus dynamic processes
Sea-ice increase and loss are driven by a combination of 'thermodynamic processes', which involve ice-atmosphere and ice-210 ocean heat fluxes, and 'dynamic processes', which involve changes to the local ice volume due to convergent or divergent ice motion (Li et al., 2014) caused by wind or ocean stress (Goosse and Fichefet, 1999). The key driver of ice motion is the 9 Figure 5. Mean sea surface salinity (in parts per thousand, ppt), for the LIG and PI. The 200-year mean over the Northern Hemisphere for each month from April-September is shown, computed as the time-average from monthly model output.
wind stress forcing (Köberle and Gerdes, 2003): present-day interannual variability in Arctic summer sea-ice extent is linked to changes to sea-ice distribution caused by wind pattern variability due, for example, to large-scale atmospheric variability (Deser and Teng, 2008;Wang et al., 2009). In order to determine to what extent the enhanced LIG sea-ice loss was caused 215 by these dynamic processes, or by thermodynamic processes related to the SW anomaly shown in Fig. 7, we examine the ice volume tendencies due to thermodynamics and due to dynamics (both of which are primary model output variables). These are shown respectively in Fig. 8 and Fig. 9. The largest difference occurred in June with 70 to 110 cm/month of ice melt during the LIG, in contrast to the 10 to 30 cm/month ice melt in the PI simulation (Fig. 8). The thickness changes caused by dynamic processes (Fig. 9) were smaller: note that different scales are used in the maps for figure clarity. While a divergent ice drift 220 reduced ice thickness during most months in the PI simulation, this was not the case during the LIG period. However, the  differences in magnitude between the two periods were generally less than 10 cm/month. This demonstrates that, during spring and summer, changes in thermodynamic rather than changes in dynamic processes were the first-order driver causing increased ice melt for the LIG simulation. Thus, the SW anomaly shown in Fig. 7 was key to the enhanced LIG ice loss. Therefore, we investigate in the next section the source of this net absorbed SW anomaly at the surface.

Melt ponds and albedo feedback
At the ice surface, 20-30% (Perovich et al., 2002a) of downwelling SW radiation is directly absorbed. Absorbed radiation may be transmitted through the ice to the ocean, which warms, leading to processes including further sea-ice melt (Perovich et al., 2002a). The ratio between downwelling and absorbed SW radiation is determined by the surface albedo. Sea-ice has a high albedo and thus tends to predominantly insulate the ocean below it by reflecting most of the downwelling SW radiation. This 230 can be seen from the similar spring SST under the PI and LIG Arctic sea-ice (Fig. 4). Therefore, the LIG-PI TOA anomaly mentioned in Section 1 led to a downwelling surface SW flux anomaly, which in turn led to the anomaly of absorbed surface SW flux that is shown in Fig. 7. This absorbed anomaly may have been amplified (or reduced) by ice-albedo feedbacks that differed between the PI and LIG (Fig. 7). In this section, we investigate downwelling surface SW radiation, and its amplification by ice-albedo feedbacks. 235 We first consider the LIG-PI anomaly of the surface albedo in the Arctic Ocean. As can be seen from Figs 10 and 11, this anomaly was small in April, but grew throughout summer. The <-5% difference above 70°N in April grew, in July, to an average of -20%, and up to -35% in some regions. Formation of melt ponds (which have lower albedos than sea-ice) contributes as the magnitude of the change in ice volume due to dynamics during the summer months is significantly smaller than the change due to thermodynamics, the colorbar scale is 1/4 of the colorbar scale shown in Fig. 8. to the albedo feedback effect (Perovich et al., 2002a;Hunke et al., 2015), so we evaluate the magnitude of the changes of pond evolution. As melt onset was not available as a model output variable, we use the mean first day each grid-cell had pond 240 fraction greater than 1% as a proxy for this. The geographical pattern of the first day of melt pond formation was similar for both periods, with ponds forming first at low latitudes and gradually spreading to higher latitudes (Fig. 12). Pond formation began on average 7.8 days earlier during the LIG compared to the PI. The difference between the two periods increases with latitude. Ponds south of 70°N began to form only 0-5 days earlier during the LIG than the PI. Ponds began to form at higher latitudes 15-20 days earlier during the LIG (Fig. 12c). As well as forming earlier in the LIG compared to the PI simulation, 245 ponds also covered a greater area fraction of the ice-covered grid-cells throughout the spring. For both periods, the long-term mean of the pond fraction of the grid-cell for each month was greater for the LIG than the PI from May-June, but smaller for the LIG than PI from July (Fig. 13). This is because melt ponds covered a larger area of the ice for the LIG than the PI in spring to early summer, but very little ice remained for the LIG simulation from July onwards, so there was less sea-ice available for ponds to cover. Fig. 14 shows the climatology of pond formation. At the end of April, the rate of pond formation was greater 250 during the LIG than the PI. From early May, the total Arctic pond area increased exponentially with rate of ∼ 0.14 day −1 , until late May when the rate of pond formation began to decelerate; maximum pond area was reached in mid-June (Guarino et al.,  2020b) which we compute as 2.67±0.37 million km 2 . For the PI, from early May, pond area increased exponentially with rate constant ∼ 0.09 day −1 , until early June when the rate of pond formation decelerated until the maximum pond area was reached at 2.54±0.45 million km 2 in mid-July. A consistently higher fraction of the LIG sea-ice was pond-covered throughout 255 the spring, with a peak value of 45.3±3.7 % (in early July) compared to peak value of 34.4±5.1 % for the PI (in late July -not shown).   In order to (i) compare downwelling SW radiation between the two time periods (ii) quantify the importance of the albedo difference (Fig. 11) and (iii) characterise to what extent melt pond formation (Fig. 13) caused this albedo difference, we compare the downwelling with the absorbed SW radiation over ice-covered grid-cells. To directly compare the same ice-260 covered geographical region between the PI and LIG simulations, for any given day of the year, all quantities are computed using only cells that were ice-covered on this day of the year in both simulations. The absorbed SW radiation was further broken down into the fraction absorbed by open water, by melt ponds, and by bare ice and snow, using their respective area fractions of the grid-cell and estimating their respective albedos, as follows. For every grid-cell for every day of model data, the mean surface albedo was calculated as described in Section 2.
3. An open-water and pond albedo respectively of 0.07 and 265 0.27, alongside the area fraction of these two components for each cell (see Section 2.3), allowed the proportion of SW flux absorbed by each of these two components to be computed. The remainder of SW flux absorbed was attributed to exposed ice and snow.
The change in surface albedo (due to changes in coverage of bare ice, snow, ponds and open water) can be compared between the LIG and PI. For both time periods, North of the Equator, from January to June, the TOA downwelling SW flux increased.

270
This increased the surface downwelling SW flux, which led to an increase of the absorbed SW flux. Fig. 15a shows these changes in surface SW flux for the LIG: over ice-covered cells, an increase of downwelling flux from 96.8 ± 11.6 Wm −2 on April 1st, to 312 ± 27 Wm −2 on June 15th, led to the absorbed flux at the surface increasing from 24.0 ± 4.1 Wm −2 to 188 ± 27.4 Wm −2 . Thus, from April 1st to June 15th, the LIG surface albedo decreased by 35% (from 75% to 40%) over ice-covered regions. Similar computations for the PI (Fig. 15b) yield an albedo decrease of only 15%. Thus, albedo changes 275 were more significant for the LIG. In May, the gradient of the downwelling SW flux decreased, but the gradient of the absorbed SW flux increased, particularly for the LIG. This increasing rate of ice melt (see also Fig 1), despite the slowing rate of change  (blue), exposed ice and snow (green), and melt ponds (black), all computed from daily model output from the first 50 years from each of the PI and LIG simulations. Further details in Section 3.1.2 of downwelling solar radiation, implies that albedo feedbacks played an especially strong role in SW absorption through May during the LIG.
Using Fig 15c, we compare downwelling SW radiation between LIG and PI, and use this to quantify to what extent albedo 280 changes (related to albedo feedbacks) modified the surface absorbed SW anomaly shown in Fig. 7. Due to the TOA LIG-PI SW anomaly outlined in Section 1, the anomaly of the surface downwelling SW flux increased from 3.18 Wm −2 on April 1st to 29.9 Wm −2 on June 15th. This increased the anomaly of the absorbed SW radiation at the surface from 2.36 to 80.8 Wm −2 , as seen from Fig. 15c. From Fig 15c, it can be seen that the surface absorbed anomaly was 2.0 times the downwelling anomaly by May 24th, with a maximum of 4.61 times the downwelling anomaly on Jun 4th: the difference in surface albedo between the 285 LIG and PI caused the difference in downwelling radiation to be amplified fourfold. As the LIG-PI albedo difference increased so significantly from April through June, we have shown that the surface absorbed SW anomaly in Fig. 7a was caused by the TOA SW anomaly, very likely significantly amplified by stronger LIG than PI albedo feedbacks (as first suggested in Guarino et al. (2020b)).
Thus, for both time periods, from spring to summer: the radiative forcing triggered the albedo feedback, and the increase 290 in the radiative forcing continued to strengthen this feedback into the summer. However, for the LIG, the stronger radiative forcing amplified the albedo feedback more significantly, so that a much greater fraction of the downwelling SW radiation was absorbed in the Arctic. This significantly changed the Arctic heat budget and ultimately resulted in a complete loss of Arctic sea-ice by August.
Ice-albedo feedbacks result from changes in ice cover and pond formation. To determine to what extent melt pond formation 295 (Fig. 13) impacted the albedo difference shown in Figs 10 and 11, and caused the surface absorbed SW anomaly, we consider the LIG-PI anomaly of the SW radiation absorbed by each surface type. From Fig. 15c: the melt pond anomaly was comparable to the open water anomaly from May to June, and much greater than the ice and snow anomaly. In particular, from May 19th up until June 23rd, just as the LIG-PI anomaly of the downwelling surface SW flux reached its peak, the magnitude of the melt pond anomaly was at least 0.5 times the open water anomaly (and at least 1.4 times the ice anomaly). The melt pond 300 anomaly from May 25th (at 15 Wm −2 ) to June 17th (at 38 Wm −2 ) was greater than both the open water, and the ice and snow, anomalies. Therefore, the role of melt ponds in decreasing the LIG albedo, thus amplifying the surface absorbed SW flux anomaly, was particularly significant as the surface downwelling SW flux anomaly grew through May and peaked in June.
This demonstrates the significant impact of melt ponds on the surface energy balance, and by extension their key role in enhancing LIG summer sea-ice loss.

Melt ponds and sea-ice predictability
Today's diminishing Arctic sea-ice has led to a new focus on seasonal forecasting of Arctic sea-ice conditions, and especially predicting the minimum sea-ice area each year (Williams et al., 2016). Spring melt ponds are a good predictor for summer seaice conditions (Schröder et al., 2014). It is of interest to see how the predictability of spring melt pond area and August-October sea-ice area varied between the PI and LIG, since this may yield insight into how predictability may change in future under 310 conditions of reduced Arctic sea-ice. Predictability is investigated here by considering interannual variability within each of the LIG and PI periods. The interannual variability in the radiative forcing at the surface was much smaller than the difference in the radiative forcing between the PI and LIG. Thus, investigating the relation between spring melt pond area and August-October sea-ice area within each of these two periods gives insight into the impact of melt ponds on the summer sea-ice area for similar TOA radiative forcing and winter ice conditions each year.

315
For the 200 years' of simulation output for both LIG and PI, Pearson's correlation coefficient was calculated between the mean August to October sea-ice area, and each of the following four variables in each of April, May and June: (1) the mean monthly melt pond fraction (2) the 'radiation-effective pond fraction', denoting the fraction of grid-cell area covered by ponds that are not covered by an ice lid and thus are expected to affect the surface albedo (DuVivier, 2018)

320
(3) the thin-ice fraction, defined as fraction of ice-covered grid-cells with ice thickness < 1.4 m, as the ice state in spring is known to affect the summer sea-ice area (Schröder et al., 2014) (4) the fraction of downwelling short-wave radiation absorbed, as this fraction over ice-covered cells accounts for albedo changes from open water as well as ponds in these cells.
We note that, for both time periods, the April through June thin-ice fraction was a statistically significant (p< 10 −5 ) predictor 325 of summer ice area, as might be expected (Schröder et al., 2014). However, of more interest here is the significance of the pond-related correlations through the spring. Whilst not significant in April for either period, by May, there was a significant negative correlation between melt pond formation and summer sea-ice area, of -0.34 for the PI and -0.42 for the LIG. This is only a slightly weaker correlation than that of the absorbed short-wave fraction, of -0.37 for the PI and -0.46 for the LIG (Fig.   16). By June both of these correlations were much stronger than the thin-ice fraction correlation: the melt pond correlation 330 was -0.55 for the PI and -0.73 for the LIG, and the short-wave correlation was -0.55 for the PI and -0.76 for the LIG. These strong correlations show that spring melt pond formation alone was nearly as reliable a predictor of September ice cover, as the presence of any water (including both ponds and open water), over ice-covered grid-cells.
Compared to the PI, the LIG spring sea-ice area above 70°N was on the order of a million km 2 lower, and the mean sea-ice thickness is about 1m less. Therefore, the correlations shown in Fig 16 imply that May-June pond-related quantities can be 335 used to make more reliable seasonal predictions of summer sea-ice for thinner (LIG) than thicker (PI) sea-ice, as thinner ice is more sensitive to pond formation and albedo changes. This is of importance for future seasonal predictions: as the climate warms and Arctic sea-ice continues to thin, the spring melt pond area is likely to become an increasingly reliable predictor of September sea-ice area.

340
Despite their significant impact on sea-ice melt, melt ponds have only recently been explicitly included in global climate models (e.g. Collins et al., 2006;Curry et al., 2001). Prior to this, models without an explicit pond scheme had the ice albedo reduced when conditions approached near-melting (Hunke et al., 2015). Until now, it has not been clear whether this melt pond parameterisation approach adequately accounts for the impact of ponds during warmer-than-present day conditions (Kwok et al., 2009). 345 Guarino et al. (2020b) found that HadGEM3 simulated an ice-free summer LIG Arctic, and suggested this was linked to HadGEM3's realistic representation of melt pond physics. Here, we have quantified the importance of melt ponds during the LIG in depth, showing that greater melt pond formation, earlier in the year, directly led to the large surface energy budget difference between LIG and PI that was first highlighted in Guarino et al. (2020b). We have demonstrated the impact of melt ponds on thermodynamic processes during the LIG, and their role in determining LIG spring albedo over ice-covered regions, 350 leading to albedo feedbacks and ultimately the extreme ice loss first presented in Guarino et al. (2020b). Additionally, our analysis of the HadGEM3 PMIP4 LIG spin-up simulation years, showing that changes to longwave radiation and ocean heat transport were not primary drivers of the observed enhanced LIG ice loss, was supported by Guarino et al. (2020b), who used the LIG production run to show negligible differences in these factors between the simulated LIG and PI.
Of the CMIP6 models that have simulated the LIG according to CMIP6-PMIP4 experimental protocol, the two to include an 355 explicit melt pond scheme are the two with the lowest minimum monthly SIA: HadGEM3, 0.2 million km 2 , and CESM2, 1.2 million km 2 (Kageyama et al., 2021). The only other model that simulated a minimum LIG SIA less than 2.0 million km 2 was NESM3, at 1.3 million km 2 . NESM3 has an unrealistic sea-ice representation for the PI period, with amplitude of the seasonal cycle and winter sea-ice significantly overestimated (Kageyama et al., 2021), thus we disregard this model in our analysis.
For the remaining models which use parameterised rather than explicitly modelled ponds, the mean LIG minimum SIA was 360 3.5 ± 1.2 million km 2 . Since only 16 models have run this simulation, it is not possible to say that implementing an explicit melt pond scheme will tend to result in lower simulated LIG summer ice area; other aspects of the model simulation set-up may also be significant in the simulation of Arctic sea-ice. For example, the IPSL-CM6L model simulates some compensating ocean circulation and cloud cover changes that contribute to preserving LIG summer sea-ice (Kageyama et al., 2021). However, despite these caveats, it is striking that the only two models to include explicitly modelled melt ponds also simulated minimum 365 LIG SIA at least 2.3 million km 2 less than (i.e. 1.9 standard deviations from) the mean of the models which use parameterised ponds, whilst simulating a realistic PI SIA annual cycle, as well as maximum LIG SIA similar to (i.e. 0.5 standard deviations from) the mean of the models which use parameterised ponds.

Conclusions
In terms of understanding mechanisms driving possible Arctic sea-ice change during the LIG, Guarino et al. (2020b) indicated 370 that albedo changes were highly influential. Here we have provided a more in-depth follow-on study which analysed the potential thermodynamic and dynamic contributors, and provided multi-model context that was not available to Guarino et al. (2020b). This shows that melt pond formation has a crucial impact on sea-ice in warm climates, potentially making the difference between ice-covered and ice-free summer conditions. Specifically, we have answered the question of how HadGEM3's melt pond scheme contributed to the simulated ice loss in the Arctic during the LIG. 375 We have identified the key thermodynamic processes that led to the ice-free summer LIG Arctic as albedo feedbacks triggered by the radiative forcing. The TOA SW positive downwelling LIG-PI radiation anomaly in the spring and summer was halved by the time it reached the surface (Guarino et al., 2020b;Kageyama et al., 2021). However, we showed that the greater albedo feedback in the LIG compared to the PI (with LIG-PI surface albedo anomaly of <5% in April, growing to 30-35% by July) amplified this small surface anomaly by the summer by a factor of four. This meant up to 80Wm −2 more short-380 wave radiation was absorbed at the surface, leading to significant differences between the LIG and PI surface heat budgets, and explaining the greatly enhanced LIG sea-ice melt compared to the PI. We further demonstrated that ponds played a key role in reducing the albedo and strengthening the feedback process: through May and June, the downwelling surface SW flux anomaly peaked, and over ice-covered cells, melt ponds and open water accounted for a similar proportion of the surface absorbed anomaly. Therefore melt ponds contributed significantly to the simulated loss of summer sea-ice.

385
Our analysis of ice volume tendencies demonstrated that the difference in HadGEM3-simulated LIG and PI summer ice melt rates was predominantly driven by thermodynamic processes. This is interesting since atmospheric and ocean circulation changes often contribute to ice loss (Li et al., 2014;Köberle and Gerdes, 2003;Kageyama et al., 2021;Goosse and Fichefet, 1999); however this was not the case here. We found that for both the PI and LIG dynamic processes, driven by wind and ocean stress, led to less than 10cm/month of Arctic ice volume change through the spring and summer months. By contrast, 390 thermodynamic processes resulted in most ice volume change for both simulations, and accounted for the enhanced LIG ice loss. Three to five times more LIG sea-ice was lost than during the PI by June, in most Arctic regions: 70-110 cm/month of LIG sea-ice was lost due to thermodynamic processes, compared to 10-30 cm/month during the PI.
Given today's new focus on seasonal forecasting of Arctic sea-ice conditions, and especially predicting the minimum seaice area each year (Williams et al., 2016), we also investigated whether melt ponds are a good predictor for summer sea-ice 395 conditions (Schröder et al., 2014). Strong correlations between May-June melt pond area and August-October sea-ice fraction showed that explicitly modelling melt ponds significantly impacts summer sea-ice area. Much stronger correlations were found for the LIG, which had thinner ice that was thus more sensitive to melt pond formation than the PI. This is of concern for future seasonal sea-ice predictions: as Arctic ice continues to thin, the spring melt pond area each year may be an increasingly important and reliable indicator of the September sea-ice area.

400
In conclusion, whilst both models with both parameterised and explicitly modelled melt ponds are relatively successful in representing present-day sea-ice behaviour (Collins et al., 2006;Curry et al., 2001;Flocco et al., 2012), we have found that they likely simulate significantly different sea-ice behaviour under forcings other than the present-day. Multi-model context, alongside our new analysis above, suggests that a better representation of the contribution of melt ponds to enhanced Arctic sea-ice melt during the Last Interglacial is important. The relatively close match of HadGEM3 surface air temperatures to those 405 derived from proxy records (Guarino et al., 2020b), and expected match of CESM2 (see figures in Otto-Bliesner et al. (2020)), suggest that explicitly modelled pond formation for the LIG period does appear to be crucial to simulate realistic areas of summer sea-ice and Arctic temperature changes in current CMIP models. This is highly relevant to future projections of seaice loss, particularly when predicting the Arctic amplification of anthropogenic forcing; this requires accurate representation of albedo feedback mechanisms (Smith et al., 2019;Stuecker et al., 2018). Thus, our study of HadGEM3 supports the idea that 410 an explicit, realistic melt pond scheme is required for both past and future sea-ice and climate projections. Competing interests. The authors declare no competing interests.