# Brief communication: Representation of heat conduction into ice in marine ice shelf melt modelling

Jonathan Wiskandt, Nicolas C Jourdain

► To cite this version:

Jonathan Wiskandt, Nicolas C Jourdain. Brief communication: Representation of heat conduction into ice in marine ice shelf melt modelling. The Cryosphere, 2025, 19 (8), pp.3253 - 3258. 10.5194/tc-19-3253-2025. hal-05361273

HAL Id: hal-05361273

https://hal.science/hal-05361273v1

Submitted on 12 Nov 2025

HAL is a multi-disciplinary open access
archive for the deposit and dissemination of sci-
entific research documents, whether they are pub-
lished or not. The documents may come from
teaching and research institutions in France or
abroad, or from public or private research centers.

L'archive ouverte pluridisciplinaire HAL, est
destinée au dépôt et à la diffusion de documents
scientifiques de niveau recherche, publiés ou non,
émanant des établissements d'enseignement et de
recherche français ou étrangers, des laboratoires
publics ou privés.

<figure><img src="image_1.png" /></figure>

Distributed under a Creative Commons Attribution 4.0 International License
<figure><img src="image_2.png" /></figure>

<figure><img src="image_3.png" /></figure>

# Brief communication: Representation of heat conduction into ice in marine ice shelf melt modelling

Jonathan Wiskandt<fnref n="1" />,<fnref n="2" /> and Nicolas C. Jourdain<fnref n="3" />

**Correspondence:** Jonathan Wiskandt (jonathan.wiskandt@misu.su.se)

Received: 17 July 2024 – Discussion started: 12 August 2024

Revised: 11 June 2025 – Accepted: 12 June 2025 – Published: 25 August 2025

**Abstract.** Basal melt of marine-terminating glaciers is a key uncertainty in predicting the future climate and the evolution of the Antarctic and Greenland ice sheets. Regional ocean circulation models rely on an approximated heat budget at the ice-ocean interface that includes turbulent heat flux from the ocean below, latent heat for phase transition, and heat conduction into the ice in order to parameterize melt rates. Here we review the formulations of Holland and Jenkins (1999) for the approximation of heat conduction into ice, utilizing three distinct observed ice temperature profiles. We show that the formulation accounting for constant vertical advection and vertical diffusion of heat into ice best represents the observed ice temperature profiles. Using idealized ocean simulations of a Greenlandic fjord and observed melt rates from Antarctica, we find a difference of up to 28% (average of 12%) in melt rates between the different approximations.

## 1 Introduction

In ocean general circulation models (OGCMs), marine-terminating glaciers are usually assumed to be static and the local melt rate is parameterized using a set of three equations consisting of a linearization of the local freezing point temperature of water as well as heat and salt budgets at the ice-ocean interface (Holland and Jenkins, 1999). The turbulent transport of heat from water to the ice-ocean interface is partly balanced by the heat conduction into the ice. Any excess heat from the combination of these two processes is used as latent heat for the solid-to-liquid phase change. Hence, a

good representation of ice shelf basal melting in ocean models requires accurate simulations of both the turbulent heat transport in the water below the ice and the heat conduction into the ice shelf.

Heat conduction into ice has often been neglected in ocean simulations resolving ice shelf cavities, with the idea that the resulting bias does not affect melt rates by more than 10% (Dinniman et al., 2016; Holland and Jenkins, 1999). Other models have adopted a crude representation of this flux by assuming a linear temperature profile between the freezing temperature at the ice shelf base and a typical annual mean air temperature at the ice shelf surface (e.g. Losch, 2008; Mathiot et al., 2017). More than 60 years ago, Wexler (1960) proposed a more elaborate formulation of this conductive heat flux by accounting for heat advection by vertical ice motion, based on the observation of a few vertical temperature profiles. The formulation was completed by Holland and Jenkins (1999), who proposed a linear approximation of Wexler’s solution so that it can be incorporated into the solution of the three-equation system used to estimate basal melt rates in ocean models.

In this brief communication, we review these approaches and assess them using, for the first time, ice temperature profiles representing multiple typical scenarios of refreezing and melting. We further discuss the approaches’ accuracies by applying them to realistic applications and estimate the influence of heat conduction on ice shelf basal melt rates around Antarctica. Additionally, we investigate the relative importance of the heat conduction in the presence of a realistic ocean circulation and possible ice-ocean feedbacks. We

