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

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 f.ETISh1.7 (ice sheet), NEMO3.6 (ocean), LIM3.6 (sea ice), COSMO5.0 (atmosphere) and CLM4.5 (land), which are here run at an 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 exchanges. 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 two-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.



Introduction
The Antarctic climate is characterized by large natural fluctuations and complex interactions between the ice sheet, ocean, sea 15 ice and atmosphere. One of the consequences of these interactions observed over the last decades lies in the Antarctic ice sheet (AIS) losing mass (Rignot et al., 2019;Edwards 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 20 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 floating 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;Holland et al., 2019). For example, oceanic processes such 25 as the intrusion of relatively warm water over the Antarctic continental 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 freshwater injection at depth, thus impacting the circulation and heat transport in the Southern Ocean (Hellmer, 2004;Jourdain 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, 30 at the circumpolar scale, no definitive conclusion 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 35 atmosphere at the high latitudes of the Southern Hemisphere. Despite their varying typical response timescales, the respective evolutions of polar Earth subcomponents are entangled and interdependent (Turner et al., 2017;Fyke et al., 2018). Coupled models are the designated tool for quantitatively investigating them, since they can incorporate feedbacks between the components. As the success of coupled model intercomparison projects (e.g. CMIP6, Eyring et al., 2016) testifies, coupled oceanatmosphere models have been a staple in climate modeling for a few decades. All state-of-the-art ocean models covering polar 40 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 development of idealized test cases for simulating ice-shelf melting (Asay-Davis et al., 2016;Zhang et al., 2017;Jordan 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;Jourdain et al., 45 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 circumpolar domains (Naughten et al., 2018;Nakayama et al., 2020;Comeau et al., 2021).

Representation of the ice-shelf cavities within the ocean model
In PARASO, Antarctic ice-shelf cavities are opened to the ocean circulation  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 immerged 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 115 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): where x ∈ {T, S} distinguishes heat and salt, Γ x is the default exchange coefficient, C isf D the ice -ocean drag coefficient, and u isf h the near-ice ocean current velocity. The ice-shelf melt parameterization also includes TEOS-10-compatible coefficients 120 for evaluating seawater freezing point temperature . In case 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 behaviour occurs in case of refreezing (heat injection, volume extraction). All constants used in the NEMO ice-shelf module are listed in Table 1. Table 1. NEMO ice-shelf-related physical constants. Here, "ice" refers to the ice sheet (not the sea ice). ΓT and ΓS are taken from Jourdain et al. (2017); ρ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 tradeoff between physical realism and numerical stability. The LIM3.6 sea-ice dynamics rely on an elastic-viscous-plastic rheology formulated on a C-grid (Bouillon et al., , 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 intercategory boundaries located around 0.5 m, 1.15 m, 2 m and 3.8 m. The choice of five categories comes from a tradeoff between computational constraints and physical realism in sea-ice mean state and variability 130 (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 135 et al. (2015) and Barthélemy et al. (2016), respectively.
2.3 Atmosphere and land model: COSMO-CLM 2 5.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 140 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 145 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), sub-gridscale turbulence (Raschendorfer, 2005) and cloud microphysics .
COSMO-CLM is coupled to the Community Land Model version 4.5 (Oleson et al., 2013) using the OASIS 3-MCT coupler . This coupling combination, which is hereinafter named CCLM 2 , forms the land-atmosphere basis for 150 PARASO, is described in Davin et al. (2011). More recently, at least two CCLM 2 Antarctic configurations have been developed (e.g., Souverijns et al., 2019;Zentek and Heinemann, 2020). In the following, we will use AEROCLOUD, 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 microphysic scheme for precipitation. The PARASO CCLM 2 version has been adapted from AEROCLOUD, and their dif-155 ferences are detailed in App. C. Moreover, a list of key PARASO-specific CCLM 2 parameterization and tuning choices is provided in Table A2.  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 being represented, the NEMO -COSMO exchanged fields are not detailed, as they have not been changed compared to existing settings.

Coupling interfaces
In this section, we describe the three novel coupling interfaces developed for PARASO. Hereinafter, we call: -Coupling window: a period between two coupling instances (intermodel exchange of information).

160
-Restart leg: 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 (timewise) CPU requirements and queuing time on supercomputers.
-Online coupling: 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).

165
-Offline coupling: 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 2 run jointly with regular exchanges of information (online coupling, see Sect. 3.1). In between NEMO -CCLM 2 restart legs, f.ETISh runs by itself with input data files coming from NEMO and CCLM 2 outputs, 170 and provides NEMO with updated ice-shelf geometry (offline coupling, see Sects. 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.
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-175 ice model. We refer to "COSMO" (atmosphere model) instead of "CCLM 2 " (atmosphere -land model), since this coupling interface does not involve land processes. COSMO and NEMO are coupled online through OASIS3-MCT_2.0 Valcke, 2013). While this coupling interface has been built from preexisting configurations involving these models (Van Pham et al., 2014;, our final setting includes several novel model developments which are described here. Figure 2 schematically represents this coupling interface.

180
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 bicubically 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 185 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; else, 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 190 same reanalysis is used for forcing the COSMO lateral boundary condition and the extra-COSMO domain NEMO surface.

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 (in PARASO's case, 3 hr). On 195 the other hand, the surface values seen by COSMO are 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.

200
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 one 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 205 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 twofold 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 210 implementation satisfies local air-sea energy conservation. Technical details on the pseudo-averaging operations are provided in App. D. 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 215 from COSMO to NEMO. As described in Rousset et al. (2015), surface fluxes injected into the sea ice are redistributed 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 nonsolar 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 charac-220 terizations (and thus drag coefficients) differ. All heat fluxes are tile-specific and thus twofold. 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 nonsolar 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 225 transfer coefficients are evaluated from a turbulent kinetic energy (TKE) derived scheme (Raschendorfer, 2005). 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 runoff to the open ocean. Over sea ice, solid precipitation contributes to snow accumulation; over the open 230 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 nonsolar 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 235 systematically computed and sent over all ocean cells using the last known (resp., initial) sea-ice surface properties in case no Rainfall rate kg m −2 s −1 sea ice is currently (resp., has ever been) present in the category. This is carried out in case sea ice appears, through advection or thermodynamics, over one grid cell within a coupling window. Coupling window n Coupling window n + 1 Figure 2. COSMO -NEMO coupling interface diagram for two coupling windows. Full blue (resp., black) arrows correspond to exchanges of tile-specific (resp., non-tiled) fields. The thick dashed vertical line delimits the transition from one coupling window to the next one.

Exchanged fields
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 (resp., outgoing) fluxes from coupling window n − 1 (resp., to coupling window n + 2) are not represented.  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 exchanges between NEMO and f.ETISh. The NEMO spinup run is used as initial conditions for the first NEMO leg only.

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): a multiple of the NEMO restart length; a divisor of the total experiment length.
NEMO always provides monthly sub-shelf 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 three months. This was chosen so that our relatively short (two-year) 260 experiments can include a few ocean -ice sheet coupling occurrences.

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).

265
-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; water column thickness.
These fields are then interpolated with a conservative method. The resulting columns are treated differently according to 270 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). After each ocean -ice-sheet coupling episode, NEMO uses the post-processed updated f.ETISh geometry as its new sub-shelf 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 290 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 smoothening effect over potentially large new openings. The initial current velocities of new cells are set to zero. On all modified columns (thinned or thickened), an 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 295 brought into (resp., taken from) the ocean upon cell opening (resp. closing). For example, when ice-sheet geometry changes lead to a NEMO cell opening (resp. closing), the corresponding volume and internal energy are added into (resp. 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 300 -ice-sheet coupling yields phantom current velocities into (resp., out of) the grounded ice-sheet upon cell opening (resp., 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 305 update. However, we have observed that enforcing the NEMO cavity geometry constraints described in Table 1 keeps critical numerical instabilities from arising, even with yearly coupling windows (not shown here).

CCLM 2 -f.ETISh interface
The exchanges between CCLM 2 and f.ETISh are one-way (from CCLM 2 to f.ETISh) and performed offline. CCLM 2 runs for one coupling time window, sending f.ETISh monthly time series of surface mass balance (SMB), which here boils down 310 to solid precipitation minus surface sublimation. The SMB is converted from [kg m −1 month −1 ] to ice equivalent thickness changes [ m a ] by using the reference ice density of ρ ice = 917kg m −3 . Interpolation between the CCLM 2 and fETISh grids are performed bilinearly. For the sake of simplicity, the coupling window partitioning for the CCLM 2 -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 2 is included within the surface boundary condition of f.ETISh. However, variations in ice-315 sheet surface elevation are not sent back to CCLM 2 , which is a limitation of PARASO, and the reason why this "coupling" procedure is referred to as "one-way". It should however be stressed out that at the time scales PARASO has been developed for (decadal at the longest), variations in Antarctic surface topography are negligible. In this subsection we detail each subcomponent of the PARASO configuration. Table 3 (Morlighem et al., 2020). This dataset, originally at a 500 m spatial resolution, has been resampled to the fETISh 8 km grid using a low-pass filter. Data was subsequently corrected to ensure grounded ice being actually grounded and floating ice being 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, 330 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 time scales: 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 335 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 (≈ 10 m draft, which is considerably smaller than typical draft values).

NEMO configuration
The NEMO PARASO set-up is derived from the GO7 configuration described in Storkey et al. (2018). The ocean grid is ePERIANT025, which includes the southernmost ice-shelf cavities . With a single lateral ocean boundary 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  with three geometrical constraints balancing physical realism and numerical stability (see Table 1): Parameters enforcing these constraints are also listed in Table 1. At the depths of the ice-shelf cavities, the vertical resolution roughly ranges from 10 m to 150 m (but we remind that partial cells are used).

CCLM 2 configuration
In PARASO, the COSMO rotated latitude-longitude grid has a horizontal resolution of 0.22 The only f.ETISh forcings are sub-shelf melt rates and SMB, which are provided by NEMO and CCLM 2 (respectively) and updated at every coupling window.

NEMO forcings
Over the ocean surface located outside of the CCLM 2 spatial coverage, the NEMO surface is forced with fluxes computed by the CORE bulk formula (Large and Yeager, 2004), using the ERA5 reanalysis as atmosphere input (Hersbach et al., 2020). The  (Stein and Stein, 1992). River runoff is adapted from Dai and Trenberth (2002) and prescribed as a climatological data set (Bourdallé-Badie and Treguier, 2006). While icebergs are not dynamically simulated, their freshwater melt input is injected as additional runoff from an interannual dataset obtained from a previous iceberg-including NEMO simulation covering 1979 to 2015 (Marsh et al., 2015;Merino et al., 2016; Jourdain the river deltas and on a stripe between 20 km 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 (N. Jourdain, personal communication, 2019). No salinity restoring is applied.

CCLM 2 forcings 385
At its lateral boundaries, COSMO is relaxed to ERA5 using a one-way interactive nesting based on Davies (1976). This consists in defining a relaxation zone where the internal COSMO-CLM solution is nudged against the external ERA5 solution. Within this zone, the variables of the driving 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 boundary relaxation inside the boundary zone has a tangent hyperbolic form (Kållberg, 1977). The width 390 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 documentation (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 interpolation program INT2LM (Schättler and Blahak, 2013 for COSMO are wind speeds (zonal and meridional), air temperature, pressure deviation from a reference pressure (1000 hPa at sea level), specific water vapour 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 2 atmosphere dynamics. A sponge layer with Rayleigh damping in the upper levels of the model domain is activated in order to avoid artificial reflections of gravity waves.
A cosine damping profile with maximum damping at the top and zero damping at the base of the sponge layer (11000 m 400 elevation) is applied.

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

405
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 coefficients (Pattyn, 2017). The method has been further optimized to enhance convergence (Bernales et al., 2017). In addition to Pollard and DeConto (2012), we also implemented a regularization term that uses a 2D Gaussian smoothing kernel to filter out high-frequency noise in the basal sliding coefficients. For the initialization, 410 f.ETISh is run with reduced complexity, using the shallow-ice approximation (SIA) for the evolution of the grounded ice sheet, thereby keeping ice shelves in the observed geometry. 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 thermomechanical 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, surface mass balance, and geothermal heat flux. The latter is based on the seismic-based dataset due to 415 Shapiro and Ritzwoller (2004). For the relaxation, ice shelves are allowed to evolve and the hybrid shallow-shelf/shallow-ice approximation (HySSA; Bueler and Brown, 2009) is used to obtain the velocity field. The applied sub-shelf melt rates are taken as the constant (timewise) 1984 -1993 mean from a NEMO stand-alone run. The PARASO setting of f.ETISh (e.g. constant calving front, comp. 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 tradeoff between gaining quasi-steady-state 420 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 2 stand-alone experiment .

NEMO initialization
The ocean is initialized from a one-year spinup run performed with NEMO stand-alone, using forcings corresponding to year 425 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 3D 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 430 to start from an upper ocean state roughly matching the CCLM 2 initialization, which is also based on ERA5, up to the limits permitted by the differences in air-sea bulk formulae.

CCLM 2 initialization
The land and ice masks used by CCLM 2 to determine the fraction of land vs. ocean and land ice vs. ice-free areas are calculated on the CCLM 2 grid based on the geometry provided by f.ETISh, so that the geometry of the fully coupled setup is consistent 435 from one subcomponent to another.
COSMO 3D 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, 440 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 −3 . The snow liquid water and ice contents for layer i = 1 to 5 are set to w liq,i = 0 kg m −2 and w ice,i = ∆z i ρ sno kg m −2 , respectively. Although a 1 m SWE snow pack is too small to fully represent the firn layer, which can extend to more than one hundred meters 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 2 as long as the adaptations 445 from van Kampenhout et al. (2017) are included to the CLM snow module (except for the snow pack depth). For future PARASO-related work, keeping a maximum snow depth of 1 m is also convenient to limit the spin-up 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 comes 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 Chapter 1 of Vertenstein et al.,450 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 assumes all ice shelves to be exactly at sea level.  to PARASO ocean (resp. AEROCLOUD ocean to PARASO land), are displayed with purple (resp. red) shadings.

Results
In this section, diagnostics from a 2-year (2000 -2001) Table   470 4. Name a : same experimental design as PARASO except for the surface flux tiling which has not been implemented in CCLM 2 stand-alone.
The PARASO biases are more generally discussed and put into perspective in Sect. 6. Computational aspects related to PARASO are briefly covered in App. B. PARASO-specific tunings used for obtaining a realistic sea-ice seasonal cycle are described in App. E.

Ice sheet and ice shelves
475 Figure 6 displays the sub-shelf melt rates for PAROCE and PARASO in different regions. Both experiments' regional melt rate patterns agree relatively well with each other (Figs. 6(b)-(e)), but the total melt of PARASO is significantly increased in the Weddell (+31%), Indian Ocean (+85%) and West Pacific (+162%) sectors (see Figs. 6(a)-(c) and 6(f)) due to the warmer and more saline Eastern Antarctic surface water in PARASO (see Figs. 6(g), 9(k,m) and 10(k,m)). However, the West-Antarctic melt rates remain of similar magnitudes from PAROCE to PARASO (−9% in the Amundsen sea and −16% in the Bellingshausen 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 490 (Schodlok et al., 2012;Rosier et al., 2018;Goldberg et al., 2019;Brisbourne et al., 2020;Goldberg et al., 2020;Wei et al., 2020).  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 (Figs. 7 (b)-(c)).

495
That thickening also occurs for PARCRYO SMB (Fig. 7 (d)) and agrees well with the differences between the forcing SMB data sets ( Fig. 7 (g)). PARACRYO SMB and PARACRYO apply the SMB from PARATMO, while PARACRYO CTRL 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. 7(f)) reveals a slightly larger thickening of the grounded ice in coastal areas for PARASO.
This is in accordance with the increased SMB for those regions provided by CCLM 2 in PARASO ( Fig. 7(i)). Figure 7(j) also 500 emphasizes the relationship between the SMB forcing and the thickness evolution of the grounded ice. While the SMB is quite constant for the first two months, the change in ice volume (VAF) can be attributed to the dynamical response of the modelled 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 SMB , which are forced by the identical SMB data set (from PARATMO), show a very similar evolution of the grounded ice ( Fig. 7(e)), while the change The forcing with sub-shelf melt, which has only been applied for PARCRYO (sub-shelf melt from PAROCE) and PARASO, has a very limited effect on the grounded ice on the short evaluation period of one year, but clearly affects the ice shelves ( Fig. 7(e)). The ice shelves show a moderate thinning of up to a few meters almost everywhere, except in the Bellingshausen  ice sheet (Fig. 7(f)) and is in accordance with the increased sub-shelf melt of PARASO, which also mainly affects the East Antarctic ice shelves (Figs. 6(a) and 7(i)).
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 2 , 525 sub-shelf 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 two years, which is a short period for the relatively slow (e.g., in comparison with the ocean and atmosphere) ice-sheet system.

Ocean and sea ice
In this section, we compare PARASO results to observations as well as to PAROCE, a NEMO stand-alone experiment forced 530 by ERA5, with the same experimental design (except the coupling, see Table 4). diagnostics from PAROCE and PARASO. Generally speaking, the PARASO Antarctic sea-ice extent seasonal cycle shares similar features compared to PAROCE (see Fig. 8(a)), 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 535 extent at ≈ 10 6 km 2 instead of 3 × 10 6 km 2 ), which is a persistent bias for most coupled models (Turner et al., 2013;Mahlstein et al., 2013;Roach et al., 2018Roach et al., , 2020. PAROCE and PARASO also suffer 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 × 10 6 km 2 instead of 18 × 10 6 km 2 ). In 540 contrast with the sea-ice extent, the PARASO maximum sea-ice volume is generally larger than PAROCE's (see Fig. 8(b)). This may be linked to the relatively stronger PARASO winds and their coastwards orientations, making sea ice more compact but less extensive (see Fig. A7). Regarding the sea-ice growth cycle (Figs. 8(c)-(e)), the PAROCE and PARASO biases are similar, with the PARASO growth being delayed by about 15 days 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 (Figs. 8(f)-(h)) occurs much earlier than 545 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 one 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 550 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). tialization, 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   December (dashed), over six circumpolar basins (same as Fig. 9). Note the nonlinear z-axis (following the NEMO vertical discretization).
For WOA18, the data is taken from a 1995 -2004 climatology; for PAROCE and PARASO, from simulation outputs for the year 2001 (second simulated year).
injection at depth, which is a novel feature for both configurations. While these general biases are present, the vertical structure 565 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 ± 4 Sv (1 Sv = 10 6 m 3 s −1 ) in PAROCE and 43 ± 4 Sv in PARASO (see also Fig. A6). This is severely underestimated compared with observations (e.g., shelf (29 Sv) counterbalancing the eastwards transport, thus leading to a degraded net transport. While both PAROCE and PARASO Drake transport 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.

Atmosphere
575 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. 11(a)) compared to the winter season ( Fig. 11(c)), 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 580 PARASO and ERA5 over the ice shelves, but this effect might contribute to the temperature signal visible over the Berkner and Roosevelt islands. 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 . Moreover, differences of comparable magnitude, 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 AEROCLOUD, a pre-existing CCLM 2 stand-alone experiment that was thoroughly evaluated using in-585 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 -AEROCLOUD 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 AEROCLOUD are probably related to deficiencies in the representation of turbulent fluxes in the boundary layer, which is known to be a considerable challenge 590 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,

595
PARASO was also compared to PARATMO (see Table 4), a CCLM 2 stand-alone run with the same tuning paramaters as for PARASO, except for the surface flux tiling which has not been implemented for CCLM 2 stand-alone experiments (instead, sea ice is assumed to fully cover a grid cell where the surface temperature is lower than −1.8 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 605 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 610 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 615 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 2 ) are compared, and the observed intermodel spread suggests considerable uncertainties even in uncoupled, atmosphere-only simulations. is primarily attributed to less precipitation in PARASO over this region (not shown). We also analysed 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 630 partly related to the representation of atmospheric dynamics, and not solely to the coupling.

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 require-635 ments 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 circumpolar 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 2 coupled configuration (i.e., excluding coupling with f.ETISh). As briefly mentioned in Sect.
5.1, even at longer time scales (40 years), the ice-sheet initialization appears 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 645 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).

650
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, 655 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 2 geometry instead of the standard data set made the PARASO setup more consistent intermodel-wise. How-660 ever, CCLM 2 parameters had been tuned to this standard topographical data set, 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 App. 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 665 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 670 to either NEMO or COSMO stand-alone. 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 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 untiled 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 680 may be of use for future studies intending at investigating the sole impact of nonbinary 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.

685
Second, coupled PARASO experiments suffer from excessive latent heat (and thus evaporation) over sea ice, as illustrated in Fig. E1(e). 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. E1(c)). 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, 690 even more snow melts at the surface. While some specific tuning (see App. 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 data sets it relied on (SMB from COSMO, ice-shelf 695 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 yr NEMO -f.ETISh coupled experiments (without CCLM 2 , see Fig. A9) have suggested that at least on these timescales, 700 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 time scales -and therefore the simulated cavity changes -are smaller than the ice geometry changes during the multi-millenial 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 705 tool. As previously mentioned, the objective of this paper is to document our tool, assess its performance as-is, and share it with the community. Taking advantage of the numerical stability of PARASO and the expertise 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 noting 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 described in Sect. 5 remains circumscribed to relatively small amplitudes, hence that PARASO does not keep on drifting at the longer term, and that its performance remains reasonably satisfactory at least in terms of Antarctic sea ice.
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: August 24th 2021). To download, the user needs to become a member of the CLM- Appendix A: Additional figures and tables Table A1. List of key NEMO parameters (excluding those related to ice-shelf cavities, which are already covered in Table 1).

Description Value
Open ocean  Table A2. List of key CCLM 2 parameterizations.

Description Parameter Value
Turbulence: 2.5 th order Mellor and Yamada (1982) TKE closure with surface thermal heterogeneity (Raschendorfer, 2005). 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 time step is 1, NEMO's is 2, the NEMO-CCLM coupling frequency is 4 and the ice sheet model coupling frequency is 8. Inbetween legs, the NEMO-CCLM system is stopped. Black boxes represent data written on the disk. Green boxes represent OASIS operations. Full arrows represent exchanges 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.
Perform conserv interp(is ice,is wet,ice thickness, ice elevation,bedrock bathy,init draft) No ice, dry → bedrock continent draft(i, j) = 0 bedrock bathy(i, j) = 0 f.ETISh keeps the front from moving but the post-processing can lead to such cases. if draft(i, j) > 0 and init draft(i, j) = 0 then draft(i, j) = 0 Force retreat advancing ice-shelf front else if draft(i, j) = 0 and init draft(i, j) > 0 then draft(i, j) = 10 m Force fill retreating ice-shelf front 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 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.  Kobayashi et al., 2015) differences in seasonally averaged 2m air temperature for the second simulated year. The used colorbar is the same as in Fig. 11. Hatching denotes areas where the (PARASO -JRA-55) and (PARASO -ERA5) (see Fig. 11) differences are of contrary signs.    Figure A9. Integrated global Antartic 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 2 ) at distinct frequencies (3 month 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). In (b), the gap around 1988 is due to a shortage in observational data.
While only one two-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 along 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:

750
-Reducing the COSMO time step from its AEROCLOUD  value of 120 s to 90 s, in order to avoid CFL-criterion errors. Such errors may probably linked to the finer PARASO vertical resolution compared with AEROCLOUD (+50%).
-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 755 requires transfer coefficients (C D , C H ) and surface properties (T S , q s ) instead of fluxes per se. Hence, transfert 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 (moistures) 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 760 identified and traced from timestep-by-timestep 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.6GHz. With these resources, the wall time was about 3.5hr for one 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 765 required for simulating one month and the parallel efficiency of PARASO, defined as: e(n proc ) = n min proc T (n min proc ) n proc T (n proc ) where n proc is the number of cores used; n min proc = 56 is the minimal number of cores required for PARASO (constrained by random access memory requirements); and T (n proc ) is the wall time required for one month of PARASO simulated time. As   Fig. B1 shows, PARASO scales reasonably well up to ≈ 10 3 cores.

770
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 inbetween NEMO -CCLM 2 legs. Its computational cost is negligible (less than 0.1%) Figure B1. (a) Wall time required for simulating one month and (b) parallel efficiency of PARASO (see Eq. B1). Note that on (a), both axes are log-scaled; on (b), the x-axis is.  compared with all other models listed in Table B1. The same applies for the f.ETISh-including offline coupling interfaces: 775 the equivalent walltime they require is negligible compared to what is shown in Table B1. This was to be expected, as these coupling interfaces only process 2D 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 remain a considerable computer science challenge in the context of climate models . Like most climate model configurations,

780
PARASO feature restart procedures, during which the content of each model's prognostic variables are 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 785 of a leg for ensuring 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, 2015). The magnitude of the differences in 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 instabilites. The non-transparency of restart procedures on one of them (Breniac) is a strong argument supporting the non-replicability across distinct architectures. It should 795 however be noted that roughly speaking, simulations 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 2 differences between AEROCLOUD and PARASO In this appendix, we list the major differences in the CCLM 2 model (tuning parameters, domain, parameterizations) between AEROCLOUD (a previous CCLM 2 stand-alone Antarctic configuration, see Souverijns et al., 2019) and PARASO.

800
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 (AEROCLOUD) and 14 PARASO levels (11 in AEROCLOUD) 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 s to 90 s.

805
The spatial domain in AEROCLOUD 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 • S 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 810 have been made in AEROCLOUD. 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.

815
More information can be found in Souverijns et al. (2019). The inclusion of the AEROCLOUD 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.

820
Compared to AEROCLOUD, 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 2 dynamics is preserved.
PARASO uses the default snow roughness length of CLM4.5 of 2.4 · 10 −3 m, which is constant in time and space. This variables are incorrect and that in the AEROCLOUD 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 subgridscale 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, The main novelties are linked to the TKE-derived surface transfer scheme (Raschendorfer, 2005). On every cell, this scheme updates:

845
κ, 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;

850
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 2m except winds which are at 10m).

