Ocean–Ice Sheet Coupling in the Totten Glacier Area, East Antarctica: Analysis of the Feedbacks and Their Response to a Sudden Ocean Warming

Abstract

We coupled together high-resolution versions of the ocean–sea ice model NEMO and the ice sheet model BISICLES configured to the Totten Glacier area and ran a series of simulations over the recent past (1995–2014) and under warming conditions (2081–2100; SSP4-4.5) with NEMO in stand-alone mode and with the coupled model to assess the effects of the coupling. During the recent past, the ocean–ice sheet coupling has increased the time-averaged value of the basal melt rate in both the Totten and Moscow University ice shelf cavities by $6.7\%$ and $14.2\%$, respectively. The relationship between the changes in ice shelf thickness and ice shelf basal melt rate suggests that the effect of the coupling is not a linear response to the melt rate but rather a more complex response, driven partly by the dynamical component of the ice sheet model. The response of the ice sheet–ocean coupling due to the ocean warming is a 10% and 3% basal melt rate decrease in the Totten and Moscow University ice shelf cavities, respectively. This indicates that the ocean–ice sheet coupling under climate warming conditions dampens the basal melt rates. Our study highlights the importance of incorporating ocean–ice sheet coupling in climate simulations, even over short time periods.

1. Introduction

The Totten Glacier area, located on the Sabrina coast in East Antarctica, holds the largest potential ice discharge from the East Antarctic ice sheet [1], with an ice volume equivalent to 3.9 m of global sea level rise [2]. With most of its drainage basin grounded well below sea level [3], the Totten Glacier is susceptible to ocean erosion and to ice sheet instability [4].

Recent observations indicate the thinning of the Totten ice shelf (TIS) during the last decade [5], with most of this thinning driven by higher ocean temperatures [6,7]. In addition to this thinning, a 1 to 3 km retreat of the ice shelf grounding line has been detected during the 1996–2013 period [2]. Observations of ice motion over the last few decades have revealed a complex pattern of changes, with the counteraction of increases and decreases in surface ice speed. The glacier exceeded its balance speed in 1989–1996, slowed down in 2000 to bring its ice flux into balance with accumulation, and then accelerated until 2007, and its speed remained constant thereafter [7,8]. As for the ice thinning, there is a strong linkage between the ice dynamics and the ocean temperature [7]. Several studies reported the presence of warm, salty modified Circumpolar Deep Water (mCDW) in the bottom layer of the continental shelf near the TIS cavity, which drives rapid basal melt [9,10,11,12,13]. Other studies have also investigated the effect of the Dalton polynya [5,9], the wind [14], and the ocean intrinsic variability [15] on the TIS basal melt rate.

The ocean can induce modifications in continental ice sheets through changes at the ocean–ice shelf interface [16]. As the ice shelf controls the distribution of normal stresses at the grounding line, changes in ice thickness and grounding line reduce the back-stress that the ice shelf exerts on the upstream ice sheet, leading to an acceleration in the grounded ice speed and to enhanced ice discharge, and thus contributes to the global sea level rise [17,18]. In turn, the induced melting results in freshwater injection at depth, known as Ice Shelf Water. This Ice Shelf Water flows upward along the ice shelf base under the influence of buoyancy, friction, and Coriolis forces in the form of a turbulent plume, continually mixing with warmer seawater [19]. The upward motion of the plume induces an inflow of possibly warmer ocean water into the ice shelf cavity, creating more melt. On the other hand, depending on the stratification of the ocean water inside the cavity, the plume may reach a level of neutral buoyancy from which it is no longer driven upward. The subsequent melting, freezing and density-driven plume circulation alter nearby seawater properties [20,21]. Moreover, ice shelf basal melting is sensitive to the ice shelf cavity shape [16,22].

Over the last few years, an increasing number of coupled ocean–ice sheet models have been developed, such as Global Climate Models (GCMs) of intermediate complexity [16,23], state-of-the-art ocean models [24,25,26] and even state-of-the-art GCMs and Regional Climate Models (RCMs) resolving all the climate components, called Earth system models [27,28]. Several studies, with regional coupled ocean–ice sheet models, have investigated the impact of the coupling on various ice shelf cavities. Seroussi et al. [29] used the coupled MITgcm–ISSM model to study the coupling effects on the Thwaites Glacier under present climate conditions. They observed a decrease in basal melt in response to the thinning of the ice shelf, which causes the ice shelf base to rise above the mCDW layers, before the ice starts increasing again (negative feedback loop). Timmermann and Goeller [24] used the coupled FESOM–RIMBAY model in the Filchner-Ronne Glacier area under present climate conditions. They found higher basal melt rates in coupled mode compared to uncoupled mode, under present-day climate conditions. They also found that an increased water column thickness, due to ice thinning, allows for more transport of warm water to the deepest part of the cavity, increases the slope of the ice base and reinforces the strongest melt rates (positive feedback loop). A coupled ocean–ice sheet study in the TIS area has already been conducted by Pelle et al. [25], using the coupled MITgcm–ISSM model, where they found a 56% increase in TIS melt in a warming climate scenario compared to a control run, and a widespread grounding line retreat along the TIS eastern grounding zone.

Nevertheless, all these studies hold some limitations, such as the absence of tides in all models, which might underestimate the ice shelf basal melt rates as described in Huot et al. [30], or such as constant temperature and salt exchange coefficients in the basal melt rate formulation in Seroussi et al. [29], which should be velocity-dependent following Jenkins et al. [31]. Moreover, these studies use ocean–ice sheet models with relatively coarse horizontal resolutions of the model ocean components (4 km in Pelle et al. [25] and from 1 km at the grounding line to 16 km at the ice shelf front in Timmermann and Goeller [24]), thus not allowing them to resolve eddies, which play a crucial role in controlling the flow across the continental shelf and in the vicinity of grounded zones, or the relatively coarse horizontal resolution of the ice sheet model components (10 km in Timmermann and Goeller [24], from 500 m to 15 km in Seroussi et al. [29] and from 1 km to 20 km in Pelle et al. [25]).

