# PARASO, a circum-Antarctic fully coupled ice-sheet-ocean-sea-ice-atmosphere-land model involving f.ETISh1.7, NEMO3.6, LIM3.6, COSMO5.0 and CLM4.5

Charles Pelletier, Thierry Fichefet, Hugues Goosse, Konstanze Haubner, Samuel Helsen, Pierre-Vincent Huot, Christoph Kittel, François Klein, Sébastien Le Clec'H, Nicole P. M. van Lipzig, et al.

► To cite this version:

Charles Pelletier, Thierry Fichefet, Hugues Goosse, Konstanze Haubner, Samuel Helsen, et al. PARASO, a circum-Antarctic fully coupled ice-sheet-ocean-sea-ice-atmosphere-land model involving f.ETISh1.7, NEMO3.6, LIM3.6, COSMO5.0 and CLM4.5. Geoscientific Model Development, 2022, 15, pp.553-594. 10.5194/gmd-15-553-2022. insu-03706399

HAL Id: insu-03706399

https://insu.hal.science/insu-03706399

Submitted on 27 Jun 2022

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>

# PARASO, a circum-Antarctic fully coupled ice-sheet-ocean-sea-ice-atmosphere-land model involving f.ETISh1.7, NEMO3.6, LIM3.6, COSMO5.0 and CLM4.5

Charles Pelletier<fnref n="1" />, Thierry Fichefet<fnref n="1" />, Hugues Goosse<fnref n="1" />, Konstanze Haubner<fnref n="2" />, Samuel Helsen<fnref n="3" />, Pierre-Vincent Huot<fnref n="1" />, Christoph Kittel<fnref n="4" />, François Klein<fnref n="1" />, Sébastien Le clec'h<fnref n="5" />, Nicole P. M. van Lipzig<fnref n="3" />, Sylvain Marchi<fnref n="3" />, François Massonnet<fnref n="1" />, Pierre Mathiot<fnref n="6" />, <fnref n="7" />, Ehsan Moravveji<fnref n="3" />, <fnref n="8" />, Eduardo Moreno-Chamarro<fnref n="9" />, Pablo Ortega<fnref n="9" />, Frank Pattyn<fnref n="2" />, Niels Souverijns<fnref n="3" />, <fnref n="10" />, Guillian Van Achter<fnref n="1" />, Sam Vanden Broucke<fnref n="3" />, Alexander Vanhulle<fnref n="5" />, Deborah Verfaillie<fnref n="1" />, and Lars Zipf<fnref n="2" />

**Correspondence:** Charles Pelletier (chrles.pelletier@protonmail.com)

Received: 13 September 2021 – Discussion started: 4 October 2021

Revised: 16 December 2021 – Accepted: 20 December 2021 – Published: 25 January 2022

**Abstract.** We introduce PARASO, a novel five-component fully coupled regional climate model over an Antarctic circumpolar domain covering the full Southern Ocean. The state-of-the-art models used are the fast Elementary Thermomechanical Ice Sheet model (f.ETISh) v1.7 (ice sheet), the Nucleus for European Modelling of the Ocean (NEMO) v3.6 (ocean), the Louvain-la-Neuve sea-ice model (LIM) v3.6 (sea ice), the COnsortium for Small-scale MOdeling (COSMO) model v5.0 (atmosphere) and its CLimate Mode (CLM) v4.5 (land), which are here run at a horizontal resolution close to 1/4°. One key feature of this tool resides in a novel two-way coupling interface for representing ocean-ice-sheet interactions, through explicitly resolved ice-shelf cavities. The impact of atmospheric processes on the Antarctic ice sheet is also conveyed through computed COSMO-CLM-f.ETISh surface mass exchange. In this technical paper, we briefly introduce each model's configuration and document the developments that were carried out in order to establish PARASO. The new offline-based NEMO-f.ETISh

coupling interface is thoroughly described. Our developments also include a new surface tiling approach to combine open-ocean and sea-ice-covered cells within COSMO, which was required to make this model relevant in the context of coupled simulations in polar regions. We present results from a 2000–2001 coupled 2-year experiment. PARASO is numerically stable and fully operational. The 2-year simulation conducted without fine tuning of the model reproduced the main expected features, although remaining systematic biases provide perspectives for further adjustment and development.

## 1 Introduction

The Antarctic climate is characterized by large natural fluctuations and complex interactions between the ice sheet, ocean, sea ice and atmosphere. One of the consequences of these in-

<footnote n="1">¹Earth and Life Institute (ELI), UCLouvain, Louvain-la-Neuve, Belgium</footnote>

<footnote n="2">²Laboratoire de Glaciologie, Université Libre de Bruxelles, Brussels, Belgium</footnote>

<footnote n="3">³Department of Earth and Environmental Sciences, KU Leuven, Leuven, Belgium</footnote>

<footnote n="4">⁴Laboratory of Climatology, Department of Geography, SPHERES, University of Liège, Liège, Belgium</footnote>

<footnote n="5">⁵Earth System Science and Department Geografie, Vrije Universiteit Brussel, Brussels, Belgium</footnote>

<footnote n="6">⁶Met Office, Exeter, United Kingdom</footnote>

<footnote n="7">⁷Université Grenoble Alpes/CNRS/IRD/G-INP, IGE, Grenoble, France</footnote>

<footnote n="8">⁸ICTS, KU Leuven, Leuven, Belgium</footnote>

<footnote n="9">⁹Barcelona Supercomputing Center (BSC), Barcelona, Spain</footnote>

<footnote n="10">¹⁰Environmental Modelling Unit, Flemish Institute for Technological Research (VITO), Mol, Belgium</footnote>
teractions observed over the last decades lies in the Antarc-
tic ice sheet (AIS) losing mass (Rignot et al., 2019; Ed-
wards et al., 2021), with the most spectacular changes located
in West Antarctica (Bell and Seroussi, 2020). Atmospheric
processes have been underlined as the main driver behind
the large Antarctic climate variability (Jones et al., 2016;
Raphael et al., 2016; Schlosser et al., 2018; Lenaerts et al.,
2019). Since precipitation is the only source of mass gain for
the AIS, moisture advection and atmospheric rivers have a
direct impact on the evolution of the AIS (van Lipzig and
van den Broeke, 2002; Gorodetskaya et al., 2014; Souverijns
et al., 2018; Wille et al., 2021). However, the main source of
AIS mass loss is the melting of the ice shelves (i.e., the float-
ing part of the ice sheet), which are in direct contact with
the Southern Ocean (Paolo et al., 2015; Jenkins et al., 2018;
Shepherd et al., 2018). At the decadal timescale, atmospheric
circulation patterns have an important indirect impact on ice-
shelf melting through their direct influence on the oceanic
circulation (Dutrieux et al., 2014; Jenkins et al., 2016; Hol-
land et al., 2019). For example, oceanic processes such as the
intrusion of relatively warm water over the Antarctic con-
tinental shelf, reaching the bottom of the ice shelves, have
been suggested as an explanation for the recent increase in
ice-shelf melt rates (Jacobs et al., 2011; Darelius et al., 2016;
Rintoul et al., 2016). In turn, ice-shelf melting leads to fresh-
water injection at depth, thus impacting the circulation and
heat transport in the Southern Ocean (Hellmer, 2004; Jour-
dain et al., 2017; Jeong et al., 2020). The impact on Antarctic
sea-ice cover of this additional freshwater injection harbors
complex phenomena which are still open to discussion. As
a consequence, at the circumpolar scale, no definitive con-
clusion on the impact of this coupling mechanism has been
drawn yet. While some preliminary studies suggested ice-
shelf melting may have played a crucial role in the 1979–
2016 Antarctic sea-ice extent (SIE) increase (Hellmer, 2004;
Bintanja et al., 2013), other investigations have suggested
that its impact is negligible (Pauling et al., 2016).

The examples listed above underline the importance of the
interactions between the ice sheet, ocean (including sea ice)
and atmosphere at the high latitudes of the Southern Hemi-
sphere. Despite their varying typical response timescales, the
respective evolutions of polar Earth subcomponents are en-
tangled and interdependent (Turner et al., 2017; Fyke et al.,
2018). Coupled models are the designated tool for quanti-
tatively investigating them, since they can incorporate feed-
backs between the components. As the success of coupled
model intercomparison projects (e.g., CMIP6, Eyring et al.,
2016) testifies, coupled ocean-atmosphere models have been
a staple in climate modeling for a few decades. All state-of-
the-art ocean models covering polar regions include a sea-
ice submodel (see Bertino and Holland, 2017, for a review).
More recent advances in ocean modeling (Dinniman et al.,
2016; Asay-Davis et al., 2017; Mathiot et al., 2017; Zhou
and Hattermann, 2020; Gladstone et al., 2021) and the de-
velopment of idealized test cases for simulating ice-shelf

melting (Asay-Davis et al., 2016; Zhang et al., 2017; Jor-
dan et al., 2018; Seroussi and Morlighem, 2018; Favier et al.,
2019) have paved the way towards realistic ocean simulations
with explicitly resolved circulation within ice-shelf cavities
on regional configurations (Donat-Magnin et al., 2017; Jour-
dain et al., 2017; Kusahara et al., 2017; Hazel and Stewart,
2020; Hausmann et al., 2020; Huot et al., 2021; Nakayama
et al., 2021; Van Achter et al., 2022), but also on circum-
polar domains (Naughten et al., 2018; Nakayama et al.,
2020; Comeau et al., 2022). Coupled ocean-ice-sheet con-
figurations, with dynamical AIS geometry, have also been
developed, first using global climate models of intermedi-
ate complexity (Swingedouw et al., 2008; Goelzer et al.,
2016; Ganopolski and Brovkin, 2017) but also state-of-the-
art ocean models over local domains (Seroussi et al., 2017;
Timmermann and Goeller, 2017; Pelle et al., 2021). More
recently, Kreuzer et al. (2020) performed global-scale, cou-
pled ocean-ice-sheet pluri-centennial simulations relying on
parameterized ice-shelf cavity circulation. Model configura-
tions including ice-sheet coupling with cavities explicitly re-
solved are emerging in a few global Earth system models
(Gierz et al., 2020; Smith et al., 2021), as well as in more
idealized settings, e.g., for constraining sea-level-rise pro-
jections (Muntjewerf et al., 2021) or performing sensitivity
studies (Richter et al., 2021).

In this paper, we introduce PARASO, a new five-
component fully coupled circumpolar regional model cover-
ing the whole Antarctic continent and Southern Ocean, de-
signed for decadal timescale applications. Compared with
current global climate models, the regional PARASO con-
figuration is run at a relatively high resolutions (≈ 20 km), at
which key Antarctic processes, such as the land–atmosphere
surface mass balance (SMB; see Mottram et al., 2021) or
ocean warm water inflow into the continental shelf (Thompson
et al., 2018), start being appropriately resolved. Moreover,
the regional nature of PARASO maintains its computational
cost reasonable: about 2 simulated years per day with
approximately 250 cores (see Appendix B).

Our tool relies on distinct models for each treated Earth
system subcomponent: f.ETISh v1.7 for the ice sheet (Pattyn,
2017; Sun et al., 2020); COSMO v5.0 for the atmospheric
circulation, coupled to CLM v4.5 as a land module (Rockel
et al., 2008); NEMO v3.6 for the ocean (Madec et al., 2017),
coupled to LIM v3.6 for the sea ice (Rousset et al., 2015).

A specificity of PARASO lies in the computation of the
subshelf melting, which is directly assessed from the ex-
plicitly resolved oceanic circulation within ice-shelf cavities.
In cavity-deprived ocean simulations, the impact of the ice-
shelf melting on the intra-cavity circulation cannot be repre-
sented and the ice-sheet geometry is static, which are simpli-
fications. Another less simplified option relies on parameter-
ized cavities (e.g., Reese et al., 2018). This is a reasonable
compromise for global simulations at relatively low resolu-
tions covering long periods (e.g., Kreuzer et al., 2020). In
the case of PARASO, explicitly resolving cavities has been
preferred, since it allows directly representing the ocean-ice-sheet boundary layer from prognostic ocean variables, and thus assessing ice-shelf melting in a more organic, low-level manner. The ice-sheet model then dynamically responds to these melt rates estimated from the cavities, as well as to the SMB from the atmosphere and land models, by providing the ocean model an evolving cavity geometry.

In this paper, we focus on the technical aspects related to the new coupling interfaces specifically developed for PARASO, evaluate their numerical soundness, and investigate their impact by comparing with stand-alone simulations from PARASO components. The paper is organized as follows. In Sect. 2, each model component used in PARASO is briefly introduced. The novel coupling interfaces specifically designed for our setup are more thoroughly discussed in Sect. 3. The configuration and the external forcings are described in Sect. 4. Results from a 2-year long fully coupled simulation are commented on in Sect. 5. A broader discussion, along with concluding remarks, is provided in Sect. 6.

# 2 Model descriptions

## 2.1 Ice-sheet model: f.ETISh v1.7

The f.ETISh (fast Elementary Thermomechanical Ice Sheet) model v1.7 is a vertically integrated hybrid finite-difference ice-sheet/ice-shelf model with vertically integrated thermomechanical coupling (see Pattyn, 2017, for a complete overview of v1.0). Differences between v1.0 and v1.7 pertain to an improved numerical stability of the different solvers, as well as improvements on the initialization procedure described in Sect. 4.3.1. These encompass the improved method due to Bernales et al. (2017) and the implementation of a regularization term to filter out high-frequency noise in the optimized basal sliding coefficients. In addition to this, v1.7 was also specifically adapted for the coupling with NEMO by allowing time-varying subshelf melt rates as boundary conditions and imposing the calving front to remain static (since the PARASO land-sea mask cannot evolve; see Sect. 3.2).

While f.ETISh is intrinsically a two-dimensional (plane) model, the transient englacial temperature field is calculated in a three-dimensional fashion using shape functions for vertical velocities. Ice dynamics are solved through the combination of the shallow-shelf approximation (SSA) for basal sliding and the shallow-ice approximation (SIA) for internal ice deformation (Bueler and Brown, 2009). The marine boundary is represented by a grounding-line flux condition according to Schoof (2007), coherent with power-law basal sliding. We employ a Weertman sliding law with exponent $m = 2$. The constant fluid densities (ice and ocean) as well as the minimal ice thickness have been set to values consistent with NEMO's (see Table 1).

## 2.2 Ocean and sea-ice model: NEMO-LIM3.6

For simulating the ocean-sea-ice system, we use the Nucleus for European Modelling of the Ocean v3.6 (NEMO3.6, Madec et al., 2017), which includes OPA (Océan PARallélisé) for the open ocean coupled with the Louvain-la-Neuve sea-ice model LIM3.6 (Vancoppenolle et al., 2012; Rousset et al., 2015). Hereinafter, this combination of ocean and sea-ice model is referred to as NEMO.

### 2.2.1 Generalities on the NEMO3.6 open-ocean model

NEMO3.6 is a primitive-equation, free-surface, finite-difference, hydrostatic ocean model relying on a C grid (Arakawa, 1966). In our setting, the ocean dynamics are based on a split-explicit formulation. A $z^*$ vertical coordinate is used with varying cell thickness over the whole column (Adcroft and Campin, 2004; Bruno et al., 2007). Vertical mixing is parameterized using a turbulent kinetic energy (TKE) scheme derived from Gaspar et al. (1990). Lateral diffusion of momentum is carried out with a bi-Laplacian viscosity (Griffies and Hallberg, 2000), convection is parameterized as enhanced vertical diffusivity, and double-diffusive mixing is included (Merryfield et al., 1999). A free-slip lateral boundary condition is applied on the coastlines, and the Beckmann and Döscher (1997) bottom boundary layer scheme is applied to both tracers and prognostic variables. The equation of state is based on a 75-term polynomial approximation of the Thermodynamic Equation of SeaWater 2010 (TEOS-10) (IOC, 2010; Roquet et al., 2015). A list of NEMO key parameters used for PARASO is provided in Table A1.

### 2.2.2 Representation of the ice-shelf cavities within the ocean model

In PARASO, Antarctic ice-shelf cavities are opened to the ocean circulation (Mathiot et al., 2017) with configuration-dependent geometrical constraints defined in Sect. 4.1 and listed in Table 1. As in Barnier et al. (2006), partial (i.e., vertically cut) cells are used to better represent the ice-shelf draft (depth of the immersed part of the ice sheet), the top cavity cell being defined as the ice-ocean boundary layer (Losch, 2008). Heat and mass fluxes associated with melt are computed from the conservative form of the three-equation formula from Hellmer and Olbers (1989); Holland and Jenkins (1999); Jenkins et al. (2010), relying on velocity-dependent heat and salt exchange coefficients, here parameterized as in Dansereau et al. (2014):

$$ \gamma_x = \sqrt{C_D^{\text{isf}} \|\mathbf{u}_h^{\text{isf}}\| \Gamma_x}, \quad (1) $$

where $x \in \{\text{T, S}\}$ distinguishes heat and salt, $\Gamma_x$ is the default exchange coefficient, $C_D^{\text{isf}}$ the ice-ocean drag coefficient, and $\|\mathbf{u}_h^{\text{isf}}\|$ the near-ice ocean current velocity. The ice-shelf melt parameterization also includes TEOS-10-compatible coefficients
ficients for evaluating seawater freezing point temperature (Jourdain et al., 2017). In the event of ice-shelf melting, the associated heat flux is extracted from the ocean-neighboring cell, and the associated mass flux is injected through a divergence flux, yielding an increase in global volume (i.e., mean sea surface height increase). The opposite behavior occurs in the event of refreezing (heat injection, volume extraction). All constants used in the NEMO ice-shelf module are listed in Table 1.

### 2.2.3 Sea-ice model: LIM3.6

The LIM3.6 sea-ice dynamics rely on an elastic-viscous-plastic rheology formulated on a C grid (Bouillon et al., 2009, 2013). Subgrid-scale sea-ice heterogeneity is rendered with a five-category ice-thickness distribution (Bitz et al., 2001); i.e., each sea-ice field is five-fold, with distinct values over different ice-thickness ranges. The sea-ice-thickness category redistribution follows Lipscomb (2001), with inter-category boundaries located around 0.5, 1.15, 2 and 3.8 m. The choice of five categories comes from a trade-off between computational constraints and physical realism in sea-ice mean state and variability (Massonnet et al., 2019; Moreno-Chamarro et al., 2020). The sea-ice thermodynamics is based on an energy conserving scheme (Bitz and Lipscomb, 1999) and includes an explicit representation of the salt content of the sea ice (Vancoppenolle et al., 2009). The snow cover is represented with one category defined on each ice-thickness category (one snow category per ice category). The sea-ice surface albedo depends on ice surface temperature, ice thickness, snow depth and cloudiness (Grenfell and Perovich, 2004; Brandt et al., 2005). More details on LIM3.6 and the open-ocean-sea-ice coupling are given in Rousset et al. (2015) and Barthélemy et al. (2016), respectively.

## 2.3 Atmosphere and land model: COSMO-CLM² v5.0

To simulate the atmosphere, we use version 5.0 of the COSMO-CLM regional climate model. The COSMO model is a limited-area model, developed by the COnsortium for Small-scale MOdeling (COSMO) and used both for numerical weather predictions and in the context of regional climate modeling. COSMO-CLM refers to the CLimate Mode (CLM) of the COSMO model that is used for the purpose of Regional Climate Modeling. The model is developed and maintained by the German Weather Service (DWD) and is now used by a broad community of scientists gathered within the CLM community (Rockel et al., 2008).

COSMO-CLM is a non-hydrostatic atmospheric regional climate model based on primitive thermo-hydrodynamical equations. The model equations are formulated on rotated geographical coordinates and a generalized terrain-following height coordinate is used (Doms and Baldauf, 2018). The COSMO-CLM model comprises several physical parameterization schemes representing a wide range of processes such

as radiative transfer (Ritter and Geleyn, 1992), subgrid-scale turbulence and cloud microphysics (Doms et al., 2018).

COSMO-CLM is coupled to the Community Land Model version 4.5 (Oleson et al., 2013) using the OASIS3 MCT Model Coupling Toolkit (MCT) coupler (Will et al., 2017). This coupling combination, which is hereinafter named CCLM², forms the land-atmosphere basis for PARASO and is described in Davin et al. (2011). More recently, at least two CCLM² Antarctic configurations have been developed (e.g., Souverijns et al., 2019; Zentek and Heinemann, 2020). In the following, we will use AERCLOUD, the circum-Antarctic configuration described in Souverijns et al. (2019), as a comparison reference. This configuration features polar-specific recalibration of the atmosphere boundary layer scheme, the implementation of a new snow scheme and the use of a two-moment microphysics scheme for precipitation. The PARASO CCLM² version has been adapted from AERCLOUD, and their differences are detailed in Appendix C. Moreover, a list of key PARASO-specific CCLM² parameterization and tuning choices is provided in Table A2.

# 3 Coupling interfaces

In this section, we describe the three novel coupling interfaces developed for PARASO. Hereinafter, we use the following terms:

- Coupling window is a period between two coupling instances (intermodel exchange of information).

- Restart leg is a simulation “chunk”, at the end of which the execution of the models is stopped and their content is written to disk, for the simulation to be continued later as if it had not been interrupted. Restarts are typically used to reduce consecutive (time-wise) CPU requirements and queuing time on supercomputers.

- Online coupling is coupling with model executables running at the same time, and exchanging information as they are running via a coupling interface (e.g., NEMO-COSMO via OASIS).

- Offline coupling is coupling with models running one after another, writing information to disk for exchanging fields (e.g., NEMO-f.ETISh). Offline coupling episodes typically occur in between restart legs, at much lower frequency than online coupling.

Generally speaking, NEMO and CCLM² run jointly with regular exchange of information (online coupling; see Sect. 3.1). In between NEMO-CCLM² restart legs, f.ETISh runs by itself, with input data files coming from NEMO and CCLM² outputs, and provides NEMO with updated ice-shelf geometry (offline coupling; see Sect. 3.2 and 3.3). Figure 1 offers a schematic view of the fully coupled system, with Fig. A1 giving more details on the timing strategy.
<table><caption>**Table 1.** NEMO ice-shelf-related physical constants. Here, “ice” refers to the ice sheet (not the sea ice). $\Gamma_T$ and $\Gamma_S$ are taken from Jourdain et al. (2017); $\rho_i$ is taken from f.ETISh, the ice-sheet model; all other melt thermodynamics parameters have been kept to their default NEMO values. All cavity geometry parameters have been set from a trade-off between physical realism and numerical stability.</caption><thead><tr><th></th><th>Name</th><th>Description</th><th>Value</th></tr></thead><tbody><tr><td rowspan="8">Melt thermodynamics</td><td>ρ<sub>i</sub></td><td>Ice density</td><td>917 kg m<sup>-3</sup></td></tr><tr><td>Δ<sub>f</sub></td><td>Latent heat of fusion of ice</td><td>3.34 × 10<sup>5</sup> J kg<sup>-1</sup></td></tr><tr><td>c<sup>w</sup><sub>p</sub></td><td>Ocean specific heat</td><td>3992 J kg<sup>-1</sup> K<sup>-1</sup></td></tr><tr><td>c<sup>i</sup><sub>p</sub></td><td>Ice specific heat</td><td>2000 J kg<sup>-1</sup> K<sup>-1</sup></td></tr><tr><td>K<sub>i</sub></td><td>Ice heat diffusivity</td><td>1.54 × 10<sup>-6</sup> m<sup>2</sup> s<sup>-1</sup></td></tr><tr><td>C<sup>isf</sup><sub>D</sub></td><td>Top ice drag coefficient</td><td>10<sup>-3</sup></td></tr><tr><td>Γ<sub>T</sub></td><td>Default heat exchange coef.</td><td>2.21 × 10<sup>-2</sup></td></tr><tr><td>Γ<sub>S</sub></td><td>Default salt exchange coef.</td><td>6.19 × 10<sup>-4</sup></td></tr><tr><td rowspan="4">Cavity geometry</td><td>d<sup>min</sup><sub>i</sub></td><td>Min. ice-shelf draft</td><td>10 m</td></tr><tr><td>h<sup>min</sup><sub>w</sub></td><td>Min. bathymetry</td><td>20 m</td></tr><tr><td>h<sup>min</sup><sub>i</sub></td><td>Min. water column thickness in cavities</td><td>50 m</td></tr><tr><td>h<sup>bl</sup><sub>i</sub></td><td>Ice-ocean boundary layer thickness</td><td>Top cell thickness</td></tr></tbody></table>