<footnote n="1">¹Department of Meteorology, Stockholm University, Stockholm, Sweden</footnote>

<footnote n="2">²Bolin Centre for Climate Research, Stockholm, Sweden</footnote>

<footnote n="3">³Univ. Grenoble Alpes/CNRS/IRD/G-INP/INRAE, IGE, Grenoble, France</footnote>
propose future directions to improve the calculation of these fluxes in ocean or Earth system models.

## 2 Assessment of the various approaches

The heat budget at the ice-ocean interface that is used to derive melt rates in ocean models can be expressed as

$$m L = \rho_w c_w u_* \Gamma_T (T_w - T_{z_d}) + \rho_i c_i \kappa_i \left(\frac{\partial T_i}{\partial z}\right)_{z_d}, \quad (1)$$

where m is the melt rate (kg m⁻² s⁻¹) (negative in the case of refreezing), ρw and ρi the densities of seawater and ice, L the melting or freezing latent heat of water, cw and ci the heat capacity of seawater and ice, u* the friction velocity calculated by the ocean model, ΓT the turbulent transfer parameter of temperature, κi the heat diffusivity of ice, Tj the ice temperature, Tw the ocean temperature at some distance from the ice-ocean interface, Tzd the local freezing point temperature, and zd the ice draft vertical position (z negative below sea level and increasing upward).

The three common approximations made in ocean models to deal with heat conduction in ice are

$$\rho_i \kappa_i \left(\frac{\partial T_i}{\partial z}\right)_{z_d} = 0 \qquad (A),$$

$$\rho_i \kappa_i \left(\frac{\partial T_i}{\partial z}\right)_{z_d} = \rho_i \kappa_i \frac{T_s - T_{z_d}}{z_s - z_d} \qquad (B), \quad (2)$$

$$\rho_i \kappa_i \left(\frac{\partial T_i}{\partial z}\right)_{z_d} = (T_s - T_{z_d}) \max\{m, 0\} \qquad (C),$$

where Ts is the ice surface temperature at height zs.

Approximation (A) in Eq. (2) would be a good approximation if there were ice shelves close to the ocean freezing point over their entire thickness, but such ice shelves are likely unstable (Morris and Vaughan, 2003). It can still be a reasonable approximation for locations where ice shelves have a layer of marine ice at their base (Fig. 1a), although the permissible nature of this type of ice makes proper calculations complicated (Wang et al., 2022).

Approximation (B) in Eq. (2) assumes no heat advected by the ice and a temperature profile that has reached a steady state between the freezing temperature at the base and a constant surface temperature. This seems to be a good approximation for the Ross Ice Shelf (RIS, Fig. 1b). The ice typically takes 1000 years to go from the grounding line to the borehole of the shown temperature profile (Mouginot et al., 2019), which corresponds to the typical time taken by heat to be conducted from the surface or base to the middle of a 400 m thick ice shelf (estimated as (H/2)²/κi, where H is the ice thickness). The ice may take less than a century to cross some of the smaller ice shelves of similar thickness (Mouginot et al., 2019), which makes it unlikely that a thermal steady state is reached for these other ice shelves. Furthermore, RIS is characterized by very low basal melt rates

(Rignot et al., 2013) and a relatively low accumulation rate at the surface (Kittel et al., 2021), which is expected to result in relatively low ice vertical velocities. Therefore, the thermal regime is very different for small ice shelves with high melt rates, such as the Pine Island Ice Shelf characterized by a strong temperature gradient at the base (Fig. 1c).

Approximation (C) in Eq. (2) is made to account for vertical heat advection and was obtained by Holland and Jenkins (1999). They considered a vertical ice velocity directly proportional to the basal melt rate. That is, they assumed that the melted ice was instantaneously compensated for by surface accumulation or ice flow convergence, which is correct for a constant ice shelf shape. The formulation of Eq. (2C) is a linearization of a more complex expression that tends to 0 for refreezing rates greater than ~0.1 m of ice per year (Holland and Jenkins, 1999). This formulation represents both a uniform ice temperature gradient throughout the ice shelf thickness associated with weak melt (Fig. 1b) and strong ice temperature gradients at the ice shelf base associated with high melt rates (Fig. 1c). Here, a steady state is still assumed, but the characteristic time to reach the steady state is imposed by the basal melt rate and not by heat conduction (see Supplement for a detailed discussion). The typical aspect ratio and horizontal ice velocities make this approximation reasonable for many ice shelves of Antarctica and Greenland (Rignot et al., 2013; Mouginot et al., 2019).