As these ice shelf–ocean interactions may lead to drastic changes in the Totten ice sheet’s stability in the future, it is essential to properly and explicitly resolve these ocean–ice stream interactions, with the evolving grounding line position and ice shelf thickness, in a coupled ocean–ice sheet model framework. Furthermore, in order to adequately resolve the ice–ocean interactions, such a framework must include a high-resolution ocean model capable of partially resolving eddies, velocity-dependent heat and salt exchange coefficients at the ice–ocean interface, an explicit representation of fast ice and the representation of grounded icebergs. The goal of the present study is twofold. Firstly, we aim to model the ocean–ice sheet interactions in the Totten Glacier area and study the effects of the coupling on the ice shelf cavities. Secondly, we aim to assess the response of the coupling to sudden warming. To do so, we couple high-resolution, regional configurations of the ocean–sea ice model NEMO3.6-LIM3 (resolution of approximately 2 km) and the ice sheet model BISICLES (500 m resolution for the fast-flowing ice). Both regional configurations cover the Totten Glacier area in East Antarctica. Then, we compare coupled and ocean stand-alone simulations over the years 1995–2014 and with EC-Earth3 anomalies (2081–2100 [32]) applied at the atmospheric and oceanic forcings of the ocean–sea ice model.

This manuscript is organised as follows. The models, regional configurations, model initialization, coupling interface and experimental design are detailed in Section 2. We analyse the effects of the coupling on the ice shelf cavities and on the ocean and sea ice properties, and the response of the coupling to a sudden warming, in Section 3. Conclusions are finally given in Section 4.

2. Materials and Methods

2.1. Ice Flow Model: BISICLES

Changes in land ice and ice shelves are simulated by the ice sheet model BISICLES [33], a finite-volume model with vertically integrated ice flow based on the hybrid approach of Schoof and Hindmarsh [34] to include membrane stresses. The BISICLES temporal and spatial adaptive mesh refinement is set between 4 km in the slow-moving regions and up to 500 m resolution around the grounding line and fast-flowing areas, where membrane stresses are important. Time steps are adapted automatically to satisfy the Courant–Friedrichs–Lewy condition. At each time step, the ice velocity, surface and base elevation and grounding line position are updated. Basal sliding is based on a Coulomb-limited sliding law after Tsai et al. [35]: ${\tau }_b=-min(\alpha N,C|u_b{|}^m)·|u_b{|}^{-1}{u}_b$ with ${\tau }_b$ the basal shear stress, $\alpha$ = 0.5, N the effective pressure at the bed, C the spatial sliding law coefficient, m = 1/3, and $u_b$ the basal velocity. Throughout the simulation, the ice front position is kept fixed, filling ice-free areas with a minimum ice thickness of 11 m to avoid changing the open ocean area in NEMO, while keeping the buttressing effect of the ice shelf itself low. Basal melt is applied solely on fully floating grid cells, and partial melting underneath grounded grid cells is avoided [36].

The BISICLES Configuration

The ice domain covers the Aurora basin, including the catchments of the Totten and Moscow University ice shelves [37] that are defined by surface slopes from Ice, Cloud and Land Elevation satellite (ICESat) data [38]. Bathymetry, bedrock and ice surface elevation are provided by Bedmachine v.2 [3,39]. For the model initialisation, we first estimate the basal friction and ice stiffness coefficients by solving an inverse problem following MacAyeal [40] and Morlighem et al. [41]. The inversion aims to minimise mismatch between modeled and observed ice surface velocities [42,43], using a conjugate gradient method. Given the estimated coefficients, the ice sheet model is run for 12 years to minimise model drift while keeping the ice thickness as close to present-day observations (Bedmachine v.2) as possible. During the 12-year relaxation simulation, the model is forced by the constant 1995–2014 mean surface mass balance from the regional atmospheric model MAR (Modèle Atmosphérique Régional) at 10 km resolution [44] and mean NEMO melt rates from Van Achter et al. [45] at approximately 2 km resolution around the grounding line. In the coupled simulations, BISICLES is forced by monthly mean ocean melt rates provided by NEMO. To have an estimate for the cells where ice retreats during a coupling window, basal melt rates from NEMO are linearly extrapolated on the 500 m grid to cover also grounding line and first grounded elements. Since the focus of this study is the effects of ice–ocean interactions, the surface mass balance is kept constant beyond the initialisation.

2.2. Ocean and Sea Ice Model: NEMO3.6-LIM3

We use a regional East Antarctic configuration of the NEMO model (Nucleus for European Modelling of the Ocean, [46]) to simulate the ocean–sea ice system. NEMO includes the ocean model OPA (Océan Parallélisé), which is coupled with the Louvain–la Neuve sea ice model (LIM3, [47,48]). This model version, hereafter referred to as NEMO, has been frequently used to study both the whole Southern Ocean, e.g., [28,49,50], and specific Antarctic sectors, e.g., [30,51,52]. The numerical discretisation relies on a C grid [53] and is based upon primitive equations, assuming hydrostatic balance and using the Boussinesq approximation [54]. The vertical discretisation relies on the so-called $z^{*}$ coordinate, with varying cell thicknesses on the full water column [55]. The vertical mixing is based on a turbulent kinetic energy (TKE) scheme [56,57,58]. In our setting, ocean convection is represented through enhancing the vertical diffusivity to 20 m${}^2$ s${}^{-1}$ in case of an unstable water column. The drag coefficient is set to 7.1 × 10${}^{-3}$ at the sea ice–ocean interface and 2 × 10${}^{-3}$ at the sea ice–atmosphere one [59]. The sea ice model relies on an elastic–viscous–plastic (EVP) rheology formulated on a C grid [60], and uses a five-category subgrid-scale distribution of sea ice thickness [61].

