Comparison of calibration methods of a PICO basal ice shelf melt module implemented in the GRISLI v2.0 ice sheet model

Uncertainties in future sea level rise are mainly due to uncertainties in Antarctic ice sheet projections. Indeed, modelling the future of the Antarctic ice sheet presents many challenges. One of them is being able to model the physical interactions between the ocean and the ice shelves. As a result of technical challenges related to computational resources, implementation, and different modelling timescales, these interactions are often parameterised rather than explicitly resolved in ice sheet models. These parameterisations vary in complexity and calibration method, eventually leading to differences in resulting sea level rise contribution of several metres. Here we present the implementation of the Potsdam Ice-shelf Cavity mOdel (PICO) basal ice shelf melt module in the GRISLI v2.0 ice sheet model. We compare six different statistical methods to calibrate PICO and assess how robust these methods are if applied at different resolutions and areas of the Antarctic ice sheet. We show that computing the mean absolute error of the bins is the best method as it allows us to match the entire distribution of melt rates retrieved from satellite data at different resolutions as well as for different Antarctic ice shelves. It also results in a smaller parameter space than the other tested methods. This method makes use of melt rate bins and minimises the differences between the values of the bins of the model and the ones of the observational target. It gives equal weight to the full distribution of melt values: low, medium, and high values. We find that, using this method, region-specific calibration of ice–ocean interactions is not needed and we can avoid using ocean temperature bias corrections. Finally, we assess the impact of the implementation of PICO in GRISLI and of the calibration choice on future projections of the Antarctic ice sheet up to the year 2300.

1 Introduction

The future evolution of the Antarctic ice sheet is the largest uncertainty in sea level rise projections for the end of the century (Edwards et al., 2021). The mass loss of the Antarctic ice sheet is primarily driven by basal melting of ice shelves (Pritchard et al., 2012). The ice shelves have a buttressing effect, slowing the ice flow towards the ocean (Dupont and Alley, 2005). Their thinning observed over the last decades (Rignot et al., 2013; Paolo et al., 2015; Adusumilli et al., 2020a) is due to increased heat provided by circumpolar deep water in the cavities directly beneath the ice shelves (Schmidtko et al., 2014; Stewart and Thompson, 2015; Jenkins et al., 2018). Subsurface melt of the ice shelves impacts the oceanic circulation in the ice cavities as well as larger-scale oceanic circulation (Bennetts et al., 2024). These feedbacks and the large range of spatiotemporal scales at play, from turbulence to large-scale ocean circulation, make the ice–ocean interaction a complex process challenging to model accurately (Bennetts et al., 2024). Additionally, on a retrograde bathymetry, such as in West Antarctica, the thinning of ice shelves and retreat of the grounding line can trigger marine ice sheet instabilities (Weertman, 1974; Schoof, 2007), leading to irreversible commitment to sea level rise. Hence, understanding and having the ability to model the ice–ocean interaction accurately is crucial to constrain uncertainties of projections of the future rise in sea level.

With our current computational resources, it is necessary to use parameterisations in ice sheet models to compute the physical interactions between the ocean and the ice. Over the last decade, several basal melt parameterisations have been developed and implemented in ice sheet models with different complexities in the simplification of the ocean circulation beneath the ice shelves and its physical interaction with the ice (Reese et al., 2018a; Lazeroms et al., 2019; Pelle et al., 2019; Jourdain et al., 2020; Lambert et al., 2023). Berends et al. (2023) demonstrated that the choice of the sub-shelf melt parameterisation has a strong impact on the Antarctic ice sheet retreat for idealised as well as realistic geometries. Moreover, the values of the parameters are poorly constrained, and in some cases the parameterisations require ocean temperature corrections up to 2 K to be able to match the basal melt rates observed (Jourdain et al., 2020; Reese et al., 2023).

Here, we present the implementation of the Potsdam Ice-shelf Cavity mOdel (PICO) (Reese et al., 2018a) in the GRISLI v2.0 ice sheet model (Quiquet et al., 2018) but also a comparison of methodologies to calibrate PICO. We aim at calibrating the module to match the whole distribution of values from the observations to the extent possible. In this research no ocean temperature corrections are added in the calibration process in order to have a more physical relationship between the forcing and the computed melt rates. The article first presents the methodology in Sect. 2, then the results in detail in Sect. 3 including comparison of calibration methods, sensitivity estimates, and future projections. Section 4 discusses the limitations and possible improvements for future studies. Finally, Sect. 5 concludes and gives some perspectives.

2 Methodology

2.1 PICO basal ice shelf melt module

PICO is a parameterisation that computes the basal melt rates under the ice shelves. It is described in detail in Reese et al. (2018a). We present here only the key concepts of PICO. It is a box model, based on the work of Olbers and Hellmer (2010). The ice shelves are divided into boxes, and the shape and number of boxes in one ice shelf depend on two variables: distance to the ice shelf front and distance to the grounding line. The number of boxes n_D in one ice shelf D is then defined by

$$\begin{array}{}\text{(1)}& {\displaystyle }{\displaystyle }{n}_D=\mathrm{1}+\text{rd}\left(\sqrt{d_{\text{GL}}\left(D\right)/d_{max}}\left(n_{max}-\mathrm{1}\right)\right),\end{array}$$

with d_GL the distance of each grid cell to the grounding line, d_max the maximum distance between the grounding line and the ice shelf front among all the ice shelves of the ice sheet, and n_max the maximum number of boxes kept here at 5.

PICO accounts for one-dimensional overturning circulation in ice shelf cavities (Lewis and Perkin, 1986). The overturning flux under the ice shelf q is driven by the density difference between the ocean box B0 (ρ₀) and the first box B1 under the ice shelf at the level of the grounding line (ρ₁).

$$\begin{array}{}\text{(2)}& {\displaystyle }{\displaystyle }q=C\left({\mathit{\rho }}_{\mathrm{0}}-{\mathit{\rho }}_{\mathrm{1}}\right)\end{array}$$

The value of q must be greater than zero. A single overturning flux value is calculated for all boxes of the same ice shelf. The constant overturning coefficient C (Sv m³ kg⁻¹) captures effects due to friction, rotation, and bottom form stress; more details are given in Olbers and Hellmer (2010). C is one of the two PICO parameters that we calibrate in the present study. To compute basal melt rates m_k in the box B_k, PICO requires two ocean inputs (ocean temperature T_k−1 and salinity S_k−1) and one ice sheet input: the ice draft to calculate the pressure at the ice–ocean interface under the ice shelf p_k using $p_k={\mathit{\rho }}_{\text{SeaWater}}\cdot g\cdot z_{\text{IceDraft}}$, with ρ_SeaWater=1033 kg m⁻³. For the box B₁, the ocean inputs are the average temperature (T₀) and salinity (S₀) at the continental shelf depth in front of the corresponding ice shelf. For the next boxes B_k, the forcing temperature and salinity depend on the overturning q and the temperature and salinity computed for the previous box T_k−1 and S_k−1. The details of the analytical derivation are given in Reese et al. (2018a) in Appendices A and B. The melt rate in the box k is then computed as follows:

$$\begin{array}{}\text{(3)}& {\displaystyle }{\displaystyle }{m}_k(x,y)=-{\displaystyle \frac{{\mathit{\gamma }}_T^{*}}{\mathit{\nu }\mathit{\lambda }}}\left(a{S}_{k-\mathrm{1}}+b-c{p}_k(x,y)-T_{k-\mathrm{1}}\right),\end{array}$$

where ${\mathit{\gamma }}_T^{*}$ is the heat exchange coefficient (m s⁻¹), the second of the two PICO parameters calibrated in this paper, $\mathit{\nu }={\mathit{\rho }}_i/{\mathit{\rho }}_w\sim \mathrm{0.89}$, $\mathit{\lambda }=L/c_p\sim \mathrm{84}$ °C. The coefficients from linearisation of the equation of state for the freezing point of seawater are as follows: a is the liquidus slope coefficient and equals −0.0572 °C PSU⁻¹, b is the liquidus intercept coefficient and equals 0.0788 °C, and c is the liquidus pressure coefficient and equals $\mathrm{7.77}\times {\mathrm{10}}^{-\mathrm{8}}$ °C Pa⁻¹.

2.2 Choices for the implementation of PICO in GRISLI

To implement PICO in GRISLI and to select the values of the two parameters (C and ${\mathit{\gamma }}_T^{*}$) we make some choices that differ from Reese et al. (2018a). We give an overview of these choices here. All the PICO implementation has been done in Fortran90 to correspond to GRISLI development language.

First, the heaviest computation part in the module is the computation of the geometry of the ice shelves. To reduce computational costs we avoid identifying individual ice shelves at each time step. Instead we decide to solve one instance of the PICO equations per drainage basin, with the basins as presented in Mouginot et al. (2017) and used in Rignot et al. (2019), with the difference that we combine the two drainage basins of the two largest ice shelves as done by Jourdain et al. (2020): Ronne and Filchner ice shelves as well as Ross East and Ross West ice shelves. Each drainage basin defines the external borders of the geometry of the ice shelves and is used to compute the oceanic forcing inputs. This implies that if two different ice shelves are in the same drainage basin, they are seen as one ice shelf for PICO. Inversely, if one ice shelf has two drainage basins, it is then seen as two separate ice shelves for PICO. The number of boxes in each drainage remains relative to the maximum distance between the ice shelf front and the grounding line of all the ice shelves, as defined in the Eq. (1). And since within the same drainage basin there are roughly similar sizes of ice shelves, the division by drainage basin does not cause discrepancies such as small ice shelves with up to five boxes or larger ice shelves with few boxes (see Fig. S1 in the Supplement). In all cases, the ice shelf front and the grounding line are defined in the same way by being neighbours to open ocean or grounded grid cells, respectively. The geometry of the ice shelves is recomputed at every time step to readjust the boxes to the changing grounding line and ice shelf front position. This simplification of solving one PICO instance per drainage basin enables us to compute the number of boxes and their corresponding temperature and salinity input faster for each ice shelf.

Reese et al. (2018a) used four selection criteria to calibrate PICO and define the values of the PICO parameters C and ${\mathit{\gamma }}_T^{*}$. The two first criteria are (1) to not have freezing dominating the melt rate values in the first ice shelf box and (2) the overall mean basal melt rates must decrease between the first and second box of the ice shelf. Criteria (3) and (4) are constraints on the average values the melt rates that should be in the cold ice cavities of Filchner–Ronne Ice Shelf (FRIS) and in the warm ice cavities of Pine Island glacier, respectively. Here, we do not follow any of these criteria; we apply different criteria. The rationale is that while the PICO equations assume that (1) and (2) are true for the one horizontal dimensional case, melt rate patterns are complex in two horizontal dimensions. Indeed, the retrieved basal melt rates from remote sensing (Adusumilli et al., 2020a; Paolo et al., 2023) show refreezing in some areas close to the grounding line but also higher melt rates close to the ice shelf front. Instead of criteria (3) and (4), the calibration methods we tested here (presented in Sect. 2.5) are designed to be able to capture the whole distribution of values, not only the average values, at an Antarctic-wide scale as well as at an ice shelf scale.

The grounding line can be defined in different ways and can therefore lead to different PICO box geometries. Here, in GRISLI-PICO, we consider any ice point to be on the grounding line if it has some neighbours that are grounded and others that are floating, and we consider ice front to be any ice point that is floating and adjacent to ocean. In contrast, in PISM-PICO Reese et al. (2018a) did not include the grounding line of ice rises and also excluded holes in ice shelves as ice shelf front when identifying PICO boxes. The grounding lines of ice rises are defined as not being directly connected to the main grounded part of the ice sheet, which is identified by the size of the connected grounded region. Thus, in PISM-PICO it is possible to have ice shelves without a grounding line connected to the main ice sheet, where PICO cannot define a box geometry. In these places, the parameterisation of Beckmann and Goosse (2002) was used in PISM to have a rough estimate of the basal melt rates.

GRISLI incorporates a dynamic calving front that advances based on a balance between the Lagrangian ice flux and local surface and basal mass balances (Quiquet et al., 2018). To evaluate potential ice shelf advance at each time step, the model must compute these mass balances even beyond the current ice extent. Unlike alternative approaches such as the level set method, which do not require mass balance information outside the ice mask, this is a necessary feature for GRISLI. Thus, in regions beyond the ice shelf, over open ocean, we apply the parameterisation of DeConto and Pollard (2016), defined as follows.

$$\begin{array}{}\text{(4)}& {\displaystyle }{\displaystyle }m={\displaystyle \frac{K_T{\mathit{\rho }}_w{C}_w}{{\mathit{\rho }}_i{L}_{\mathrm{f}}}}\left|T_{\mathrm{o}}-T_{\mathrm{f}}\right|\left(T_{\mathrm{o}}-T_{\mathrm{f}}\right)\end{array}$$

The main difference is the inclusion of a quadratic dependence between the melt rate and the difference of the temperature between the ocean T_o and the ocean freezing point at the ice base T_f. This quadratic relation enables us to limit the growth of the ice shelves towards the ocean. The combined factor $\frac{K_T{\mathit{\rho }}_w{C}_w}{{\mathit{\rho }}_i{L}_{\mathrm{f}}}$ equals 0.224 $\mathrm{m}{\mathrm{yr}}^{-\mathrm{1}}\mathrm{°}{\mathrm{C}}^{-\mathrm{2}}$ as we keep the same K_T value of 15.77 $\mathrm{m}{\mathrm{yr}}^{-\mathrm{1}}\mathrm{°}{\mathrm{C}}^{-\mathrm{1}}$ from DeConto and Pollard (2016) and Pollard and Deconto (2012). The temperature input T₀ is computed the same way as for PICO. The parameterisation chosen for the open ocean does not impact the calibration results but does impact the transient ice sheet simulations.

Finally, in GRISLI v2.0 the iceberg calving is defined by a simple ice thickness threshold criterion (Quiquet et al., 2018). The threshold value varies in space and time as it is dependent on the depth of the bathymetry at the location of the ice shelf front.

2.3 Calibration ensemble