When considering the variety of ice temperature profiles found in Antarctica (Fig. 1a-c), it appears that approximation (A) only works in the presence of refreezing (Fig. 1a), while approximation (B) only works for large ice shelves with very weak basal melt rates (Fig. 1b). Approximation (C) appears to be the only one able to represent the variety of temperature profiles. In the case of refreezing (m < 0 in Eq. 2C), no heat is conducted into the ice (Fig. 2a), which is consistent with the vertical temperature profiles observed in the marine ice layer beneath the Amery Ice Shelf (Fig. 1a). In the case of actual melting (m > 0 in Eq. 2C), the conductive heat flux is proportional to the melt rates, consistent with the strong temperature gradients observed at the base of Pine Island Ice Shelf. The range of very weak melt rates (-0.05 < m < 0.05 m yr⁻¹) is the only one where the linear temperature profile is a better approximation than the linearization of Holland and Jenkins (1999) (Fig. 2a).

If we assume positive melt rates and use Eq. (1), we can estimate the relative difference between approximation (A) and approximation (C):

$$\frac{\Delta m}{m} = \frac{c_i (T_s - T_{z_d})}{L}. \qquad (3)$$

Considering the parameter values used in Holland and Jenkins (1999), melt rates calculated without considering heat conduction in the ice are overestimated by 11% for a surface temperature of -20 °C and a freezing point temperature of -2 °C at the ice base. The overestimation is often
<figure><img src="image_4.png" /><figcaption>**Figure 1.** Ice temperature profiles measured by thermistor strings in boreholes at (a) two sites on the Amery Ice Shelf (69.4° S, 71.4° E and 69.9° S, 70.3° E, referred to as AM01 and AM04, respectively), (b) one site in the middle of the Ross Ice Shelf (82.4° S, 168.6° W, referred to as J9), and (c) three sites on Pine Island Ice Shelf (75.1° S, 100.5° W). The black dotted and dashed lines in panels (b) and (c) show temperature profiles approximated from Eqs. (2B) and (2C) using observed melt rates and the annually averaged surface temperature. A layer of marine ice is present underneath the ice shelf at sites AM01 and AM04 (more than 100 m thick) and at site J9 (6 m thick).</figcaption></figure>

similar when using the linear temperature profile rather than nothing, except in the absence of basal melting.

Applying Eq. (3) to the observed melt rates in Antarctica gives an idea of where differences in the heat flux estimation have the largest effect on melt (Fig. 3). We use the 40-year (1980–2019) average snow surface temperatures as $T_s$ and calculate the pressure-dependent freezing point temperature $T_{zd}$ using a typical salinity of 34.6 kg m$^{-1}$ everywhere. Ice shelves in the Amundsen embayment and the large ice shelves show a difference of around 12%, with even higher values seen at the Amery Ice Shelf (AmIS) and the eastern RIS (both above 27%). Under most of RIS this leads to only small differences in absolute numbers (below 0.1 m yr$^{-1}$, Fig. 3b), since melt rates there are very low (Rignot et al., 2013). In crucial areas, such as close to the grounding lines of the Filchner–Ronne Ice Shelf (FRIS), most of AmIS, and

the smaller ice shelves draining towards the Amundsen Sea, the difference is several metres per year, surpassing even $10 \text{ m yr}^{-1}$ at Pine Island Glacier and AmIS. At the calving fronts of FRIS and RIS, melt rate differences are around $1 \text{ m yr}^{-1}$.

### 3 Application to ocean simulations of a small cavity