The Antarctic ice shelf cavities are opened to the ocean circulation, with explicit ocean–ice shelf interactions [62]. The heat and mass fluxes associated with the melt or freezing are determined by the three-equation formulation of Jenkins [63]. The transfer coefficients for both heat (${\gamma }_T$) and salt (${\gamma }_S$) between the ocean and ice shelves are velocity-dependent [64]: ${\gamma }_{T,S}={\Gamma }_{T,S}\times u_{*}$. The friction velocity is given by $u_{*}=C_d\times \sqrt{u_{TML}^2}$ and constant values of ${\Gamma }_T$ and ${\Gamma }_S$ taken from Jourdain et al. [51] are used (${\Gamma }_T=2.21\times {10}^{-2}$ and ${\Gamma }_S=6.19\times {10}^{-4}$ for temperature and salinity, respectively), $C_d$ = $4\times {10}^{-3}$ is the ocean–ice sheet drag coefficient and $u_{TML}$ is the ocean velocity in the top mixed layer, which is either the top 30 m of the water column or the top model layer (if thicker than 30 m) [65].

The Totten24 Model Configuration

The NEMO Totten24 setup is similar to the configuration used in Van Achter et al. [52] and Van Achter et al. [45]. The ocean grid is a 1/24${}^{\circ }$ refinement (less than 2 km grid spacing) of the eORCA1 tripolar grid, centered on the continental shelf in front of the TIS, East Antarctica, and covering an area between 108–129${}^{\circ }$ E and 63–68${}^{\circ }$ S (Figure 1). The vertical discretisation includes 75 levels, with increasing thickness from ∼1 m at the surface to ∼200 m at depth. Partial cells are used for a better representation of the seafloor and the ice shelf base [66]. The OPA and LIM time steps are 150 s and 900 s, respectively. The ocean layer directly underneath the ice shelf base varies between 30 m near the cavity front and 80 m in the center of the cavity. The bathymetry dataset comes from the NASA MEaSUREs program (BedMachine Antarctica v2, [3]), which contains a bathymetry map of Antarctica based on mass conservation, streamline diffusion and other methods [3]. The initial ice shelf thickness for the ocean model is the relaxed geometry obtained after a 12-year relaxation simulation as described in Section 2.1. The Antarctic surface continental mask is constant time-wise.

At the lateral open-ocean boundaries, NEMO is driven using outputs from a 1979–2014 simulation with an eORCA025 (1/4${}^{\circ }$) NEMO configuration covering the whole Southern Ocean [28] (hereafter referred to as PARASO). PARASO was chosen over reanalyses such as GLORYS12v1 due to its much better representation of ocean temperatures over the continental shelf, which compares well with observations from Rintoul et al. [10]. In order to correct a fresh bias present in PARASO, the ocean salinity lateral boundary condition is uniformly increased by 0.25 g kg${}^{-1}$. The model is forced daily at the west, north and east boundaries, with SSH, temperatures, salinities, ocean velocities, sea ice concentrations, sea ice thicknesses and snow thicknesses over sea ice. A flow relaxation scheme [67] is applied at the lateral boundaries to the three-dimensional ocean variables and two-dimensional sea ice variables. To allow gravity waves generated internally to exit the model boundaries, a Flather radiation scheme [68] is used at the lateral boundary conditions for barotropic velocities and sea surface elevation (SSH). Furthermore, the SSH and barotropic velocities from the FES2014 tide model [69] are added to the boundary for the tide components K1, K2, M2, P1, O1, S2, 2N2, Mm, M4, Mf, Mtm, MU2, N2, NU2, Q1, S1, L2 and T2, as in Maraldi et al. [70], Jourdain et al. [71] and Huot et al. [30]. Surface boundary fluxes are computed using the CORE bulk formula [72], with atmosphere input coming from ERA5 [73]. No surface salinity restoring is applied.

While the iceberg calving is not dynamically simulated by the model, grounded icebergs are extracted from the remotely sensed mosaic ‘RAMP AMM-1 SAR Image Mosaic of Antarctica, Version 2’ [74] and incorporated into the model as follows: a model cell is imposed as a grounded iceberg (bathymetry value is set to zero) when the cumulative area of all observed icebergs within the cell reaches 2 km${}^2$. Then, to avoid creating iceberg walls and to ensure regular oceanic circulation, we filter the icebergs to avoid prescribing neighbouring iceberg grid cells. Finally, only the icebergs located in oceanic areas shallower than 450 m are kept. There is no freshwater injection related to the icebergs melting in the model. A fast ice representation is included through the combination of the effect of grounded icebergs with a sea ice tensile strength parameterisation [52]. The EVP rheology of the sea ice is modified to include both isotropic and uniaxial tensile strength based on the tensile strength parameterisation of Lemieux et al. [75]. The eccentricity of the elliptical yield curve and the isotropic tensile strength parameter are set to 1.2 and 0.2, respectively. The number of sub-iterations of the EVP scheme is set to 720 to solve the momentum and the stress equations.

Figure 1. Model bathymetry and domain. The contour interval is 50 m up to 500 m depth and 500 m up to 4500 m depth. Ice shelf cavities are surrounded by a thick black line. The 0.75 fast ice observed frequency from [76] is shown by the shaded gray areas. The inset (top left) shows Antarctica with the model domain denoted by the grey area (ice sheet model domain) and by the red rectangle (ocean model domain).