To calibrate the two PICO parameters C and ${\mathit{\gamma }}_T^{*}$ we run an ensemble of 169 members of PICO implemented in GRISLI corresponding to all possible combinations between 13 values for the parameter C (ranging from 0.01 to 15.00 Sv m³ kg⁻¹) and 13 values for the parameter ${\mathit{\gamma }}_T^{*}$ (from $\mathrm{0.01}\times {\mathrm{10}}^{-\mathrm{5}}$ to $\mathrm{15.00}\times {\mathrm{10}}^{-\mathrm{5}}$ m s⁻¹). The range of values for the two parameters has been chosen based on the literature (Reese et al., 2018a, 2023; Burgard et al., 2022) and adjustments in such a way that the best values are not at one of the extremes of the range of values. The geometry of the ice sheet and the ice shelves is kept fixed to remove the influence of ice shelf geometry changes on the computed basal melt rate. The fixed geometry corresponds to Bedmap2 (Fretwell et al., 2013) with a 30-year relaxation with GRISLI. This relaxation is needed because with the basal drag coefficient inversion methodology used for ISMIP6, we compute the ice sheet internal thermal equilibrium with a long (60 kyr) experiment with fixed observed geometry. Thus, to avoid any artificial drift when releasing this constraint we run a 30-year relaxation experiment with the same boundary conditions as for the control experiment from ISMIP6 (Seroussi et al., 2024). Doing the calibration of PICO in a coupled GRISLI-PICO with fixed geometry enables us to facilitate the transition between the PICO calibration and the GRISLI-PICO transient experiments (see Sect. 2.6) without impacting the results. With this ensemble of 169 simulations, we apply six different methods, presented in Sect. 2.5, to evaluate which members of the ensemble provide the best fit with respect to the observational dataset.

2.4 Data: ocean forcing and basal melt rate target

The oceanic forcing we use for the calibration ensemble presented above in Sect. 2.3 is the dataset produced by Jourdain et al. (2020). This is a present-day estimate of three-dimensional fields of temperature and salinity of the ocean surrounding the Antarctic ice sheet. Jourdain et al. (2020) computed this estimate by using the following datasets: a pre-release of the NOAA World Ocean Atlas 2018 covering the period 1995–2017 (Locarnini et al., 2018; Zweng et al., 2019), the Met Office EN4 subsurface ocean profiles for the period 1995–2017 (Good et al., 2013), and Marine Mammals Exploring Oceans from Pole to Pole for the period 2004–2018 (Treasure et al., 2017). The final dataset created by Jourdain et al. (2020) includes extrapolation of the ocean properties into the ice shelf cavities where observations are not available. The end product is on a polar stereographic grid with a resolution of 8 km horizontally and 60 m vertically.

Our target is the average basal melt rates retrieved by Adusumilli et al. (2020a), a dataset that used CryoSat-2 altimetry to create an average value estimate of basal melt rates of the ice shelves of Antarctica for the period 2010–2018 at a resolution of 500 m (Adusumilli et al., 2020b).

These ocean forcings and basal melt rates differ from Reese et al. (2018a) where they used Schmidtko et al. (2014) for the ocean forcing and Rignot et al. (2013) as the target for the basal melt rates. The selected datasets in the present study are more up to date and have good overlap with the time period of data retrieval between the forcings and the target. For all the methods and analysis, all the datasets, forcings, and the observational target are scaled up to the same resolution as the ice sheet model (16 or 40 km) using Climate Data Operator (CDO) bilinear interpolation.

2.5 Six statistical methods of calibration of the two PICO parameters: C and ${\mathit{\gamma }}_T^{*}$

The six statistical methods of calibration compared in this study are explained here. The overview of the methods is given in Table 1, including names, equations, and a short description. We present first the three methods that do not use binning of melt rates, then how we process the binning, and finally the three methods that use binning. For all the methods, the model data and the observational data are expressed in m yr⁻¹ (m of ice equivalent per year). The analysis methods that do not use binning are the following.

• Absolute difference of averages (ADA): we compute the average value of each ensemble member and of the target and compute the absolute difference between each ensemble member average and the average of the target.

• Two-dimensional root mean square error (2D RMSE): we compute the RMSE cell to cell with the same geographical location between each ensemble member and the target.

• Two-dimensional mean absolute error (2D MAE): we compute the MAE cell to cell with the same geographical location between each ensemble member and the target. Since no squaring is used in the error computation of the MAE, the MAE is less sensitive to outliers than the RMSE.

Table 1 Statistical methods applied to rank the ensemble members compared to the target (observations): names, equations, and short descriptions. Here, n is the number of grid cells existing in both datasets, x_member,i (x_target,i) refers to any single grid cell in one ensemble member (grid cell in the target), m is the number of bins (10 in this study), and ${\overline{B}}_{\text{member},j}$ is each bin of a single ensemble member. The same nomenclature is applied for ensemble members and the observations.

By applying the first three ranking methods, the ranking metrics do not enable us to systematically pick the ensemble members with the best fit to the distribution of values of the observational dataset. To improve that, we decide to bin each dataset: the ensemble members as well as the target. The aim of the binning is to be able to force the method to pick ensemble members that better fit the target distribution, including the higher and lower tails of the distribution. We proceed with the binning of the melt rates as shown in the schematic in Fig. 1. Each bin ($\overline{B}$) is the average value of 10 % of the total number of ordered data points. The data points must be ordered to make each bin representative of a specific share of the dataset. For instance, if an ensemble member has 200 data points, the 20 points with the lowest melt rate values are averaged and become one bin value. We proceed similarly for all the following bins and end up with 10 bins for each of the 169 members of the ensemble ${\overline{B}}_{\text{member},j}$ and 10 bins for the target ${\overline{B}}_{\text{target},j}$. Once the binning is done, we apply the following statistical analysis methods. They are similar to the three presented above but applied to the 10 bin values rather than 2D data fields.

• Average difference of averages of the bins (ADA of bins): we compute the average value of the bins of each ensemble member and of bins of the target and compute the absolute difference between the two. This method does not lead to exactly the same results as the ADA without binning because some bins might contain a different number of grid cells than others, meaning that they cover slightly different areas but still get exactly the same weighting in the binning approach. As we intend here to compare methodologies, we consider it relevant to also test this method, even if it gives results very close to ADA without binning, to quantify how much they differ.

• Root mean square error of the bins (RMSE of bins): we compute the RMSE value between the bins of each ensemble member and the bins of the target.

• Mean absolute error of the bins (MAE of bins): we compute the MAE value between the bins of each ensemble member and the bins of the target.

Figure 1 Schematic showing how the binning of melt rates is done. ${\overline{B}}_x$ stands for the value of each bin, which is the average value of 10 % of the ordered values. The darker and lighter shading represents each 10 % of the corresponding dataset. Panel (a) is the distribution of the values of the target, (b) is the distribution of the values of the first ensemble member, and (c) is the one of the last ensemble member. Once the bins for the target and all the ensemble members are calculated we apply the statistical methods to rank the ensemble members.