For a practical implementation of approximation (C) in an ocean model, the criterion on $m$ to bound the heat conduction to zero is not very convenient, as this would require some iterations of the three-equation system because $m$ is a solution that is not known a priori. A more practical solution is to check the sign of the thermal forcing ($T_w - T_{zd}$ in Eq. 1) to decide whether or not heat conduction is set to zero. This is a common procedure for implementing the three equations
<figure><img src="image_5.png" /><figcaption>**Figure 2.** (a) Conductive heat flux at the ice shelf base as a function of the basal melt rate expressed in metres of ice per year. The inset shows the same lines but with a broader range of melt rates. The orange line is the full analytical solution provided by Holland and Jenkins (1999), the dotted blue line is the linear simplification that they suggest (approximation (C) in Eq. 2), and the green line is derived from the assumption of a linear temperature profile across the ice shelf thickness (approximation (B) in Eq. 2). (b) Simulated average melt rate from identical MITgcm ocean simulations of an idealized, 2D Greenland fjord-like domain (Wiskandt et al., 2023) using the three different approximations of Eq. (2) as a function of bottom-layer ocean temperature. (c) Melt rate difference with respect to approximation (C).</figcaption></figure>

<figure><img src="image_6.png" /><figcaption>**Figure 3.** Map of Antarctica showing relative melt rate differences (Eq. 3) between approximations (A) and (C) and respective absolute values when applied to melt rate data for floating ice shelves (b). White areas indicate regions of refreezing.</figcaption></figure>
in ocean modelling. We test this implementation in the MITgcm ocean model (Adcroft et al., 2004) and evaluate whether the relative importance of heat conduction, as derived from theory, is still valid in the presence of a realistic ocean circulation and possible ice-ocean feedbacks.

To show the differences between the three approaches in the melt parameterization in a regional circulation model, we set up a suite of idealized, 2D, and non-rotational simulations (Wiskandt et al., 2023) in the MITgcm with varying ocean temperatures comparing the three different approaches introduced in Eq. (2). The geometry is representative of a long and narrow fjord typically found around Greenland. The domain is 30 km long and 1 km deep with an ice shelf covering the first 20 km of the domain, a grounding line depth of 950 m, and a 50 m deep vertical front. Note that none of the present simulations exhibits refreezing. For all details of the simulation, please refer to Wiskandt et al. (2023).

Given the constant small difference between the conductive heat fluxes from approximations (A) and (B) (Fig. 2a), the simulated melt rates are very similar in these two approximations (Fig. 2b and c). The differences in absolute melt rate are barely discernible to the eye, and the relative differences are < 1 % for the simulated temperature range, except for the coldest experiment (Fig. 2c). The difference between approximation (C) and either (A) or (B) is between 10 % and 15 % in most of the experiments. For example, those with temperatures ≥ 0 °C show a difference of 14 %–15 %.

## 4 Conclusions

Previous ocean circulation studies at both Greenlandic and Antarctic marine-terminating glaciers implemented approximation (A) (e.g. Gwyther et al., 2020; Comeau et al., 2022) or approximation (B) (e.g. De Rydt and Naughten, 2024; Wiskandt et al., 2023; Holland et al., 2008; Losch, 2008) in Eq. (2) for heat conduction into ice. Combining for the first time a suite of ice temperature profiles representing a variety of melting conditions, we have shown that only approximation (C) estimates the temperature gradient in the ice accurately across the full range of possible melting and refreezing scenarios. We provide a novel map of the local impact of heat conduction parameterization on melt rates in Antarctica and find that approximation (A) overestimates melt rates by up to 25 % (12 % average) compared to approximation (C). Based on Fig. 3a, we can expect a similar overestimation of melt from approximation (B).

We have shown that an inaccurate representation of heat conduction into ice induces biases in melt rates with a high spatial variability. This cannot be corrected by tuning a constant drag or heat transfer coefficient as often done in ocean models. Accurate local melt rates are also crucial for estimating the stability of the ice shelf, e.g. close to the grounding line, where changes in the ice shelf geometry can trigger rapid ice loss (De Rydt and Naughten, 2024). In ocean cir-

culation models, a melt-dependent local temperature gradient is relatively straightforward to implement and, in fact, is available in the MITgcm in the SHELFICE package (Losch, 2008). The Finite-volumE Sea ice-Ocean Model (FESOM2) (Danilov et al., 2017) also uses an approach similar to approximation (C).