<figure><img src="image_4.png" /><figcaption>**Figure 1.** Fully coupled system submodel hierarchy. Each box represents an Earth system component and its corresponding model. The arrow colors distinguish the different coupling methods used, with legend given in the bottom right. “SST” stands for “sea surface temperature”, “SIC” – “sea-ice concentration”, “IST” – “ice surface temperature”, “alb.” – “albedo” and “SMB” – “surface mass balance”. Above, “NEMO-OPA” represents the open-ocean part of NEMO (i.e., NEMO without LIM). Although they are represented, the NEMO–COSMO exchanged fields are not detailed, as they have not been changed compared to existing settings.</figcaption></figure>

## 3.1 COSMO–NEMO coupling interface

### 3.1.1 General description

In this section, we introduce the coupling interface between COSMO, the atmosphere model, and NEMO, the ocean–sea-ice model. We refer to “COSMO” (atmosphere model) instead of “CCLM²” (atmosphere–land model), since this coupling

interface does not involve land processes. COSMO and NEMO are coupled online through OASIS3-MCT_2.0 (Valcke et al., 2013; Valcke, 2013). While this coupling interface has been built from pre-existing configurations involving these models (Van Pham et al., 2014; Will et al., 2017), our final setting includes several novel model developments which are described here. Figure 2 schematically represents this coupling interface.

Overall, NEMO sends COSMO surface values which are then used by COSMO to compute air–sea fluxes, eventually sent back to NEMO. The exchanged fields are detailed in Sect. 3.1.4. Spatial interpolations from one grid to the other are all performed bilinearly, except wind stresses which are interpolated bi-cubically to preserve curl. Despite being derived from the same initial f.ETISh geometry, the COSMO and NEMO continental masks may differ due to rounding effects at the coastline. In such cases, the nearest ocean neighbor is used to prevent COSMO from computing fluxes from land surface temperature with an ocean surface scheme, or NEMO from receiving fluxes that had been computed from the CLM land model.

The COSMO domain is smaller than NEMO’s; hence, the ocean surface boundary condition is coupled to COSMO if possible; otherwise, outside of the COSMO domain, it is forced with atmospheric reanalyses (see Sect. 4.2 for more details). In order to avoid discontinuities and potential instabilities between both regimes, a linear transition over 20 NEMO grid points (the full grid being ≈ 600 grid points) is applied southwards from the limit of the COSMO domain. For the sake of consistency, the same reanalysis is used for forcing the COSMO lateral boundary condition and the extra-COSMO domain NEMO surface.
### 3.1.2 Timing strategy