For all of the methods above, the ranking of the best ensemble members is given by the lowest values obtained with the equations given in Table 1. The results of the analysis applied to all ice shelves of the Antarctic ice sheet are given in Sect. 3.1. We also apply the above methods at a more local scale than Antarctica to test whether a local calibration is needed. To accomplish this we apply the above methods to two additional areas of the Antarctic ice sheet: the Filchner–Ronne Ice Shelf (FRIS) and the sector of the Bellingshausen and Amundsen seas (BA seas). The results are given in Sect. 3.2. To test the robustness and sensitivity of the methods under more conditions we also test them with different resolutions of the ice sheet model, ocean temperature and salinity forcings, and targets of basal melt rates. We discuss these results in Sect. 4.

2.6 Future applications: ISMIP 2300

To make a preliminary assessment of the relevance of this implementation and of the calibration methods, we run a small ensemble of future scenarios. For these simulations, the PICO parameter values are consistent with the results of the analysis of the calibration ensemble defined in Sect. 2.5 and presented in Fig. 2. We make use of the ISMIP6 2300 protocol in which GRISLI and other ice sheet models using PICO participated (Seroussi et al., 2024). The basal melt parameterisation used in the submitted GRISLI simulations for the ISMIP 2300 was the quadratic non-local melting parameterisation from Jourdain et al. (2020), which, in the following, will be referred to as QuadNL. For these simulations with PICO we do not re-calibrate GRISLI mechanic parameters and use the same initial state as for QuadNL. Here we repeat the experiment “expAE05” from ISMIP 2300, which corresponds to climate forcing computed by the UK Earth System Model (UKESM1-0-LL) for the scenario SSP5-8.5 (Seroussi et al., 2024). With our new simulations we will assess (i) the difference in sensitivity between PICO and QuadNL with the same model and same forcings, (ii) the importance of calibration choices for model results, and (iii) how the response of GRISLI-PICO differs from other ISMIP participating models that also used PICO. The results are presented in Sect. 3.7.

Figure 2 Comparison of the six methods. Panels (a–c) and (g–i) are the results of the distribution of the best five members according to each method. Panels (d–f) and (j–l) show the ranking of all the ensemble members according to each method. The red square shows the best member. The black dots with numbers show the best five members for each method, the same members as in the distribution panels.

3 Results

3.1 Can we capture the spatial or binned distribution of melt rates using any of the six calibration methods?

Here we present the results for the six calibration methods explained in Sect. 2 applied to all ice shelves of Antarctica at a 16 km resolution. The main calibration results for each method are presented in Fig. 2, and more detailed results are shown in Sect. S3 in the Supplement. In the Fig. 2, we show in panels (a–c) and (g–i) the distribution of the values of the five best ensemble members according to each method that we can compare directly with the distribution of the observations from Adusumilli et al. (2020a). In panels (d–f) and (j–l) we show on the heatmaps the ranking of all the ensemble members for each corresponding method. We also highlight the five best members with the black dots showing their member number.

Panels (a) and (d) show the results using the ADA method. We see that the best members using the ADA method cover a large range of values from $\mathrm{0.1}\times {\mathrm{10}}^{-\mathrm{5}}$ to $\mathrm{3.0}\times {\mathrm{10}}^{-\mathrm{5}}$ m s⁻¹ for ${\mathit{\gamma }}_T^{*}$ and from 0.1 to 5.0 Sv m³ kg⁻¹ for C. Also, the matching with the target distribution of the top five members is in some cases good (members 36 and 35) and in others not (members 133, 56, and 44). The large spread of the top five members demonstrates that with this method the best parameters will not systematically be of the same order of magnitude of values and can depend heavily on the sampling ensemble. The results for the 2D RMSE and 2D MAE methods are shown in panels (b), (e), (c), and (f), respectively. Both methods gives quite similar responses, and the top five members for both methods are side by side. In comparison to the ADA methods, these two methods enable us to have a narrower range of PICO parameter values, in particular for the parameter ${\mathit{\gamma }}_T^{*}$. However, none of the selected top five members, for both methods, match the distribution of the observations (panels b and c).

The following panels, from (g) to (l), correspond to the results of the methods using binning. The first of the three, the ADA of bins (panels g and j), selects almost the same members as the ADA without binning. Therefore, using binning before computing the ADA is not enough to have a small set of best members fitting the distribution.

The next methods, RMSE of bins and MAE of bins, minimise the differences between the bins of the ensemble members and the bins of the target systematically. Because of this we obtain a selection of the best ensemble members that systematically fit the distribution of the target (panels h and i). But also the range of parameters corresponding to the best five members is small, with all the top members side by side (panels k and l). This gives us confidence that these two last methods consistently give the same range of parameter values for different ensemble sampling that would also best match the distribution of the target.

The inability of ADA of bins to give a small range of parameter values can be explained by the fact that it allows compensation between bins. Figure 3a, c, and e show the values of the 10 bins of the five best members according to the three methods using bins. We see a compensating effect between the lower and higher parts of the values of the bins, in particular for members 133, 56, and 120. In other words, these members score well with this method because, in this case, the positive differences from the target in the lower bins are compensated for by the negative difference from the target in the higher bins, whereas the RMSE of bins and the MAE of bins have systematically smaller anomaly values and do not allow for a compensating effect.

Figure 3 Anomaly of the values of the bins for the five best members according to each method: (a) ADA of bins, (c) RMSE of bins, and (e) MAE of bins, with regard to the Adusumilli et al. (2020a)target. The closer the value is to 0, the closer the value of the bin of the ensemble member is to the value of the bin of the target. Panels (b, d) and (f) show the error in percentages between the model result and the observations for the lowest bin (b), highest bin (d), and the sum of the absolute error of the 10 bins (f). The values of the bins are the same for all three methods using binning, and therefore results shown in (b, d) and (f) are relevant for the three methodologies. The black dots with numbers in panel (f) show the five best members according to the MAE of bins, and the red square shows the best member for this same method.

The largest discrepancies between the binned values of modelled melt rates and those from observational datasets occur at the extremes, which are the bins with the lowest 10 % and highest 10 % of the values from the distribution. We therefore further analyse these two bins for all the PICO configurations in panels (b) and (d) of Fig. 3, respectively. These panels reveal distinct sensitivities to the two PICO parameters. For example, at fixed values of C, increasing ${\mathit{\gamma }}_T^{*}$ consistently raises the bin values in both the lowest (panel b) and highest (panel d) deciles. This implies a reduction in error when the model underestimates melt (negative error, in blue) or an amplification of error when it overestimates melt (positive error, in red). In contrast, at fixed ${\mathit{\gamma }}_T^{*}$, varying C can produce divergent effects between the lowest and highest bins. For instance, at ${\mathit{\gamma }}_T^{*}=\mathrm{2.0}\times {\mathrm{10}}^{-\mathrm{5}}$ m s⁻¹, increasing C leads to a decrease in bin error values in the lower 10 % bin (b) and an increase in bin error values in the upper 10 % bin (d). Finally, by computing the sum of the absolute values of the errors for all 10 bins we can find the combinations of PICO parameters that minimise this error the most. Panel (f) of Fig. 3 shows the results with the best members according to the methods using MAE of bins shown in panel (e) superposed. We find that the two metrics, MAE of bins and sum of absolute errors of the bins, lead to a similar selection of the best ensemble members.