We argue here that the parameterization using a linear dependency of the heat flux on the melt rate is the most accurate one. It is nonetheless based on several approximations (Holland and Jenkins, 1999). To accurately model the ice-ocean interactions, a coupled ice sheet-ocean model should be used, with a representation of both vertical and lateral heat conduction and heat advection in the ice (Sergienko et al., 2013). A large number of the current ice sheet models do not however explicitly simulate the ice sheet temperatures and sometimes use approximations of the Stokes equation that prevent a good representation of vertical heat advection in ice. Hence, we suggest that coupled models currently under development should consider calculating the ice sheet temperature gradients or heat fluxes and transfer them to the ocean model to accurately calculate the melt rate using the three-equation system. For non-coupled models, the authors strongly suggest implementing approach (C) of Eq. (2) when parameterizing marine basal melt, as it more accurately represents the heat conduction into the ice and comes at very little additional computational cost.

**Code and data availability.** The configuration files necessary to reproduce the simulations used in the present study employed the MITgcm (downloaded in 2020 from https://github.com/MITgcm/MITgcm/releases/tag/checkpoint67s, last access: 4 July 2023, Adcroft et al., 2004) and are made available through the Bolin Research Centre's Data Centre (https://git.bolin.su.se/bolin/wiskandt-2024-ice-shelf-heat-conduction, last access: 14 August 2025, Wiskandt and Jourdain, 2025). Time-averaged fields and key diagnostics from the simulations and code needed to reproduce the shown figures are also available. See the Supplement for a detailed account of all of the data used.

**Supplement.** The supplement related to this article is available online at https://doi.org/10.5194/tc-19-3253-2025-supplement.

**Author contributions.** JW and NCJ initiated this work and reviewed the theory together. JW ran the MITgcm simulations and analysed the results. All of the authors discussed the results and contributed to the manuscript.

**Competing interests.** At least one of the (co-)authors is a member of the editorial board of *The Cryosphere*. The peer-review process was guided by an independent editor, and the authors also have no other competing interests to declare.
*Disclaimer.* Publisher's note: Copernicus Publications remains neutral with regard to jurisdictional claims made in the text, published maps, institutional affiliations, or any other geographical representation in this paper. While Copernicus Publications makes every effort to include appropriate place names, the final responsibility lies with the authors.

*Acknowledgements.* The MITgcm simulations were enabled by resources provided by the National Academic Infrastructure for Supercomputing in Sweden (NAISS) at the National Supercomputer Centre (NSC) and partially funded by the Swedish Research Council through grant no. 2022-06725. The authors would furthermore like to thank the four reviewers for their insightful contributions to improving this paper.

*Financial support.* This research has been supported by the Stockholms Universitet (grant no. SUFV-1.2.1-0124-17), the European Horizon's 2020 programme (grant no. 101003536), and the Agence Nationale de la Recherche (grant no. ANR-22-EXTR-0010).

The publication of this article was funded by the Swedish Research Council, Forte, Formas, and Vinnova.

*Review statement.* This paper was edited by Ed Blockley and reviewed by four anonymous referees.

## References

Adcroft, A. J., Hill, C., Campin, J. M., Marshall, J., and Heimbach, P.: Overview of the formulation and numerics of the MIT GCM, in: Proceedings of the ECMWF seminar series on Numerical Methods, Recent developments in numerical methods for atmosphere and ocean modelling, ECMWF, https://www.ecmwf.int/en/elibrary/7642-overview-formulation-and-numerics-mit-gcm (last access: 4 July 2023), 139–149, 2004.

Comeau, D., Asay-Davis, X. S., Begeman, C. B., Hoffman, M. J., Lin, W., Petersen, M. R., Price, S. F., Roberts, A. F., Van Roekel, L. P., Veneziani, M., Wolfe, J. D., Fyke, J. G., Ryngler, T. D., and Tuner, A. T.: The DOE E3SM v1.2 cryosphere configuration: description and simulated Antarctic ice-shelf basal melting, J. Adv. Model. Earth Sy., 14, e2021MS002468, https://doi.org/10.1029/2021MS002468, 2022.

Danilov, S., Sidorenko, D., Wang, Q., and Jung, T.: The Finite-volumE Sea ice–Ocean Model (FESOM2), Geosci. Model Dev., 10, 765–789, https://doi.org/10.5194/gmd-10-765-2017, 2017.

De Rydt, J. and Naughten, K.: Geometric amplification and suppression of ice-shelf basal melt in West Antarctica, The Cryosphere, 18, 1863–1888, https://doi.org/10.5194/tc-18-1863-2024, 2024.

Dinniman, M. S., Asay-Davis, X. S., Galton-Fenzi, B. K., Holland, P. R., Jenkins, A., and Timmermann, R.: Modeling ice shelf/ocean interaction in Antarctica: a review, Oceanography, 29, 144–153, 2016.

Gwyther, D. E., Kusahara, K., Asay-Davis, X. S., Dinniman, M. S., and Galton-Fenzi, B. K.: Vertical processes and resolution impact ice shelf basal melting: a multi-model study, Ocean Model., 147, 101569, https://doi.org/10.1016/j.ocemod.2020.101569, 2020.

Holland, D. M. and Jenkins, A.: Modeling thermodynamic ice-ocean interactions at the base of an ice shelf, J. Phys. Oceanogr., 29, 1787–1800, https://doi.org/10.1175/1520-0485(1999)029<1787:mtioia>2.0.co;2, 1999.

Holland, P. R., Jenkins, A., and Holland, D. M.: The response of Ice shelf basal melting to variations in ocean temperature, J. Climate, 21, 2558–2572, https://doi.org/10.1175/2007JCLI1909.1, 2008.

Kittel, C., Amory, C., Agosta, C., Jourdain, N. C., Hofer, S., Delhasse, A., Doutreloup, S., Huot, P.-V., Lang, C., Fichefet, T., and Fettweis, X.: Diverging future surface mass balance between the Antarctic ice shelves and grounded ice sheet, The Cryosphere, 15, 1215–1236, https://doi.org/10.5194/tc-15-1215-2021, 2021.

Losch, M.: Modeling ice shelf cavities in a z coordinate ocean general circulation model, J. Geophys. Res.-Oceans, 113, C08043, https://doi.org/10.1029/2007JC004368, 2008.

Mathiot, P., Jenkins, A., Harris, C., and Madec, G.: Explicit representation and parametrised impacts of under ice shelf seas in the z\* coordinate ocean model NEMO 3.6, Geosci. Model Dev., 10, 2849–2874, https://doi.org/10.5194/gmd-10-2849-2017, 2017.

Morris, E. M. and Vaughan, D. G.: Spatial and Temporal Variation of Surface Temperature on the Antarctic Peninsula And The Limit of Viability of Ice Shelves. In Antarctic Peninsula Climate Variability: Historical and Paleoenvironmental Perspectives, edited by: Domack, E., Levente, A., Burnet, A., Bindschadler, R., Convey P., and Kirby, M., https://doi.org/10.1029/AR079p0061, 2003.

Mouginot, J., Rignot, E., and Scheuchl, B.: Continent-wide, interferometric SAR phase, mapping of Antarctic ice velocity, Geophys. Res. Lett., 46, 9710–9718, 2019.

Rignot, E., Jacobs, S., Mouginot, J., and Scheuchl, B.: Ice-shelf melting around Antarctica, Science, 341, 266–270, 2013.

Sergienko, O. V., Goldberg, D. N., and Little, C. M.: Alternative ice shelf equilibria determined by ocean environment, J. Geophys. Res.-Earth, 118, 970–981, https://doi.org/10.1002/jgrf.20054, 2013.

Wang, Y., Zhao, C., Gladstone, R., Galton-Fenzi, B., and Warner, R.: Thermal structure of the Amery Ice Shelf from borehole observations and simulations, The Cryosphere, 16, 1221–1245, https://doi.org/10.5194/tc-16-1221-2022, 2022.

Wexler, H.: Heating and melting of floating ice shelves, J. Glaciol., 3, 626–645, 1960.

Wiskandt, J. and Jourdain, N.: Code for a study of heat conduction into the ice in marine ice shelf melt modelling, Software version 1.1.0, Bolin Centre Code Repository [code], https://doi.org/10.57669/wiskandt-2024-ice-shelf-heat-conduction-1.1.0, 2025.

Wiskandt, J., Koszalka, I. M., and Nilsson, J.: Basal melt rates and ocean circulation under the Ryder Glacier ice tongue and their response to climate warming: a high-resolution modelling study, The Cryosphere, 17, 2755–2777, https://doi.org/10.5194/tc-17-2755-2023, 2023.