Figure 1. Model bathymetry and domain. The contour interval is 50 m up to 500 m depth and 500 m up to 4500 m depth. Ice shelf cavities are surrounded by a thick black line. The 0.75 fast ice observed frequency from [76] is shown by the shaded gray areas. The inset (top left) shows Antarctica with the model domain denoted by the grey area (ice sheet model domain) and by the red rectangle (ocean model domain).

2.3. Initialisation and Coupling Interface

For the ocean–ice sheet coupled and ocean stand-alone simulations, the NEMO model is initialised from a 2-year spin-up run performed with NEMO stand-alone, using forcings corresponding to the years 1993 and 1994. This spin-up is started from a 20-year climatology of a 1995–2014 NEMO stand-alone run, which is long enough to reach equilibrium in the ice shelf cavities [45]. The 20-year climatology of the 1995–2014 NEMO stand-alone run is used to provide an average basal melt rate distribution that is later applied in a BISICLES stand-alone relaxed run (described in Section 2.1). This run is performed to generate the initial geometry that will later serve as an initial condition for NEMO for coupled and stand-alone simulations.

The coupling interface between NEMO and BISICLES is based on the coupling interface between NEMO and the ice sheet model f.ETISh used in PARASO [28], and following Smith et al. [27] for ice sheet geometry changes within NEMO. Both the model initialisation and the coupling procedure, which is an off-line coupling, are summarised in Figure 2. NEMO runs first for a coupling time window (here 3 months), sending monthly time series of ice shelf melt fluxes to BISICLES. While NEMO waits, BISICLES runs for the same time window, forced with the received melt rates as a boundary condition. Afterwards, BISICLES provides NEMO with a new ice shelf geometry. NEMO can then start the next 3-month period with its updated ice shelf geometry. These exchanges between both models are repeated over each coupling window until the end of the coupled simulations. The data are bilinearly interpolated from one grid to the other.

When NEMO receives the updated ice shelf geometry from BISICLES and before NEMO starts running again, ocean cells may be opened or closed due to the change in ice shelf geometry. To account for these potential changes in geometry, we follow the methodology of Smith et al. [27], as described in Pelletier et al. [28]. For any subshelf ocean cells to open (or close), the following constraints must be met; otherwise, the ocean cell remains closed (or opened). (1) The ice shelf thickness must be greater than 11 m to keep the ice shelf extent constant and (2) the water column thickness must be greater than 50 m. When an ocean cell is opened, the SSH, the temperature and the salinity from neighbouring cells are extrapolated into the newly opened cell [77]. The initial ocean current velocities of the new opened cells are set to zero and a horizontal divergence correction is applied in order to avoid spurious vertical ocean velocities. This approach does not comply with energy and mass conservation, but the amounts are so small that they are negligible for the overall outcome.

When BISICLES receives monthly mean subshelf melt fluxes, mask discrepancies are addressed with a bilinear extrapolation to cover BISICLES ice shelf cavity columns, which correspond to closed columns in NEMO. Ice shelf cavity columns in NEMO, which correspond to grounded ice or open ocean in BISICLES, are disregarded.

2.4. Experimental Design

Our experimental design consists of four 20-year-long simulations (

Table 1. Names and descriptions of the simulations used in this study.

	Model	Forcing and Lateral Boundary Conditions
REF	NEMO stand-alone	Recent past (ERA5, PARASO, 1995–2014)
REF3mcpl	BISICLES–NEMO coupled	Recent past (ERA5, PARASO, 1995–2014)
WARM	NEMO stand-alone	REF + anomalies derived from EC-Earth3 (2081–2100)
WARM3mcpl	BISICLES–NEMO coupled	REF3mcpl + anomalies derived from EC-Earth3 (2081–2100)

). The first set of two simulations, devoted to the study of the ocean–ice sheet coupling and its impact on the ice–ocean interactions, is composed of the REF and REF3mcpl simulations, which cover the 1995–2014 period. REF is a NEMO stand-alone simulation and REF3mcpl is a coupled NEMO–BISICLES simulation with a coupling time step of 3 months. Both simulations have the same initial ice shelf thickness, initial ocean conditions and atmospheric and oceanic forcings.

The second set of two simulations, which aims to investigate the response of the ocean–ice sheet coupling to sudden warming, is formed by the WARM and WARM3mcpl simulations. Both simulations correspond to climate conditions that might prevail during the 2081–2100 period. As described in Van Achter et al. [45], annual cycles of the EC-Earth3 climate anomalies from the SSP4-4.5 scenario [32], computed as the differences between 2081–2100 and 1995–2014, have been added to all the fields of the atmospheric and oceanic forcings used for the 1995–2014 period in REF and REF3mcpl. The 2081–2100 simulations use the same initial ice shelf geometry as the 1995–2014 simulations. This experiment is designed to estimate the effects of a present-day versus warmer climate scenario to an ice–ocean coupling setup on a regional scale. Where keeping the same initial ice geometry facilitates the comparison of simulation results, we do not attempt to project the ocean, sea ice, ice sheet and ice shelf characteristics that might prevail by the end of the 21st century. Rather, WARM3mcpl aims at evaluating the response of the NEMO–BISICLES coupling to the sudden warming of ice shelf cavities triggered by changes in atmospheric forcings and oceanic conditions at the boundaries of the domain.

A similar simulation to REF was conducted by Van Achter et al. [52] and evaluated against observations (sea ice concentration, fast ice, sea ice production, sea ice thickness, polynya locations and ocean temperature and salinity distributions) and both REF and WARM simulations are compared in detailed in Van Achter et al. [45].

3. Results

3.1. Ocean–Ice Sheet Coupling Effects on the Ice Shelf Cavities