Nonetheless, being able to match the distribution can also mean spatial compensation between different locations. Therefore we look at the spatial distribution of the values in Fig. 4. It shows the single best member of each method corresponding to the red square on the heatmaps of Fig. 2. We see that four methods (ADA, ADA of bins, 2D RMSE, and 2D MAE) result in a spatial distribution with little contrast between higher and lower values; they do not even have values more negative than −1 m yr⁻¹ in blue (Fig. 4a–c). It could be because this selection gives low ${\mathit{\gamma }}_T^{*}$ values of $\mathrm{0.1}\times {\mathrm{10}}^{-\mathrm{5}}$ and $\mathrm{0.25}\times {\mathrm{10}}^{-\mathrm{5}}$ m s⁻¹, whereas the best single member following the RMSE of bins or the MAE of bins has a lot more contrast (Fig. 4d and e), which corresponds better to what is seen in the observations (Fig. 4f). These two methods result in higher ${\mathit{\gamma }}_T^{*}$ values: $\mathrm{1.5}\times {\mathrm{10}}^{-\mathrm{5}}$ and $\mathrm{2.0}\times {\mathrm{10}}^{-\mathrm{5}}$ m s⁻¹.

Figure 4 Spatial distribution of the basal melt rate values for the six calibration methods tested (a–e). Panel (a) represents two methods (ADA and ADA of bins) as they give the same best ensemble member. Panel (f) is the spatial distribution of the observations from Adusumilli et al. (2020a). In (f), the rectangles show the two chosen areas to test the calibration at two smaller scales than Antarctic-wide as presented in Sect. 3.2; the red further to the left is the area of the Bellingshausen and Amundsen seas (BA seas) and the blue further to the right is the Filchner–Ronne Ice Shelf (FRIS) area.

Overall, we see that calculating the average, with or without binning, can lead to very different optimal PICO parameter values, which can be explained by the possibility that there is compensation between negative and positive values. The four other methods show more systematic results where the best member points are all side by side. However, between the methods without (2D RMSE and 2D MAE) and with binning (RMSE of bins and MAE of bins) the selected members are in different parameter spaces. By observing the spread of the rankings combined with the distributions, we can consider the methods RMSE of bins and MAE of bins to be the best ones among the six tested here. This is justified by the fact that the selected members (i) are better able to match the distribution curve from the target (Fig. 2a–c, and g–i), (ii) systematically give the best values in the same small range of values (Fig. 2d–f and j–l), and (iii) the magnitude of the spatial patterns is similar to the target (Fig. 4).

3.2 Do we need to calibrate the PICO parameters locally?

Here we show how different the results would be if we calibrated PICO for a specific domain of the Antarctic ice sheet, and we assess whether an Antarctic-wide calibration is suited to domain-wide applications. Figure 5 presents a selection of the analysis similar to the previous section but applied to the two domains: BA seas on the left side and FRIS on the right side. The results shown here correspond only to the method MAE of bins considered to be one of the two best methods; results for the other methods can be seen in Sects. S4 and S5 in the Supplement. First, with panels (a) and (b) we see that it is more challenging to match the observation distribution in the BA seas than for the FRIS. But also, in these two panels we show the best Antarctic-wide calibration following the MAE of bins. In the BA seas, the top member, 34, is the same as for the Antarctic-wide selection. For FRIS, the best Antarctic-wide selection is also among the top five members. Second, with panels (c) and (d), we see a strong difference between the two sectors in the sensitivity of the average basal melt to a change of the PICO parameter values. A change of the overturning coefficient of $+\mathrm{0.4}\times {\mathrm{10}}^{-\mathrm{5}}$ m s⁻¹ would lead to an average basal melt value higher than 4 m yr⁻¹ above the average of the target in the case of the BA seas, whereas it would barely make any difference for the FRIS sector. Thus, the best Antarctic-wide calibration using the MAE of bins is similar to the best possible local calibration, with negligible differences in less sensitive areas.

Figure 5 Comparison of the calibration between the Bellingshausen/Amundsen seas (a, c) and Filchner–Ronne Ice Shelf sector (b, d). Panels (a) and (b) are the distribution of the five best members for both sectors using methods with the MAE of bins, with the distribution of the best calibration applied Antarctic-wide (in purple). Panels (c) and (d) show the absolute difference between the average of each of the 169 ensemble members and the average of the observations. However, the five best members shown with the numbered black dots correspond to the best members with the method using MAE of bins (as in a and b). The best member for the local and Antarctic-wide calibration using the MAE of bins is given by the red square and the purple hexagon, respectively.

3.3 How robust are the methods if applied at different resolutions and areas of the ice sheet?

To test how robust the results are, we run an additional calibration ensemble of 169 members with a resolution of 40 km. We use the same forcing from Jourdain et al. (2020) and target from Adusumilli et al. (2020a) regridded at 40 km. With this new calibration ensemble we repeat the analysis Antarctica-wide as well as for the BA seas and FRIS. We therefore test each statistical method with six different cases: two resolutions × three areas of interest. The detailed analyses are shown in the Supplement for all six conditions in Sects. S3–S8 in the Supplement. We summarise all these conditions by aggregating the top five members in all six conditions for each statistical method and plot them in Fig. 6. It enables us to visualise the spread of the top PICO parameters aggregated over the different cases. We can see that the methods RMSE of bins and MAE of bins are the two methods that give consistent optimal PICO parameters under all different conditions. But also they give more systematically the same members (see Appendix A1), suggesting the use of the same order of magnitude of the parameters for all the tested conditions. This consistency can matter when intercomparing (i) different parameterisations, (ii) when the same parameterisation is used in different ice sheet models, and (iii) results at different resolutions. No clear trend can be seen between the best members at 16 and 40 km resolution (see Appendix A2). It is particularly relevant to have robust methods over different resolutions to make the comparisons more systematic such as in the scope of the Ice Sheet Model Intercomparison Project (ISMIP) (Seroussi et al., 2024), the Marine Ice Sheet–Ocean Model Intercomparison Project (MISOMIP) (Rydt et al., 2024), or comparison between paleo-ice-sheet and future ice sheet behaviours (Golledge et al., 2021).

Figure 6 Average ranking of each calibration method tested under two different resolutions (40 and 16 km) and applied to three different sectors (Antarctic-wide, BA seas, and FRIS). The top five members of all six conditions are shown with the black and white hexagons.

3.4 How sensitive are the methods to the forcings?