The coupling frequency is set to the LIM time step (5400 s). Exchanged fields are time averaged over the past coupling window by OASIS. The coupling is asynchronous: NEMO receives surface fluxes from COSMO with a one-coupling-window delay, which is less than the typical forcing frequency used for running either model in stand-alone mode (in PARASO's case, 3 h). On the other hand, the surface values seen by COSMO are in real time, albeit originating from NEMO which had been constrained by time-staggered fluxes. At the end of a restart leg, the content of the coupler is written to disk in order to guarantee smooth restartability. During the first coupling window of the first restart leg, NEMO receives fluxes obtained from monthly averages of a previous coupled run.

### 3.1.3 Tiling representation in COSMO

The main coupling-related development lies in the implementation of tiling in COSMO to represent heterogeneous ocean surface cells, containing both open ocean and sea ice. Such methods have first been developed in land models (e.g., Koster and Suarez, 1992), then generalized to mixed surfaces (e.g., Best et al., 2004) and subsequently implemented in state-of-the-art surface modules (e.g., Voldoire et al., 2017). The standard COSMO configuration requires a single set of surface values over each cell, which is then used to compute vertical diffusion. The general strategy is to have COSMO compute a set of two tile-specific surface fluxes, from two distinct sets of surface properties corresponding to the open ocean and the sea ice. This is performed on all surface cells, regardless of the sea-ice concentration. The tile-specific fluxes are then sent to NEMO. The surface properties perceived by COSMO are degenerated from their tile-specific 2-fold values to a single set of values through pseudo-averaging. The pseudo-averaging consists in re-evaluating surface properties so that a given COSMO vertical column receives the concentration-wise average of the tile-specific fluxes it had computed. In that regard, the COSMO tiling implementation satisfies local air-sea energy conservation. Technical details on the pseudo-averaging operations are provided in Appendix D.

### 3.1.4 Exchanged fields

Table 2 lists the exchanged fields between NEMO and COSMO. Sea-ice surface properties are averaged over thickness categories so that only one set of sea-ice surface properties is read by COSMO, and only one set of air-sea-ice fluxes is sent from COSMO to NEMO. As described in Rousset et al. (2015), surface fluxes injected into the sea ice are re-distributed so that LIM perceives distinct fluxes over each thickness category. With this redistribution, the net solar heat flux is the same as the one that would have arisen from a

category-specific coupling and the net non-solar heat flux is a first-order development (w.r.t. surface temperature) of the category-specific case. While the wind stress computations do not take into account surface ocean currents or sea-ice velocities, their values above open ocean and sea ice are distinct as the surface roughness characterizations (and thus drag coefficients) differ. All heat fluxes are tile-specific and thus two-fold. Over each tile, the solar heat flux relies on its tile-specific albedo, with the sea-ice albedo taken as the average over ice-thickness categories from a parameterization depending on ice and snow thicknesses (Grenfell and Perovich, 2004; Brandt et al., 2005). The tiled non-solar heat fluxes contain net longwave (including direct downward atmosphere and upward surface contributions), sensible and latent contributions. The turbulent heat fluxes (sensible and latent heat) depend on surface properties, and tile-specific turbulent heat transfer coefficients are evaluated from a turbulent kinetic energy (TKE)-derived scheme (Doms et al., 2018). The latent heat fluxes are used for the diagnostic of sea-ice sublimation and open-ocean evaporation. Sea-ice sublimation is permitted but solid condensation is not: in other words, over sea ice, downward latent heat can only be negative (i.e., it is set to zero if the surface scheme assumes it is positive). Rainfall is injected as freshwater at sea surface temperature; over sea ice, it is assumed to immediately run off to the open ocean. Over sea ice, solid precipitation contributes to snow accumulation; over the open ocean, it is assumed to immediately melt, and the corresponding latent heat is removed from the local ocean surface cell. The computation of the surface temperature over sea ice requires the non-solar heat flux sensitivity (i.e., its derivative w.r.t. surface temperature) to ensure numerical stability, since the sea-ice surface temperature may significantly evolve within a coupling window (Iovino et al., 2013). The sensitivity estimate is local both time- and space-wise: it includes the temperature's influence on surface moisture but assumes the heat transfer coefficient to be constant (which is a simplification). Fluxes over sea ice are systematically computed and sent over all ocean cells using the last known (initial) sea-ice surface properties in the event that no sea ice is currently (has ever been) present in the category. This is carried out in the event that sea ice appears, through advection or thermodynamics, over one grid cell within a coupling window.

## 3.2 NEMO-f.ETISh coupling interface

### 3.2.1 General procedure

Figure 3 schematically represents the NEMO-f.ETISh coupling interface. Like most idealized and realistic settings (e.g., Timmermann and Goeller, 2017; Favier et al., 2019), the coupling between NEMO and f.ETISh is offline. NEMO runs first for a coupling time window, sending f.ETISh monthly time series of ice-shelf melt fluxes (both melt and refreezing are permitted). While NEMO computes both mass
<table><caption>**Table 2.** NEMO–COSMO list of exchanged variables.</caption><thead><tr><th></th><th>Name</th><th>Unit</th></tr></thead><tbody><tr><td rowspan="4">NEMO → COSMO</td><td>Sea surface temperature</td><td>K</td></tr><tr><td>Sea-ice surface temperature</td><td>K</td></tr><tr><td>Sea-ice albedo</td><td>–</td></tr><tr><td>Sea-ice concentrations</td><td>–</td></tr><tr><td rowspan="12">COSMO → NEMO</td><td>Wind stress over ocean</td><td>N m<sup>-2</sup></td></tr><tr><td>Wind stress over sea ice</td><td>N m<sup>-2</sup></td></tr><tr><td>Net solar heat flux over ocean</td><td>W m<sup>-2</sup></td></tr><tr><td>Net solar heat flux over sea ice</td><td>W m<sup>-2</sup></td></tr><tr><td>Net non-solar heat flux over ocean</td><td>W m<sup>-2</sup></td></tr><tr><td>Net non-solar heat flux over sea ice</td><td>W m<sup>-2</sup></td></tr><tr><td>Non-solar heat flux sensitivity over sea ice</td><td>W m<sup>-2</sup> K<sup>-1</sup></td></tr><tr><td>Snow and/or ice sublimation rate</td><td>kg m<sup>-2</sup> s<sup>-1</sup></td></tr><tr><td>Total evaporation and/or sublimation rate</td><td>kg m<sup>-2</sup> s<sup>-1</sup></td></tr><tr><td>Snowfall rate</td><td>kg m<sup>-2</sup> s<sup>-1</sup></td></tr><tr><td>Rainfall rate</td><td>kg m<sup>-2</sup> s<sup>-1</sup></td></tr></tbody></table>

<figure><img src="image_5.png" /><figcaption>**Figure 2.** COSMO–NEMO coupling interface diagram for two coupling windows. Full blue (black) arrows correspond to exchange of tile-specific (non-tiled) fields. The thick dashed vertical line delimits the transition from one coupling window to the next one. NEMO sends tile-specific surface properties to the COSMO surface module, which are used to compute tile-specific fluxes. These fluxes are then (i) sent back to NEMO on the next coupling window (hence NEMO getting delayed fluxes); (ii) pseudo-averaged to be injected into the COSMO dynamics. In addition, non-tiled fluxes (e.g., precipitation) are also sent from COSMO to NEMO. For the sake of readability, the incoming (outgoing) fluxes from coupling window $n-1$ (to coupling window $n+1$) are not represented.</figcaption></figure>

and latent heat fluxes related to the ice-shelf melting and incorporate them as ocean boundary conditions at the ice-shelf base, only the mass flux is sent to f.ETISh. The grounding line (including cavity opening and closing) and ice thickness are evolving freely with the minimum bound on ice-shelf thickness of 11 m (very thin compared to typical values) to keep the ice-shelf extent constant, so that the PARASO land–

<figure><img src="image_6.png" /><figcaption>**Figure 3.** NEMO–f.ETISh coupling procedure. The box represents the coupling loop (one iteration per NEMO–f.ETISh coupling window). Colored arrows represent the data exchange between NEMO and f.ETISh. The NEMO spinup run is used as initial conditions for the first NEMO leg only.</figcaption></figure>

sea mask does not evolve. Afterwards, NEMO waits while f.ETISh runs for the same coupling window, including the received melt rates as a boundary condition. After this, f.ETISh provides NEMO an updated ice-sheet geometry. NEMO then uses this new ice-shelf draft dataset starting from the beginning of the next coupling window. This process is iteratively repeated over each coupling window until the end of the simulation.
### 3.2.2 Timing strategy

Our implementation is using a collaborative job submission script manager for NEMO called Coral. The coupling frequency can be changed, as long as it is (all conditions must be met)

- an integer number of months;

- a multiple of the NEMO restart length; and

- a divisor of the total experiment length.

NEMO always provides monthly subshelf melt rates to f.ETISh to include seasonality regardless of the coupling frequency and receives the simulated ice-shelf draft from f.ETISh at the end of the coupling step. For the results presented in Sect. 5, the NEMO-f.ETISh coupling frequency has been set to 3 months. This was chosen so that our relatively short (2-year) experiments can include a few ocean-ice-sheet coupling occurrences.

### 3.2.3 Exchanged fields

f.ETISh sends NEMO

- `is_ice`, a Boolean variable distinguishing ice-sheet columns (true for ice shelves and grounded ice) from ocean ones (true for open-ocean);

- `is_dry`, a Boolean variable distinguishing dry columns (true for grounded ice) from wet ones (true for ice-shelf cavities and open-ocean);

- ice surface elevation; and

- water column thickness.

These fields are then interpolated with a conservative method. The resulting columns are treated differently according to the interpolated Boolean variable values with a 50 % threshold, as described in Fig. A2. With this procedure, grounding-line displacement is theoretically bounded by the NEMO grid southern limits, which, for PARASO, covers most of Antarctica (see the NEMO grid imprint in Fig. 4).

NEMO sends f.ETISh monthly mean subshelf melt rates as freshwater mass fluxes [kg m⁻² s⁻¹], which are converted to ice equivalent melt rates [m yr⁻¹] by applying an ice density of ρ_ice = 917 kg m⁻³. The mask discrepancies between NEMO and f.ETISh (e.g., arising from the NEMO-specific criteria for opening columns described in Sect. 4.1.2) are dealt with on a cavity-by-cavity basis. If a f.ETISh cavity is at least partly covered by NEMO, the NEMO melt rates are interpolated bilinearly (applied for about 84 % of all ice-shelf grid cells in f.ETISh) with potential nearest-neighbor extrapolation to cover f.ETISh cavity columns which correspond to closed columns in NEMO (applied for about 15 % of all ice-shelf grid cells). Only a few very small and dynamically insignificant ice shelves are not represented in NEMO,

summing up about 1 % of all ice-shelf grid cells in f.ETISh. Melt rates of cavity columns in NEMO which correspond to grounded ice or open ocean in f.ETISh are disregarded (applicable for about 2.5 % of the ice-shelf area in NEMO). The distribution between those procedures (84 % interpolation, 15 % extrapolation and 1 % not represented in NEMO) is dominated by the different grids of NEMO and f.ETISh and stays almost constant even for longer runs (a few decades; not shown here); the coupling frequency only has a very minor influence on this phenomenon.

### 3.2.4 Procedure for opening and closing cavity cells in NEMO

After each ocean-ice-sheet coupling episode, NEMO uses the post-processed updated f.ETISh geometry as its new subshelf cavity geometry. In PARASO, the post-ice-sheet-coupling NEMO procedure of Smith et al. (2021) is employed. For the sake of completeness, its main steps and characteristics are briefly described below.

The geometry constraints described in Table 1 are enforced at each coupling time step. The sea surface height, temperature and salinity of a newly opened ocean column are extrapolated from their neighboring cells (Favier et al., 2019). The extrapolation procedure is repeated 100 times, which yields a smoothing effect over potentially large new openings. The initial current velocities of new cells are set to zero. On all modified columns (thinned or thickened), a horizontal divergence correction is applied in order to avoid spurious vertical velocities.

Two minor caveats keep this post-coupling procedure from being conservative. First, artificial mass and heat fluxes are brought into (taken from) the ocean upon cell opening (closing). For example, when ice-sheet geometry changes lead to a NEMO cell opening (closing), the corresponding volume and internal energy are added into (extracted from) the ocean. These fluxes are not related to ice-shelf melting and refreezing, whose associated mass and heat fluxes are already included as ocean boundary conditions at the shelf base and at each NEMO time step by its cavity module, regardless of whether ocean-ice-sheet coupling is activated (see Sect. 2.2.2). Second, the divergence correction applied after each ocean-ice-sheet coupling yields phantom current velocities into (out of) the grounded ice sheet upon cell opening (closing), which implies minor artificial ocean mass and internal energy changes. Since PARASO is a regional configuration with desired simulation times of a few decades at best, this lack of conservation has been deemed acceptable.

Geometry variations lead to spurious barotropic currents dissipating over the first few days following a mesh update. Practically, this implies that the coupling numerical stability is conditioned by the amplitude of the geometry variations through one update. However, we have observed that enforcing the NEMO cavity geometry constraints described in Ta-
ble 1 keeps critical numerical instabilities from arising, even with yearly coupling windows (not shown here).

## 3.3 CCLM²–f.ETISh interface

The exchange between CCLM² and f.ETISh runs one way (from CCLM² to f.ETISh) and performed offline. CCLM² runs for one coupling time window, sending f.ETISh monthly time series of surface mass balance (SMB), which here boils down to solid precipitation minus surface sublimation. The SMB is converted from $m yr^{-1}$ by using the reference ice density of $\rho_{ice} = 917 \text{ kg m}^{-3}$. Interpolations between the CCLM² and f.ETISh grids are performed bilinearly. For the sake of simplicity, the coupling window partitioning for the CCLM²–f.ETISh interface is the same as the NEMO–f.ETISh one, and the workflow is managed with the same job submission tool (Coral).

The SMB received from CCLM² is included within the surface boundary condition of f.ETISh. However, variations in ice-sheet surface elevation are not sent back to CCLM², which is a limitation of PARASO, and the reason why this “coupling” procedure is referred to as “one way”. It should however be stressed that at the timescales PARASO has been developed for (decadal at the longest), variations in Antarctic surface topography are negligible.

# 4 The PARASO configuration

## 4.1 General description and geometry

In this subsection, we detail each subcomponent of the PARASO configuration. Table 3 contains a brief, cross-model overview, and Fig. 4 shows and compares each model’s configuration geometry.

### 4.1.1 f.ETISh configuration

The f.ETISh model is run on a regular polar stereographic grid of 701 × 701 grid points centered around the South Pole, with a horizontal resolution of 8 km. The f.ETISh time step is 1/60 yr. The domain encompasses the continental shelf and shelf break, which limits the maximum possible extension of the ice sheet. Bedrock topography, bathymetry, observed ice thickness and grounding-line position are taken from Bed-Machine Antarctica (Morlighem et al., 2020). This dataset, originally at a 500 m spatial resolution, has been resampled to the f.ETISh 8 km grid using a low-pass filter. Data were subsequently corrected to ensure grounded ice being actually grounded and floating ice actually floating. The initial ice-sheet geometry for the PARASO setting is derived from the initialization and relaxation described in Sect. 4.3. For the PARASO configuration, iceberg calving has been omitted (but they are part of the NEMO forcing; see Sect. 4.1.2). The land-ice extent is kept constant as the NEMO-COSMO interface has not been designed to deal with an evolving land–sea

mask. This is one clear limitation which prevents PARASO from being relevant at multidecadal and longer timescales: in that case, the rapidly melting ice-shelf cavities would be carved from below and reach an unrealistically thin shape. Growing ice shelves are cut off at the calving front derived from BedMachine Antarctica (Morlighem et al., 2020), but the equivalent iceberg melt flux is not injected into NEMO (an external iceberg melt dataset is used; see Sect. 4.2.2). Moreover, in accordance with the constraints listed in Table 1, shrinking ice shelves cannot become thinner than 11 m ($\approx$ 10 m draft, which is considerably smaller than typical draft values).

### 4.1.2 NEMO configuration

The NEMO PARASO setup is derived from the GO7 configuration described in Storkey et al. (2018). The ocean grid is ePERIANT025, which includes the southernmost ice-shelf cavities (Mathiot et al., 2017). With a single lateral ocean boundary at 30° S, NEMO has the widest spatial coverage of all PARASO subcomponents. The effective ocean model resolution decreases from 24 km at the 30° S boundary to 14 km at 60° S, eventually reaching 3.8 km in its southernmost cells (in the Ross cavity at approximately 86° S). The vertical discretization includes 75 levels with thickness increasing from 1 m at the surface to ~ 200 m at depth. The NEMO and LIM time steps are 900 and 5400 s, respectively. The Antarctic continental shelf bedrock bathymetry is taken from Bed-Machine Antarctica (Morlighem et al., 2020) with a linear transition to ETOPO1 (Amante and Eakins, 2009) to cover the remaining PARASO domain (northwards from roughly 63° S). For numerical stability, a minimal ocean depth of 20 m is imposed. The Antarctic surface continental mask is constant time-wise (calving is not permitted) and taken from the ice-sheet geometry obtained from the initial f.ETISh run described in Sect. 4.3. Ice-shelf cavities are opened to ocean circulation (Mathiot et al., 2017) with three geometrical constraints balancing physical realism and numerical stability (see Table 1):

1. the minimal ice-shelf draft and water column height are set to 10 and 50 m, respectively;

2. columns within ice-shelf cavities must have at least two vertical levels to allow overturning circulation;

3. subglacial lakes (i.e., floating continental ice surrounded by grounded ice and separated from the ocean), ice-shelf crevasses and “chimneys” (vertical ocean segments surrounded by ice in all directions and connected to the cavity) are removed by NEMO (i.e., they are artificially filled with ice from their nearest-neighbor ice-shelf column).

Parameters enforcing these constraints are also listed in Table 1. At the depths of the ice-shelf cavities, the vertical res-
<table><caption>**Table 3.** Summary of PARASO model setups. “res.” stands for “resolution”.</caption><thead><tr><th>Model</th><th>f.ETISh</th><th>NEMO</th><th>CCLM<sup>2</sup></th></tr></thead><tbody><tr><td>Earth component</td><td>Ice sheet</td><td>Ocean and sea ice</td><td>Atmosphere and soil</td></tr><tr><td>Domain boundary</td><td>Ice-shelf front (REMA*)</td><td>30° S</td><td>Between 50 and 40° S</td></tr><tr><td>Grid type</td><td>Stereographic projection</td><td>Quasi-isotropic bipolar</td><td>Rotated lat-long</td></tr><tr><td>Horizontal res.</td><td>8 km</td><td>0.25°</td><td>0.22°</td></tr><tr><td>Vertical res.</td><td>—</td><td>75 levels (1 m → 200 m)</td><td>60 levels (10 m → 1000 m)</td></tr><tr><td>Time steps</td><td>1/60 yr (525 600 s)</td><td>900 s (NEMO) &amp; 5400 s (LIM)</td><td>90 s (COSMO) &amp; 3600 s (CLM)</td></tr><tr><td>Top forcing</td><td>—</td><td>ERA5 (outside of CCLM<sup>2</sup> domain)</td><td>ERA5</td></tr><tr><td>Lateral forcing</td><td>—</td><td>ORAS5</td><td>ERA5</td></tr></tbody></table>

<figure><img src="image_7.png" /><figcaption>**Figure 4.** PARASO configuration geometry over (a) the full NEMO domain (cut at 30° S) and (b) Antarctica. The two color maps indicate the bedrock bathymetry and the ice surface elevation (Morlighem et al., 2020) over Antarctica (note the nonlinear scale). The lines represent the f.ETISh, NEMO and COSMO grids (for the sake of readability, not all cells are drawn).</figcaption></figure>

olution roughly ranges from 10 to 150 m (but we remind the reader that partial cells are used).

### 4.1.3 CCLM² configuration

In PARASO, the COSMO rotated latitude–longitude grid has a horizontal resolution of 0.22° (approximately 25 km at the Antarctic coastline), covering the whole Antarctic ice sheet as well as a significant part of the Southern Ocean, with the northern boundary located between 50 and 40° S. The terrain-following vertical discretization has 60 levels. The lowest model level is at 5 m height and cell thicknesses span from 10 m at the bottom to 1000 m at the top. The COSMO

and CLM time steps are 90 and 3600 s, respectively. The COSMO-CLM coupling is performed at every CLM time step (3600 s frequency).

## 4.2 External forcing

### 4.2.1 f.ETISh forcings

The only f.ETISh forcings are subshelf melt rates and SMB, which are provided by NEMO and CCLM² (respectively) and updated at every coupling window.

<footnote>* Reference elevation model of Antarctica (Howat et al., 2019).</footnote>
### 4.2.2 NEMO forcings

Over the ocean surface located outside of the CCLM² spa-
tial coverage, the NEMO surface is forced with fluxes com-
puted by the CORE bulk formula (Large and Yeager, 2004),
using the ERA5 reanalysis as atmosphere input (Hersbach
et al., 2020). The resulting air-sea fluxes are different from
the ones that the COSMO surface scheme (prescribed over
the coupled part of the domain) would have provided. This
limiting trade-off is related to the nature of the COSMO
surface scheme, requiring model-specific atmospheric input
fields which are typically not distributed in reanalyses (e.g.,
the TKE at the lowest two atmospheric levels). At its 30° S
lateral open-ocean boundary, NEMO is constrained with the
NEMO-based ORAS5 ocean reanalysis (Zuo et al., 2019)
at monthly frequency. The lateral boundary condition is en-
forced with a radiation scheme (see Eq. 9 in Flather, 1994)
for 2-D dynamics (sea surface height and barotropic veloci-
ties) and a flow relaxation scheme for the baroclinic veloci-
ties, temperature and salinity. A constant $86.4 \times 10^{-3} \text{ W m}^{-2}$
geothermal flux is imposed at the bottom of the ocean (Stein
and Stein, 1992). River runoff is adapted from Dai and
Trenberth (2002) and prescribed as a climatological dataset
(Bourdallé-Badie and Treguier, 2006). While icebergs are
not dynamically simulated, their freshwater melt input is in-
jected as additional runoff from an interannual dataset ob-
tained from a previous iceberg-including NEMO simulation
covering 1979 to 2015 (Marsh et al., 2015; Merino et al.,
2016; Jourdain et al., 2019). Buoyant plume mixing result-
ing from runoff is parameterized with enhanced mixing in
the shallowest 10 m, over the river deltas and on a stripe be-
tween 20 and 200 km from the Antarctic coast (where most
of the iceberg melting occurs). Enhanced mixing is excluded
in the direct vicinity of the coast to avoid interferences with
ice-shelf meltwater (Nicolas C. Jourdain, personal commu-
nication, 2019). No salinity restoration is applied.

### 4.2.3 CCLM² forcings

At its lateral boundaries, COSMO is relaxed to ERA5 us-
ing a one-way interactive nesting based on Davies (1976).
This consists in defining a relaxation zone where the inter-
nal COSMO-CLM solution is nudged against the external
ERA5 solution. Within this zone, the variables of the driv-
ing model are gradually combined with their corresponding
variables in COSMO-CLM by adding a relaxation forcing
term to the tendency equations that govern their evolution.
The attenuation function which specifies the lateral bound-
ary relaxation inside the boundary zone has a tangent hy-
perbolic form (Kållberg, 1977). The width of the relaxation
layer is set to 220 km, which corresponds to approximately
10 times the typical cell size. More information about the
one-way nesting can be found in the COSMO documenta-
tion (Doms and Baldauf, 2018). The boundary forcings (as
well as the initial conditions) for the COSMO model have

been prepared at discrete time intervals using the interpo-
lation program INT2LM (Schättler and Blahak, 2013). The
interval between two consecutive sets of boundary data (frequ-
quency of the lateral forcing) is 3 h. Within this time interval,
boundary values are interpolated linearly in time. The 3-D at-
mospheric variables for COSMO are wind speeds (zonal and
meridional), air temperature, pressure deviation from a ref-
erence pressure (1000 hPa at sea level), specific water vapor
content, specific cloud water content and specific cloud ice
content. No spectral nudging has been applied to the upper
atmosphere model levels to preserve the CCLM² atmosphere
dynamics. A sponge layer with Rayleigh damping in the up-
per levels of the model domain is activated in order to avoid
artificial reflections of gravity waves. A cosine damping pro-
file with maximum damping at the top and zero damping at
the base of the sponge layer (11 000 m elevation) is applied.

## 4.3 Initialization

Prior to running the fully coupled system, a f.ETISh stand-
alone initialization and relaxation run (see Sect. 4.3.1) is per-
formed in order to generate a steady Antarctic ice-sheet state
used as an initial geometry for all three subcomponents (see
Fig. 4).

### 4.3.1 f.ETISh initialization

The f.ETISh model initialization is derived from an adapted
iterative procedure based on Pollard and DeConto (2012) to
fit the model as close as possible to present-day observed
thickness (BedMachine Antarctica; Morlighem et al., 2020)
and flow field by iteratively updating local basal friction co-
efficients (Pattyn, 2017). The method has been further opti-
mized to enhance convergence (Bernales et al., 2017). In ad-
dition to Pollard and DeConto (2012), we also implemented a
regularization term that uses a 2-D Gaussian smoothing ker-
nel to filter out high-frequency noise in the basal sliding co-
efficients. For the initialization, f.ETISh is run with reduced
complexity, using the SIA for the evolution of the grounded
ice sheet, thereby keeping ice shelves in the observed ge-
ometry. Optimized basal sliding coefficients and steady-state
modeled geometry for the Antarctic ice sheet were obtained
after a forward-in-time integration of 60 000 yr with thermo-
mechanical coupling. The initial temperature field (prior to
that forward run) is derived from a steady-state temperature
solution based on the present-day surface temperature, sur-
face mass balance and geothermal heat flux. The latter is
based on the seismic-based dataset due to Shapiro and Ritz-
woller (2004). For the relaxation, ice shelves are allowed
to evolve and the hybrid shallow-shelf-shallow-ice approx-
imation (HySSA; Bueler and Brown, 2009) is used to obtain
the velocity field. The applied subshelf melt rates are taken
as the constant (time-wise) 1984–1993 mean from a NEMO
stand-alone run. The PARASO setting of f.ETISh (e.g., con-
stant calving front; compare Sect. 4) is run for 20 years to
gain the relaxed geometry (see Fig. A8), representing the initial geometry for the PARASO configuration. The duration of 20 years is the result of a trade-off between gaining quasi-steady-state ice-flow velocities and reducing the ice-shelf thinning during the relaxation. For the initialization and the relaxation, f.ETISh is forced by a constant surface temperature and SMB climatology provided by AEROCLOUD, a previous CCLM² stand-alone experiment (Souverijns et al., 2019).

### 4.3.2 NEMO initialization

The ocean is initialized from a 1-year spinup run performed with NEMO stand-alone, using forcings corresponding to the year 1999, to be as consistent as possible with ones used in PARASO: ORAS5 for the ocean lateral boundary condition and ERA5 for the ocean-atmosphere surface. As already mentioned in Sect. 4.2.2, the resulting fluxes are not identical to the coupled ones since the NEMO coupled (COSMO surface scheme) and stand-alone (CORE) bulk formulae differ. For the spinup, the 3-D initial conditions for temperature and salinity are taken from the final 10-year climatology of a cavity-including 40-year NEMO stand-alone run, so that the spinup is in equilibrium, in terms of both dynamics and thermodynamics. This allows PARASO to start from an upper-ocean state roughly matching the CCLM² initialization, which is also based on ERA5, up to the limits permitted by the differences in air-sea bulk formulae.

### 4.3.3 CCLM² initialization

The land and ice masks used by CCLM² to determine the fraction of land vs. ocean and land ice vs. ice-free areas are calculated on the CCLM² grid based on the geometry provided by f.ETISh, so that the geometry of the fully coupled setup is consistent from one subcomponent to another.

COSMO 3-D fields are initialized from the ERA5 atmospheric reanalysis (Hersbach et al., 2020). The CLM model starts from hard-coded initial conditions (Oleson et al., 2013); hence, the atmosphere-land system does not start from equilibrium. The visible and near-infrared albedos for glacier ice are respectively set to 0.80 and 0.55. Glacier temperatures are initialized at -23.15 °C. The overlying snow-pack is modeled with up to five layers, depending on the total snow depth. At the initialization, five layers are used to represent a snow water equivalent (SWE) of 1 m (the maximum snow depth in SWE), which corresponds to a 4 m thick snow layer for a bulk snow density of 250 kg m⁻³. The snow liquid water and ice content for layer *i* = 1 to 5 is set to *w*<sub>liq,i</sub> = 0 kg m⁻² and *w*<sub>ice,i</sub> = Δ*z*<sub>*i*</sub>ρ<sub>sno</sub> kg m⁻², respectively. Although a 1 m SWE snowpack is too small to fully represent the firn layer, which can extend to more than 100 m in Antarctica (van den Broeke et al., 2009), Souverijns et al. (2019) showed that this choice can lead to a decent SMB representation in CCLM² as long as the adaptations

from van Kampenhout et al. (2017) are included in the CLM snow module (except for the snowpack depth). For future PARASO-related work, keeping a maximum snow depth of 1 m is also convenient to limit the spinup phase of the snow-pack to a few years at most. Multiple surface datasets are also required to initialize CLM. They are listed in Table 2.5 of Oleson et al. (2013) and come from a variety of sources. Depending on the domain (location, resolution), datasets can be downloaded and regridded on a case-by-case basis using the CESM tools (the procedure is described in chap. 1 of Vertstein et al., 2013). Note that one specificity of PARASO lies in the modifications of the land mask and mean elevation, carried out to be consistent with the f.ETISh relaxation run (see Fig. 5). In particular, the ice-shelf surface elevation is represented, while the original CESM data simply assume all ice shelves to be exactly at sea level.

## 5 Results

In this section, diagnostics from a 2-year (2000–2001) PARASO experiment are presented and evaluated. Unless explicitly specified otherwise, the results shown are taken from the second simulated year (2001) in order to limit the influence of the initialization.

PARATMO, a CCLM² stand-alone experiment whose atmospheric experimental design is identical to PARASO's, serves as a comparison basis for atmospheric results. PAROCE is an analogous twin ocean experiment. In Sect. 5.3, PARASO is also compared to AEROCLOUD (a previous CCLM² stand-alone experiment detailed in Souverijns et al., 2019), which does not share the same experimental design. Moreover, three f.ETISh stand-alone experiments are performed. These experiments only differ in their forcing, which correspond to coupled fields in PARASO, and have been specifically designed to gradually investigate the impact of the two coupling interfaces involving f.ETISh. PARCRYO<sup>CTRL</sup> is forced with the constant SMB used during the initialization (provided by AEROCLOUD; see Sect. 4.3.1), and no subshelf melt is applied. It serves as a comparison basis since the CCLM² inputs to f.ETISh are from AEROCLOUD rather than any PARASO-like experiments, and no input whatsoever from NEMO is included. PARCRYO<sup>SMB</sup> also does not apply subshelf melt but is forced with monthly SMB from PARATMO, to isolate the effects of the CCLM²–f.ETISh coupling interface from the NEMO–f.ETISh one. Finally, PARCRYO is forced with monthly SMB from PARATMO and monthly subshelf melt rates from PAROCE. It includes inputs from both CCLM² and NEMO but does not account for ocean-ice-sheet feedback as PARASO does through the NEMO–f.ETISh coupling interface. A summary of all experiments designed for this study and used for diagnostics is provided in Table 4.

The PARASO biases are more generally discussed and put into perspective in Sect. 6. Computational aspects related
<figure><img src="image_8.png" /><figcaption>**Figure 5.** Surface elevation for (a) a previous CCLM² stand-alone model configuration (AEROCLOUD), (b) PARASO and (c) their absolute difference (PARASO-AEROCLOUD). Note that all color scales are nonlinear. In panel (c), nodes that change nature, from AEROCLOUD land to PARASO ocean (AEROCLOUD ocean to PARASO land), are displayed with purple (red) shadings.</figcaption></figure>

<table><caption>**Table 4.** List of experiments from which the presented diagnostics are derived. Besides the presence (or lack thereof) of model components and coupling interfaces, all experiments share the same design and cover the 2000–2001 period.</caption><thead><tr><th>Name</th><th>CCLM<sup>2</sup></th><th>NEMO</th><th>f.ETISh</th></tr></thead><tbody><tr><td>PARATMO*</td><td>Y</td><td>N</td><td>N</td></tr><tr><td>PAROCE</td><td>N</td><td>Y</td><td>N</td></tr><tr><td>PARCRYO<sup>CTRL</sup></td><td>N</td><td>N</td><td>Y</td></tr><tr><td>PARCRYO<sup>SMB</sup></td><td>N</td><td>N</td><td>Y</td></tr><tr><td>PARCRYO</td><td>N</td><td>N</td><td>Y</td></tr><tr><td>PARASO</td><td>Y</td><td>Y</td><td>Y</td></tr></tbody></table>

to PARASO are briefly covered in Appendix B. PARASO-specific tunings used for obtaining a realistic sea-ice seasonal cycle are described in Appendix E.

## 5.1 Ice sheet and ice shelves

Figure 6 displays the subshelf melt rates for PAROCE and PARASO in different regions. Both experiments' regional melt rate patterns agree relatively well with each other (Fig. 6b–e), but the total melt of PARASO is significantly increased in the Weddell (+31 %), Indian Ocean (+85 %) and West Pacific (+162 %) sectors (see Fig. 6a–c and f) due to the warmer and more saline eastern Antarctic surface water in PARASO (see Figs. 6g, 9k, m and 10k, m). However, the West Antarctic melt rates remain of similar magnitude from PAROCE to PARASO (−9 % in the Amundsen Sea and −16 % in the Bellingshausen Sea; see also Fig. 6a and d–e). This adds up to a 2001 average circumpolar PARASO melt rate of 1103 Gt yr⁻¹, +33 % compared with PAROCE's 829 Gt yr⁻¹. Compared to observed ice-shelf melt rates from

Rignot et al. (2013) and Adusumilli et al. (2020), both experiments underestimate subshelf melt. While the simulated melt rates for the Ross and Weddell seas match observations, the melt for the West Pacific sector, Bellingshausen Sea and Amundsen Sea in particular, is considerably lower (Fig. 6a). The underestimated melt rates likely stem from the initial PARASO ice-shelf geometry arising from the f.ETISh initialization (see Sect. 4.3 and Fig. A8), which differs more from the observed geometry in western Antarctica compared to the other sectors. A similar NEMO stand-alone run (not shown here) with an observation-based ice-sheet geometry (from BedMachine Antarctica, Morlighem et al., 2020) provides significantly larger melt rates in those regions, in better agreement with Rignot et al. (2013) and Adusumilli et al. (2020). These results support that the geometry is a main driver of the observed melt biases. The importance of cavity geometry on simulated melt rates has already been emphasized by previous studies (Schodlok et al., 2012; Rosier et al., 2018; Goldberg et al., 2019; Brisbourne et al., 2020; Goldberg et al., 2020; Wei et al., 2020).

Figure 7 compares PARASO and the PARCRYO experiments, corresponding uncoupled stand-alone f.ETISh simulations for different forcings (see Table 4). Overall, PARCRYO and PARASO display similar behavior across most coastal regions, and in particular the Antarctic Peninsula, with a moderate grounded ice-thickness increase up to a few meters (Fig. 7b–c). That thickening also occurs for PARCRYO<sup>SMB</sup> (Fig. 7d) and agrees well with the differences between the forcing SMB datasets (Fig. 7g). PARACRYO<sup>SMB</sup> and PARACRYO apply the SMB from PARATMO, while PARACRYO<sup>CTRL</sup> is forced by the SMB based on AEROCLOUD that has been used for the f.ETISh initialization (see Sect. 4.3.1). The direct comparison of PARCRYO and PARASO (see Fig. 7f) reveals a slightly larger thickening of the grounded ice in coastal areas for PARASO. This is in accordance with the increased SMB for

<footnote>* Same experimental design as PARASO except for the surface flux tiling which has not been implemented in the CCLM² stand-alone experiment.</footnote>
<figure><img src="image_9.png" /><figcaption>**Figure 6.** (a) Total ice-shelf melt rate over all sectors (see initials on indicative map) from observations (Rignot et al., 2013, and Adusumilli et al., 2020), PAROCE and PARASO. For observations, error bars correspond to measures uncertainties; for PAROCE and PARASO, to intra-annual STD. (b–c) Zoom-in on the Filchner–Ronne cavity for (b) PARASO (negative values indicate ice-shelf refreezing) and (c) PARASO–PAROCE differences (positive values indicate stronger PARASO melting). (d–e) Zoom-in for the western Antarctic sector for (d) PARASO (note the nonlinear color map scale) and (e) PARASO–PAROCE differences. (f–g) Zoom-in on the Amery Ice Shelf for (f) melt rate and (g) top conservative temperature PARASO–PAROCE differences. All simulation diagnostics are averaged over 2001 (second simulated year).</figcaption></figure>

those regions provided by CCLM² in PARASO (Fig. 7i). Figure 7j also emphasizes the relationship between the SMB forcing and the thickness evolution of the grounded ice. While the SMB is quite constant for the first 2 months, the change in ice volume (VAF) can be attributed to the dynamical response of the modeled ice sheet, as it has probably not reached a steady state. Furthermore, the presented mass changes (ΔVAF) are the mass changes due to the forcing as the VAF of the control run has been subtracted. PARCRYO and PARCRYO<sup>SMB</sup>, which are forced by the identical SMB dataset (from PARATMO), show a very similar evolution of

the grounded ice (Fig. 7e), while the change of the volume above floatation (VAF) for PARASO clearly follows the enhanced SMB. A more detailed discussion of the differences between AEROCLOUD, PARATMO and PARASO can be found in Sect. 5.3. Other ice-thickness changes occur at the grounding line (Fig. 7a–c), as the relatively coarse resolution and numerical uncertainties of the grounding-line flux induce small oscillations between neighboring grid cells in the grounding-line position. Furthermore, the experiments show a thickening of many narrow ice streams (e.g., Pine Island Glacier, Byrd Glacier, Rutford Ice Stream), and only
<figure><img src="image_10.png" /><figcaption>**Figure 7.** (a–c) Ice-thickness change during the year 2001 for (a) PARCRYO<sup>CTRL</sup>, (b) PARCRYO and (c) PARASO. (d–f) Difference of the ice-thickness change between (d) PARCRYO<sup>SMB</sup> and PARCRYO<sup>CTRL</sup>, (e) PARCRYO and PARCRYO<sup>SMB</sup> and (f) PARASO and PARCRYO. (g) Difference between the mean SMB applied for PARCRYO (SMB based on PARATMO) and for PARCRYO<sup>CTRL</sup> (SMB from initialization run based on AEROCLOUD). (h) Mean subshelf melt rates (from PAROCE) applied for PARCRYO (equivalent to Δ(PARCRYO–PARCRYO<sup>CTRL</sup>) since no subshelf melt is applied in PARCRYO<sup>CTRL</sup>). Out of consistency with panel (g), negative (positive) values indicate melting (refreezing). (i) Difference of the total mass balance between the coupled PARASO and the stand-alone PARCRYO experiment. (j) Change of the volume above floatation (VAF) and the integrated SMB and (k) change of the cavity volume and the integrated subshelf melt. (j–k) PARCRYO<sup>SMB</sup>, PARCRYO and PARASO are plotted as anomalies compared to PARCRYO<sup>CTRL</sup>. In panel (j), the PARCRYO<sup>SMB</sup> SMB anomalies are not plotted as they are identical with that of PARCRYO, by design. In panel (k), the PARCRYO<sup>SMB</sup> integrated subshelf melt anomalies are not plotted as no melt is applied for that experiment. All results are from 2001 (second simulated year).</figcaption></figure>
Thwaites Glacier shows significant grounded ice mass loss (Fig. 7a-c). The suspicious thickening is due to both the limited resolution and the initialization (see Sect. 4.3.1) and causes a certain model drift.

The forcing with subshelf melt, which has only been applied for PARCRYO (subshelf melt from PAROCE) and PARASO, has a very limited effect on the grounded ice on the short evaluation period of 1 year but clearly affects the ice shelves (Fig. 7e). The ice shelves show a moderate thinning of up to a few meters almost everywhere, except in the Bellingshausen Sea, where substantial ice-shelf thickening occurs (Fig. 7b-c), as the increase of the SMB compared to PARCRYO<sup>CTRL</sup> exceeds the applied subshelf melt (Fig. 7g and h), while for the other regions the subshelf melt is the dominating forcing for the ice shelves. Figure 7k illustrates the influence of the different SMB forcing on the ice shelves. The enhanced SMB of PARCRYO<sup>SMB</sup> compared to PARCRYO<sup>CTRL</sup> leads to a decrease of the cavity volume and hence ice-shelf thickening. Figure 7k also exhibits enhanced ice-shelf thinning (stronger cavity volume increase) of PARASO compared to PARCRYO, following the increased subshelf melt of PARASO. The increased ice-shelf thinning of PARASO mainly occurs across the East Antarctic ice sheet (Fig. 7f) and is in accordance with the increased subshelf melt of PARASO, which also mainly affects the East Antarctic ice shelves (Figs. 6a and 7i).

In terms of ice-ocean coupling, the ice sheet behaves very similarly for both coupled and uncoupled experiments (PARASO and PARCRYO, respectively), and their differences are explained by discrepancies in ice-sheet forcing (SMB from CCLM², subshelf melt rates from NEMO). No biases or additional noise due to the coupling itself have been found for the ice-sheet model. It should be kept in mind, however, that both simulations cover only 2 years, which is a short period for the relatively slow (e.g., in comparison with the ocean and atmosphere) ice-sheet system.

## 5.2 Ocean and sea ice

In this section, we compare PARASO results to observations as well as to PAROCE, a NEMO stand-alone experiment forced by ERA5, with the same experimental design (except the coupling; see Table 4).

Figure 8 displays sea-ice observations (sea-ice index from the National Snow & Ice Data Center, Fetterer et al., 2017) and diagnostics from PAROCE and PARASO. Generally speaking, the PARASO Antarctic sea-ice extent seasonal cycle shares similar features compared to PAROCE (see Fig. 8a), with larger biases. The chronology of the first-year cycle (i.e., dates of maximum and minimum) is well simulated, still with both configurations retaining too little sea ice in the summer (minimum extent at $\approx 10^6$ km² instead of $3 \times 10^6$ km²), which is a persistent bias for most coupled models (Turner et al., 2013; Mahlstein et al., 2013; Roach et al., 2018, 2020). PAROCE and PARASO also suf-

fer from another well-known bias in forced NEMO simulations, displaying too large maximal Antarctic sea-ice extent and too rapid subsequent melting (Vancoppenolle et al., 2009; Rousset et al., 2015; Barthélemy et al., 2018). The second-year PARASO seasonal cycle is considerably degraded compared with the first year, with a more significant low bias in maximum sea-ice extent ($16 \times 10^6$ km² instead of $18 \times 10^6$ km²). In contrast with the sea-ice extent, the PARASO maximum sea-ice volume is generally larger than PAROCE's (see Fig. 8b). This may be linked to the relatively stronger PARASO winds and their coastward orientations, making sea ice more compact but less extensive (see Fig. A7). Regarding the sea-ice growth cycle (Fig. 8c-e), the PAROCE and PARASO biases are similar, with the PARASO growth being delayed by about 15 d compared to PAROCE, which then translates into smaller maximum sea-ice extent. On the other hand, in many regions, the PARASO sea-ice decay cycle (Fig. 8f-h) occurs much earlier than PAROCE. The clearest examples are the Indian and western Pacific sectors, whose decay cycles start from a too small cover, and where the formation of coastal polynyas linked to strong katabatic winds (see Fig. A7) contributes to leaving only a thin sea-ice strip from mid-December on. In western Antarctica, virtually no sea ice is present in PARASO from mid-December on, whereas in observations it can persist year-round (and in PAROCE, melt 1 month later). In the Ross and Bellingshausen seas, most of the sea-ice melt is occurring more than a month earlier in PARASO compared with PAROCE, which is in better agreement with observations. In comparison with the aforementioned regions, the PARASO sea-ice pack in the Weddell sector melts later, but still leaves smaller persisting sea-ice cover compared to PAROCE. The only exception to this rapid PARASO sea-ice decay lies in front of the Ross Ice Shelf, where strong winds blowing towards the coast lead to the formation of a large and thick persisting sea-ice pack (see Fig. A7).

Figures 9 and 10 show in situ temperature and practical salinity diagnostics from the regridded World Ocean Atlas 2018 (Boyer et al., 2018), as well as those from PAROCE and PARASO. The PARASO temperature and salinity biases are generally similar to PAROCE's, even close to the air-sea interface (Figs. 9a-b and 10a-b). This can be explained by the NEMO initialization, which starts from a biased mean state corresponding to this specific configuration's equilibrium. It should however be noted that in PARASO, the eastern Antarctic subsurface is perceivably warmer and less saline than in PAROCE, matching the increased ice-shelf melt rates observed in Sect. 5.1. This is also related to the very rapid sea-ice decay in that area: from early spring, the ocean surface can undergo warming and the enhanced sea-ice melting leads to additional early freshwater release. An examination of the vertical profiles drawn in Figs. 9c-n and 10 confirms that the main origins of PARASO biases are related to the initialization, since they are similar to PAROCE. In general, deep waters are too cold ($\Delta T \approx -1$ °C) and fresh
<figure><img src="image_11.png" /><figcaption>**Figure 8.** (a) Daily sea-ice extent from observations (NSIDC-G02315), PAROCE and PARASO, over 2000–2001. (b) Same for sea-ice volume (no observations available). (c–e) Sea-ice growth progression over the 2000 freezing season from (c) observations (NSIDC-0051), (d) PAROCE and (e) PARASO. The shadings indicate the earliest time of the year at which sea ice appears after the yearly minimum (in February 2000). (f–h) Same for the sea-ice decay progression over the 2000–2001 melting season. The shadings indicate the earliest time of the year at which sea ice disappears after the yearly maximum (early September 2000). On both color bars, “pers” indicates persisting sea ice throughout each full freezing or melting process. Sea-ice presence is defined by a 15% concentration threshold.</figcaption></figure>

($\Delta P_s \approx 0.25$ psu) in both PAROCE and PARASO. One possible cause for this bias is the presence of ice-shelf meltwater injection at depth, which is a novel feature for both configurations. While these general biases are present, the vertical structure of the water column is represented well within PARASO, and interseasonal changes are of the same magnitude as in WOA18.

In 2001, the Antarctic Circumpolar Current (ACC) transport through the Drake Passage is $47 \pm 4$ Sv (1 Sv = $10^6 \text{ m}^3 \text{ s}^{-1}$) in PAROCE and $43 \pm 4$ Sv in PARASO (see also Fig. A6). This is severely underestimated compared with observations (e.g., $173.3 \pm 10.7$ Sv in Donohue et al., 2016) or other simulations (e.g., 117 Sv in Mathiot et al., 2011). The only directly identified cause of Drake ACC weakening in PARASO is the presence of spurious westwards currents along the Antarctic continental shelf (29 Sv) counterbalancing the eastwards transport, thus leading to a degraded net transport. While both PAROCE and PARASO Drake trans-

port values display large biases, which will be a tuning focus for the NEMO configuration, the fact that the values are not significantly different from each other suggests that the coupling itself does not bear a significant impact on ACC transport.

## 5.3 Atmosphere

Figure 11 shows seasonally averaged 2 m air temperature differences between PARASO and ERA5. Large differences (up to 15 °C in absolute value) are found over the ice shelves and the ocean close to Antarctica. The differences are much smaller in summer (Fig. 11a) compared to the winter season (Fig. 11c), when more sea ice is present and the atmosphere over the ice shelves is very stable. Over the Ross and Filchner–Ronne ice shelves, the 2 m air temperature is systematically lower in PARASO than in ERA5. This systematic difference is too large to be solely attributed to differences in elevation between PARASO and ERA5 over
<figure><img src="image_12.png" /></figure>

**Figure 9.** (a–b) In situ temperature biases for (a) PAROCE and (b) PARASO, with respect to the World Ocean Atlas 2018 (WOA18), averaged on the top 100 m of the ocean and the full seasonal cycle. (c–n) In situ temperature vertical profiles for WOA18, PAROCE and PARASO, in August (full) and December (dashed), over six circumpolar basins drawn on indicative maps. Note the nonlinear z axis. On the indicative maps, bright red shadings indicate the continental shelf (defined by the 1200 m isobath; see panels c, e, g, i, k, m); dark blue shadings indicate the “near-shelf” area (see panels d, f, h, j, l, n), defined as the area outside of the continental shelf and south of 60° S. On the continental shelf vertical profile panels (c, e, g, i, k, m), the purple shadings represent each sector’s vertical distribution of ice-shelf cavity presence (i.e., sector-specific, horizontally integrated area of ice-shelf-cavity-open-ocean interfaces, in normalized units). For WOA18, the data are taken from a 1995–2004 climatology; for PAROCE and PARASO, they are taken from simulation outputs for the year 2001 (second simulated year).

the ice shelves, but this effect might contribute to the tem-
perature signal visible over the Berkner and Roosevelt is-
lands. Due to the scarcity of in situ observations over the
ice shelves, a fraction of the difference between PARASO
and ERA5 may be due to uncertainties in ERA5 (Gossart

et al., 2019). Moreover, differences of comparable magni-
tude, but of opposite sign, are observed when comparing
PARASO with the Japanese 55-year Reanalysis (JRA-55; see
Fig. A3). PARASO was also compared to AERCLOUD, a
pre-existing CCLM² stand-alone experiment that was thor-
<figure><img src="image_13.png" /><figcaption>**Figure 10.** (a–b) Practical salinity biases for (a) PAROCE and (b) PARASO, with respect to WOA18, averaged on the top 100 m of the ocean and the full seasonal cycle. (c–n) Practical salinity vertical profiles for WOA18, PAROCE and PARASO, in August (full) and December (dashed), over six circumpolar basins (same as Fig. 9). Note the nonlinear z axis (following the NEMO vertical discretization). For WOA18, the data are taken from a 1995–2004 climatology; for PAROCE and PARASO, they are taken from simulation outputs for the year 2001 (second simulated year).</figcaption></figure>

oughly evaluated using in situ observational data (not shown; Souverijns et al., 2019). While it was found that both experiments behave similarly on the Filchner–Ronne Ice Shelf, the PARASO–AERCLOUD differences over the Ross Ice Shelf were similar to the PARASO–ERA5 ones, hinting that

the bias development in that sector is likely related to the coupling. The boundary layer over the ice shelves is very stable and the large differences between ERA5, JRA-55, PARASO and AERCLOUD are probably related to deficiencies in the representation of turbulent fluxes in the boundary layer,
which is known to be a considerable challenge in polar regions (e.g., Cerenzia et al., 2014; Vignon et al., 2018, and references therein). Due to the lack of data over the ice shelves, it is difficult to definitively attribute this deficiency to one of the products mentioned above.

Over the ocean, the PARASO-ERA5 differences are largest in the East Antarctic sector, where the 2 m air temperature is systematically higher in PARASO than ERA5. This is consistent with the increased ice-shelf melt rates (see Sect. 5.1) and near-surface ocean warming (see Sect. 5.2) previously highlighted in this sector. In addition to ERA5 and AEROCLOUD, PARASO was also compared to PARATMO (see Table 4), a CCLM² stand-alone run with the same tuning parameters as for PARASO, except for the surface flux tiling which has not been implemented for CCLM² stand-alone experiments (instead, sea ice is assumed to fully cover a grid cell where the surface temperature is lower than -1.8 °C). The large East Antarctic warm difference is less prevalent in PARATMO (see Fig. A4), thus suggesting an origin related to coupling. Moreover, in PARASO, large positive temperature anomalies are simulated over a restricted area in the Ross Sea. They are associated with the development of an open-ocean polynya, where the excessive vertical mixing induces the release of a substantial amount of energy. Some deep convection and open-ocean polynya formation have occurred during the past decades in the Southern Ocean, but many atmosphere-ocean-sea-ice coupled climate models that participated in the latest CMIP6 simulate them too frequently and at locations where they have not been observed (Heuzé, 2021), such as in the Ross Sea in PARASO.

Figure 12 displays the relative differences in seasonally averaged 10 m horizontal wind speed between PARASO and ERA5. Major differences are located in the Antarctic Peninsula and the Transantarctic Mountains. There, the magnitude of the PARASO winds can be more than twice the value in ERA5, with little interseasonal differences. More generally, the regions of highest wind speed differences match the regions where the elevation differs the most. The effect of subgrid-scale orography is accounted for through an effective roughness length (Lott and Miller, 1997), which is likely lower in PARASO compared to ERA5. Over the continent, outside of mountainous regions, the surface wind speed is generally lower in PARASO than ERA5. Gossart et al. (2019) showed that ERA5 underestimates the 10 m wind speed in coastal areas and the interior of Antarctica. Similar conclusions apply to PARASO.

The differences in PARASO and ERA5 SMB, calculated as total precipitation minus surface sublimation, are displayed on Fig. 13. Precipitation is the dominant source of the Antarctic ice-sheet SMB. Accordingly, the differences in SMB between PARASO and ERA5 primarily result from differences in the simulated precipitation. As surface melt runoff (only significant for areas below 1000 m elevation) is not included and snowdrift processes are not represented, the estimated SMB may depart from observations. However,

the scarcity of in situ observations prevents us from accurately estimating the AIS SMB, and we chose to compare our SMB estimate to ERA5, whose quality is discussed in Gossart et al. (2019). Moreover, in Mottram et al. (2021), SMB estimates from different Antarctic regional climate models (including CCLM²) are compared, and the observed inter-model spread suggests considerable uncertainties even in uncoupled, atmosphere-only simulations.

Compared to ERA5, PARASO provides higher SMB values over East Antarctica and the Weddell Peninsula (Fig. 13). For East Antarctica, the largest differences are located in the coastal areas, where the SMB in PARASO can be more than twice the SMB in ERA5. Those differences may partly result from a too low SMB estimate in ERA5, particularly below 500 m elevation (Gossart et al., 2019). Souverijns et al. (2019) found, when investigating the spatial pattern of the simulated SMB in AEROCLOUD with the reconstruction based on ice cores and ERA-Interim of Medley and Thomas (2018), a significant underestimation of the SMB for most of the coastal sites including the Antarctic Peninsula. A comparison of PARASO to AEROCLOUD reveals higher SMB values in those regions, suggesting a better agreement with observations. Irrespective of the comparison product used (AEROCLOUD or ERA5), PARASO simulates lower SMB values over the Ross Ice Shelf. This is primarily attributed to less precipitation in PARASO over this region (not shown). We also analyzed the SMB provided by PARATMO (see Fig. A5), for which similar results to PARASO are found. This argues that the model behavior can at least be partly related to the representation of atmospheric dynamics and not solely to the coupling.

# 6 Discussion and conclusions

In this technical paper, we have introduced PARASO, a new five-component coupled configuration for simulating the Earth system in the high latitudes of the Southern Hemisphere. Aside from the novel coupling interfaces, establishing this new tool required substantial model developments in the COSMO atmosphere model, in order to adapt its surface scheme to requirements from NEMO (mostly related to flux tiling), the ocean-sea-ice model it is coupled to in PARASO. To our knowledge, PARASO is the first publicized circum-polar Antarctic ocean-atmosphere coupled regional configuration that includes an ice-sheet model, with ice-shelf cavities explicitly resolved. The ocean-ice-sheet offline coupling interface has also been thoroughly described in this paper. However, our results suggest that at the short timescales investigated for this technical paper, the practical impact of this particular coupling interface is minor, and that the main features of PARASO would be reproduced with a similar NEMO-CCLM² coupled configuration (i.e., excluding coupling with f.ETISh). As briefly mentioned in Sect. 5.1, even at longer timescales (40 years), the ice-sheet initialization ap-
<figure><img src="image_14.png" /><figcaption>**Figure 11.** (PARASO-ERA5) differences in seasonally averaged 2 m air temperature for the second simulated year.</figcaption></figure>

<figure><img src="image_15.png" /><figcaption>**Figure 12.** Relative differences in seasonally averaged 10 m horizontal wind speed norm between ERA5 and PARASO (i.e., (PARASO-ERA5)/ERA5). Note the nonlinear color scale.</figcaption></figure>

pears to be more determining than the presence of an ocean-ice-sheet coupling interface.

In addition to the functioning of the coupling interfaces, the major PARASO achievement lies in its numerical stability. Yet, the specific biases observed in PARASO (e.g., warm and moist near-surface bias in eastern Antarctica, generalized excessive latent heat over sea ice, the representation of stable atmospheric boundary layers above ice shelves) are a significant drawback, which prevents PARASO from being considerably more skillful than some global coupled configurations. Global climate models suffer from similar biases in the high latitudes of the Southern Hemisphere (Wang et al., 2014; Schneider and Reusch, 2016) and the PARASO biases are of comparable magnitude to the differences between distinct reanalyses (see Figs. 11 and A3).

Overall, PARASO is a novel tool and further calibration could reduce the biases. For this first evaluation of our tool, the objective was to check whether the biases were affected by the coupling interfaces themselves (which they are not), rather than to correct the more general issue related to each model's biases. Extra tuning for limiting biases introduced by our choice of model combination and data input changes is however of high priority. While recognizing that clearly distinguishing each PARASO component's contribution to the biases from those purely due to their coupling interfaces is far from straightforward, we consider the former to be beyond the scope of this study.

Besides the coupling itself, the PARASO biases may be imputable to changes in each model's configuration that led to apply each component in conditions different from the

ones in which they were developed and calibrated. Yet, our choice has been to retain these changes, since they are relevant. For example, taking the Antarctic surface topography from f.ETISh as initial CCLM² geometry instead of the standard dataset made the PARASO setup more consistent intermodel-wise. However, CCLM² parameters had been tuned to this standard topographical dataset, and this might have impacted the results. In other words, we have preferred fundamental consistency over the precise agreement with current observations. Incorporating our novelties while limiting biases also identifies clear paths for model improvements. Below, we provide three potential bias sources and perspectives related to their corrections or limitations. We also refer to Appendix E for more details on tuning experiments which have led to PARASO, as the “out-of-the-box” coupling between PARASO components led to even stronger biases.

First, in the vicinity of their common interface, the ocean and atmosphere are highly sensitive to their boundary conditions (Miller et al., 1992; Large et al., 1997; Torres et al., 2019a, b). For PARASO, the interface boundary condition has been altered on both sides for enforcing COSMO-NEMO intermodel compatibility. Besides being assessed from a dynamical model instead of an external dataset (such as a reanalysis), ocean-atmosphere fluxes are computed differently in PARASO compared to either NEMO or COSMO stand-alone experiments. On the NEMO side, the ocean receives turbulent fluxes from the COSMO surface scheme instead of an ocean-calibrated bulk formula (e.g., CORE Large and Yeager, 2004). Performing NEMO stand-alone
<figure><img src="image_16.png" /><figcaption>**Figure 13.** Surface mass balance reconstruction for 2001, the second simulated year. (a) PARASO; (b) PARASO absolute difference w.r.t. ERA5; (c) PARASO relative difference w.r.t. ERA5. Note the nonlinear color scales. Contours denote elevation at 500, 2000 and 3000 m above mean sea level.</figcaption></figure>

experiments with a COSMO-derived bulk formulation is a considerable challenge, as the COSMO surface scheme requires several near-surface atmosphere properties (e.g., TKE at bottom levels or laminar transfer coefficients) as input. On the COSMO side, the atmosphere receives tiled fluxes which harbor subgrid-scale heterogeneity in surface properties. Since the COSMO surface scheme has been tuned for non-tiled fluxes (it uses instead a binary approach), the new tiling method, which is a requirement for coupling with NEMO, may lead to undesired systematic biases. Performing COSMO stand-alone experiments with tiled fluxes has not been pursued because this would have required further adapting the COSMO sources and forcing datasets. In a recent study, Heinemann et al. (2021) present a nonlinear tiling approach for COSMO (similar to that of PARASO) as an enhancement of a non-conserving method (Gutjahr et al., 2016), which also relies on sea-ice thickness (not considered in PARASO). This may be of use for future studies intending at investigating the sole impact of non-binary tiling. Alternative approaches are possible to locally reduce the biases in the air-sea fluxes, in particular by nudging them to reanalysis data, as done for instance in Ho-Hagemann et al. (2020) in a recent study coupling NEMO and COSMO. This has the advantage of providing model results closer to observations but at the cost of perturbing the atmosphere-ocean feedbacks and adding artificial energy sources and sinks at the surface, thus limiting the applications of the model outside of present-day conditions.

Second, coupled PARASO experiments suffer from excessive latent heat (and thus evaporation) over sea ice, as illustrated in Fig. E1e. Since coupled and forced snowfall rates are roughly similar (not shown), this evaporation bias is the main driver beneath the observed decrease in snow thickness over sea ice (see Fig. E1c). This leads to a positive feedback further inhibiting the presence of snow over sea ice and thus the radiative balance at the sea-ice surface: the excess of

evaporation degrades the snow cover, the surface albedo is reduced, more solar radiation is absorbed by the sea-ice-snow system; hence, even more snow melts at the surface. While some specific tuning (see Appendix E) helped reducing this bias, its positive feedback has not been fully controlled, and this aspect represents one of the main challenges for future PARASO developments.

Third, another potentially large source of biases is the updated Antarctic continent geometry used in PARASO. As described in Sect. 4.3, the initial Antarctic topography is derived from a f.ETISh relaxation run, which displays clear differences from observations, related to this model's specific biases and to the forcing datasets it relied on (SMB from COSMO, ice-shelf melt rates from NEMO). Out of inter-model consistency, all PARASO initial geometries have been directly derived from this f.ETISh relaxation run. In the atmosphere, the signature of the new topography can be seen on the differences between AEROCLOUD and PARASO (see Fig. 5 and Sect. 5.3). In the ocean, the near-Antarctic bathymetry is crucial to the warm water intrusion into the cavities, and thus to the ice-shelf melt rates (see Goldberg et al., 2019 for an overview). Results from 30-year NEMO-f.ETISh coupled experiments (without CCLM²; see Fig. A9) have suggested that at least on these timescales, the initial cavity geometry plays a much bigger role on ice-shelf melt rates than the ice-sheet-ocean coupling. This is to be expected, since the melt rates induced by the coupling on decadal timescales – and therefore the simulated cavity changes – are smaller than the ice geometry changes during the multi-millennial f.ETISh initialization (see Sect. 4.3).

The three sources of biases listed above have been mostly identified from comparison with the different stand-alone versions of each used model. However, other bias sources enhanced by the coupled nature of PARASO may also be present within this tool. As previously mentioned, the objective of this paper is to document our tool, assess its perfor-
mance as is and share it with the community. Taking advan-
tage of the numerical stability of PARASO and the exper-
tise gained in its development process, additional tuning and
calibration experiments aiming for further bias reduction are
currently at the designing stage. Finally, it may be worth not-
ing that in comparison with the short experiments of Sect. 5,
the sea-ice extent biases remain stable or are even reduced
with time in a longer fully coupled PARASO experiment (see
Fig. A10). This suggests that the ocean surface warming de-
scribed in Sect. 5 remains circumscribed to relatively small
amplitudes; hence, PARASO does not keep on drifting in the
longer term, and its performance remains reasonably satis-
factory at least in terms of Antarctic sea ice.

**Appendix A: Additional figures and tables**

<table><caption>**Table A1.** List of key NEMO parameters (excluding those related to ice-shelf cavities, which are already covered in Table 1).</caption><thead><tr><th></th><th>Description</th><th>Value</th></tr></thead><tbody><tr><td rowspan="9">Open ocean</td><td>Reference seawater density</td><td>1026 kg m<sup>-3</sup></td></tr><tr><td>Surface albedo</td><td>0.088</td></tr><tr><td>Horizontal bilaplacian eddy viscosity</td><td>-1.5 &times; 10<sup>11</sup> m<sup>4</sup> s<sup>-1</sup></td></tr><tr><td>Lateral eddy diffusivity</td><td>150 m<sup>2</sup> s<sup>-1</sup></td></tr><tr><td>Background vertical eddy viscosity</td><td>1.2 &times; 10<sup>-4</sup> m<sup>2</sup> s<sup>-1</sup></td></tr><tr><td>Background vertical eddy diffusivity</td><td>1.2 &times; 10<sup>-5</sup> m<sup>2</sup> s<sup>-1</sup></td></tr><tr><td>Minimal* TKE value</td><td>10<sup>-4</sup> m<sup>2</sup> s<sup>-2</sup></td></tr><tr><td>Bottom drag coefficient</td><td>10<sup>-3</sup></td></tr><tr><td rowspan="6">Sea ice</td><td>P*(ice strength; see Hibler, 1979)</td><td>4 &times; 10<sup>4</sup></td></tr><tr><td>Ocean–sea-ice drag coefficient</td><td>7.5 &times; 10<sup>-3</sup></td></tr><tr><td>Horizontal diffusivity</td><td>350 m<sup>2</sup> s<sup>-1</sup></td></tr><tr><td>Dry-snow albedo</td><td>0.87</td></tr><tr><td>Melting-snow albedo</td><td>0.82</td></tr><tr><td>Dry-ice albedo</td><td>0.65</td></tr></tbody></table>

<footnote>* The TKE is brought back to this threshold in the event that computations lead to lower values.</footnote>
<table><caption>**Table A2.** List of key CCLM² parameterizations.</caption><thead><tr><th>Description</th><th>Parameter</th><th>Value</th></tr></thead><tbody><tr><td rowspan="2">Turbulence: 2.5th-order Mellor and Yamada (1982) TKE closure with surface thermal heterogeneity (Doms et al., 2018)</td><td>Minimum turbulent diffusivity</td><td>0.05 m<sup>2</sup> s<sup>-1</sup></td></tr><tr><td>Length scale for subgrid thermal surface patterns</td><td>20 m</td></tr><tr><td>Clouds: prognostic ice-cloud variable</td><td>Default parameters from Doms et al. (2018)</td><td></td></tr><tr><td>Radiation scheme from Ritter and Geleyn (1992)</td><td>Subgrid-scale variability factor for ice-cloud optical thickness</td><td>0.3</td></tr><tr><td rowspan="3">Fresh snow density param. (Liston et al., 2007) with drifting-induced compaction (Vionnet et al., 2012) and capping (van Kampenhout et al., 2017)</td><td>Snow roughness length</td><td>10<sup>-5</sup> m</td></tr><tr><td>Maximum snow thickness</td><td>1 m</td></tr><tr><td>Snow density</td><td>[50; 350] kg m<sup>-3</sup></td></tr></tbody></table>

<figure><img src="image_17.png" /><figcaption>**Figure A1.** Fully coupled system timing scheme. Here, two restart legs are represented, separated by the thick horizontal line. The top and bottom timelines represent the NEMO and CCLM physical clocks, respectively. They are deliberately staggered to account for NEMO receiving delayed fluxes. For the sake of readability, above the CCLM, the time step is 1, NEMO's is 2, the NEMO-CCLM coupling frequency is 4 and the ice-sheet model coupling frequency is 8. In between legs, the NEMO-CCLM system is stopped. Black boxes represent data written on the disk. Green boxes represent OASIS operations. Full arrows represent exchange of data between models and OASIS. Dashed lines represent reading or writing to disk. “SST” stands for sea surface temperature; “IST”, ice surface temperature; “alb.”, albedo; “SIC”, sea-ice concentration; “SMB”, surface mass balance.</figcaption></figure>
<figure><img src="image_18.png" /></figure>

**Figure A2.** Pseudocode for the ice-sheet geometry post-processing. “**conserv_interp**” is a conservative interpolation from f.ETISh’s grid to NEMO’s. On the f.ETISH grid, “**is_ice**” is 0 for the open-ocean column, else (cavity or grounded ice) it is 1; “**is_wet**” is 0 for grounded ice, else (cavity or ocean) it is 1. “**draft**” is the ice-shelf draft; like the bedrock bathymetry, it is positive downwards and set to zero at the ocean surface. “**ice_elevation**” is the ice surface elevation compared to floating point level, defined positive upwards. “**init draft**”, the initial NEMO ice-shelf draft, is used to ensure that the procedure has not led to ice-shelf front displacement.

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

**Figure A3.** (PARASO-JRA-55) (Kobayashi et al., 2015) differences in seasonally averaged 2 m air temperature for the second simulated year. The used color bar is the same as that in Fig. 11. Hatching denotes areas where the PARASO-JRA-55 and PARASO-ERA5 (see Fig. 11) differences are of contrary signs.

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

**Figure A4.** The 2001 average of 2 m air temperature for (a) PARATMO, (b) PARASO and (c) the PARASO-PARATMO difference. Note that in panels (a)–(b), distinct color scales are used for the ocean (“oce”) and continent (“cont”).
<figure><img src="image_21.png" /></figure>

**Figure A5.** 2001 average of SMB for (a) PARATMO, (b) PARASO-PARATMO absolute difference and (c) PARASO-PARATMO relative difference. Note the nonlinear color scales. Contours denote elevation at 500, 2000 and 3000 m above mean sea level.

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

**Figure A6.** Barotropic streamfunction for (a) PAROCE and (b) PARASO, averaged over 2001 (second simulated year), in Sv (1 Sv = $10^6 \text{ m}^3 \text{ s}^{-1}$). The streamfunction has been set to zero over Antarctica. Positive values indicate clockwise circulation. Negative values are not represented.

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

**Figure A7.** Sea-ice thickness (color map) and air-sea-ice wind stress (arrows) in (a, b) August and (c, d) December of 2001 (second simulated year), for (a, c) PAROCE and (b, d) PARASO. Sea-ice thicknesses are only drawn on areas with sea-ice concentrations exceeding 15 %. While the wind stress scale is linear, note that the sea-ice thickness scale is not.
<figure><img src="image_24.png" /></figure>

**Figure A8.** Differences between PARASO initial draft geometry taken from the f.ETISh initialization run (see Sect. 4.3) and BedMachine (Morlighem et al., 2020) (PARASO-BedMachine) for (a) the full Antarctic and (b) a western Antarctica zoom-in. Blue nodes indicate grounded columns in the BedMachine dataset that become cavities in the PARASO initial geometry. Red nodes indicate the opposite (cavity columns turning into grounded ones).

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

**Figure A9.** Integrated global Antarctic subshelf melt rates from different experiments (lines) and observational estimates (bars). “PAROCE” is presented in Sect. 5.2. The two “NEMO–f.ETISh” experiments share the same initial bathymetry derived from the initialization described in Sect. 4.3, and feature ocean–ice-sheet coupling (no CCLM²) at distinct frequencies (3 months and 1 year). “NEMO, IBSCO” and “NEMO, BedMachine” are obtained from NEMO stand-alone experiments, with bathymetry and static cavity geometry taken from the International Bathymetric Chart of the Southern Ocean (Arndt et al., 2013) and BedMachine v2 (Morlighem et al., 2020), respectively. Observational datasets are taken from Rignot et al. (2013) and the steady state of Adusumilli et al. (2020).

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

**Figure A10.** Antarctic sea-ice extent from observations (NSIDC-G02135, Fetterer et al., 2017) and from an additional fully coupled PARASO experiment (not presented here) covering the 1985–1994 period: (a) seasonal cycle and (b) bias time series (PARASO-NSIDC-G02135). In panel (b), the gap around 1988 is due to a shortage in observational data.
# Appendix B: Computational aspects

## B1 Numerical stability

While only one 2-year simulation is presented in this paper, a considerable amount of sensitivity and tuning experiments have been performed for this study. In total, over 50 years of equivalent fully simulated years have been performed with the final PARASO source code, mostly for limiting the large bias developments. With these “latest” sources (distributed throughout this paper), no numerical instability has been observed yet, which suggests that PARASO can be considered as numerically stable. In our case, this was achieved by

- reducing the COSMO time step from its AEROCLOUD (Souverijns et al., 2019) value of 120 s to 90 s, in order to avoid CFL-criterion errors – such errors may probably be linked to the finer PARASO vertical resolution compared with AEROCLOUD (+50 %); and

- adapting the land “transfer coefficient limiter” implemented in COSMO by Davin et al. (2011). Over land, CLM computes all fluxes (including turbulent ones) and sends them to COSMO. However, the COSMO vertical physics scheme requires transfer coefficients ($C_D$, $C_H$) and surface properties ($T_s$, $q_s$) instead of fluxes *per se*. Hence, transferred coefficients are re-evaluated from CLM-originating fluxes through dividing by the near-surface temperature and moisture gradients (for sensible and latent heat, respectively). In cases where the surface temperatures (moisture) are very close to near-surface ones, this may lead to floating-point instabilities (small flux divided by small gradient), which can generate very large heat fluxes and subsequent crashes. In some of our tuning simulations, such crashes have been formally identified and traced from time-step-by-time-step COSMO outputs. Hence, in our simulation, $|C_H|$ is limited to $10^{-2}$.

## B2 Performance and load balancing

The PARASO experiments presented here have been performed on the Skylake nodes of the Flemish Tier-1 Breniac cluster. We asked for 224 cores running at 2.6 GHz. With these resources, the wall time was about 3.5 h for 1 month of simulated time. This includes all models and their I/O, as well as the offline f.ETISh coupling interfaces. Figure B1 displays the wall time required for simulating 1 month and the parallel efficiency of PARASO, defined as

$$e(n_{\text{proc}}) = \frac{n_{\text{proc}}^{\min} T(n_{\text{proc}}^{\min})}{n_{\text{proc}} T(n_{\text{proc}})}, \quad (B1)$$

where $n_{\text{proc}}$ is the number of cores used; $n_{\text{proc}}^{\min} = 56$ is the minimal number of cores required for PARASO (constrained by random access memory requirements); and $T(n_{\text{proc}})$ is the

<table><caption>**Table B1.** Intermodel load balancing. XIOS is the NEMO I/O program, which is run as a separate executable. f.ETISh is not listed here, as it is run in between NEMO–CCLM² legs.</caption><thead><tr><th rowspan="2">Model</th><th colspan="2">Core</th></tr><tr><th>No.</th><th>%</th></tr></thead><tbody><tr><td>COSMO</td><td>121</td><td>54.0</td></tr><tr><td>CLM</td><td>49</td><td>21.9</td></tr><tr><td>NEMO</td><td>52</td><td>23.2</td></tr><tr><td>XIOS</td><td>2</td><td>0.9</td></tr></tbody></table>

wall time required for 1 month of PARASO simulated time. As Fig. B1 shows, PARASO scales reasonably well up to $\approx 10^3$ cores.

Message Passing Interface (MPI)-based domain decomposition parallelism is implemented in NEMO, COSMO and CLM. The CPU balancing between these models, which is shown in Table B1, has been empirically tuned. As expected, COSMO is the biggest CPU consumer, amounting for about twice as much as NEMO or CLM, which use roughly equivalent resources. f.ETISh is not shown in Table B1, since it is run sequentially in between NEMO–CCLM² legs. Its computational cost is negligible (less than 0.1 %) compared with all other models listed in Table B1. The same applies for the f.ETISh—including offline coupling interfaces: the equivalent wall time they require is negligible compared to what is shown in Table B1. This was to be expected, as these coupling interfaces only process 2-D fields at relatively low temporal frequencies (monthly at most frequent).

## B3 Replicability and restartability

Replicability and restartability are highly desirable properties for scientific models, but achieving them remains a considerable computer science challenge in the context of climate models (Massonnet et al., 2020). Like most climate model configurations, PARASO features restart procedures, during which the content of each model’s prognostic variables is written to disk in order to interrupt the execution and limit the wall-time requirements. In PARASO, the restart procedure is also needed for the f.ETISh coupling, which is performed offline in between restarts. In other words, even with unlimited CPU time, PARASO would require restarts.

Restart procedures have been implemented in all used models, plus OASIS, which writes the content of the coupler at the end of a leg to ensure smooth restartability. Hence, PARASO restarts should theoretically be fully transparent (i.e., not have an impact on model trajectory), but empirically, it has been observed that this was not the case. Further experiments have identified COSMO as a potential source for non-restartability. Floating-point rounding errors (from double to single precision) employed in restart I/O may lead to potentially diverging trajectories (Corden and Kreitzer,
<figure><img src="image_27.png" /><figcaption>**Figure B1.** (a) Wall time required for simulating 1 month and (b) parallel efficiency of PARASO (see Eq. B1). Note that in panel (a), both axes are log scaled; in panel (b), the x axis is.</figcaption></figure>

2015). The magnitude of the differences in NEMO trajectories induced by a restart is insignificant in comparison with that induced by a COSMO restart. This was to be expected, since the atmosphere is physically much more chaotic than the ocean.

While the PARASO experiment presented above has been performed on Breniac (Tier-1 cluster from the Vlaams Supercomputer Centre), the fully coupled setup has also been successfully run on Lemaitre3 (Tier-2 cluster from the French-speaking Belgian *Consortium des Équipements de Calcul Intensif*), without numerical instabilities. The non-transparency of restart procedures on one of them (Breniac) is a strong argument supporting the non-replicability across distinct architectures. It should however be noted that, roughly speaking, simulation results from both machines were similar to each other. Moreover, outputs from identically launched experiments on one given machine–architecture combination are identical, down to bit precision.

## Appendix C: CCLM² differences between AERCLOUD and PARASO

In this appendix, we list the major differences in the CCLM² model (tuning parameters, domain, parameterizations) between AERCLOUD (a previous CCLM² stand-alone Antarctic configuration; see Souverijns et al., 2019) and PARASO.

The number of vertical levels in PARASO has been increased from 40 to 60. It allows a better representation of the stable boundary layer over the Antarctic ice sheet (Genthon et al., 2013; Handorf et al., 1999; King et al., 2006). The lowest full model level is located at 5 m (10 m) in PARASO (AERCLOUD) and 14 PARASO levels (11 in AERCLOUD) are located on the bottom 1000 m of the atmosphere. To prevent numerical instabilities related to this new vertical resolution from arising, the COSMO time step has been decreased from 120 to 90 s.

The spatial domain in AERCLOUD is limited to the Antarctic ice sheet. To study the Antarctic climate variability

at the decadal timescale, the PARASO domain has been extended to the whole Southern Ocean to explicitly simulate the interactions between the atmosphere and the ocean. The boundaries of the PARASO domain are thus located between 50 and 40° S. The horizontal resolution is unchanged (0.22°).

In order to improve the representation of perennial snow and land surface processes, some Antarctic-specific adjustments have been made in AERCLOUD. These adjustments concern

1. the representation of the snowpack in the Community Land Model following van Kampenhout et al. (2017);

2. the roughness length of snow to have a correct representation of the katabatic winds at the coastal margins of the Antarctic ice sheet and the ice shelves;

3. the turbulence scheme to account for the strong stable conditions.

More information can be found in Souverijns et al. (2019). The inclusion of the AERCLOUD two-moment cloud-microphysics scheme is disregarded in PARASO, because over this extended domain and with the increased number of vertical levels, only small improvements were found in the representation of clouds and solid precipitation. The computational resources necessary are much larger than those needed for the default parameterization scheme based on the one-moment Kessler scheme used to calculate the effects of grid-scale clouds and precipitation.

Compared to AERCLOUD, the spectral nudging has been switched off. A prime objective of future work with PARASO is to study the mesoscale interactions and their potential impact on larger-scale atmospheric and oceanographic structures. Accordingly, we did not impose any relaxation in the inner domain of the atmospheric model towards the large-scale driving model (ERA5). In this way, the CCLM² dynamics is preserved.

PARASO uses the default snow roughness length of CLM4.5 of $2.4 \times 10^{-3}$ m, which is constant in time and space. This contrasts with AERCLOUD, for which the
snow roughness length had been set to $10^{-5}$ m. Moreover, PARASO uses the standard snow and ice module of CLM4.5 (Oleson et al., 2013) and does not include the modified scheme implemented in AEROCLLOUD. Future sensitivity studies of the modified snow and ice scheme and of the value for the surface roughness length are planned.

In COSMO, the turbulence scheme uses a 1-D closure with a prognostic equation for TKE. To better represent very stable conditions over the Antarctic ice sheet, a reduction of the limit on the vertical diffusion coefficients has been considered in both AEROCLLOUD and PARASO, based on Cerenzia et al. (2014). Practically, the minimal diffusion coefficients for vertical scalar (heat) transport is set to $0.05 \text{ m}^2 \text{ s}^{-1}$ and the effective length scale of subscale surface patterns over land (used to compute the energy transfer from subgrid-scale coherent eddies to turbulent scale, and referred to as “thermal circulation term” in Cerenzia et al., 2014) is set to 20 m. Note that the values provided in Souverijns et al. (2019) for those two variables are incorrect, and that in the AEROCLLOUD model integrations the values mentioned above were used.

## Appendix D: CCLM tiling

Here, we describe the pseudo-averaging of ocean CCLM surface properties in order to conserve local energy in the subgrid-scale tiling. The pseudo-averaging operation has an impact on the CCLM-perceived surface temperature, surface moisture, turbulent transfer coefficients, laminar transfer coefficients, lowest-level TKE, lowest-level heat and momentum diffusivity, albedo and radiative surface temperature.

Throughout this appendix, the index $j = 0, 1$ distinguishes both tiles, with $j = 0$ corresponding to the open ocean and $j = 1$ to the sea ice. ($f_0, f_1 \in [0; 1]^2$ are the tile concentrations and satisfy $f_0 + f_1 = 1$.

The main novelties are linked to the TKE-derived surface transfer scheme (Doms et al., 2018). On every cell, this scheme updates

- $\kappa$, the TKE on the lowest level (defined on cell vertices), also used as input;

- $K_m$ and $K_h$, the momentum and heat diffusivities, on the lowest level (defined on cell centers), also used as input;

- $C_D$ and $C_H$, the turbulent momentum and heat transfer coefficients (in CCLM, the moisture transfer coefficient is $C_H$ as well);

- $z_0$, the surface roughness length;

- $f_h$, the scalar laminar transfer coefficient;

- $f_v$, the laminar reduction factor for evaporation;

- near-surface diagnostics: temperature, specific water vapor content, dew-point temperature, relative humidity,

winds (all values are assessed at 2 m except winds which are at 10 m).

# D1 Turbulent transfer coefficients and surface properties

The pseudo-averaging operation combines surface scheme output so that the three turbulent fluxes (momentum, latent heat and sensible heat) received are the average of all tiles weighted by their respective concentrations. The pseudo-averaged heat transfer coefficient is a simple weighted average:

$$C_H^{\text{eq}} = \overline{C_H^{j}}^{f_j}, \qquad (\text{D1})$$

where $\overline{x}_j^{f_j} = f_0 x_0 + f_1 x_1$, and ($C_H^j$) is the tile-specific heat transfer coefficients. The pseudo-averaged surface moisture is assessed by enforcing mass-flux conservation, i.e., by determining $q_s^{\text{eq}}$ so that the mass flux absorbed by CCLM is the tile-concentration average of tile-specific mass fluxes. This yields

$$q_s^{\text{eq}} = \frac{\widetilde{C}_H^{\text{eq}} q^{\text{low}} - \overline{F}_q^{j} \overline{T}_s^{j}}{\widetilde{C}_H^{\text{eq}} + \left(\frac{\mathcal{R}_v}{\mathcal{R}_d} - 1\right) \overline{F}_q^{j} \overline{T}_s^{j}}, \qquad (\text{D2})$$

$$\text{where } \widetilde{C}_H^{\text{eq}} = \frac{C_H^{\text{eq}} \|\boldsymbol{u}_h^{\text{low}}\| p_s}{\mathcal{R}_d}. \qquad (\text{D3})$$

$q_s^{\text{low}}$ is the lowest-level moisture, ($T_s^j$) the tile-specific surface temperatures, ($F_q^j$) the tile-specific mass fluxes linked to latent heat, $p_s$ the surface pressure, $\|\boldsymbol{u}_h^{\text{low}}\|$ the lowest-level horizontal wind norm, and $\mathcal{R}_d$ and $\mathcal{R}_v$ are the gas constants for dry air and water vapor, respectively. In the event that Eq. (D3) leads to $q_s^{\text{eq}} < 0$ (which in practice happens in less than once in $10^8$ calculations), we set $q_s^{\text{eq}} = 0$ and accordingly increase $C_H^{\text{eq}}$ and $\widetilde{C}_H^{\text{eq}}$ to enforce mass flux conservation.

Similarly to $q_s^{\text{eq}}$, the pseudo-averaged surface temperature $T_s^{\text{eq}}$ is assessed by enforcing sensible heat conservation, which yields

$$T_s^{\text{eq}} = \theta^{\text{low}} - \frac{\overline{Q}_h^{j} \overline{T}_s^{j}}{c_p^{\text{a}} \widetilde{C}_H^{\text{eq}}} \left(1 + \left(\frac{\mathcal{R}_v}{\mathcal{R}_d} - 1\right) q_s^{\text{eq}}\right), \qquad (\text{D4})$$

where $\theta^{\text{low}}$ is the lowest-level potential temperature and $c_p^{\text{a}}$ the dry air heat capacity. Using Eqs. (D1)–(D4) ensures mass flux, sensible heat and latent heat quasi-conservation, with minor discrepancies being imputable to their reliance on the tile-averaged surface temperature $\overline{T}_s^{j}$ instead of $T_s^{\text{eq}}$ for evaluating a small surface pressure correction term.

The equivalent momentum transfer coefficient is

$$C_D^{\text{eq}} = \frac{\|\boldsymbol{\tau}_j\|^{f_j} \mathcal{R}_d \theta_{s,v}^{\text{eq}}}{\|\boldsymbol{u}_h^{\text{low}}\|^2 p_s}, \qquad (\text{D5})$$
where $\theta_{s,v}^{eq} = T_s^{eq}(1 + (\frac{R_v}{R_d} - 1)q_s^{eq})$ is the surface virtual potential temperature and $(\tau_j)$ the tile-specific wind stresses. It should be underlined that since surface currents or ice velocities are neglected in wind stress computations, all $(\tau_j)$ values are co-aligned with $\mathbf{u}_h^{low}$; hence, averaging $\|\tau_j\|$ is transparent. Equation (D5) ensures momentum flux conservation.

The pseudo-averaged surface roughness length is assessed as a geometric average, tile-concentration-wise:

$$z_0^{eq} = \Pi_{j=0}^{1} (z_{0,j})^{f_j}. \qquad (D6)$$

A geometric average is used here because $\ln z_0^{eq}$ is the quantity of interest in terms of energy transfer.

## D2 Atmosphere boundary layer

The general procedure for getting all the other surface scheme output is then to calculate them from the turbulent transfer coefficients had a non-TKE based transfer scheme been used (e.g., Louis, 1979). In the event that such diagnostics are not readily applicable, simple algebraic averaging is performed. $\kappa_{bot}$, the bottom-level TKE velocity, is thus assessed from $C_D^{eq}$:

$$\kappa_{bot} = c_{tke} C_D^{eq} \| \mathbf{u}_h^{low} \|, \qquad (D7)$$

where $c_{tke} = \sqrt[3]{16.6}$. The bottom momentum and heat diffusivity are then updated as

$$\kappa_{m}^{\text{bot}} = \max \{ \| \mathbf{u}_{h}^{\text{low}} \|, u_{\min} \} d_{\text{bot}} C_{D}^{\text{eq}} \qquad (D8)$$

$$\kappa_{h}^{\text{bot}} = \max \{ \| \mathbf{u}_{h}^{\text{low}} \|, u_{\min} \} d_{\text{bot}} C_{H}^{\text{eq}} \qquad (D9)$$

$$\text{where } d_{\text{bot}} = z_0^{\text{eq}} \ln \left( 1 + \frac{\delta z_{\text{bot}}}{2z_0^{\text{eq}}} \right). \qquad (\text{D10})$$

In Eqs. (D8)–(D9), $u_{min} = 10^{-2} \text{ m s}^{-1}$; in Eq. (D10), $\delta z_{bot} \approx 10 \text{ m}$ is the lowest-level cell thickness. The heat laminar transfer coefficient is then reevaluated as

$$f_h = \frac{1}{1 + f_{\text{lam}}} \qquad (\text{D11})$$

$$\text{where } f_{\text{lam}} = \frac{\lambda_{\text{lam}} \kappa_{h}^{\text{bot}}}{d_{\text{bot}} \nu_{h,0}} \qquad (\text{D12})$$

$$\text{and } \lambda_{\text{lam}} = z_0^{\text{eq}} \sqrt{\frac{\nu_{m,0}}{\kappa_{m}^{\text{bot}}}}. \qquad (\text{D13})$$

In Eqs. (D12)–(D13), $\nu_{h,0} = 2.2 \times 10^{-5} \text{ m}^2 \text{ s}^{-1}$ is the scalar conductivity and $\nu_{m,0} = 1.5 \times 10^{-5} \text{ m}^2 \text{ s}^{-1}$ is the kinematic viscosity. $f_v$, the laminar reduction factor for evaporation, is evaluated as a tile-concentration-wise algebraic average.

All other near-surface diagnostics (temperature, dew-point temperature, moisture content, winds, specific humidity), which are not used in the model's dynamics, are assessed as tile-concentration-wise algebraic averages.

## D3 Radiative adjustments

Minor adjustments are also implemented in the radiative scheme transfer in order to ensure surface balance of the net shortwave and longwave radiative fluxes. The equivalent cell albedo is taken as a tile-concentration average:

$$\alpha^{\text{eq}} = f_0 \alpha_0 + f_1 \alpha_1, \qquad (\text{D14})$$

which ensures net shortwave radiation conservation. The surface temperature used in the radiative scheme is distinct from Eq. (D4). Instead, it is evaluated so that the upward longwave radiation is the same as the tile-concentration-averaged one, i.e.,

$$T_{s,\text{rad}}^{\text{eq}} = \frac{1}{\epsilon_d} \left( \epsilon_0 f_0 (T_S^0)^4 + \epsilon_1 f_1 (T_S^1)^4 \right)^{1/4}, \qquad (\text{D15})$$

where $(T_S^j)$ is the tile-specific surface temperature, $\epsilon_d = 0.996$ is the CCLM default emissivity, $\epsilon_0 = 1$ is the NEMO ocean emissivity, and $\epsilon_1 = 0.95$ is the NEMO sea-ice emissivity. Setting the radiative surface temperature as specified by Eq. (D15) ensures net longwave radiation conservation.

## Appendix E: Specific PARASO tunings

Our first fully coupled simulations suffered from much larger biases than the ones presented in Sect. 5. Therefore, several sensitivity and tuning experiments (not all shown here) have been performed in order to counterbalance and reduce the coupling-induced biases. This appendix documents the PARASO parameters that were changed from stand-alone CCLM² and NEMO simulations (see Table E1), and presents results from some of our tuning experiments (see Fig. E1). Here we present the following two experiments:

- v1: no tuning, default setup;

- v2: adjustments in the COSMO turbulent scheme; the COSMO cloud scheme; the NEMO albedo, sea-ice rheology and schemes.

More details on the differences between these two versions are given in Table E1. f.ETISh is not included here, since including it yielded no significant changes in terms of model bias. The results presented as the PARASO runs in Sect. 5 were obtained with v2.

One of the most significant problems that has been encountered when coupling NEMO and CCLM², compared to NEMO stand-alone experiments, is the degraded snow cover over sea ice. While the increase in surface albedo does slightly counterbalance the degraded snow cover, this aspect remains problematic in PARASO. Compared with the CORE bulk formula used in NEMO, systematic biases in surface humidity have been identified (dry bias for sea ice, moist bias for open ocean; see Fig. E2). Theoretically, the dry bias over
<table><caption>**Table E1.** Details on PARASO tunings. For each parameter, one of two values or methods is given: v1 (raw, non-tuned) or v2 (tuned). The “New?” column specifies whether this tuning parameter was introduced for PARASO (Y) or already existed in the model (N).</caption><thead><tr><th rowspan="2">Model</th><th rowspan="2">Comp.</th><th rowspan="2">Parameter</th><th rowspan="2">Description</th><th colspan="2">Value</th><th rowspan="2">Unit</th><th rowspan="2">New?</th></tr><tr><th>v1</th><th>v2</th></tr></thead><tbody><tr><td rowspan="3">COSMO</td><td rowspan="3">Surface turb.</td><td>Q<sup>oce</sup><sub>s</sub> (T<sup>oce</sup><sub>s</sub>)</td><td>Open-ocean surface humidity</td><td>–<sup>a</sup></td><td rowspan="3">CORE<sup>b</sup></td><td>kg kg<sup>-1</sup></td><td>Y</td></tr><tr><td>Q<sup>ice</sup><sub>s</sub> (T<sup>ice</sup><sub>s</sub>)</td><td>Sea-ice surface humidity</td><td>–<sup>a</sup></td><td>kg kg<sup>-1</sup></td><td>Y</td></tr><tr><td>C<sup>oce</sup><sub>H</sub></td><td>Open-ocean heat transfer coeff.</td><td>–<sup>a</sup> ×1.3<sup>c</sup></td><td>–</td><td>Y</td></tr><tr><td>Cloud</td><td></td><td>radqi_fact</td><td>Subgrid-scale variability factor for ice-cloud optical thickness</td><td>0.5</td><td>0.3</td><td>–</td><td>N<sup>d</sup></td></tr><tr><td rowspan="5">LIM</td><td rowspan="4">Albedo</td><td>α<sup>dry</sup><sub>sn</sub></td><td>Dry-snow albedo</td><td>85</td><td>87</td><td>%</td><td>N</td></tr><tr><td>α<sup>mlt</sup><sub>sn</sub></td><td>Melting snow albedo</td><td>75</td><td>82</td><td>%</td><td>N</td></tr><tr><td>α<sup>dry</sup><sub>ice</sub></td><td>Dry-ice albedo</td><td>60</td><td>65</td><td>%</td><td>N</td></tr><tr><td>α<sup>mlt</sup><sub>ice</sub></td><td>Bare-puddled ice albedo</td><td>50</td><td>58</td><td>%</td><td>N</td></tr><tr><td rowspan="2">Dynamics</td><td>C<sup>oce-ice</sup><sub>D</sub></td><td>Sea-ice-ocean drag coefficient</td><td>5</td><td>7.5</td><td>10<sup>-3</sup></td><td>N</td></tr><tr><td>P<sup>*</sup><sub>ice</sub></td><td>Sea-ice strength param.</td><td>2</td><td>4</td><td>10<sup>4</sup> N m<sup>-2</sup></td><td>N</td></tr><tr><td>NEMO</td><td>Albedo</td><td>α<sup>oce</sup></td><td>Open-ocean albedo</td><td>6.6</td><td>8.8</td><td>%</td><td>Y</td></tr></tbody></table>

<figure><img src="image_28.png" /><figcaption>**Figure E1.** Integrated diagnostics from several PARASO tuning experiments: (a) sea-ice extent (also featuring observations from Fetterer et al., 2017); (b) sea-ice volume; (c) average snow thickness over sea ice (only sea-ice-covered cells counted); (d) net downward surface solar radiation over sea ice; (e) downward latent heat flux over sea ice; (f) average SSTs (only sea-ice-free cells counted). Panels (c), (d) and (e) are averaged over all sea-ice-covered cells. Panel (f) is averaged over all sea-ice-free cells south of 55°S. Sea-ice coverage has been defined with a 15% concentration threshold. Results shown above are all from 1 year (2000) of experiments.</figcaption></figure>

<footnote>“Comp.”: component; “Surface turb.”: surface turbulence. <sup>a</sup> Regular COSMO scheme. <sup>b</sup> Surface humidity diagnostics from CORE bulk formula (Large and Yeager, 2004). <sup>c</sup> Regular COSMO scheme value increased by 30%. <sup>d</sup> Adapted from Muskatel et al. (2021).</footnote>
<figure><img src="image_29.png" /><figcaption>**Figure E2.** Comparison in surface moisture diagnostics (i.e., $T_s \rightarrow Q_s(T_s)$) over (a) sea ice and (b) open ocean for the COSMO and the CORE (Large and Yeager, 2004) parameterizations.</figcaption></figure>

sea ice may partly be responsible for the enhanced evaporation, yet correcting it by using the CORE humidity diagnostics only yielded slightly perceivable improvements. The COSMO heat and moisture transfer coefficient (same coefficient) has also been slightly increased to further prevent biases related to latent heat, but once again, this only yielded slight improvements. This was to be expected, as the excess in latent heat is probably related to a moist bias in the near-surface atmosphere: hence, tuning the moisture transfer coefficients only affects the speed at which the evaporation bias develops at the ocean surface, without altering the longer-term bias from developing.

All surface albedo values were slightly increased in order to limit the ocean surface warming (and sea-ice decline) that was observed in PARASO. While we admit that this was solely done to counterbalance the PARASO warm biases, the values used are all within the observational range. Moreover, default LIM3 sea-ice albedo parameters had mostly been tuned for Arctic simulations, suggesting that some Southern Ocean adaptations may indeed be required.

Finally, sea-ice dynamical properties were also slightly retuned to limit sea-ice compaction related to the relatively strong PARASO winds. One signature of this phenomenon is the relatively small bias in sea-ice volume (see Fig. 8b) compared to extent: the PARASO sea-ice pack is small in area but significant in volume; hence, it is quite thick (see Fig. A7). While the increase of these sea-ice dynamical parameters is quite significant (+50 %), the uncertainty of their actual value is also quite high (Tsamados et al., 2014; Ungermann et al., 2017).

*Code and data availability.* The full PARASO sources are available to CLM-Community members on their RedC (https://redc.clm-community.eu/ → COSMO-CLM → “Downloads”). In addition to this tarball, all PARASO sources, at the exception of their COSMO parts, are publicly available for didactic purposes at https://doi.org/10.5281/zenodo.5576201 (Pelletier et al., 2021a). All the input data for running a 3-month long PARASO experiment are publicly available at https://doi.org/10.5281/zenodo.5588468 and https://doi.org/10.5281/zenodo.5590053 (Pelletier et al., 2021b;

Pelletier and Helsen, 2021). Scripts and data for generating all figures contained herein are publicly available at https://doi.org/10.5281/zenodo.5763819 (Pelletier, 2021).

The COSMO-CLM model is free of charge for all research applications; however, access is license restricted (see http://www.cosmo-model.org/content/consortium/licencing.htm, last access: 24 August 2021). To download, the user needs to become a member of the CLM-Community or the respective institute needs to hold an institutional license. The documentation of the COSMO model is maintained by the COSMO-Consortium at https://www.dwd.de/EN/ourservices/cosmo_documentation/cosmo_documentation.html (last access: 24 August 2021).

NEMO, LIM and XIOS (a NEMO-compatible I/O library) are developed by the NEMO consortium and distributed under the CeCILL license. The NEMO-LIM version used in PARASO has been built from the standard 3.6 version, with the ice shelf following two modifications: (i) an undocumented lateral sea-ice melt scheme (Jonathan Raulier, UCLouvain); (ii) the ice-shelf coupling module from revision 11248 of the dev_isf_remapping_UKESM_GO6package_r9314 NEMO development branch, available from https://forge.ipsl.jussieu.fr/nemo browser/branches/UKMO (Mathiot and Storkey, 2018). Complete NEMO documentation is available from the NEMO consortium website: http://www.nemo-ocean.eu (last access: 21 January 2022).

f.ETISh (Fast Elementary Thermomechanical Ice Sheet model of intermediate complexity) v1.7 has been developed by Frank Pattyn and co-workers (Pattyn, 2017). This program is a free software: you can redistribute it and/or modify it under the terms of the GNU General Public License as published by the Free Software Foundation, either version 3 of the license or (at your option) any later version.

The OASIS-MCT coupler is developed by the CERFACS (Toulouse, France) and the CNRS (Paris, France). It is distributed under the Lesser GNU General Public License (LGPL). OASIS-MCT_2.0 can be downloaded from the its official website: https://portal.enes.org/oasis (Valcke, 2013).

The Coral job submission tool is developed by the CISM in the UCLouvain (Louvain-la-Neuve, Belgium) and is distributed under the Creative Commons CC0 1.0 Universal license.

The ocean-ice-sheet coupling interface also relies on the Climate Data Operators (CDO, Schulzweida, 2019) version 1.9.5, developed at the Max-Planck-Institute for Meteorology (see https://code.mpimet.mpg.de/projects/cdo/, last access: 24 August 2021), and released under the terms of the GNU General Public License v2 (GPL). NetCDF Operators (NCO; see Zender, 2008) version
4.7.6 was developed at the University of California, Irvine, and distributed under the GNU GPL v2 license (see http://nco.sourceforge.net, last access: 24 August 2021).

**Author contributions.** CP ran the numerical experiments and coordinated the PARASO developments and the writing of this manuscript. CP, SVB and SH adapted COSMO from AERO-CLOUD into PARATMO. PM provided expertise on the NEMO ice-shelf cavity module and its inclusion within PARASO. KH, CP and LZ developed the NEMO-f.ETISh offline coupling interface. CP developed the CCLM²-f.ETISh offline coupling interface. FK and EM provided technical support for the compilation and job submission tool used for PARASO. SH, SM, CP and LZ wrote the manuscript, with inputs from all co-authors.

**Competing interests.** The contact author has declared that neither they nor their co-authors have any competing interests.

**Disclaimer.** Publisher's note: Copernicus Publications remains neutral with regard to jurisdictional claims in published maps and institutional affiliations.

**Acknowledgements.** This work was supported by the PARAMOUR project, "Decadal predictability and variability of polar climate: the role of atmosphere-ocean-cryosphere multiscale interactions", supported by the Fonds de la Recherche Scientifique - FNRS and the FWO under the Excellence of Science (EOS) program (grant no. O0100718F, EOS ID no. 30454083).

The computational resources and services used in this work were provided by the VSC (Flemish Supercomputer Center), funded by the Research Foundation - Flanders (FWO) and the Flemish Government, the supercomputing facilities of the Université catholique de Louvain (CISM/UCL) and the Consortium des Équipements de Calcul Intensif en Fédération Wallonie Bruxelles (CÉCI), funded by the Fond de la Recherche Scientifique de Belgique (F.R.S.-FNRS) under convention 2.5020.11 and by the Walloon region.

Furthermore, the authors would like to thank Nicolas C. Jourdain (CNRS/IGE, France) for fostering fruitful discussions and providing the iceberg melt rate dataset; the JWCRP Joint Marine Modelling Programme for providing support and access to GO7 model configuration and output; Ulrich Blahak (German Meteorological Service) for his help with tuning COSMO; Pierre-Yves Barriat and François Damien (CISM/UCL, Belgium) for their help with Coral; and Olivier Marti for editing this paper, as well as the two anonymous referees who provided insightful comments.

The ERA5 data (Hersbach et al., 2018) were downloaded on 1 September 2019 from the Copernicus Climate Change Service (C3S) Climate Data Store. The results contain modified Copernicus Climate Change Service information from 2020. Neither the European Commission nor ECMWF is responsible for any use that may be made of the Copernicus information or data it contains.

The BedMachine Antarctica data (Morlighem et al., 2020) were downloaded on 1 February 2020 from the National Snow & Ice Data Center.

**Financial support.** This research has been supported by the Fonds De La Recherche Scientifique - FNRS (grant no. O0100718F).

**Review statement.** This paper was edited by Olivier Marti and reviewed by two anonymous referees.

## References

Adcroft, A. and Campin, J.-M.: Rescaled height coordinates for accurate representation of free-surface flows in ocean circulation models, Ocean Model., 7, 269–284, https://doi.org/10.1016/j.ocemod.2003.09.003, 2004.

Adusumilli, S., Fricker, H. A., Medley, B., Padman, L., and Siegfried, M. R.: Interannual variations in meltwater input to the Southern Ocean from Antarctic ice shelves, Nat. Geosci., 13, 616–620, https://doi.org/10.1038/s41561-020-0616-z, 2020.

Amante, C. and Eakins, B. W.: ETOPO1 arc-minute global relief model: procedures, data sources and analysis, Tech. rep., National Center for Atmospheric Research [data set], https://doi.org/10.7289/V5C8276M, 2009.

Arakawa, A.: Computational design for long-term numerical integration of the equations of fluid motion: Two-dimensional incompressible flow, Part I, J. Comput. Phys., 1, 119–143, https://doi.org/10.1016/0021-9991(66)90015-5, 1966.

Arndt, J. E., Schenke, H. W., Jakobsson, M., Nitsche, F. O., Buys, G., Goleby, B., Rebesco, M., Bohoyo, F., Hong, J., Black, J., Greku, R., Udintrsev, G., Barrios, F., Reynoso-Peralta, W., Taisei, M., and Wigley, R.: The International Bathymetric Chart of the Southern Ocean (IBCSO) Version 1.0 – A new bathymetric compilation covering circum-Antarctic waters, Geophys. Res. Lett., 40, 3111–3117, https://doi.org/10.1002/grl.50413, 2013.

Asay-Davis, X. S., Cornford, S. L., Durand, G., Galton-Fenzi, B. K., Gladstone, R. M., Gudmundsson, G. H., Hattermann, T., Holland, D. M., Holland, D., Holland, P. R., Martin, D. F., Mathiot, P., Pattyn, F., and Seroussi, H.: Experimental design for three interrelated marine ice sheet and ocean model intercomparison projects: MISMIP v. 3 (MISMIP +), ISOMIP v. 2 (ISOMIP +) and MISOMIP v. 1 (MISOMIP1), Geosci. Model Dev., 9, 2471–2497, https://doi.org/10.5194/gmd-9-2471-2016, 2016.

Asay-Davis, X. S., Jourdain, N. C., and Nakayama, Y.: Developments in Simulating and Parameterizing Interactions Between the Southern Ocean and the Antarctic Ice Sheet, Curr. Clim. Change Rep., 3, 316–329, https://doi.org/10.1007/s40641-017-0071-0, 2017.

Barnier, B., Madec, G., Penduff, T., Molines, J.-M., Treguier, A.-M., Le Sommer, J., Beckmann, A., Biastoch, A., Böning, C., Dengg, J., Derval, C., Durand, E., Gulev, S., Remy, E., Talandier, C., Theetten, S., Maltrud, M., McClean, J., and De Cuevas, B.: Impact of partial steps and momentum advection schemes in a global ocean circulation model at eddy-permitting resolution, Ocean Dynam., 56, 543–567, 2006.

Barthélemy, A., Fichefet, T., and Goosse, H.: Spatial heterogeneity of ocean surface boundary conditions under sea ice, Ocean Model., 102, 82–98, https://doi.org/10.1016/j.ocemod.2016.05.003, 2016.

Barthélemy, A., Goosse, H., Fichefet, T., and Lecomte, O.: On the sensitivity of Antarctic sea ice model biases to atmo-
spheric forcing uncertainties, Clim. Dynam., 51, 1585–1603,
https://doi.org/10.1007/s00382-017-3972-7, 2018.

Beckmann, A. and Döscher, R.: A Method for Improved Representation of Dense Water Spreading over Topography in Geopotential-Coordinate Models, J. Phys. Oceanogr., 27, 581–591, https://doi.org/10.1175/1520-0485(1997)027<0581:AMFIRO>2.0.CO;2, 1997.

Bell, R. E. and Seroussi, H.: History, mass loss, structure, and dynamic behavior of the Antarctic Ice Sheet, Science, 367, 1321–1325, https://doi.org/10.1126/science.aaz5489, 2020.

Bernales, J., Rogozhina, I., Greve, R., and Thomas, M.: Comparison of hybrid schemes for the combination of shallow approximations in numerical simulations of the Antarctic Ice Sheet, The Cryosphere, 11, 247–265, https://doi.org/10.5194/tc-11-247-2017, 2017.

Bertino, L. and Holland, M. M.: Coupled ice-ocean modeling and predictions, J. Mar. Res., 75, 839–875, https://doi.org/10.1357/002224017823524017, 2017.

Best, M. J., Beljaars, A., Polcher, J., and Viterbo, P.: A Proposed Structure for Coupling Tiled Surfaces with the Planetary Boundary Layer, J. Hydrometeor., 5, 1271–1278, https://doi.org/10.1175/JHM-382.1, 2004.

Bintanja, R., van Oldenborgh, G. J., Drijfhout, S. S., Wouters, B., and Katsman, C. A.: Important role for ocean warming and increased ice-shelf melt in Antarctic sea-ice expansion, Nat. Geosci., 6, 376–379, https://doi.org/10.1038/ngeo1767, 2013.

Bitz, C. M. and Lipscomb, W. H.: An energy-conserving thermodynamic model of sea ice, J. Geophys. Res.-Oceans, 104, 15669–15677, https://doi.org/10.1029/1999JC900100, 1999.

Bitz, C. M., Holland, M. M., Weaver, A. J., and Eby, M.: Simulating the ice-thickness distribution in a coupled climate model, J. Geophys. Res.-Oceans, 106, 2441–2463, https://doi.org/10.1029/1999JC000113, 2001.

Bouillon, S., Morales Maqueda, M. Á., Legat, V., and Fichefet, T.: An elastic–viscous–plastic sea ice model formulated on Arakawa B and C grids, Ocean Model., 27, 174–184, https://doi.org/10.1016/j.ocemod.2009.01.004, 2009.

Bouillon, S., Fichefet, T., Legat, V., and Madec, G.: The elastic–viscous–plastic method revisited, Ocean Model., 71, 2–12, https://doi.org/10.1016/j.ocemod.2013.05.013, 2013.

Bourdallé-Badie, R. and Treguier, A. M.: A climatology of runoff for the global ocean-ice model ORCA025, Tech. rep., Mercator-Ocean, available at: https://www.drakkar-ocean.eu/publications/reports/runoff-mercator-06.pdf (last access: 21 January 2022), 2006.

Boyer, T. P., Garcia, H. E., Locarnini, R. A., Zweng, M. M., MisHonov, A. V., Reagan, J. R., Weathers, K. A., Baranova, O. K., Seidov, D., and Smolyar, I. V.: World Ocean Atlas 2018, Statistical means of temperature and salinity on 0.25° grid, available at: https://accession.nodc.noaa.gov/NCEI-WOA18 (last access: 5 May 2021), 2018.

Brandt, R. E., Warren, S. G., Worby, A. P., and Grenfell, T. C.: Surface Albedo of the Antarctic Sea Ice Zone, J. Climate, 18, 3606–3622, https://doi.org/10.1175/JCLI3489.1, 2005.

Brisbourne, A., Kulessa, B., Hudson, T., Harrison, L., Holland, P., Luckman, A., Bevan, S., Ashmore, D., Hubbard, B., Pearce, E., White, J., Booth, A., Nicholls, K., and Smith, A.: An updated seabed bathymetry beneath Larsen C Ice

Shelf, Antarctic Peninsula, Earth Syst. Sci. Data, 12, 887–896,
https://doi.org/10.5194/essd-12-887-2020, 2020.

Bruno, L., Tréguier, A.-M., Madec, G., and Garnier, V.: Free surface and variable volume in the NEMO code, Tech. rep., Marine EnviRonment and Security for the European Area, https://doi.org/10.5281/zenodo.3244182, 2007.

Bueler, E. and Brown, J.: Shallow shelf approximation as a “sliding law” in a thermomechanically coupled ice sheet model, J. Geophys. Res.-Earth, 114, F03008, https://doi.org/10.1029/2008JF001179, 2009.

Cerenzia, I., Tampieri, F., and Tesini, M.: Diagnosis of Turbulence Schema in Stable Atmospheric Conditions and Sensitivity Tests, in: COSMO News Letter No. 14, available at: http://www.cosmo-model.org/content/model/documentation/newsLetters/newsLetter14/default.htm (last access: January 21 2022), 2014.

Comeau, D., Asay-Davis, X., Begeman, C., Hoffman, M., Lin, W., Petersen, M., Price, S., Roberts, A., Van Roekel, L., Veneziani, M., Wolfe, J., Fyke, J., Ringler, T., and Turner, A.: Ice-shelf basal melt rates from a global Earth system model, J. Adv. Model. Earth Sy., 14, e2021MS002468, https://doi.org/10.1029/2021MS002468, 2022.

Corden, M. J. and Kreitzer, M.: Consistency of Floating-Point Results using the Intel Compiler or Why doesn't my application always give the same answer?, Tech. rep., Intel, available at: https://software.intel.com/en-us/articles/consistency-of-floating-point-results-using-the-intel-compiler, (last access: 13 July 2021), 2015.

Dai, A. and Trenberth, K. E.: Estimates of Freshwater Discharge from Continents: Latitudinal and Seasonal Variations, J. Hydrometeorol., 3, 660–687, 2002.

Dansereau, V., Heimbach, P., and Losch, M.: Simulation of subice shelf melt rates in a general circulation model: Velocity-dependent transfer and the role of friction, J. Geophys. Res.-Oceans, 119, 1765–1790, https://doi.org/10.1002/2013JC008846, 2014.

Darelius, E., Fer, I., and Nicholls, K. W.: Observed vulnerability of Filchner-Ronne Ice Shelf to wind-driven inflow of warm deep water, Nat. Commun., 7, 12300, https://doi.org/10.1038/ncomms12300, 2016.

Davies, H. C.: A lateral boundary formulation for multi-level prediction models, Q. J. Roy. Meteor. Soc., 102, 405–418, https://doi.org/10.1002/qj.49710243210, 1976.

Davin, E. L., Stöckli, R., Jaeger, E. B., Levis, S., and Senevi-ratne, S. I.: COSMO-CLM2: A new version of the COSMO-CLM model coupled to the Community Land Model, Clim. Dynam., 37, 1889–1907, https://doi.org/10.1007/s00382-011-1019-z, 2011.

Dinniman, M., Asay-Davis, X., Galton-Fenzi, B., Holland, P., Jenkins, A., and Timmermann, R.: Modeling Ice Shelf/Ocean Interaction in Antarctica: A Review, Oceanography, 29, 144–153, https://doi.org/10.5670/oceanog.2016.106, 2016.

Doms, G. and Baldauf, M.: COSMO-Model Version 5.05: A Description of the Nonhydrostatic Regional COSMO-Model – Part II: Physical Parameterizations, Tech. rep., Consortium for Small-Scale Modelling, https://doi.org/10.5676/DWD_PUB/NWV/COSMO-DOC_5.05_I, 2018.
Doms, G., Förstner, J., Heise, E., Herzog, H.-J., Mironov, D., Raschendorfer, M., Renhardt, T., Ritter, B., Schrodin, R., Schulz, J.-P., and Vogel, G.: COSMO-Model Version 5.05: A Description of the Nonhydrostatic Regional COSMO-Model – Part I: Dynamics and Numerics, Tech. rep., Consortium for Small-Scale Modelling, https://doi.org/10.5676/DWD_PUB/NWV/COSMO-DOC_5.05_II, 2018.

Donat-Magnin, M., Jourdain, N. C., Spence, P., Sommer, J. L., Galée, H., and Durand, G.: Ice-Shelf Melt Response to Changing Winds and Glacier Dynamics in the Amundsen Sea Sector, Antarctica, J. Geophys. Res.-Oceans, 122, 10206–10224, https://doi.org/10.1002/2017JC013059, 2017.

Donohue, K. A., Tracey, K. L., Watts, D. R., Chidichimo, M. P., and Chereskin, T. K.: Mean Antarctic Circumpolar Current transport measured in Drake Passage, Geophys. Res. Lett., 43, 11760–11767, https://doi.org/10.1002/2016GL070319, 2016.

Dutrieux, P., Rydt, J. D., Jenkins, A., Holland, P. R., Ha, H. K., Lee, S. H., Steig, E. J., Ding, Q., Abrahamsen, E. P., and Schröder, M.: Strong Sensitivity of Pine Island Ice-Shelf Melting to Climatic Variability, Science, 343, 174–178, https://doi.org/10.1126/science.1244341, 2014.

Edwards, T. L., Nowicki, S., Marzeion, B., Hock, R., Goelzer, H., Seroussi, H., Jourdain, N. C., Slater, D. A., Turner, F. E., Smith, C. J., McKenna, C. M., Simon, E., Abe-Ouchi, A., Gregory, J. M., Larour, E., Lipscomb, W. H., Payne, A. J., Shepherd, A., Agosta, C., Alexander, P., Albrecht, T., Anderson, B., Asay-Davis, X., Aschwanden, A., Barthel, A., Bliss, A., Calov, R., Chambers, C., Champollion, N., Choi, Y., Cullather, R., Cuzzone, J., Dumas, C., Felikson, D., Fettweis, X., Fujita, K., Galton-Fenzi, B. K., Gladstone, R., Golledge, N. R., Greve, R., Hattermann, T., Hoffman, M. J., Humbert, A., Huss, M., Huybrechts, P., Immerzeel, W., Kleiner, T., Kraaijenbrink, P., clec'h, S. L., Lee, V., Leguy, G. R., Little, C. M., Lowry, D. P., Malles, J.-H., Martin, D. F., Maussion, F., Morlighem, M., O'Neill, J. F., Nias, I., Pattyn, F., Pelle, T., Price, S. F., Quiquet, A., Radić, V., Reese, R., Rounce, D. R., Rückamp, M., Sakai, A., Shafer, C., Schlegel, N.-J., Shannon, S., Smith, R. S., Straneo, F., Sun, S., Tarasov, L., Trusel, L. D., Breedam, J. V., van de Wal, R., van den Broeke, M., Winkelmann, R., Zekollari, H., Zhao, C., Zhang, T., and Zwinger, T.: Projected land ice contributions to twenty-first-century sea level rise, Nature, 593, 74–82, https://doi.org/10.1038/s41586-021-03302-y, 2021.

Eyring, V., Bony, S., Meehl, G. A., Senior, C. A., Stevens, B., Stouffer, R. J., and Taylor, K. E.: Overview of the Coupled Model Intercomparison Project Phase 6 (CMIP6) experimental design and organization, Geosci. Model Dev., 9, 1937–1958, https://doi.org/10.5194/gmd-9-1937-2016, 2016.

Favier, L., Jourdain, N. C., Jenkins, A., Merino, N., Durand, G., Gagliardini, O., Gillet-Chaulet, F., and Mathiot, P.: Assessment of sub-shelf melting parameterisations using the ocean-ice-sheet coupled model NEMO(v3.6)–Elmer/Ice(v8.3), Geosci. Model Dev., 12, 2255–2283, https://doi.org/10.5194/gmd-12-2255-2019, 2019.

Fetterer, F., Knowles, K., Meier, W. N., Savoie, M., and Windnagel, A. K.: Sea Ice Index, Version 3. Daily Antarctic [data set], https://doi.org/10.7265/N5K072F8, 2017.

Flather, R. A.: A Storm Surge Prediction Model for the Northern Bay of Bengal with Application to the Cyclone Disaster in April 1991, J. Phys.

Oceangr., 24, 172–190, https://doi.org/10.1175/1520-0485(1994)024<0172:ASSPMF>2.0.CO;2, 1994.

Fyke, J., Sergienko, O., Löfverström, M., Price, S., and Lenaerts, J. T. M.: An Overview of Interactions and Feedbacks Between Ice Sheets and the Earth System, Rev. Geophys., 56, 361–408, https://doi.org/10.1029/2018RG000600, 2018.

Ganopolski, A. and Brovkin, V.: Simulation of climate, ice sheets and CO2 evolution during the last four glacial cycles with an Earth system model of intermediate complexity, Clim. Past, 13, 1695–1716, https://doi.org/10.5194/cp-13-1695-2017, 2017.

Gaspar, P., Grégoris, Y., and Lefevre, J.-M.: A simple eddy kinetic energy model for simulations of the oceanic vertical mixing: Tests at station Papa and long-term upper ocean study site, J. Geophys. Res.-Oceans, 95, 16179–16193, https://doi.org/10.1029/JC095iC09p16179, 1990.

Genthon, C., Six, D., Gallée, H., Grigioni, P., and Pellegrini, A.: Two years of atmospheric boundary layer observations on a 45-m tower at Dome C on the Antarctic plateau, J. Geophys. Res.-Atmos., 118, 3218–3232, https://doi.org/10.1002/jgrd.50128, 2013.

Gierz, P., Ackermann, L., Rodehacke, C. B., Krebs-Kanzow, U., Stepanek, C., Barbi, D., and Lohmann, G.: Simulating interactive ice sheets in the multi-resolution AWI-ESM 1.2: A case study using SCOPE 1.0, Geosci. Model Dev. Discuss. [preprint], https://doi.org/10.5194/gmd-2020-159, 2020.

Gladstone, R., Galton-Fenzi, B., Gwyther, D., Zhou, Q., Hattermann, T., Zhao, C., Jong, L., Xia, Y., Guo, X., Petrakopoulos, K., Zwinger, T., Shapero, D., and Moore, J.: The Framework For Ice Sheet-Ocean Coupling (FISOC) V1.1, Geosci. Model Dev., 14, 889–905, https://doi.org/10.5194/fisoc-14-889-2021, 2021.

Goelzer, H., Huybrechts, P., Loutre, M.-F., and Fichefet, T.: Last Interglacial climate and sea-level evolution from a coupled ice sheet–climate model, Clim. Past, 12, 2195–2213, https://doi.org/10.5194/cp-12-2195-2016, 2016.

Goldberg, D. N., Gourmelen, N., Kimura, S., Millan, R., and Snow, K.: How Accurately Should We Model Ice Shelf Melt Rates?, Geophys. Res. Lett., 46, 189–199, https://doi.org/10.1029/2018GL080383, 2019.

Goldberg, D. N., Smith, T. A., Narayanan, S. H. K., Himbach, P., and Morlighem, M.: Bathymetric Influences on Antarctic Ice-Shelf Melt Rates, J. Geophys. Res.-Oceans, 125, e2020JC016370, https://doi.org/10.1029/2020JC016370, 2020.

Gorodetskaya, I. V., Tsukernik, M., Claes, K., Ralph, M. F., Neff, W. D., and van Lipzig, N. P. M.: The role of atmospheric rivers in anomalous snow accumulation in East Antarctica, Geophys. Res. Lett., 41, 6199–6206, https://doi.org/10.1002/2014GL060881, 2014.

Gossart, A., Helsen, S., Lenaerts, J. T. M., Broucke, S. V., van Lipzig, N. P. M., and Souverijns, N.: An Evaluation of Surface Climatology in State-of-the-Art Reanalyses over the Antarctic Ice Sheet, J. Climate, 32, 6899–6915, https://doi.org/10.1175/jcli-d-19-0030.1, 2019.

Grenfell, T. C. and Perovich, D. K.: Seasonal and spatial evolution of albedo in a snow-ice-land-ocean environment, J. Geophys. Res.-Oceans, 109, 8044, https://doi.org/10.1029/2003JC001866, 2004.

Griffies, S. M. and Hallberg, R. W.: Biharmonic Friction with a Smagorinsky-Like Viscosity for Use in Large-Scale Eddy-Permitting Ocean Models, Mon. Weather
Rev., 128, 2935-2946, https://doi.org/10.1175/1520-0493(2000)128<2935:BFWASL>2.0.CO;2, 2000.

Gutjahr, O., Heinemann, G., Preußer, A., Willmes, S., and Drüe, C.: Quantification of ice production in Laptev Sea polynyas and its sensitivity to thin-ice parameterizations in a regional climate model, The Cryosphere, 10, 2999-3019, https://doi.org/10.5194/tc-10-2999-2016, 2016.

Handorf, D., Foken, T., and Kottmeier, C.: The Stable Atmospheric Boundary Layer over an Antarctic ice Sheet, Bound.-Lay. Meteorol., 91, 165-189, https://doi.org/10.1023/a:1001889423449, 1999.

Hausmann, U., Sallée, J.-B., Jourdain, N. C., Mathiot, P., Rousset, C., Madec, G., Deshayes, J., and Hattermann, T.: The Role of Tides in Ocean-Ice Shelf Interactions in the Southwestern Weddell Sea, J. Geophys. Res.-Oceans, 125, e2019JC015847, https://doi.org/10.1029/2019JC015847, 2020.

Hazel, J. E. and Stewart, A. L.: Bistability of the Filchner-Ronne Ice Shelf Cavity Circulation and Basal Melt, J. Geophys. Res.-Oceans, 125, e2019JC015848, https://doi.org/10.1029/2019JC015848, 2020.

Heinemann, G., Willmes, S., Schefczyk, L., Makshtas, A., Kustov, V., and Makhotina, I.: Observations and Simulations of Meteorological Conditions over Arctic Thick Sea Ice in Late Winter during the Transarktika 2019 Expedition, Atmosphere, 12, 174, https://doi.org/10.3390/atmos12020174, 2021.

Hellmer, H. and Olbers, D.: A two-dimensional model for the thermohaline circulation under an ice shelf, Antarctic Sci., 1, 325-336, https://doi.org/10.1017/S0954102089000490, 1989.

Hellmer, H. H.: Impact of Antarctic ice shelf basal melting on sea ice and deep ocean properties, Geophys. Res. Lett., 31, L10307, https://doi.org/10.1029/2004GL019506, 2004.

Hersbach, H., Bell, B., Berrisford, P., Biavati, G., Horányi, A., Muñoz Sabater, J., Nicolas, J., Peubey, C., Radu, R., Rozum, I., Schepers, D., Simmons, A., Soci, C., Dee, D., and Thépaut, J.-N.: ERA5 hourly data on single levels from 1979 to present, Climate Data Store (Copernicus) [data set], https://doi.org/10.24381/cds.adbb2d47, 2018.

Hersbach, H., Bell, B., Berrisford, P., Hirahara, S., Horányi, A., Muñoz-Sabater, J., Nicolas, J., Peubey, C., Radu, R., Schepers, D., Simmons, A., Soci, C., Abdalla, S., Abellan, X., Balsamo, G., Bechtold, P., Biavati, G., Bidlot, J., Bonavita, M., De Chiara, G., Dahlgren, P., Dee, D., Diamantakis, M., Dragani, R., Fleming, J., Forbes, R., Fuentes, M., Geer, A., Haimberger, L., Healy, S., Hogan, R. J., Hólm, E., Janisková, M., Keeley, S., Laloyaux, P., Lopez, P., Lupu, C., Radnoti, G., de Rosnay, P., Rozum, I., Vamborg, F., Villaume, S., and Thépaut, J.-N.: The ERA5 global reanalysis, Q. J. Roy. Meteor. Soc., 146, 1999-2049, https://doi.org/10.1002/qj.3803, 2020.

Heuzé, C.: Antarctic Bottom Water and North Atlantic Deep Water in CMIP6 models, Ocean Sci., 17, 59-90, https://doi.org/10.5194/os-17-59-2021, 2021.

Hibler, W. D.: A Dynamic Thermodynamic Sea Ice Model, J. Phys. Oceanogr., 9, 815-846, https://doi.org/10.1175/1520-0485(1979)009<0815:ADTSIM>2.0.CO;2, 1979.

Ho-Hagemann, H. T. M., Hagemann, S., Grayek, S., Petrik, R., Rockel, B., Staneva, J., Feser, F., and Schrum, C.: Internal Model Variability of the Regional Coupled System Model GCOAST-AHOI, Atmosphere, 11, 227, https://doi.org/10.3390/atmos11030227, 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., Bracegirdle, T. J., Dutrieux, P., Jenkins, A., and Steig, E. J.: West Antarctic ice loss influenced by internal climate variability and anthropogenic forcing, Nat. Geosci., 12, 718-724, https://doi.org/10.1038/s41561-019-0420-9, 2019.

Howat, I. M., Porter, C., Smith, B. E., Noh, M.-J., and Morin, P.: The Reference Elevation Model of Antarctica, The Cryosphere, 13, 665-674, https://doi.org/10.5194/tc-13-665-2019, 2019.

Huot, P.-V., Fichefet, T., Jourdain, N. C., Mathiot, P., Rousset, C., Kittel, C., and Fettweis, X.: Influence of ocean tides and ice shelves on ocean-ice interactions and dense shelf water formation in the D'Urville Sea, Antarctica, Ocean Model., 162, 101794, https://doi.org/10.1016/j.ocemod.2021.101794, 2021.

IOC: The International thermodynamic equation of seawater: calculation and use of thermodynamic properties, Tech. rep., UN-ESCO, IOC/2010/MG/56 available at: https://unesdoc.unesco.org/ark:/48223/pf0000188170 (last access: 21 January 2022), 2010.

Iovino, D., Vancoppenolle, M. and Fichefet, T.: Implementation of LIM Sea Ice Model in the CMCC Global Ocean High-Resolution Configuration, CMCC Research Paper No. 209, https://doi.org/10.2139/ssrn.2639669, 2013.

Jacobs, S. S., Jenkins, A., Giulivi, C. F., and Dutrieux, P.: Stronger ocean circulation and increased melting under Pine Island Glacier ice shelf, Nat. Geosci., 4, 519-523, https://doi.org/10.1038/ngeo1188, 2011.

Jenkins, A., Nicholls, K. W., and Corr, H. F. J.: Observation and Parameterization of Ablation at the Base of Ronne Ice Shelf, Antarctica, J. Phys. Oceanogr., 40, 2298-2312, https://doi.org/10.1175/2010JPO4317.1, 2010.

Jenkins, A., Dutrieux, P., Jacobs, S., Steig, E., Gudmundsson, H., Smith, J., and Heywood, K.: Decadal Ocean Forcing and Antarctic Ice Sheet Response: Lessons from the Amundsen Sea, Oceanography, 29, 106-117, https://doi.org/10.5670/oceanog.2016.103, 2016.

Jenkins, A., Shoosmith, D., Dutrieux, P., Jacobs, S., Kim, T. W., Lee, S. H., Ha, H. K., and Stammerjohn, S.: West Antarctic Ice Sheet retreat in the Amundsen Sea driven by decadal oceanic variability, Nat. Geosci., 11, 733-738, https://doi.org/10.1038/s41561-018-0207-4, 2018.

Jeong, H., Asay-Davis, X. S., Turner, A. K., Comeau, D. S., Price, S. F., Abernathey, R. P., Veneziani, M., Petersen, M. R., Hoffman, M. J., Mazloff, M. R., and Ringler, T. D.: Impacts of Ice-Shelf Melting on Water-Mass Transformation in the Southern Ocean from E3SM Simulations, J. Climate, 33, 5787-5807, https://doi.org/10.1175/JCLI-D-19-0683.1, 2020.

Jones, J. M., Gille, S. T., Goosse, H., Abram, N. J., Canziani, P. O., Charman, D. J., Clem, K. R., Crosta, X., De Lavergne, C., Eisenman, I., England, M. H., Fogt, R. L., Frankcombe, L. M., Marshall, G. J., Masson-Delmotte, V., Morrison, A. K., Orsi, A. J., Raphael, M. N., Renwick, J. A., Schneider, D. P., Simpkins, G. R., Steig, E. J., Stenni, B., Swingedouw, D., and Vance, T. R.: Assessing recent trends in high-latitude Southern Hemisphere surface climate, Nat. Clim. Change, 6, 917-926, https://doi.org/10.1038/nclimate3103, 2016.
Jordan, J. R., Holland, P. R., Goldberg, D., Snow, K., Arth-
ern, R., Campin, J.-M., Heimbach, P., and Jenkins, A.: Ocean-Forced Ice-Shelf Thinning in a Synchronously Cou-
pled Ice-Ocean Model, J. Geophys. Res.-Oceans, 123, 864–882,
https://doi.org/10.1002/2017JC013251, 2018.

Jourdain, N. C., Mathiot, P., Merino, N., Durand, G., Sommer, J. L., Spence, P., Dutrieux, P., and Madec, G.: Ocean circulation and sea-ice thinning induced by melting ice shelves in the Amundsen Sea, J. Geophys. Res.-Oceans, 122, 2550–2573, https://doi.org/10.1002/2016JC012509, 2017.

Jourdain, N. C., Merino, N., Le Sommer, J., Durand, G., and Mathiot, P.: Interannual iceberg meltwater fluxes over the Southern Ocean, Zenodo [data set] https://doi.org/10.5281/zenodo.3514728, 2019.

Källberg, P.: Test of a lateral boundary relaxation scheme in a barotropic model, Internal Report 3, ECMWF, Shinfield Park, Reading, available at: https://www.ecmwf.int/node/10587, 1977.

King, J. C., Argentini, S. A., and Anderson, P. S.: Contrasts between the summertime surface energy balance and boundary layer structure at Dome C and Halley stations, Antarctica, J. Geophys. Res., 111, D02105, https://doi.org/10.1029/2005jd006130, 2006.

Kobayashi, S., Ota, Y., Harada, Y., Ebita, A., Moriya, M., Onoda, H., Onogi, K., Kamahori, H., Kobayashi, C., Endo, H., Miyaoka, K., and Takahashi, K.: The JRA-55 Reanalysis: General Specifications and Basic Characteristics, J. Meteorol. Soc. Jpn. Ser. II, 93, 5–48, https://doi.org/10.2151/jmsj.2015-001, 2015.

Koster, R. D. and Suarez, M. J.: Modeling the land surface boundary in climate models as a composite of independent vegetation stands, J. Geophys. Res., 97, 2697, https://doi.org/10.1029/91JD01696, 1992.

Kreuzer, M., Reese, R., Huiskamp, W. N., Petri, S., Albrecht, T., Feulner, G., and Winkelmann, R.: Coupling framework (1.0) for the PISM (1.1.4) ice sheet model and the MOM5 (5.1.0) ocean model via the PICO ice shelf cavity model in an Antarctic domain, Geosci. Model Dev., 14, 3697–3714, https://doi.org/10.5194/gmd-14-3697-2021, 2021.

Kusahara, K., Hasumi, H., Fraser, A. D., Aoki, S., Shimada, K., Williams, G. D., Massom, R., and Tamura, T.: Modeling Ocean–Cryosphere Interactions off Adélie and George V Land, East Antarctica, J. Climate, 30, 163–188, https://doi.org/10.1175/JCLI-D-15-0808.1, 2017.

Large, W. G. and Yeager, S. G.: Diurnal to Decadal Global Forcing For Ocean and Sea-Ice Models: The Data Sets and Flux Clima-
tologies, Tech. rep., National Center for Atmospheric Research,
https://doi.org/10.5065/D6KK98Q6, 2004.

Large, W. G., Danabasoglu, G., Doney, S. C., and McWilliams, J. C.: Sensitivity to Surface Forcing and Boundary Layer Mixing in a Global Ocean Model: Annual-Mean Climatology, J. Phys. Oceanogr., 27, 2418–2447, https://doi.org/10.1175/1520-0485(1997)027<2418:STSFAB>2.0.CO;2, 1997.

Lenaerts, J. T. M., Medley, B., van den Broeke, M. R., and Wouters, B.: Observing and Modeling Ice Sheet Surface Mass Balance, Rev. Geophys., 57, 376–420, https://doi.org/10.1029/2018RG000622, 2019.

Lipscomb, W. H.: Remapping the thickness distribution in sea ice models, J. Geophys. Res.-Oceans, 106, 13989–14000,
https://doi.org/10.1029/2000JC000518, 2001.

Liston, G. E., Haehnel, R. B., Sturm, M., Hiemstra, C. A., Bere-
zovskaya, S., and Tabler, R. D.: Simulating complex snow distri-
butions in windy environments using SnowTran-3D, J. Glaciol.,
53, 241–256, https://doi.org/10.3189/172756507782202865,
2007.

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

Lott, F. and Miller, M. J.: A new subgrid-scale orographic drag parametrization: Its formulation and testing, Q. J. Roy. Meteor. Soc., 123, 101–127, https://doi.org/10.1002/qj.49712353704, 1997.

Louis, J.-F.: A parametric model of vertical eddy fluxes in the atmosphere, Bound.-Lay. Meteorol., 17, 187–202, https://doi.org/10.1007/BF00117978, 1979.

Madec, G., Bourdallé-Badie, R., Bouttier, P.-A., Bricaud, C., Bruciaferri, D., Calvert, D.,nut, J., Clementi, E., Cow-
ard, A., Delrosso, D., Ethé, C., Flavoni, S., Graham, T.,
Harle, J., Iovino, D., Lea, D., Lévy, C., Lovato, T., Martin,
N., Masson, S., Mocavero, S., Paul, J., Rousset, C., Storkey,
D., Storto, A., and Vancoppenolle, M.: NEMO ocean en-
gine, Tech. rep., Insitut Pierre-Simon Laplace, Zenodo [code],
https://doi.org/10.5281/zenodo.3248739, 2017.

Mahlstein, I., Gent, P. R., and Solomon, S.: Historical Antarc-
tic mean sea ice area, sea ice trends, and winds in
CMIP5 simulations, J. Geophys. Res.-Atmos., 118, 5105–5110,
https://doi.org/10.1002/jgrd.50443, 2013.

Marsh, R., Ivchenko, V. O., Skliris, N., Alderson, S., Bigg, G. R., Madec, G., Blaker, A. T., Aksenov, Y., Sinha, B., Coward, A. C., Le Sommer, J., Merino, N., and Zalesny, V. B.: NEMO-ICB (v1.0): interactive icebergs in the NEMO ocean model globally configured at eddy-permitting resolution, Geosci. Model Dev., 8, 1547–1562, https://doi.org/10.5194/gmd-8-1547-2015, 2015.

Massonnet, F., Barthélemy, A., Worou, K., Fichefet, T., Vancop-
penolle, M., Rousset, C., and Moreno-Chamarro, E.: On the dis-
cretization of the ice thickness distribution in the NEMO3.6-
LIM3 global ocean–sea ice model, Geosci. Model Dev., 12,
3745–3758, https://doi.org/10.5194/gmd-12-3745-2019, 2019.

Massonnet, F., Ménégoz, M., Acosta, M., Yepes-Arbós, X., Exarchou, E., and Doblas-Reyes, F. J.: Replicability of the EC-Earth3 Earth system model under a change in computing environment, Geosci. Model Dev., 13, 1165–1178, https://doi.org/10.5194/gmd-13-1165-2020, 2020.

Mathiot, P. and Storkey, D.: NEMO model code, MetOffice (UK) branch dev_isf_remapping_UKESM_GO6package_r9314, revision 11248, MetOffice [code], available at: https://forge.ipsl.jussieu.fr/nemo browser/branches/UKMO/dev_isf_remapping_UKESM_GO6package_r9314?rev=15667 (last access 21 January 2022), 2018.

Mathiot, P., Goosse, H., Fichefet, T., Barnier, B., and Gallée, H.: Modelling the seasonal variability of the Antarctic Slope Current, Ocean Sci., 7, 455–470, https://doi.org/10.5194/os-7-455-2011, 2011.

Mathiot, P., Jenkins, A., Harris, C., and Madec, G.: Explicit repre-
sentation 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.

Medley, B. and Thomas, E. R.: Increased snowfall over the Antarc-
tic Ice Sheet mitigated twentieth-century sea-level rise, Nat.
Clim. Change, 9, 34–39, https://doi.org/10.1038/s41558-018-0356-x, 2018.

Mellor, G. L. and Yamada, T.: Development of a turbulence closure model for geophysical fluid problems, Rev. Geophys., 20, 851, https://doi.org/10.1029/rg020i004p00851, 1982.

Merino, N., Le Sommer, J., Durand, G., Jourdain, N. C., Madec, G., Mathiot, P., and Tournadre, J.: Antarctic icebergs melt over the Southern Ocean: Climatology and impact on sea ice, Ocean Model., 104, 99–110, https://doi.org/10.1016/j.ocemod.2016.05.001, 2016.

Merryfield, W. J., Holloway, G., and Gargett, A. E.: A Global Ocean Model with Double-Diffusive Mixing, J. Phys. Oceanogr., 29, 1124–1142, https://doi.org/10.1175/1520-0485(1999)029<1124:AGOMWD>2.0.CO;2, 1999.

Miller, M. J., Beljaars, A. C. M., and Palmer, T. N.: The Sensitivity of the ECMWF Model to the Parameterization of Evaporation from the Tropical Oceans, J. Climate, 5, 418–434, https://doi.org/10.1175/1520-0442(1992)005<0418:TSOTEM>2.0.CO;2, 1992.

Moreno-Chamarro, E., Ortega, P., and Massonnet, F.: Impact of the ice thickness distribution discretization on the sea ice concentration variability in the NEMO3.6–LIM3 global ocean–sea ice model, Geosci. Model Dev., 13, 4773–4787, https://doi.org/10.5194/gmd-13-4773-2020, 2020.

Morlighem, M., Rignot, E., Binder, T., Blankenship, D., Drews, R., Eagles, G., Eisen, O., Ferraccioli, F., Forsberg, R., Fretwell, P., Goel, V., Greenbaum, J. S., Gudmundsson, H., Guo, J., Helm, V., Hofstede, C., Howat, I., Humbert, A., Jokat, W., Karlsson, N. B., Lee, W. S., Matsuoka, K., Millan, R., Mouginot, J., Paden, J., Pattyn, F., Roberts, J., Rosier, S., Ruppel, A., Seroussi, H., Smith, E. C., Steinhage, D., Sun, B., van den Broeke, M. R., van Ommen, T. D., van Wessem, M., and Young, D. A.: Deep glacial troughs and stabilizing ridges unveiled beneath the margins of the Antarctic ice sheet, Nat. Geosci., 13, 132–137, 2020.

Mottram, R., Hansen, N., Kittel, C., van Wessem, J. M., Agosta, C., Amory, C., Boberg, F., van de Berg, W. J., Fettweis, X., Gossart, A., van Lipzig, N. P. M., van Meijgaard, E., Orr, A., Phillips, T., Webster, S., Simonsen, S. B., and Souverijns, N.: What is the surface mass balance of Antarctica? An intercomparison of regional climate model estimates, The Cryosphere, 15, 3751–3784, https://doi.org/10.5194/tc-15-3751-2021, 2021.

Muntjewerf, L., Sacks, W. J., Lofverstrom, M., Fyke, J., Lipscomb, W. H., da Silva, C. E., Vizcaino, M., Thayer-Calder, K., Lenaerts, J. T. M., and Sellevold, R.: Description and Demonstration of the Coupled Community Earth System Model v2 – Community Ice Sheet Model v2 (CESM2-CISM2), J. Adv. Model. Earth Sy., 13, e2020MS002356, https://doi.org/10.1029/2020MS002356, 2021.

Muskatel, H. B., Blahak, U., Khain, P., Levi, Y., and Fu, Q.: Parametrizations of Liquid and Ice Clouds’ Optical Properties in Operational Numerical Weather Prediction Models, Atmosphere, 12, 89, https://doi.org/10.3390/atmos12010089, 2021.

Nakayama, Y., Timmermann, R., and H. Hellmer, H.: Impact of West Antarctic ice shelf melting on Southern Ocean hydrography, Cryosphere, 14, 2205–2216, https://doi.org/10.5194/tc-14-2205-2020, 2020.

Nakayama, Y., Cai, C., and Seroussi, H.: Impact of Subglacial Freshwater Discharge on Pine Island Ice

Shelf, Geophys. Res. Lett., 48, e2021GL093923, https://doi.org/10.1029/2021GL093923, 2021.

Naughten, K. A., Meissner, K. J., Galton-Fenzi, B. K., England, M. H., Timmermann, R., Hellmer, H. H., Hattermann, T., and Debernard, J. B.: Intercomparison of Antarctic ice-shelf, ocean, and sea-ice interactions simulated by MetROMS-iceshelf and FESOM 1.4, Geosci. Model Dev., 11, 1257–1292, https://doi.org/10.5194/gmd-11-1257-2018, 2018.

Oleson, K. W., Lawrence, D. M., Bonan, G. B., Drewniak, B., Huang, M., Koven, C. D., Levis, S., Li, F., Riley, W. J., Subin, Z. M., Swenson, S. C., Thornton, P. E., Bozbiyik, A., Fisher, R., Heald, C. L., Kluzek, E., Lamarque, J.-F., Lawrence, P. J., Leung, L. R., Lipscomb, W., Muszala, S., Ricciuto, D. M., Sacks, W., Sun, Y., Tang, J., and Yang, Z.-L.: Technical Description of version 4.5 of the Community Land Model (CLM), Tech. Rep. July, NCAR, available at: http://www.cesm.ucar.edu/models/cesm1.2/clm/CLM45_Tech_Note.pdf (last access: 21 January 2022), 2013.

Paolo, F. S., Fricker, H. A., and Padman, L.: Volume loss from Antarctic ice shelves is accelerating, Science, 348, 327–331, https://doi.org/10.1126/science.aaa0940, 2015.

Pattyn, F.: Sea-level response to melting of Antarctic ice shelves on multi-centennial timescales with the fast Elementary Thermomechanical Ice Sheet model (f.ETISH v1.0), The Cryosphere, 11, 1851–1878, https://doi.org/10.5194/tc-11-1851-2017, 2017.

Pauling, A. G., Bitz, C. M., Smith, I. J., and Langhorne, P. J.: The Response of the Southern Ocean and Antarctic Sea Ice to Freshwater from Ice Shelves in an Earth System Model, J. Climate, 29, 1655–1672, https://doi.org/10.1175/JCLI-D-15-0501.1, 2016.

Pelle, T., Morlighem, M., Nakayama, Y., and Seroussi, H.: Widespread Grounding Line Retreat of Totten Glacier, East Antarctica, Over the 21st Century, Geophys. Res. Lett., 48, e2021GL093213, https://doi.org/10.1029/2021GL093213, 2021.

Pelletier, C.: Scripts and data for the PARASO Geoscientific Model Development paper figures, Zenodo [code], https://doi.org/10.5281/zenodo.5763819, 2021.

Pelletier, C. and Helsen, S.: PARASO ERA5 forcings, Zenodo [data set], https://doi.org/10.5281/zenodo.5590053, 2021.

Pelletier, C., Klein, F., Zipf, L., Haubner, K., Mathiot, P., Pattyn, F., Moravveji, E., and Vanden Broucke, S.: PARASO source code (no COSMO), Zenodo [code], https://doi.org/10.5281/zenodo.5576201, 2021a.

Pelletier, C., Klein, F., Zipf, L., Vanden Broucke, S., Haubner, K., and Helsen, S.: Input data for PARASO, a circum-Antarctic fully-coupled 5-component model, Zenodo [data set], https://doi.org/10.5281/zenodo.5588468, 2021b.

Pollard, D. and DeConto, R. M.: A simple inverse method for the distribution of basal sliding coefficients under ice sheets, applied to Antarctica, The Cryosphere, 6, 953–971, https://doi.org/10.5194/tc-6-953-2012, 2012.

Raphael, M. N., Marshall, G. J., Turner, J., Fogt, R. L., Schneider, D., Dixon, D. A., Hosking, J. S., Jones, J. M., and Hobbs, W. R.: The Amundsen Sea Low: Variability, Change, and Impact on Antarctic Climate, B. Am. Meteorol. Soc., 97, 111–121, https://doi.org/10.1175/BAMS-D-14-00018.1, 2016.

Reese, R., Albrecht, T., Mengel, M., Asay-Davis, X., and Winkelmann, R.: Antarctic sub-shelf melt rates via PICO, The Cryosphere, 12, 1969–1985, https://doi.org/10.5194/tc-12-1969-2018, 2018.
Richter, O., Gwyther, D. E., Galton-Fenzi, B. K., and Naughten, K. A.: The Whole Antarctic Ocean Model (WAOM v1.0): Development and Evaluation, Geosci. Model Dev. Discuss. [preprint], https://doi.org/10.5194/gmd-2020-164, in review, 2020.

Rignot, E., Jacobs, S., Mouginot, J., and Scheuchl, B.: Ice-Shelf Melting Around Antarctica, Science, 341, 266–270, https://doi.org/10.1126/science.1235798, 2013.

Rignot, E., Mouginot, J., Scheuchl, B., Broeke, M. V. D., Wessem, M. J. V., and Morlighem, M.: Four decades of Antarctic Ice Sheet mass balance from 1979–2017, P. Natl. Acad. Sci. USA, 116, 1095–1103, 2019.

Rintoul, S. R., Silvano, A., Pena-Molino, B., Wijk, E. V., Rosenberg, M., Greenbaum, J. S., and Blankenship, D. D.: Ocean heat drives rapid basal melt of the Totten Ice Shelf, Sci. Adv., 2, e1601610, https://doi.org/10.1126/sciadv.1601610, 2016.

Ritter, B. and Geleyn, J.: A comprehensive radiation scheme for numerical weather prediction models with potential applications in climate simulations, Mon. Weather Rev., 120.2, 303–325, https://doi.org/10.1175/1520-0493(1992)120<0303:ACRSFN>2.0.CO;2, 1992.

Roach, L. A., Dean, S. M., and Renwick, J. A.: Consistent biases in Antarctic sea ice concentration simulated by climate models, The Cryosphere, 12, 365–383, https://doi.org/10.5194/tc-12-365-2018, 2018.

Roach, L. A., Dörr, J., Holmes, C. R., Massonnet, F., Blockley, E. W., Notz, D., Rackow, T., Raphael, M. N., O’Farrell, S. P., Bailey, D. A., and Bitz, C. M.: Antarctic Sea Ice Area in CMIP6, Geophys. Res. Lett., 47, e2019GL086729, https://doi.org/10.1029/2019GL086729, 2020.

Rockel, B., Will, A., and Hense, A.: The regional climate model COSMO-CLM (CCLM), Meteorol. Z., 17, 347–348, https://doi.org/10.1127/0941-2948/2008/0309, 2008.

Roquet, F., Madec, G., McDougall, T., and Barker, P.: Accurate polynomial expressions for the density and specific volume of seawater using the TEOS-10 standard, Ocean Model., 90, 29–43, https://doi.org/10.1016/j.ocemod.2015.04.002, 2015.

Rosier, S. H. R., Hofstede, C., Brisbourne, A. M., Hattermann, T., Nicholls, K. W., Davis, P. E. D., Anker, P. G. D., Hillenbrand, C.-D., Smith, A. M., and Corr, H. F. J.: A New Bathymetry for the Southeastern Filchner-Ronne Ice Shelf: Implications for Modern Oceanographic Processes and Glacial History, J. Geophys. Res.-Oceans, 123, 4610–4623, https://doi.org/10.1029/2018JC013982, 2018.

Rousset, C., Vancoppenolle, M., Madec, G., Fichefet, T., Flavoni, S., Barthélemy, A., Benshila, R., Chanut, J., Levy, C., Masson, S., and Vivier, F.: The Louvain-La-Neuve sea ice model LIM3.6: global and regional capabilities, Geosci. Model Dev., 8, 2991–3005, https://doi.org/10.5194/gmd-8-2991-2015, 2015.

Schättler, U. and Blahak, U.: COSMO-Model Version 5.00: A Description of the Nonhydrostatic Regional COSMO-Model – Part V: Initial and Boundary Data for the COSMO-Model, Tech. rep., Deutscher Wetteringenst, https://doi.org/10.5676/DWD_PUB/NWV/COSMO-DOC_5.00_V, 2013.

Schlosser, E., Haumann, F. A., and Raphael, M. N.: Atmospheric influences on the anomalous 2016 Antarctic sea ice decay, The Cryosphere, 12, 1103–1119, https://doi.org/10.5194/tc-12-1103-2018, 2018.

Schneider, D. P. and Reusch, D. B.: Antarctic and Southern Ocean Surface Temperatures in CMIP5 Models in the Context of the Surface Energy Budget, J. Climate, 29, 1689–1716, https://doi.org/10.1175/JCLI-D-15-0429.1, 2016.

Schodlok, M. P., Menemenlis, D., Rignot, E., and Studinger, M.: Sensitivity of the ice-shelf/ocean system to the sub-ice-shelf cavity shape measured by NASA IceBridge in Pine Island Glacier, West Antarctica, Ann. Glaciol., 53, 156–162, https://doi.org/10.3189/2012AoG60A073, 2012.

Schoof, C.: Ice sheet grounding line dynamics: Steady states, stability, and hysteresis, J. Geophys. Res.-Earth Surf., 112, F03S28, https://doi.org/10.1029/2006JF000664, 2007.

Schulzweida, U.: CDO User Guide, Zenodo [code], https://doi.org/10.5281/zenodo.3539275, 2019.

Seroussi, H. and Morlighem, M.: Representation of basal melting at the grounding line in ice flow models, The Cryosphere, 12, 3085–3096, https://doi.org/10.5194/tc-12-3085-2018, 2018.

Seroussi, H., Nakayama, Y., Larour, E., Menemenlis, D., Morlighem, M., Rignot, E., and Khazendar, A.: Continued retreat of Thwaites Glacier, West Antarctica, controlled by bed topography and ocean circulation, Geophys. Res. Lett., 44, 6191–6199, https://doi.org/10.1002/2017GL072910, 2017.

Shapiro, N. M. and Ritzwoller, M. H.: Inferring surface heat flux distributions guided by a global seismic model: particular application to Antarctica, Earth Planet. Sc. Lett., 223, 213–224, https://doi.org/10.1016/j.epsl.2004.04.011, 2004.

Shepherd, A., Ivins, E., Rignot, E., Smith, B., van den Broeke, M., Velicogna, I., Whitehouse, P., Briggs, K., Joughin, I., Krinner, G., Nowicki, S., Payne, T., Scambos, T., Schlegel, N., A. G., Agosta, C., Ahlström, A., Babonis, G., Barletta, V., Blazquez, A., Bonin, J., Csatho, B., Cullather, R., Felikson, D., Fettweis, X., Forsberg, R., Gallee, H., Gardner, A., Gilbert, L., Groh, A., Gunter, B., Hanna, E., Harig, C., Helm, V., Horvath, A., Horwath, M., Khan, S., Kjeldsen, K. K., Konrad, H., Langen, P., Lecavalier, B., Loomis, B., Luthcke, S., McMillan, M., Melini, D., Mernild, S., Mohajerani, Y., Moore, P., Mouginot, J., Moyano, G., Muir, A., Nagler, T., Nield, G., Nilsson, J., Noel, B., Otosaka, I., Pattle, M. E., Peltier, W. R., Pie, N., Rietbroek, R., Rott, H., Sandberg-Sørensen, L., Sasgen, I., Save, H., Scheuchl, B., Schrama, E., Schröder, L., Seo, K.-W., Simonsen, S., Slater, T., Spada, G., Sutterley, T., Talpe, M., Tarasov, L., van de Berg, W. J., van der Wal, W., van Wessem, M., Vishwakarma, B. D., Wiese, D., Wouters, B., and The IMBIE team: Mass balance of the Antarctic Ice Sheet from 1992 to 2017, Nature, 558, 219–222, https://doi.org/10.1038/s41586-018-0179-y, 2018.

Smith, R. S., Mathiot, P., Siahaan, A., Lee, V., Cornford, S. L., Gregory, J. M., Payne, A. J., Jenkins, A., Holland, P. R., Riedley, J. K., and Jones, C. G.: Coupling the U.K. Earth System Model to Dynamic Models of the Greenland and Antarctic Ice Sheets, J. Adv. Model. Earth Sy., 13, e2021MS002520, https://doi.org/10.1029/2021MS002520, 2021.

Souverijns, N., Gossart, A., Gorodetskaya, I. V., Lhermitte, S., Mangold, A., Laffineur, Q., Delcloo, A., and van Lipzig, N. P. M.: How does the ice sheet surface mass balance relate to snowfall? Insights from a ground-based precipitation radar in East Antarctica, The Cryosphere, 12, 1987–2003, https://doi.org/10.5194/tc-12-1987-2018, 2018.

Souverijns, N., Gossart, A., Demuzere, M., Lenaerts, J. T. M., Medley, B., Gorodetskaya, I. V., Vanden Broucke, S., and van
Lipzig, N. P. M.: A New Regional Climate Model for POLAR-CORDEX: Evaluation of a 30-Year Hindcast with COSMO-CLM2 Over Antarctica, J. Geophys. Res.-Atmos., 124, 1405–1427, https://doi.org/10.1029/2018JD028862, 2019.

Stein, C. A. and Stein, S.: A model for the global variation in oceanic depth and heat flow with lithospheric age, Nature, 359, 123–129, 1992.

Storkey, D., Blaker, A. T., Mathiot, P., Megann, A., Aksenov, Y., Blockley, E. W., Calvert, D., Graham, T., Hewitt, H. T., Hyder, P., Kuhlbrodt, T., Rae, J. G. L., and Sinha, B.: UK Global Ocean GO6 and GO7: a traceable hierarchy of model resolutions, Geosci. Model Dev., 11, 3187–3213, https://doi.org/10.5194/gmd-11-3187-2018, 2018.

Sun, S., Pattyn, F., Simon, E. G., Albrecht, T., Cornford, S., Calov, R., Dumas, C., Gillet-Chaulet, F., Goelzer, H., Golledge, N. R, Greve, R., Hoffman, M. J., Humbert, A., Kazmierczak, E., Kleiner, T., Leguy, G. R., Lipscomb, W. H., Martin, D., Morlighem, M., Nowicki, S., Pollard, D., Price, S., Quiquet, A., Seroussi, H., Schlemm, T., Sutter, J., van de Wal, R. S. W., Winkelmann, R., and Zhang, T.: Antarctic ice sheet response to sudden and sustained ice-shelf collapse (ABUMIP), J. Glaciol., 66, 891–904, https://doi.org/10.1017/jog.2020.67, 2020.

Swingedouw, D., Fichefet, T., Huybrechts, P., Goosse, H., Driesschaert, E., and Loutre, M.-F.: Antarctic ice-sheet melting provides negative feedbacks on future climate warming, Geophys. Res. Lett., 35, L17705, https://doi.org/10.1029/2008GL034410, 2008.

Thompson, A. F., Stewart, A. L., Spence, P., and Heywood, K. J.: The Antarctic Slope Current in a Changing Climate, Rev. Geophys., 56, 741–770, https://doi.org/10.1029/2018RG000624, 2018.

Timmermann, R. and Goeller, S.: Response to Filchner–Ronne Ice Shelf cavity warming in a coupled ocean–ice sheet model – Part 1: The ocean perspective, Ocean Sci., 13, 765–776, https://doi.org/10.5194/os-13-765-2017, 2017.

Torres, O., Braconnot, P., Hourdin, F., Roehrig, R., Marti, O., Belamari, S., and Lefebvre, M.-P.: Competition Between Atmospheric and Surface Parameterizations for the Control of Air-Sea Latent Heat Fluxes in Two Single-Column Models, Geophys. Res. Lett., 46, 7780–7789, https://doi.org/10.1029/2019GL082720, 2019a.

Torres, O., Braconnot, P., Marti, O., and Gential, L.: Impact of air-sea drag coefficient for latent heat flux on large scale climate in coupled and atmosphere stand-alone simulations, Clim. Dynam., 52, 2125–2144, https://doi.org/10.1007/s00382-018-4236-x, 2019b.

Tsamados, M., Feltham, D. L., Schroeder, D., Flocco, D., Farrell, S. L., Kurtz, N., Laxon, S. W., and Bacon, S.: Impact of Variable Atmospheric and Oceanic Form Drag on Simulations of Arctic Sea Ice, J. Phys. Oceanogr., 44, 1329–1353, https://doi.org/10.1175/JPO-D-13-0215.1, 2014.

Turner, J., Bracegirdle, T. J., Phillips, T., Marshall, G. J., and Hosking, J. S.: An Initial Assessment of Antarctic Sea Ice Extent in the CMIP5 Models, J. Climate, 26, 1473–1484, https://doi.org/10.1175/JCLI-D-12-00068.1, 2013.

Turner, J., Orr, A., Gudmundsson, G. H., Jenkins, A., Bingham, R. G., Hillenbrand, C.-D., and Bracegirdle, T. J.: Atmosphere-ocean-ice interactions in the Amundsen Sea Em-

bayment, West Antarctica, Rev. Geophys., 55, 235–276, https://doi.org/10.1002/2016RG000532, 2017.

Ungermann, M., Tremblay, L. B., Martin, T., and Losch, M.: Impact of the ice strength formulation on the performance of a sea ice thickness distribution model in the Arctic, J. Geophys. Res.-Oceans, 122, 2090–2107, https://doi.org/10.1002/2016JC012128, 2017.

Valcke, S.: The OASIS3 coupler: a European climate modelling community software, Geosci. Model Dev., 6, 373–388, https://doi.org/10.5194/gmd-6-373-2013, 2013.

Valcke, S., Craig, T., and Coquart, L.: OASIS3-MCT user guide, Tech. rep., CERFACS, available at: https://oasis.cerfacs.fr/wp-content/uploads/sites/114/2021/02/GLOBC-TR-oasis3mct_UserGuide3.0_052015.pdf (last access: 21 January 2022), 2013.

Van Achter, G., Fichefet, T., Goosse, H., Pelletier, C., Sterlin, J., Huot, P.-V., Lemieux, J.-F., Fraser, A. D., Haubner, K., and Porter-Smith, R.: Modelling landfast sea ice and its influence on ocean–ice interactions in the area of the Totten Glacier, East Antarctica, Ocean Model., 169, 101920, https://doi.org/10.1016/j.ocemod.2021.101920, 2022.

Vancoppenolle, M., Fichefet, T., Goosse, H., Bouillon, S., Madec, G., and Maqueda, M. A. M.: Simulating the mass balance and salinity of Arctic and Antarctic sea ice. 1. Model description and validation, Ocean Model., 27, 33–53, https://doi.org/10.1016/j.ocemod.2008.10.005, 2009.

Vancoppenolle, M., Bouillon, S., Fichefet, T., Goosse, H., Lecomte, O., Morales Maqueda, M. A., and Madec, G.: LIM The Louvain-la-Neuve sea Ice Model, Tech. Rep. 31, Note du Pôle de Modélisation de l'Institut Pierre-Simon Laplace No. 31, ISSN No 1288-1619, available at: https://cmc.ipsl.fr/images/publications/scientific_notes/lim3_book.pdf (last access 21 January 2022), 2012.

van den Broeke, M., Bamber, J., Ettema, J., Rignot, E., Schrama, E., van de Berg, W. J., van Meijgaard, E., Velicogna, I., and Wouters, B.: Partitioning Recent Greenland Mass Loss, Science, 326, 984–986, https://doi.org/10.1126/science.1178176, 2009.

van Kampenhout, L., Lenaerts, J. T. M., Lipscomb, W. H., Sacks, W. J., Lawrence, D. M., Slater, A. G., and van den Broeke, M. R.: Improving the Representation of Polar Snow and Firn in the Community Earth System Model, J. Adv. Model. Earth Sy., 9, 2583–2600, https://doi.org/10.1002/2017ms000988, 2017.

van Lipzig, N. P. M. and van den Broeke, M. R.: A model study on the relation between atmospheric boundary-layer dynamics and poleward atmospheric moisture transport in Antarctica, Tellus A, 54, 497–511, https://doi.org/10.1034/j.1600-0870.2002.201404.x, 2002.

Van Pham, T., Brauch, J., Dieterich, C., Frueh, B., and Ahrens, B.: New coupled atmosphere-ocean-ice system COSMO-CLM/NEMO: assessing air temperature sensitivity over the North and Baltic Seas, Oceanologia, 56, 167–189, https://doi.org/10.5697/oc.56-2.167, 2014.

Vertenstein, M., Bertini, A., Craig, T., Edwards, J., Levy, M., Mai, A., and Schollenberger, J.: CESM user's guide (CESM1.2 release series user's guide), Tech. rep., National Center for Atmospheric Research, available at: https://www.cesm.ucar.edu/models/cesm1.2/cesm/doc/usersguide/book1.html (last access: July 2021), 2013.
Vignon, E., Hourdin, F., Genthon, C., de Wiel, B. J. H. V., Gallée, H., Madeleine, J.-B., and Beaumet, J.: Modeling the Dynamics of the Atmospheric Boundary Layer Over the Antarctic Plateau With a General Circulation Model, J. Adv. Model. Earth Sy., 10, 98–125, https://doi.org/10.1002/2017ms001184, 2018.

Vionnet, V., Brun, E., Morin, S., Boone, A., Faroux, S., Le Moigne, P., Martin, E., and Willemet, J.-M.: The detailed snow-pack scheme Crocus and its implementation in SURFEX v7.2, Geosci. Model Dev., 5, 773–791, https://doi.org/10.5194/gmd-5-773-2012, 2012.

Voldoire, A., Decharme, B., Pianezze, J., Lebeaupin Brossier, C., Sevault, F., Seyfried, L., Garnier, V., Bielli, S., Valcke, S., Alias, A., Accensi, M., Ardhuin, F., Bouin, M.-N., Ducrocq, V., Faroux, S., Giordani, H., Léger, F., Marsaleix, P., Rainaud, R., Redelsperger, J.-L., Richard, E., and Riette, S.: SURFEX v8.0 interface with OASIS3-MCT to couple atmosphere with hydrology, ocean, waves and sea-ice models, from coastal to global scales, Geosci. Model Dev., 10, 4207–4227, https://doi.org/10.5194/gmd-10-4207-2017, 2017.

Wang, C., Zhang, L., Lee, S.-K., Wu, L., and Mechoso, C. R.: A global perspective on CMIP5 climate model biases, Nat. Clim. Change, 4, 201–205, https://doi.org/10.1038/nclimate2118, 2014.

Wei, W., Blankenship, D. D., Greenbaum, J. S., Gourmelen, N., Dow, C. F., Richter, T. G., Greene, C. A., Young, D. A., Lee, S., Kim, T.-W., Lee, W. S., and Assmann, K. M.: Getz Ice Shelf melt enhanced by freshwater discharge from beneath the West Antarctic Ice Sheet, The Cryosphere, 14, 1399–1408, https://doi.org/10.5194/tc-14-1399-2020, 2020.

Will, A., Akhtar, N., Brauch, J., Breil, M., Davin, E., Hohagemann, H. T. M., Maisonnave, E., Thürkow, M., and Weiher, S.: The COSMO-CLM 4.8 regional climate model coupled to regional ocean, land surface and global earth system models using OASIS3-MCT: description and performance, Geosci. Model Dev., 10, 1549–1586, https://doi.org/10.5194/gmd-10-1549-2017, 2017.

Wille, J. D., Favier, V., Gorodetskaya, I. V., Agosta, C., Kittel, C., Beeman, J. C., Jourdain, N. C., Lenaerts, J. T. M., and Cordron, F.: Antarctic Atmospheric River Climatology and Precipitation Impacts, J. Geophys. Res.-Atmos., 126, e2020JD033788, https://doi.org/10.1029/2020JD033788, 2021.

Zender, C. S.: Analysis of Self-describing Gridded Geoscience Data with netCDF Operators (NCO), Environ. Modell. Softw., 23, 1338–1342, https://doi.org/10.1016/j.envsoft.2008.03.004, 2008.

Zentek, R. and Heinemann, G.: Verification of the regional atmospheric model CCLM v5.0 with conventional data and lidar measurements in Antarctica, Geosci. Model Dev., 13, 1809–1825, https://doi.org/10.5194/gmd-13-1809-2020, 2020.

Zhang, T., Price, S., Ju, L., Leng, W., Brondex, J., Durand, G., and Gagliardini, O.: A comparison of two Stokes ice sheet models applied to the Marine Ice Sheet Model Intercomparison Project for plan view models (MISMIP3d), The Cryosphere, 11, 179–190, https://doi.org/10.5194/tc-11-179-2017, 2017.

Zhou, Q. and Hattermann, T.: Modeling ice shelf cavities in the unstructured-grid, Finite Volume Community Ocean Model: Implementation and effects of resolving small-scale topography, Ocean Model., 146, 101536, https://doi.org/10.1016/j.ocemod.2019.101536, 2020.

Zuo, H., Balmaseda, M. A., Tietsche, S., Mogensen, K., and Mayer, M.: The ECMWF operational ensemble reanalysis–analysis system for ocean and sea ice: a description of the system and assessment, Ocean Sci., 15, 779–808, https://doi.org/10.5194/os-15-779-2019, 2019.