Time series of area-averaged ice shelf basal melt rates in REF and REF3mcpl (continuous and dotted blue lines in Figure 3) indicate that the NEMO–BISICLES coupling over the years 1995–2014 tends to slightly increase the basal melt rate in the TIS cavity and its neighbouring Moscow University ice shelf (MUIS) cavity. In both cavities, the coupled simulation shows higher basal melt rates with a clear increase after the fifth year of the simulation. The TIS basal melt rates averaged over the years 1995–2014 are 10 ± 1.88 and 10.68 ± 1.95 m year${}^{-1}$ for REF and REF3mcpl, respectively. Likewise, the MUIS basal melt rates are 6.2 ± 2.38 and 7.09 ± 2.58 m year${}^{-1}$ for REF and REF3mcpl, respectively. These simulated basal melt rates are in good agreement with the estimates of Rignot et al. [78], namely 10.74 ± 0.7 m year${}^{-1}$ for TIS and 4.7 ± 0.8 m year${}^{-1}$ for MUIS. The basal melt rates in REF3mcpl are increased by 6.7% in TIS and 14.2% in MUIS, and the differences between the two experiments are much smaller than the simulated interannual melt rate variability (SD in REF: 3.56 m year${}^{-1}$ in TIS and 5.68 m year${}^{-1}$ in MUIS), which is partly driven by the oceanic lateral boundary conditions (correlation of 0.5 between the ocean temperature at the ocean boundaries and the MUIS basal melt rate). The effect of the coupling on the basal melt rate is not constant over time and fluctuates between +0.5 and +2 m year${}^{-1}$ and is related to the interannual variability of the basal melt rate: a high melting year is related to a stronger effect of the coupling on the basal melt rate and vice versa.

The spatial distributions of the mean basal melt rates for REF and REF3mcpl are displayed in Figure 4a,b, respectively. Maximum melt rates in both simulations are around 60 m year${}^{-1}$ and occur in the deepest part of the cavity near the grounding line, where the ice base reaches 1100 and 1300 m below sea level, and the freezing point is approximately 1 ${}^{\circ }$C below the surface freezing point. Differences between REF3mcpl and REF (Figure 4c) show that the melt rates below the TIS in the coupled run are decreased along the eastern flank of the cavity and increased along the western flank. Moreover, in both the TIS and MUIS cavities, high values of the basal melt rate in REF are found in both the middle and the bottom of the cavity in REF, while, in REF3mcpl, they are located only in the bottom of the cavity. Furthermore, part of the basal melt rate heterogeneity in REF is smoothed out in REF3mcpl, which exhibits a relatively more homogeneous spatial distribution of the basal melt rate compared to REF.

The impacts of the coupling on the ocean temperature and velocity at the ice–ocean interface (first ocean layer beneath the ice shelves) are depicted in Figure 5 and Figure 6, respectively. The ocean temperature at the ice–ocean interface is decreased in the coupled simulation in all cavities near the cavity fronts, down to −0.2 ${}^{\circ }$C. On the other hand, the ocean temperature is increased by more than 0.1 ${}^{\circ }$C near the bottom of the cavities. The TIS shows a higher mean ocean temperature than the MUIS, with values reaching −0.9 ${}^{\circ }$C at the bottom of the TIS cavity and −1.8 ${}^{\circ }$C in the MUIS cavity (see Figure 5a). Despite their differences in mean ocean temperature, both TIS and MUIS cavities exhibit the same ocean warming at the ice–ocean interface in the coupled simulation compared to the fixed ice shelf geometry simulation (see Figure 5b). The cooling, however, is much larger at the TIS cavity front than at the MUIS cavity front. The differences in cooling between both fronts of the TIS and MUIS cavities might partly be related to the larger amount of cold melt water released by the MUIS cavity, which cools down the ocean (even at the surface, as shown in Figure 5b) on the western side of the MUIS in REF3mcpl compared to REF. Cold Ice Shelf Water (water colder than −1.92 ${}^{\circ }$C) is found in the MUIS cavity but does not exit the cavity, as it is warmed up along its way towards the ice shelf front.

The ocean current speeds at the ice–ocean interface are intensified in the coupled simulation compared to the fixed ice shelf geometry simulation. The mean ocean velocities in REF at the ice–ocean interface are characterised in both the TIS and MUIS cavities by a clockwise circulation, with a slow entering current along the eastern side of the cavities and a faster (up to 0.2 m s${}^{-1}$) outgoing current along the western side (Figure 6a). In the coupled run, the ocean speeds at the ice–ocean interface are enhanced up to 0.1 m s${}^{-1}$, mostly in areas where the ice shelf thickness thinned during the coupling (as on the western side of the TIS cavity). Other regions display a slowdown of the ocean currents due to the change in ice shelf thickness, as on the eastern side of the TIS cavity where the ice shelf has thickened.

Vertical sections of ocean temperature, salinity and meridional velocity across the TIS cavity are shown in Figure A3 and Figure A4 (in Appendix A). The bottom of the cavity is filled with colder (−0.1 ${}^{\circ }$C) and fresher waters (−0.01 g kg${}^{-1}$) in REF3mcpl, which is the opposite of what is found at the ice–ocean interface (warmer and saltier in REF3mcpl). The ocean stratification seems to be similar in both REF3mcpl and REF. The temperature and salinity effects on the stratification at the ice–ocean interface are opposite (higher temperature and higher salinity in REF3mcpl) and may explain the lack of changes in stratification inside the cavity. The meridional circulation inside the TIS cavity shows, at the ice–ocean interface, a steady rising current from the bottom to the front of the cavity. This plume circulation exhibits the same intensity in both simulations, except where the ice shelf thickness is modified in REF3mcpl, where the plume circulation is enhanced.