In all presented results until now we have used the forcing from Jourdain et al. (2020). Here, we assess how sensitive the results are to this forcing. We run three additional calibration ensembles of 169 members at 16 km resolution using the same forcing by Jourdain et al. (2020) but with a different temperature correction on the top of it. The temperature corrections are (1) +1 K, (2) +dT from Reese et al. (2023) (see their Table S1 in the Supplement), and (3) +dT from Jourdain et al. (2020) (see their Fig. 5a for the quadratic non-local parameterisation). Reese et al. (2023) and Jourdain et al. (2020) apply a temperature correction that differs per drainage region. To remain concise, in Fig. 7 we only present the results of the ranking for the method using MAE of bins (for additional analysis see Sect. S9–S11 in the Supplement). We can observe that with a correction of +1 K both PICO parameters shift to slightly lower values. The opposite happens with the correction from Reese et al. (2023), which is expected since the temperatures are almost all negative values, reaching up to −2 K. The correction from Jourdain et al. (2020) makes only minor differences as the values are rather low with a maximum absolute value of +1.07 K. These results suggest that using warmer forcing for the calibration will lead to a calibrated PICO less sensitive to temperature changes, and vice versa with a colder forcing. Finally, we can state that, even after combining all the ranges of values suggested by the four different forcings, forcing uncertainties lead to a smaller range of PICO parameter values than using different statistical calibration methods. In other words, the choice of the calibration method is more important than the choice of the forcing. Since the best parameters do not vary much (about −0.5 for the ${\mathit{\gamma }}_T^{*}$ and −0.05 for C) between the +0 K and the +1 K forcings, we analyse in the next subsection the sensitivity of PICO to better understand what the response of PICO would be to warmer than present-day conditions.

Figure 7 Ranking using the MAE of bins for all Antarctic ice shelves with different forcings. Panel (a) is without additional temperature correction, (b) is with +1 K, (c) is with dT defined by Reese et al. (2023), and (d) is with dT defined by Jourdain et al. (2020).

3.5 What is the melt rate sensitivity of PICO to changes in ocean temperature?

Thanks to the ensembles with +0 and +1 K ocean forcings from the previous subsection we can determine the melt rate sensitivity of PICO to changes of ocean temperature for all the ensemble members. We assume here a linear sensitivity and therefore compute it as the difference of the +1 K experiment minus the +0 K experiment. The results are shown in Fig. 8 where we also differentiate the sensitivity in the three areas defined for the analysis shown in Sect. 3.2. First of all, we see that in most cases the sensitivity of PICO increases when the value of either of the two parameters is increased. Second, we see that the sensitivity varies between areas. For instance, 1 K of warming with the same PICO parameters as used by PISM-PICO (Seroussi et al., 2024) (squares in Fig. 8) would lead to an increase of 1.5 $\mathrm{m}{\mathrm{yr}}^{-\mathrm{1}}{\mathrm{K}}^{-\mathrm{1}}$ in FRIS, whereas the increase would be 8.4 $\mathrm{m}{\mathrm{yr}}^{-\mathrm{1}}{\mathrm{K}}^{-\mathrm{1}}$ in the BA sea ice shelves. We also see that the range of possible sensitivities (panel d) is about 4 times larger in the BA seas than in the FRIS. These results quantify how much the sensitivity of the basal melt rate would change, globally and regionally, by changing the values of the two PICO parameters. In all the cases, the best calibration with the MAE of bins (black hexagons) is in the lower range of all the tested combinations of parameter values. These values are also lower than the range of Antarctic ice shelf sensitivity estimates from some previous studies (Levermann et al., 2020; van der Linden et al., 2023) but closer to the PICO sensitivity obtained by Reese et al. (2023) and Lambert and Burgard (2025) when optimising parameters for present-day melting. The methodologies in the assessments of the sensitivities are, however, different in each study. Nonetheless, based on the results shown in Fig. 8 we can expect a low to moderate response of ice shelves in this calibrated version of PICO to future projections scenarios.

Figure 8 Linear sensitivity of PICO for the three areas of interest: (a) all ice shelves of Antarctica, (b) Filchner–Ronne Ice Shelf, and (c) ice shelves in the Bellingshausen and Amundsen seas. Panel (d) compares the ranges of sensitivity of the three areas. In all the panels, the hexagon represents the calibration done in this study with the MAE of bins, and the square represents the calibration used in PISM-PICO for ISMIP 2300 (Seroussi et al., 2024).

3.6 How sensitive are the methods to the target?

As presented in Sect. 2.4, we took as the target the basal melt rates retrieved from Adusumilli et al. (2020a). However, the retrieval of basal melt rates from satellite observations is poorly constrained. Hence, we could make an argument for choosing a different target for the calibration of the basal melt rate parameterisation. Therefore, to assess the uncertainty due to the choice of the melt rate target, we ran the same robustness analysis as in Sect. 3.3 but with basal melt rates retrieved by Paolo et al. (2023) as the target the instead of Adusumilli et al. (2020a). The detailed results are shown in Appendix B. Overall, they are similar to the ones obtained with Adusumilli et al. (2020a) as the target, with a slight shift towards higher ${\mathit{\gamma }}_T^{*}$ values. We also make the bin analysis in Fig. 3 with the datasets from Paolo et al. (2023) and we observe similar results (see Fig. S29 in the Supplement). This is expected as the distribution and the main statistics are very similar (see Fig. 9a, with a higher standard deviation for the Paolo et al., 2023, dataset). Figure 9b shows the spatial differences between the two datasets. The magnitude of these differences is of the same order of magnitude as the values of basal melt rates themselves, reaching values below −4 m yr⁻¹ and above 4 m yr⁻¹. This observation is made for both resolutions 16 and 40 km (see Sect. S13 in the Supplement). The standard deviation of the difference between the two datasets is 1.97 m yr⁻¹, which is about 3 times larger the mean melt rate of both datasets. Hence, despite the important spatial differences, we conclude that the calibration methods are not significantly sensitive to the choice of the target dataset as they both give a similar selection of the best ensemble members for each method.

Figure 9 (a) Histogram of values of the two datasets and main statistics. (b) Spatial distribution of the difference of basal melt rates between Adusumilli et al. (2020a) and Paolo et al. (2023). Results in this figure are shown for the mesh grid resolution of 16 km; more details are given for mesh grid resolutions of 16 and 40 km in Sect. S13. The difference of the means ($\mathrm{0.66}-\mathrm{0.68}=-\mathrm{0.02}$ m yr⁻¹) is different from the mean of the difference (0.16 m yr⁻¹) because in the mean of the difference only the grid cells with values in both datasets are taken into account.

3.7 How much does the calibration method matter for future projections?

Following the previous analysis, we select seven cases to make a first-order assessment of the impact of the calibration method choice on future projections of the Antarctic ice sheet until 2300 (Fig. 10a). Five of the seven cases correspond to the best members according to the six methods applied Antarctic-wide. We get five cases rather than six because ADA and ADA of bins give the same best member. The two additional members correspond to the parameter values chosen by the ice sheet models PISM/Elmer-ice (101) and Kori (102) that use PICO for the ISMIP 2300 experiment in Seroussi et al. (2024). Figure 10b shows the distribution of the basal melt rate at the very start of the simulations for the ISMIP 2300 simulation for the year 2015 compared to the observations from Adusumilli et al. (2020a). For comparison, we include GRISLI with the QuadNL parameterisation (Jourdain et al., 2020).

Figure 10 (a) Selection of PICO parameter values for ISMIP2300 applications in this study based on the different calibration methods and previous PICO calibrations. (b) Comparison of the basal melt rate between Adusumilli et al. (2020a) and the start of the simulation (t=2015) for each calibration.