D1 Turbulent transfer coefficients and surface properties 855
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: where x j fj = f 0 x 0 + f 1 x 1 , and (C j h ) are the tile-specific heat transfer coefficients. The pseudo-averaged surface moisture is 860 assessed by enforcing mass-flux conservation, i.e., by determining q eq s so that the mass flux absorbed by CCLM is the tileconcentration average of tile-specific mass fluxes. This yields: q low is the lowest-level moisture, (T j s ) the tile-specific surface temperatures, (F j q ) the tile-specific mass fluxes linked to latent 865 heat, p s the surface pressure, u low h the lowest-level horizontal wind norm, R d and R v the gas constants for dry air and water vapor, respectively. In case Eq. (D3) leads to q eq s < 0 (which in practice happens in less than once in 10 8 calculations), we set q eq s = 0 and accordingly increase C eq h and C eq h to enforce mass flux conservation. Similarly to q eq s , the pseudo-averaged surface temperature T eq s is assessed by enforcing sensible heat conservation, which yields: where θ low is the lowest-level potential temperature and c a p 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 relying on the tile-averaged surface temperature T j s fj instead of T eq s for evaluating a small surface pressure correction term. The equivalent momentum transfer coefficient is: where θ eq s,v = T eq s 1 + Rv R d − 1 q eq s is the surface virtual potential temperature and (τ j ) the tile-specific wind stresses.
It should be underlined that since surface currents or ice velocities are neglected in wind stress computations, all (τ j ) are coaligned with u low h , hence averaging τ j is transparent. Equation (D5) ensures momentum flux conservation. The pseudo-averaged surface roughness length is assessed as a geometric average, tile-concentration wise: A geometric average is used here because ln z eq 0 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 case such diagnostics are not readily 885 applicable, a simple algebraic averaging is performed. κ bot , the bottom-level TKE velocity, is thus assessed from C eq d : where c tke = 3 √ 16.6. The bottom momentum and heat diffusivity are then updated as: In Eqs. (D12)-(D13), ν h,0 = 2.2 × 10 −5 m 2 s −1 is the scalar conductivity and ν m,0 = 1.5 × 10 −5 m 2 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, dewpoint 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 adjustements 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: α eq = f 0 α 0 + f 1 α 1 (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.: where (T j S ) are the tile-specific surface temperatures, d = 0.996 is the CCLM default emissivity, 0 = 1 is the NEMO ocean emissivity, and 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 couplinginduced biases. This appendix documents the PARASO parameters that were changed from stand-alone CCLM 2 and NEMO simulations (see Table E1), and presents results from some of our tuning experiments (see Fig. E1). Here we present the 915 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 920 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 2 , compared to NEMO stand-alone, 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, 925 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 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 930 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 Table E1. Details on PARASO tunings. For each parameter, two values or methods are given, v1 (raw, nontuned) or v2 (tuned). The "New?" column specifies whether this tuning parameter was introduced for PARASO (Y) or already preexisting in the model (N are all within the observational range. Moreover, default LIM3 sea-ice albedo parameters had mostly been tuned for Arctic 935 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. 8(b)) 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 are quite significant (+50%), the uncertainty of their actual value are also quite 940 high (Tsamados et al., 2014;Ungermann et al., 2017). The ERA5 data (Hersbach et al., 2018) was downloaded on 01-SEP-2019 from the Copernicus Climate Change Service (C3S) Climate Data Store. The results contain modified Copernicus Climate Change Service information 2020. Neither the European Commission nor 960 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) was downloaded on 01-FEB-2020 from the National Snow and Ice Data Center.