The changes in the ice shelf draft (ice shelf thickness below sea level) and grounding line position after 20 years of the coupled simulation are illustrated in Figure 7a,b, respectively. The TIS experiences a thickening near the front of the cavity, with an ice gain reaching more than 150 m. Deeper into the cavity, the ice thickness is decreased by more than 200 m. As illustrated in Figure A2, both the thinning and thickening of the ice shelf have smoothed the ice–ocean interface. The MUIS exhibits a more heterogeneous change in thickness near its cavity front, and a large ice thickness increase near the bottom of its cavity. There is no change in the grounding line position in the two main ice shelf cavities. Nonetheless, there is a grounding line migration around 119${}^{\circ }$ E, which closes half of a small ice shelf cavity between TIS and MUIS, and is explained by the lower ocean temperature at the ice–ocean interface (−0.1 ${}^{\circ }$C) due to the increased basal melt rate in its neighbouring MUIS cavity in REF3mcpl.

A recent study by Li et al. [79] observed a grounding line retreat over the 1996–2020 period of 9.37–13.85 km at the eastern and western flanks of the MUIS cavity and of 5.95 and 3.51 km at the TIS cavity. Even the maximum grounding line retreat reaches less than 1 km year${}^{-1}$, which may explain why our 2 km ocean horizontal resolution does not reproduce this grounding line retreat in REF3mcpl. Moreover, the threshold value of a minimum of 50 m water column thickness for a cell to become ungrounded may also limit the grounding line retreat, and the basal melt rates along the grounded zone are relatively low (less than 15 m year${}^{-1}$) and may result in underestimated grounding line retreat in REF3mcpl.

Cross-sections of the ice shelf base elevation for both the TIS and MUIS are presented in Figure A2a,c (in Appendix A). In the TIS zone, the ice shelf base is unchanged in REF3mcpl from the grounding line until 67.1${}^{\circ }$ S, where the changes in ice shelf thickness increase the slope of the ice shelf base, which may reinforce the plume circulation and locally enhance the basal melt. Further north, around 67${}^{\circ }$ S, the slope of the ice shelf base in REF3mcpl is strongly decreased compared to the one in REF, which decreases the basal melt rates in REF3mcpl around this latitude, as shown in Figure 4. Near the ice shelf front, the ice shelf is thickened in REF3mcpl, which should enhance the basal melt due to the increased freezing temperature. On the other hand, the slope of the ice shelf base is decreased and may explain the locally slower plume circulation and subsequent lower basal melt rates in REF3mcpl near the TIS front. The year-to-year changes in ice shelf base depth (each line is a different year in Figure A2a,c) indicate that the thinning of TIS has reached its equilibrium after 5 years (changes in ice shelf base depth between year 5 and 20 are very small). On the other hand, the ice thickening at the TIS front is not yet at equilibrium (ice shelf base depths between year 15 and 20 are not yet overlapping). The changes in ice shelf base elevation for the MUIS are less significant.

The next step is to investigate the origin of the changes in basal melt rates between REF3mcpl and REF. Figure 8a indicates that the link between the changes in ocean temperature and the changes in basal melt rate is weak. Their linear regression shows a relationship between a higher ocean temperature at the ice–ocean interface and higher basal melt and vice versa, but the R value is low (0.129). On the other hand, there is a clear link between the changes in ocean current speed at the ice–ocean interface and the changes in basal melt rate (Figure 8b). The linear regression reveals, with an R value of 0.573, that higher basal melt rates in REF3mcpl compared to REF are generally related to increased ocean current speeds in REF3mcpl compared to REF, and the same is true for lower basal melt rates and decreased ocean current speeds. The relationship between changes in ice shelf base slope and basal melt rates, shown in Figure A1, indicates that the basal melt rates in REF3mcpl are related to steeper slopes compared to the basal melt rates in REF.

Areas of changes in mean basal melt rate in the TIS cavity related to changes in ocean temperature and current speed are highlighted in Figure 9. Regions of increased basal melt rate linked to higher ocean temperatures are few, and most of them are limited to the smaller cavity located on the eastern side of the main cavity. On the other hand, the areas of increased basal melt rate related to a larger ocean current speed cover a large part of the cavity—in particular, the western side of the cavity, where the higher basal melt rate increase occurs (see Figure 4c). Regions of decreased basal melt rate in REF3mcpl associated with a lower ocean temperature are also related to slower ocean currents, and are located near the cavity front (see Figure 9c,d). The decreased basal melt rates on the eastern flank of the cavity (see Figure 4c) are partly driven by decreased ocean current speeds (Figure 9d).

Under stable ocean vertical stratification, a reduced ice shelf thickness makes the in situ freezing temperature increase, thus implying a decrease in melt at constant inflow towards the ice shelf cavity. Therefore, with this assumption, the expected coupling response is that a reducing (increasing) ice thickness would eventually reduce (increase) the basal melt rate. Nevertheless, a more complex response is found in REF3mcpl, as revealed by Figure 10a, which shows the link between the 20-year basal melt rates and changes in ice shelf thickness. The scatter plot indicates that there are indeed areas with large values of basal melt rate related to larger ice shelf thicknesses, but there are also regions of lower ice shelf thicknesses caused by larger values of the basal melt rate. Moreover, most areas with a basal melt rate greater than 300 m are mostly related to a decreased ice shelf thickness.

Figure 10c suggests no clear link between changes in ocean current speed and changes in ice shelf thickness, with a small R value of 0.101. Areas of accelerated ocean current are found in ice shelf areas that are related to ice thickening as well as ice thinning (see Figure 10b,c). Nevertheless, in areas of high basal melt rates, there is a strong relationship between an increased ocean current speed, increased basal melt rates and reduced ice shelf thickness (Figure A5 in Appendix A). This positive feedback is explained by Timmermann and Goeller [24]. The reduced ice shelf thickness, which means a thicker water column, allows for more transport of warm water to the deepest parts of the cavity and allows for faster ocean currents at the ice–ocean interface. These faster ocean currents then enhance the basal melt rates.