The main results of this small ensemble of ISMIP 2300 members are shown in Fig. 11. In Fig. 11a we see the total basal mass balance flux (BMB flux) over time. We can see very different behaviours between ice sheet models and parameterisations. Over the whole simulation the highest values are obtained with the QuadNL parameterisation with GRISLI and Kori. Then in the medium range we have the calibrations 101 and 102 of PICO with different ice sheet models (Elmer-ice, PISM, Kori, and GRISLI). In the upper part of the figure, the group of simulations with the lowest basal melt rates corresponds to all the GRISLI-PICO simulations with the different calibration methods applied above. In Fig. 11b, the floating ice area representing the size of the ice shelves shows significant differences between the different simulations, including the GRISLI-PICO with different calibrations. The growth of the floating ice of the GRISLI-PICO calibrations 101 and 102 is something also observed in the Kori-PICO simulation. The calibration 34, considered the best one in the calibration process, shows a small decreasing trend in floating ice area over the course of the simulation compared with most other cases. Finally, Fig. 11c shows the contribution to sea level rise of all the simulations. We can see that simulations with PICO with the same PICO parameters lead to lower values of sea level contribution by 2300 for both GRISLI and Kori compared to results from the same ice sheet models but using the QuadNL parameterisation. Elmer-ice with PICO even suggests a negative sea level contribution from Antarctica throughout the simulation. Except compared to Elmer-ice with PICO, all the simulations of GRISLI-PICO result in lower values of sea level contribution by 2300 than all the other cases.

Figure 11 ISMIP 2300 applications with different PICO calibrations. Panel (a) shows the evolution of the total basal melt balance (BMB) beneath floating ice over time. Panel (b) shows the evolution of the floating ice area. Panel (c) shows the contribution to sea level rise.

4 Discussion

Here, we develop the paper further by discussing three questions: what can be the right observational target? What can we do to better understand the sensitivity of the Antarctic ice shelves to warming oceans? And why does PICO lead to lower sea level contribution estimates than QuadNL in GRISLI?

4.1 What could be the right observational target?

The observation of the disagreement between the two target datasets shown in Fig. 9 is important to justifying the usage of bins suggested in the present study. Indeed, by using the binning methods we give spatial freedom to the datasets and constrain them by their values. Since calibration methods that minimise cell-to-cell differences between modelled and observed melt rates often fail to capture the overall distribution of observed values, and given the spatial inconsistencies among observational datasets, we prioritise reproducing the correct distribution of basal melt rates over minimising spatial mismatches. We consider having a good statistical representation of the melt rates to be potentially more important for the dynamic of the ice sheet. For instance, the highest melt rate values are observed in the Amundsen Sea area, where due to the retrograde slope the West Antarctic ice sheet is exposed to the marine ice sheet instability process (Weertman, 1974; Joughin et al., 2014). Therefore, capturing these high melt rates values is potentially important for future projections. Moreover, even if we might not have the right values at the right locations within a given ice shelf, we have seen in Sect. 3.2 that the calibration method MAE of bins enables us to have the values close to the distribution of the target in local areas. Prioritising values over spatial correspondence within an ice shelf is in agreement with Joughin et al. (2021), who argue that the ocean-induced melt volume, regardless of the spatial distribution, directly paces the ice loss. However, other studies suggest that localised sub-ice-shelf melt can have a strong impact on the buttressing or that in more strongly buttressed areas sub-ice-shelf melt would have an outsized effect (Gudmundsson, 2013; Reese et al., 2018b). Additionally, one limitation here is that by scaling up the resolution of the observational datasets to the ice sheet resolution, 16 or 40 km here, we are losing most of the data points with melt rate values above 6 m yr⁻¹ (see Sect. S13). Consequently, certain ice shelves in the West Antarctic ice sheet such as Thwaites or Pine Island have very few grid cells at 40 km resolution. This can be an important limitation, yet we consider it important to calibrate and test for this coarse resolution because it is an option for paleo-ice-sheet simulations. This also means that a calibration by using the same method but for higher-resolution models or with irregular grids might have different values of the PICO parameters than the ones found here. At a higher resolution we could also consider computing more than 10 bins, but this has not been explored in this study.

Nonetheless, with this analysis using the observational dataset from Paolo et al. (2023) we can justify the fact that the method using MAE of bins is more robust than the RMSE of bins in a case where we want to use one set of parameters for all the Antarctic ice shelves. Indeed, as shown in Fig. 12, the range of values of the parameters suggested by the top five members is smaller for the MAE of bins (C ranging from 0.01 to 0.25 Sv m³ kg⁻¹ and ${\mathit{\gamma }}_T^{*}$ from $\mathrm{1.00}\times {\mathrm{10}}^{-\mathrm{5}}$ to $\mathrm{5.00}\times {\mathrm{10}}^{-\mathrm{5}}$ m s⁻¹) than the RMSE of bins (C ranging from 0.01 to 0.50 Sv m³ kg⁻¹ and ${\mathit{\gamma }}_T^{*}$ from $\mathrm{1.00}\times {\mathrm{10}}^{-\mathrm{5}}$ to $\mathrm{15.00}\times {\mathrm{10}}^{-\mathrm{5}}$ m s⁻¹).

Figure 12 Average ranking of two calibration methods, RMSE of bins (a) and MAE of bins (b), tested under two different resolutions (40 and 16 km) and applied to three different sectors (Antarctic wide, BA seas, and FRIS) with the observational dataset from Paolo et al. (2023) as the target. The top five members of all six conditions are shown with the black and white hexagons.

4.2 Towards a better understanding of the sensitivity of Antarctic ice shelves to warming oceans

Understanding the sensitivity of the Antarctic ice shelves to ocean warming is key to be able to make future projections. In Sect. 3.5 we quantified the Antarctic ice shelf sensitivity of our PICO version based on a highly simplified approach. More advanced methods have been developed; for instance, Lambert and Burgard (2025) apply both salinity and temperature anomalies that compensate for each other to maintain a present-day-like density profile. Previous studies report a wide range of sensitivity estimates, spanning up to an order of magnitude (Levermann et al., 2020; van der Linden et al., 2023; Lambert and Burgard, 2025); however, they differ in their methodology to compute the sensitivity. These differences include the choice of forcing, the sub-shelf melt parameterisation, model resolution, and whether a linear or quadratic sensitivity is assumed. While this goes beyond the present study, developing a standardised framework for quantifying the sensitivity of the ice shelves to ocean warming would be useful to facilitate comparisons across studies and models. In parallel, continuing efforts in Earth observations to extend and refine records of ocean properties and sub-shelf melt rate remain critical to improving the robustness of sensitivity estimates. Additionally, we observe a large range of possible melt rate sensitivities to ocean properties with changes in the values of the parameters (Fig. 8). Therefore, we could argue for calibrating basal shelf melt rate parameterisations directly to a target value of sensitivity. This has been done by Reese et al. (2023), leading to higher PICO sensitivity than the previous calibration (Reese et al., 2018a) and used in the future projections until 2300 (Seroussi et al., 2024).

4.3 Why does PICO lead to lower sea level contribution estimates than QuadNL in GRISLI?

Overall, including PICO in ice sheet models leads to lower sea level contribution in the simulations up to 2300 shown in Fig. 11. Indeed, Kori with the PICO parameterisations produces 1.5 m sea level equivalent less than with the QuadNL parameterisation. Similarly, all the PICO calibrations with GRISLI have a lower sea level contribution than with the QuadNL parameterisation. Elmer-ice with PICO even suggests a negative sea level contribution from Antarctica throughout the simulation. Of course, there are many other factors influencing future sea level predictions in ice sheet models that are not related to the parameterisation of the ice–ocean interactions, and it is definitely possible to have a larger sea level contribution with PICO as shown in the simulation with PISM. But here we want to provide some hypotheses that could explain this overall pattern of lower sea level contribution from simulations with PICO. We outline five possibilities.

• The overturning circulation under the ice shelves, which tends to reduce the basal melt rate, is computed differently in the parameterisations PICO and QuadNL. In PICO the overturning fluxes are computed with the overturning circulation coefficient C and the difference in densities (see Eq. 2), whereas in the QuadNL it is in the product involving the thermal forcings, which results in stronger overturning from warmer conditions (see Eq. 1 in Jourdain et al., 2020).

• PICO includes a linear relationship between the ocean temperature and the basal melt rate for high temperatures, whereas QuadNL has a quadratic one. It means that for high projections, if both start off with a similar basal melt rate, the QuadNL will project significantly higher melt rate values with increasing temperatures.

• The QuadNL parameterisation takes as input 3D fields of oceanic forcings, whereas PICO takes one value per ice shelf (T₀ and S₀), which is an average of values over the continental shelf at the depth of the continental shelf in front of the concerned ice shelf.

• PICO tends to have more smoothed-out melt rates and does not show significantly higher melt rates at grounding lines, as seen in satellite-derived fields. If the basal melt rates at the grounding line are a major factor for ice loss, it could explain a less sensitive response.

• The sensitivity to ocean warming in our calibrated version of PICO lies below some previously reported ranges of Antarctic ice shelf sensitivity (Levermann et al., 2020; van der Linden et al., 2023) and is more consistent with the PICO sensitivity range estimated by Lambert and Burgard (2025). In their study, the PICO sensitivity is lower than that obtained with other sub-shelf melt rate parameterisations, although it is not the lowest overall. By contrast, the QuadNL parameterisation lies on the higher end of the sensitivity spectrum compared to other parameterisations (Burgard et al., 2023; Lambert and Burgard, 2025). So far, neither PICO nor QuadNL sensitivity ranges, which also depend on their calibration, can be ruled out as we do not know what the right sensitivity is.

In addition, the results of the future simulations show that the ice shelves can have different behaviours depending on the calibration method and the choices of the values of the parameters. Thus, we advocate here for a calibration methodology that best fits the full distribution of the observational datasets; it is a more physical calibration of the process modelled than simply matching the average value for instance. This methodology can potentially be applied to modules in other models that benefit from existing observational datasets. However, regardless of the quality of the calibration, all parameterisations are simplifications of processes. Therefore, one should be aware of the processes represented or not in a model for the interpretation of the outputs and when possible still run an ensemble to cover the range of uncertainties.

5 Conclusions and perspectives

The Antarctic ice sheet retreat is driven by ice–ocean interaction and differences in the ice–ocean parameterisations can lead to major differences in simulations of future dynamics of the ice sheet. We presented the implementation of the PICO basal ice shelf melt module (Reese et al., 2018a) in the GRISLIv.2 ice sheet model (Quiquet et al., 2018). Then we compared six statistical calibration methods to find the best set of two PICO parameters, C and ${\mathit{\gamma }}_T^{*}$. We demonstrated that the only two methods, the RMSE of bins and the MAE of bins, which also fit the low and high extremes in the target histogram, provide a robust constraint of both parameters in a narrower range of values. They give more consistent results, making them more reliable and less dependent upon the ensemble sampling. The results from these two methods also better reproduce the range of special variability, if not the details of spatial patterns observed in the chosen target (Adusumilli et al., 2020a). By using these two methods that closely fit the entire distribution of the target for all Antarctic ice shelves combined, we also show that region-specific calibration is not necessary. According to this research, the best values of the PICO parameters in our specific setup are ${\mathit{\gamma }}_T^{*}=\mathrm{2.0}\times {\mathrm{10}}^{-\mathrm{5}}$ m s⁻¹ and C=0.1 Sv m³ kg⁻¹.

We did future simulations without re-calibration of the mechanical parameters of GRISLI, such as ice flow and drag, to have a preliminary assessment of the impact of the choice of the calibration method applied to PICO. On the one hand, we see that the re-calibration of the PICO parameters can lead to major differences in sea level contribution compared to simulations using parameter values from previous studies. We showed that only our calibration fits the whole distribution of sub-shelf melt rates from the observations. On the other hand, the choice of the calibration method does not have a major direct impact on the sea level contribution, yet it does have a significant impact the ice shelf extent. This remains coherent as the relationship between ice shelf changes and sea level contribution is still uncertain and ice-sheet-model-dependent (Sun et al., 2020). Thus, our analysis of future simulations is not enough to gain confidence in one specific projection, and further work is needed to constrain the sensitivity of the Antarctic ice sheet using paleo-records and coupled ice sheet–ocean models.

Finally, we thoroughly tested the statistical methods by assessing how robust the results are by applying them to additional cases such as different resolutions, regions of Antarctica, forcings, and targets. This assessment give us confidence in our results, confirming that the methods using RMSE of bins or the MAE of bins are the most robust ones and could avoid the need for modellers to use temperature corrections on top of the parameterisation as well as give more confidence in paleo-ice-sheet applications at lower resolutions using present-day data for the calibration (Quiquet and Roche, 2024). The principle of using bins is justified by observing the magnitude of the spatial disagreement between the observational datasets (Adusumilli et al., 2020a; Paolo et al., 2023). As the MAE of bins gives a smaller parameter space under the six conditions and different targets tested, we recommend using this method.

To progress further, we invite ice–ocean interaction modellers to test the method using MAE of bins in their own setup of ice–ocean parameterisation, ice sheet model, and initial state. Another possible improvement would be to target the sensitivity of the basal melt rate to ocean forcing changes rather than targeting a given basal melt rate for a given ocean forcing. Reese et al. (2023) calibrated PICO to a sensitivity to temperature, but it required the use of temperature corrections. We suggest that using the calibration method of MAE of bins could enable calibrating to sensitivity without additional temperature corrections. Finally, when possible, we encourage modellers to quantify the sensitivity of their sub-shelf melt parameterisation.

Appendix A: Standard deviation of the rankings for each method applied to six different conditions (two resolutions × three areas of the Antarctic ice sheet)

Figure A1 Same as Fig. 6 but showing the standard deviation of the rankings under the six different conditions instead of the average. It gives additional information about the degree of confidence and robustness of the methods.

Figure A2 Same as Fig. 6 but differentiating the two resolutions.

Appendix B: Average and standard deviation of the rankings for each method applied to six different conditions (two resolutions × three areas of the Antarctic ice sheet) with Paolo et al. (2023) as the target instead of Adusumilli et al. (2020a)

Figure B1 Same as Fig. 6 but with Paolo et al. (2023) as the target.

Figure B2 Same as Fig. 6 but the standard deviation with Paolo et al. (2023) as the target.