To account for the ice dynamics (up to 2 km year${}^{-1}$ in the fast-flowing parts of the TIS; see Figure 11), the relationship between the divergence of the ice shelf thickness flux and the changes in basal melt rate is shown in Figure 11. The divergence of the ice shelf thickness flux expresses the loss of ice volume related to basal melting while taking into account the loss of ice volume related to horizontal ice transport. The ice shelf thickness flux divergence and the basal melt rate exhibit a strong linear relationship, with higher melt rates related to more ice volume loss (the effects of the coupling on the ice sheet model will be discussed in another study). Figure 10 and Figure 11 suggest that the changes in ice shelf thickness are not the simple integrated response of the changes in basal melt rate (thinner ice implies less melt), but are more complex. On the one hand, a reduced ice thickness may locally raise the basal melt rate. On the other hand, the changes in ice shelf thickness are also driven by the dynamics of the ice sheet. This result highlights the importance and the added value of the ocean–ice sheet coupling.

3.2. Ocean–Ice Sheet Coupling Impact on the Ocean Properties and Sea Ice over the Continental Shelf

The impact of the coupling outside of the ice shelf cavities is small. The ocean temperature and salinity properties over the continental shelf in front of the TIS are illustrated in Figure 12a. Due to the overall enhanced basal melt rate in both the TIS and MUIS cavities, the larger amount of cold melt water in REF3mcpl cools down and reduces the salt content of the water masses over the continental shelf. However, the changes are rather small, with temperature differences less than 0.05 ${}^{\circ }$C and salinity differences less than 0.01 g kg${}^{-1}$.

The fast ice frequency in REF and its changes compared to REF3mcpl are shown in Figure 12b,c. The frequency differences are close to zero everywhere, except on the western side of the TIS cavity, where the differences in fast ice frequency reach 0.2. This indicates that the lower ocean temperatures in front of the TIS cavity, due to the coupling, lead to up to two additional months of fast ice per year. Nevertheless, these results must be interpreted with caution as the close location of the oceanic boundary conditions must certainly reduce the effects of the coupling on the open ocean.

3.3. Sensitivity of the Ocean–Ice Sheet Coupling Response to Sudden Warming

The ocean and sea ice changes in WARM compared to REF are similar to the ones documented in Van Achter et al. [45]. Here, we focus on the comparison between WARM and WARM3mcpl. The area-averaged basal melt rates for both WARM and WARM3mcpl are presented in Figure 13. The differences in basal melt rate are lower between WARM3mcpl and WARM than in the REF experiments. WARM3mcpl has a lower mean basal melt rate in the TIS cavity by 3.6% and a higher one by 11.6% in the MUIS cavity compared to WARM.

REF3mcpl and WARM3mcpl do not start from a balanced coupled state but from the same initial state as REF and WARM for an easier comparison between them. As a consequence, a small drift is seen in REF3mcpl compared to REF due to the coupling. The difference between WARM3mcpl and WARM includes also such a drift, in addition to the different responses of the coupling to the warming. If we assume that the drift in basal melt rate is independent of the initial conditions, it can be evaluated from the difference between REF3mcpl and REF (+6.7% in TIS and +14.2% in MUIS). We can then subtract this drift from the basal melt rate differences between WARM3mcpl and WARM to isolate the changes in the basal melt rate response to the warming associated with the coupling. This leads to a basal melt rate response to the warming of −10% in TIS and −3% in MUIS. Under this hypothesis, the ice sheet coupling reduces the effect of the ocean warming on the basal melt rates in both cavities and can be seen as a negative feedback.

The spatial distributions of the basal melt rate in WARM and WARM3mcpl (see Figure 14) show at the first order the same changes as REF and REF3mcpl (which supports our condition-independent drift hypothesis). The high values of basal melt rate, found in the middle and center of the TIS cavity in WARM, are limited to the bottom of the cavity in WARM3mcpl. Furthermore, the basal melt rates are increased in the western side of the TIS cavity and decreased in its eastern side. Moreover, as in REF3mcpl, the melt rates in WARM3mcpl are smoothed out and differences between WARM3mcpl and WARM can reach up to 30 m year${}^{-1}$.

The changes in ice shelf thickness throughout the WARM3mcpl simulation are shown in Figure 15. The MUIS cavity exhibits the same changes in ice shelf thickness as in REF3mcpl (see Figure 7a), which could be explained by the limited effect of the warmer oceanic conditions in the MUIS cavity [45]. In the TIS cavity, most of the areas of ice thickening are the same as in REF3mcpl, with the same magnitude of change (up to 150 m of ice), except the one at the bottom of the cavity, which is no longer thickening in WARM3mcpl. On the other hand, the areas of ice thinning in WARM3mcpl in TIS are almost identical to those in REF3mcpl but experience larger ice losses (up to 350 m of ice loss). Cross-sections of the ice shelf base elevation for both TIS and MUIS are shown in Figure A2b,d (in Appendix A). The changes in ice shelf thickness for the TIS increase the slope of the base of the ice shelf, from the grounding line until 67.1${}^{\circ }$ S. In the middle of the TIS cavity, the ice shelf base slope is strongly decreased, which may explain the lower basal melt rates (in addition to the thinned ice shelf). The year-to-year changes in ice shelf base elevation indicate that the ice thinning of the TIS has not yet reached its equilibrium at the end of WARM3mcpl, as the the ice shelf thickness still shows changes after 20 years of simulation. On the other hand, the ice thickening near the cavity front seems to have reached its equilibrium, as years 6 to 20 are overlapping. As in REF3mcpl, the changes in ice shelf thickness in MUIS are small.

The changes in grounding line position after 20 years of simulation exhibit a similar behaviour in WARM3mcpl as in REF3mcpl (see Figure 15b), except on the eastern side of the TIS cavity, where the grounding line has retreated. This retreat, which occurs after 19 years of coupling, is found around the same location where Pelle et al. [25] simulated grounding line retreat in the TIS cavity under warm climate conditions.

The negative feedback associated with the ice sheet coupling, which reduces the effect of the ocean warming on the basal melt rate, can be explained by the changes in basal melt rate present in the center of the cavity. In this zone, where high melt rates are observed in WARM, melt rates are decreased in WARM3mcpl. As illustrated in Figure A8 (in Appendix A), at these locations, the ice shelf starts to thin at the beginning of WARM3mcpl (from 800 to 500 m of ice) but stabilises after 3 years. As described in Seroussi et al. [29], a thinner ice shelf decreases the pressure at the ice–ocean interface, which results in a rise in the in situ freezing temperature that, at some point, becomes larger than the ocean temperature and decreases the basal melt rates (the changes in the in situ ocean salinity are too small to induce a significant effect on the freezing point temperature in our experiments).

4. Discussion and Conclusions

The first aim of this study was to analyse the impact of the ocean–ice sheet coupling in the Totten Glacier area. To reach this goal, we coupled together high-resolution versions of the ocean–sea ice model NEMO and the ice sheet model BISICLES configured to this region and ran two simulations over the 1995–2014 period: one with NEMO in stand-alone mode and another one with the coupled NEMO-BISICLES model. These simulations revealed that the basal melt rates are increased with the coupling by 6.7% in the TIS cavity and by 14.2% in the MUIS cavity. In addition, the coupling leads to the smoothing of the basal melt rate’s spatial distribution. The changes in basal melt rate due to the coupling show a strong linear relationship with the changes in ocean current speed, with faster ocean currents generating more basal melt and slower ocean currents inducing less melt. The changes in basal melt rate are not linearly related to the changes in ocean temperature. The increased basal melt rates in the TIS are largely caused by increased ocean current speeds at the ice–ocean interface, both in the center and near the bottom of the cavity. On the other hand, the decreased basal melt rates in the TIS cavity, mostly near the cavity front, originate from both lower ocean temperatures and slower ocean currents. The relationship between changes in ocean current speed at the ice–ocean interface and the changes in ice shelf thickness is not clear, as faster ocean currents are found in both areas of thicker and thinner ice shelves. However, increased ocean velocities are often found in areas of high basal melt rates and a thinned ice shelf. This is explained by the thicker water column and by the steeper slope at the ice shelf base, allowing faster ocean currents and enhancing the basal melting (positive feedback).

A nonlinear link is found between the total basal melt rate and changes in ice shelf thickness, with a large number of positive values of the total basal melt rate related to a thicker ice shelf. Nevertheless, most areas with a total basal melt rate higher than 300 m are related to a thinner ice shelf. This suggests that the changes in ice shelf thickness are not a linear response to the melt rate but rather a more complex response, driven by both the basal melt rates and by the dynamics of the ice sheet.

An effect of the ocean–ice sheet coupling outside of the ice shelf cavities is detected in the ocean properties over the continental shelf, with slightly cooler and fresher water masses near the cavities’ outflow. There is also a fast ice frequency increase (up to 2 more months of fast ice per year in some areas), but these changes are small and limited to the western side of the main cavities.

The second goal of this study was to investigate the sensitivity of the response of the ocean–ice sheet coupling to sudden warming. To do so, we ran two additional simulations with EC-Earth3 anomalies (2081–2100) added to the oceanic and atmospheric forcings of the NEMO model. The area-averaged basal melt rates in the simulation conducted with the coupled model are decreased by 3.6% in the TIS cavity and increased by 11.6% in the MUIS cavity compared to the simulation carried out with NEMO in stand-alone mode. Considering that the ice shelf melt rate drift due to the coupling in the 2081–2100 simulation is the same as the one in the 1995–2014 simulation, we can remove the effect of this drift and obtain the response of the ocean–ice sheet coupling due to the ocean warming, which corresponds to a 10% and 3% basal melt rate decrease in the TIS and MUIS cavities, respectively. This suggests that the coupling tends to attenuate the basal melt rates under warmer ocean conditions. This negative feedback is explained by the ice shelf thinning that raises the in situ freezing temperature at the ice–ocean interface and decreases the basal melt rates in the ocean–ice sheet coupled simulation compared to the ocean stand-alone run.

One limitation of our work is the short period of our simulations. As described in Goldberg et al. [16], the ocean circulation usually takes a couple of years to adjust to changes in ice shelf cavity and the glacier flow takes decades to adapt to changes in melt rates. Thus, longer simulations, as in Siahaan et al. [80] or Timmermann and Goeller [24], would allow us to better understand which part of the ice shelf thickness change is due to the initialisation or to the dynamical component of the ice sheet model, which has a slower response than the ocean. Another limitation of our model is its difficulty in representing the grounding line movements, which are underestimated in REF3mcpl in the main ice shelf cavities. Additional experiments should be conducted to study how the grounding line retreats could be better represented. Another perspective is to analyse the effect of tidal representation and high horizontal resolution on ice sheet–ocean coupling, which would help in determining the model developments that require prioritisation for ice sheet–ocean coupling in GCMs.

Although the effect of the coupling on the basal melt rate is smaller than the basal melt rate’s interannual variability, this study showed the importance of ocean–ice sheet coupling, even over a short time period such as 20 years. The coupling impact outside of the ice shelf cavities is small, but other larger ice shelves may have a greater influence on the ocean properties in coupled mode. Finally, our simulations with the coupled model indicate that ocean stand-alone simulations with fixed ice shelf geometry may overestimate the future changes in basal melt rates, at least in some regions. This highlights the importance of using coupled models for century-scale simulations of the Antarctic coastal climate or for global sea level projections.