The Framework For Ice Sheet–Ocean Coupling (FISOC) V1.1

A number of important questions concern processes at the margins of ice sheets where multiple components of the Earth system, most crucially ice sheets and oceans, interact. Such processes include thermodynamic interaction at the ice–ocean interface, the impact of meltwater on ice shelf cavity circulation, the impact of basal melting of ice shelves on grounded ice dynamics and ocean controls on iceberg calving. These include fundamentally coupled processes in which feedback mechanisms between ice and ocean play an important role. Some of these mechanisms have major implications for humanity, most notably the impact of retreating marine ice sheets on the global sea level. In order to better quantify these mechanisms using computer models, feedbacks need to be incorporated into the modelling system. To achieve this, ocean and ice dynamic models must be coupled, allowing runtime information sharing between components. We have developed a flexible coupling framework based on existing Earth system coupling technologies. The open-source Framework for Ice Sheet–Ocean Coupling (FISOC) provides a modular approach to coupling, facilitating switching between different ice dynamic and ocean components. FISOC allows fully synchronous coupling, in which both ice and ocean run on the same time step, or semi-synchronous coupling in which the ice dynamic model uses a longer time step. Multiple regridding options are available, and there are multiple methods for coupling the sub-ice-shelf cavity geometry. Thermodynamic coupling may also be activated. We present idealized simulations using FISOC with a Stokes flow ice dynamic model coupled to a regional ocean model. We demonstrate the modularity of FISOC by switching between two different regional ocean models and presenting outputs for both. We demonstrate conservation of mass and other verification steps during evolution of an idealized coupled ice–ocean system, both with and without grounding line movement.

1 Introduction

The Antarctic and Greenland ice sheets have the potential to provide the greatest contributions to global sea level rise on century timescales (Church et al., 2013; Moore et al., 2013), with the greatest uncertainty in projections being due to marine ice sheet instability (MISI; Mercer, 1978; Schoof, 2007; Robel et al., 2019). Ice dynamic behaviour is strongly sensitive to ocean currents, in particular the transport of warmer waters across the continental shelf, causing high basal melt rates under ice shelves (Hellmer et al., 2012; Thoma et al., 2015). For Antarctica's Pine Island Glacier, which is likely undergoing unstable retreat due to MISI, ocean-induced basal melting has been established as a trigger for MISI through both observational evidence (Christianson et al., 2016) and model studies (Gladstone et al., 2012; De Rydt et al., 2014; Favier et al., 2014). While MISI is fundamentally a geometrically controlled phenomenon, its onset and the resulting rate of ice mass loss are strongly dependent on tight coupling between ice dynamic behaviour and ocean processes. Importantly, ocean-driven basal melt rates respond to the evolving geometry of ice shelf cavities (Timmermann and Goeller, 2017; Mueller et al., 2018), and the grounded-ice dynamic behaviour responds to the evolving basal melt rates through their impact on the buttressing force provided by ice shelves to the grounded ice. While most model-based ice sheet studies use relatively simple parameterizations for calculating basal melt rates beneath ice shelves, recent studies have highlighted limitations of this approach (De Rydt and Gudmundsson, 2016; Favier et al., 2019). In particular, melt parameterizations as a function of depth or thermal driving do not impose conservation of heat in the system, and none of the parameterizations fully capture the impact of evolving ice geometry on cavity circulation.

Several projects to couple ice sheet and ocean models are underway, and most (including the current study) will contribute to the Marine Ice Sheet–Ocean Model Intercomparison Project first phase (MISOMIP1) and its child projects: the Marine Ice Sheet Model Intercomparison Project third phase (MISMIP+) and the Ice Shelf Ocean Model Intercomparison Project second phase (ISOMIP+; Asay-Davis et al., 2016).

Coupling projects take different approaches to handling the different timescales of ice and ocean processes. An ice sheet flow line model coupled to a five-box ocean model allows large ensemble simulations to be carried out but is limited in terms of implementation of physical processes (Gladstone et al., 2012). A temporally synchronous approach allows the cavity geometry to evolve on the ocean time step as a function of the melt rates calculated by the ocean model and the ice dynamics calculated by the ice model (Goldberg et al., 2018). Asynchronous approaches incorporate a longer time step for ice than ocean and sometimes involve coupling through file exchange and with restarts for the ocean model (Seroussi et al., 2017; De Rydt and Gudmundsson, 2016; Thoma et al., 2015).

Here, we present a new, flexible Framework for Ice Sheet–Ocean Coupling (FISOC; Sect. 2). FISOC allows runtime coupling in which ice and ocean components are compiled as runtime libraries and run through one executable. FISOC provides the user a selection of synchronicity options. Adopting Earth System Modeling Framework terminology (ESMF; Sect. 2), we refer to an ocean model coupled through FISOC as an “ocean component” and an ice sheet or ice dynamic model coupled through FISOC as an “ice component”. We use FISOC to couple two different 3-D ocean models to an ice dynamic model and present idealized simulations demonstrating mass conservation and consistent grounding line behaviour (Sect. 3). FISOC is also currently being used to contribute to ISOMIP+ and MISOMIP1.

2 Methodology

FISOC is an open-source coupling framework built using the ESMF (Hill et al., 2004; Collins et al., 2005). FISOC aims to provide seamless runtime coupling between an existing ice sheet model and an existing ocean model for application to Antarctic ice sheet–ocean systems. In its current form, FISOC assumes that the important ice sheet–ocean interactions occur at the underside of a floating ice shelf and that the lower surface of the ice shelf can be projected on to the horizontal plane.

FISOC aims to provide flexibility and computational efficiency through the following key features.

• Flexible modular architecture (Sect. 2.1) facilitates swapping between different ice components or between different ocean components according to purpose (Sect. 2.2).

• Access to ESMF tools allows multiple regridding and interpolation options, including between regular grids and unstructured meshes (Sect. 2.3).

• Multiple options for handling differing ice and ocean timescales include fully synchronous coupling, passing rates of change, time averaging of variables (Sects. 2.4 and 2.5).

• Flexible runtime control over the exchange of variables allows specific coupling modes to be (de)activated as required, e.g. geometric coupling and thermodynamic coupling.

• Grounding line movement (Sect. 2.8) is implemented using geometry change rates and a modified wet–dry scheme in the ocean component, with multiple options available for updating cavity geometry (Sect. 2.5).

• Flexibility for parallelization options is provided; currently sequential coupling is implemented, but any combination of sequential and concurrent parallelization is possible with minimal coding effort (see also Sect. 2.1.1).

• ESMF compatibility means that FISOC can be embedded within any ESMF-based modelling system, e.g. as a regional model within a global model.

• ESMF compatibility also means that additional ESMF components (e.g. an atmosphere model) could easily be added to the coupled system.

These features are described further in the following sections and in the FISOC manual, which can be found in the FISOC repository (see the code availability section at the end of this paper).

2.1 Software design

While coupled models in Earth system science have been in existence for decades and such coupled models are often viewed as single entities (ocean–atmosphere general circulation models, for example), the field of coupled ice sheet–ocean modelling is relatively young. FISOC is intended as a framework for coupling independent models rather than as a coupled model in itself. Building and running a coupled ice sheet–ocean model is currently more complex than building and running both an ice and an ocean model independently. FISOC aims to minimize the additional complexity.

The ice and ocean components may use their standard runtime input files, and their paths are set in a FISOC runtime configuration file, along with information about time stepping and variables to be exchanged.

FISOC adopts the hierarchical modular structure of the Earth System Modelling Framework. The FISOC code structures are summarized in Fig. 1. A top-level executable is called a FISOC parent module (this could in principle also be embedded within a larger coupled model framework). The parent module coordinates calling of the ice, ocean and regridding components. Regridding is one of the reasons to make use of ESMF, described further in Sect. 2.3. The ice and ocean components are independent models that are not included in the FISOC code repository and compiled as libraries to be called by FISOC during runtime. On each side (ice and ocean) of the coupling is a model-specific wrapper, whose main runtime functions are as follows:

• call the component's initialize, run, and finalize routines as required;

• convert the component's grid or mesh to ESMF format using ESMF data structures;

• read from or write to the component's required state variables, converting between the component's native data structures and ESMF data structures.

Further processing of variables (such as calculating rates of change) is implemented by the ice and ocean generic code modules.

Incorporating a new ice or ocean component into FISOC can be straightforward, depending on the existing level of ESMF compatibility of the new component. Models able to provide mesh information and variables in ESMF data structures can be very easily built in to FISOC. The only coding required for a new component is a new model-specific wrapper in the FISOC repository. Copying an existing wrapper can be a viable starting point.

Figure 1 Overview of FISOC code structures. OM and ISM are short for ocean model and ice sheet model (or component), respectively. ImpSt and ExpSt are short for import state and export state, respectively.

Figure 2 FISOC workflow. The black arrow indicates the direction of simulated time. The yellow arrows indicate the order of events during a FISOC simulation. The green boxes indicate when regridding and passing of variables between components occurs. The length of the blue boxes in the vertical indicates the relative size of time steps and coupling interval (this is illustrative; in practice there will be many more OM time steps per ISM time step and the ISM time step size will usually equal the coupling interval).

Table 1 Ice and ocean components currently coupled through FISOC.

Table 2 Model choices and input parameters used in verification experiment 1 (VE1, Sect. 3.1) and verification experiment 2 (VE2, Sect. 3.2) comprising four simulations in total: VE1_ER, VE1_EF, VE2_ER and VE2_EF. Component abbreviations in these simulation names are E (Elmer/Ice), R (ROMS), and F (FVCOM). Semi-structured refers to a mesh that is in principle unstructured but where in practice a structure can be seen (see Fig. 3 middle and lower panels). STOD stands for source to destination.

Figure 3 Unstructured meshes used in the current study. The first 15 km is shown. The ocean model in the ER simulations uses a structured grid.

2.1.1 Sequential parallelism

FISOC currently adopts a sequential parallelism paradigm. Each component runs on the full set of available persistent execution threads (PETs). PET is an ESMF abstraction catering for multiple parallelism options. FISOC has so far used only the message passing interface (MPI), in which one PET wraps one MPI process.

The sequential workflow is illustrated in Fig. 2. The order of events during time stepping is as follows. The ocean component is called for the full number of ocean time steps required to complete one coupling interval. Ocean outputs are then regridded and passed to the ice component, which also runs for as many time steps as are required to complete one coupling interval. The ice component outputs are then regridded and passed to the ocean component. The ice component time step size is equal to the coupling interval for all simulations in the current study.

The initialization is not shown in Fig. 2, but we note that this is similar to the runtime event order: the ocean component is initialized first, followed by regridding and then the ice component. There are two initialization phases for each component, allowing for the possibility that variables may be needed to be passed from ice component to ocean component in order to finalize initialization.

This ordering of events imposes a lag in the system. While the ice component receives ocean variables for the current coupling interval, the ocean component only receives ice variables for the previous coupling interval. This could be reversed (running the ice component before the ocean component) or could be modified such that both components receive variables from the other component for the previous coupling interval.

While FISOC implements sequential parallelism, ESMF also supports concurrent parallelism. Concurrent parallelism allows different components to run at the same time on different subsets of the available PETs. This approach is beneficial when different components have very different computational costs and parallel scaling: a cheap component that scales poorly is more effectively run on a subset of the available PETs, and concurrent parallelism allows this to be implemented more efficiently than sequential parallelism. This could easily be implemented in FISOC if it becomes necessary, as the components, which utilize MPI, are assigned a distinct MPI communicator during initialization. This communicator could be made to represent a subset of the available PETs. In principle, concurrent parallelism also offers sub-time-step coupling: it is possible to exchange variables between components during convergence of numerical schemes. Such coupling is unlikely to be implemented within FISOC as the timescales for ice and ocean components are so different. While sequential coupling imposes a lag between components (described above), concurrent coupling implemented in FISOC would impose a lag in both components: exchange of variables in both directions would occur at the end of the coupling interval.

2.1.2 Error handling

The ESMF adopts a defensive strategy to error handling: all errors are logged and passed back up the call stack. The calling routine has the option of attempting to continue running in the event of errors occurring. As the call structure between FISOC and ESMF is one way (FISOC routines may call ESMF routines but not vice versa), all such errors are eventually returned to FISOC.

FISOC adopts a fail-fast approach. Errors are generally considered to be fatal, in which case FISOC will log error information and finalize both ice and ocean components and ESMF. FISOC also aims to provide consistency checks, most of which are considered fatal if not passed. For example, ice and ocean input files might both contain time-stepping information, potentially duplicating information in the FISOC runtime configuration file, and these can be checked for consistency in the model-specific wrappers. The general intention is to stop running if something unexpected happens and provide a meaningful message to the user about why.

There are a few cases where ESMF errors can be handled at runtime. Details can be found in the FISOC manual, which can be accessed from the FISOC repository (see the code availability section at the end of this paper).

2.2 Components

FISOC is designed to facilitate swapping between different ocean or ice components. Currently two different ocean components and one ice component are available through FISOC. Table 1 summarizes components currently coupled into FISOC. In some cases, a non-standard build of the component is required for FISOC compatibility, and these cases are described in the FISOC manual, which can be obtained through the FISOC repository (Sect. 2.1).

The ice component Elmer/Ice (Gagliardini et al., 2013) is a powerful, flexible, state-of-the-art ice dynamic model.

The Regional Ocean Modeling System (ROMS; Shchepetkin and McWilliams, 2005) is a 3-D terrain-following, sigma coordinate ocean model that has already been adapted to use in ice shelf cavities (Galton-Fenzi et al., 2012). The module for ice shelf cavities implemented in the Finite Volume Community Ocean Model (FVCOM; Chen et al., 2003) provides non-hydrostatic options and a horizontally unstructured mesh that lends itself to refinement and may be more suited to small-scale processes such as ice shelf channels (Zhou and Hattermann, 2020).

2.3 Regridding

As stated above, FISOC provides coupling on a horizontal plane onto which the lower surface of an ice shelf can be projected. It is this plane on which ice and ocean properties are exchanged through the FISOC. Adapting the FISOC code to handle a vertical ice cliff is expected to be straightforward and would be desirable for application to the Greenland ice sheet. More complex 3-D ice–ocean interface geometries are challenging not only for FISOC but also for the current generation of ice sheet and ocean models.

FISOC has access to all the runtime regridding options provided by ESMF. These include nearest-neighbour options, conservative options, patch recovery and bilinear regridding. These options are available for structured grids and unstructured meshes. FISOC requires that both ice and ocean components define their grid or mesh on the same coordinate system and that both components use the same projection. All FISOC simulations to date have used a Cartesian coordinate system (i.e. all components have so far used Cartesian coordinates).

Our current FISOC setup does not meet the requirements for all forms of ESMF regridding. Specifically, the conservative methods, when an unstructured mesh is involved, require that field values are defined on elements and not on nodes. Elmer, by default, provides field values on nodes but can also provide element-wide values or values on integration points within elements. We will need to either map nodal values to element values or utilize element-type variables in order to use conservative regridding, and this is intended as a future development.

When using FISOC to couple Elmer/Ice to ROMS, the ROMS grid extends beyond the Elmer/Ice mesh. This is due to ROMS using a staggered grid (Arakawa C-grid) and ghost cells extending beyond the active domain. This necessitates the use of extrapolation. ESMF regridding methods provide options for extrapolation, which are used here. Simulations in the current study use either nearest “source to destination” (STOD, a form of nearest neighbour) regridding or use bilinear interpolation (in which case nearest STOD is used only for destination points that lie outside the source domain).

We use subscripts with square brackets, _[X], where X is either O (ocean component) or I (ice component), to denote a variable that exists in both ice and ocean components with the same physical meaning but potentially different values due to being represented on different grids or meshes.

2.4 Coupling timescales

The timescales for sub-shelf cavity circulation behaviour are in general much shorter than the timescales for ice flow and geometry evolution (typically minutes to days instead of years to centuries). Typical time step sizes are correspondingly smaller for ocean models (seconds to minutes) than for ice sheet models (days to months). A single ice sheet model time step, if the Stokes equations are solved in full, will typically require orders of magnitude more computational time than a single ocean time step. Due to the combination of these two reasons the ice and ocean components of FISOC will in general use different time steps, with the ice time step size being much larger. We define relevant terminology for coupling timescales below.

• Fully synchronous coupling. The ice and ocean components have the same time step size, and they exchange variables every time step.

• Semi-synchronous coupling. The ice component has a larger time step than the ocean component, but the ocean component's cavity geometry and grounding line position are allowed to evolve on the ocean time step (e.g. by using ice velocities from a previous ice time step or rates of change based on the two most recent time steps).

• Asynchronous coupling. The ice component has a larger time step than the ocean component. Cavity geometry is updated on the ice component time step or less frequently.

• Coupling interval. The time interval at which the ice and ocean components exchange variables.

In the current study, FISOC sets the coupling interval as equal to the ice component time step size. This is an exact multiple of the ocean model time step size. More generally (for potential future experiments), FISOC calls each component for a fixed time period and allows the component to determine its own time stepping within that period. In principle, adaptive time stepping could be implemented within this framework as long as each component runs for the required amount of simulated time. FISOC does not currently provide an option to vary the coupling interval during a simulation, but this could be implemented if needed.

FISOC is flexible with regard to time processing of ocean or ice variables. It is possible to cumulate variables, calculate averages, or use snapshots. In the current study, the ocean components (both ROMS and FVCOM) calculate averaged basal melt rates over the coupling interval and pass these averages through FISOC to the ice component. In the current study, as the ice component time step size is equal to the coupling interval for all simulations, no time processing of ice component variables is needed.

In principle, FISOC supports all three synchronicity options, though fully synchronous coupling is not practical to achieve when solving the Stokes equations for the ice component. The experiments carried out for this paper use semi-synchronous coupling with cavity geometry evolution as described in Sect. 2.5.

Goldberg et al. (2018) and Snow et al. (2017) implement fully synchronous coupling, whereas Seroussi et al. (2017) and Favier et al. (2019) implement asynchronous coupling with ocean restarts.

2.5 Handling cavity evolution

The evolution of cavity geometry under the ice shelf, defined by a reference ice draft, z_d (positive upward), and grounding line location, is calculated by the ice component forced by the melt rates passed from the ocean component. We refer to z_d as a “reference” ice draft because the ocean component may further modify the ice draft according to the dynamic pressure field. The ocean component's “free surface” variable, ζ, represents the height of the upper surface of the ocean domain relative to a mean sea level for the open ocean. Under the ice shelf, ζ represents the deviation of the upper surface of the ocean domain relative to the reference ice draft z_d (similar to Goldberg et al., 2018). To summarize the meaning of the key variables, z_d[I] is the reference ice draft computed by the ice component, z_d[O] is the same but regridded for the ocean component and (z_d[O]+ζ) is the actual ice draft according to the ocean component.

Given the potential for non-synchronicity of the ice and ocean component time stepping, several methods are implemented in FISOC for the ocean to update its representation of z_d. All the processing options described below are applied on the ocean grid after the ice component representation of ice geometry has been regridded (i.e. z_d[I] regridded to z_d[O]).

Most recent ice. The simplest option is that the ocean component uses the ice draft directly from the most recent ice component time step. If fully synchronous coupling is used, this option should be chosen. The main disadvantage of this approach for semi-synchronous or asynchronous coupling is that due to the much longer time step of the ice component the ocean component will experience large, occasional changes in ice draft instead of smoothly evolving ice draft. This could be both physically unrealistic and potentially numerically challenging for the ocean component.

Rate. The vertical rate of change of ice draft, $\frac{\mathrm{d}{z}_{\mathrm{d}}}{\mathrm{d}t}$, is calculated by FISOC after each ice component time step using the two most recent ice component time steps. If we assume that the ice component completes a time step at time t, the rate at this time is given by

$$\begin{array}{}\text{(1)}& {\displaystyle }{\displaystyle \frac{\mathrm{d}{z}_{\mathrm{d}[O,t]}}{\mathrm{d}t}}={\displaystyle \frac{z_{\mathrm{d}[I,t]}-z_{\mathrm{d}[I,t-\mathrm{\Delta }{t}_I]}}{\mathrm{\Delta }{t}_I}},\end{array}$$

where z_d[O,t] is the ocean component's reference ice draft at time t, z_d[I,t] is the ice component's reference ice draft at time t, $z_{\mathrm{d}[I,t-\mathrm{\Delta }{t}_I]}$ is the ice component's reference ice draft at time t−Δt_I and Δt_I is the ice component time step size. This rate of change is used by the ocean component to update the cavity geometry until the next ice component time step completes. In this sense the ocean component lags the ice component as mentioned in Sect. 2.1.1. This approach provides temporally smooth changes to the ocean representation of the ice draft but has the potential for the ice and ocean representations to diverge over time as a result of regridding artefacts.

Corrected rate. This is the same as above, except that a drift correction is applied to ensure that ice and ocean representations of cavity geometry do not diverge:

$$\begin{array}{}\text{(2)}& {\displaystyle }{\displaystyle \frac{\mathrm{d}{z}_{\mathrm{d}[O,t]}}{\mathrm{d}t}}={\displaystyle \frac{z_{\mathrm{d}[I,t]}-z_{\mathrm{d}[I,t-\mathrm{\Delta }{t}_I]}+f_{\text{cav}}\left(z_{\mathrm{d}[I,t]}-z_{\mathrm{d}[O,t]}\right)}{\mathrm{\Delta }{t}_I}},\end{array}$$

where f_cav is a cavity correction factor between 0 and 1. Equation (2) is applied at coupling time steps, and the calculated rate of cavity change is then held constant during ocean component evolution until the coupling interval completes. Conceptually, this option prioritizes ice–ocean geometry consistency over mass conservation.

Linear interpolation. The ocean representation of the ice draft is given by temporal linear interpolation between the two most recent ice sheet time steps. This imposes additional lag of the ocean component behind the ice component.

The above options are all implemented in FISOC, but only the “rate” and “corrected rate” approaches are used in the current study.

The cavity geometry may be initialized independently by ice and ocean components. In this case, the user must ensure consistency. It is also possible for the cavity geometry from the ice component to be imposed on the ocean component during FISOC initialization. This ensures consistency.

Handling cavity evolution is a little more complicated in the case of an evolving grounding line, as discussed in Sect. 2.8 below.

2.6 Thermodynamics at the ice–ocean interface

Exchange of heat at the ice–ocean interface is handled within the ocean model. Like many ocean models, FVCOM and ROMS adopt the three-equation formulation for thermodynamic exchange (Hellmer and Olbers, 1989; Holland and Jenkins, 1999; Jenkins et al., 2010). This parameterization assumes that the interface is at the in situ pressure freezing point and that there is a heat balance and salt balance at the interface. Both ROMS and FVCOM assume constant turbulent transfer coefficients for scaling the heat and salt fluxes through the interface, with thermal and saline exchange velocities calculated as the product of these coefficients with friction velocity. Further details of the ROMS- and FVCOM-specific implementations of the three-equation formulation are given by Galton-Fenzi et al. (2012) and Zhou and Hattermann (2020), respectively. An ablation or melt rate is calculated for each ocean model grid cell, which is then passed to FISOC as a boundary condition for the lower surface of the ice model at the coupling time interval.

Internally, both ocean models account for the thermodynamic effect of basal melting by imposing virtual heat and salt fluxes within a fixed geometry at each ocean model time step to mimic the effects of basal melting, rather than employing an explicit volume flux at the ice–ocean interface. Independent of this, a geometry change is passed back from the ice model through FISOC at after each coupling interval (including the effect of melting and freezing, as well as any ice dynamical response), which is used to update the ocean component cavity shape (Sect. 2.5).

For some applications, conductive heat fluxes into the ice shelf due to vertical temperature gradients in the ice at the ice–ocean interface are required by the three-equation parameterization to calculate the flux balance at the ice ocean interface. While ice–ocean thermodynamic parameterizations in ocean-only models must make an assumption about this temperature gradient, FISOC can pass the temperature gradient from the ice component directly to the ocean component. This feature is not demonstrated in the current study but will be properly tested in future studies.

Non-zero basal melt rates may be calculated by the ocean component in regions that are defined as grounded by the ice model. This could occur due to isolated patches of ungrounding upstream of the grounding line or discrepancies between the ice and ocean component's representation of the grounded region. Basal melt rates are masked using the ice component's grounded mask before being applied within the ice component. This has the potential to impact on mass conservation in the coupled system. Future studies utilizing conservative regridding will ensure that passing masked field variables between components remains conservative.

2.7 Interface pressure

Aside from the geometry evolution, an ocean boundary condition for pressure at the ice–ocean interface, P_interface, must be provided to the ocean component. FISOC can pass pressure directly from ice to ocean components. However, using actual ice overburden directly as an upper ocean boundary condition results in higher horizontal pressure gradients at the grounding line (for dry cells, see Sect. 2.8) than ocean models can typically handle (Goldberg et al., 2018). In the current study, the ocean component uses the reference ice draft (see Sect. 2.5) to estimate floatation pressure. ROMS assumes a constant reference ocean density:

$$\begin{array}{}\text{(3)}& {\displaystyle }{\displaystyle }{P}_{\text{interface}}=-g{\mathit{\rho }}_{\text{or}}{z}_{\mathrm{d}\left[O\right]},\end{array}$$

where g is acceleration due to gravity, ρ_or is a reference ocean density and z_d[O] is the ocean representation of ice draft (positive upward). For the current study, all simulations with ROMS use ${\mathit{\rho }}_{\text{or}}=\mathrm{1027}\mathrm{kg}{\mathrm{m}}^{-\mathrm{3}}$. FVCOM assumes a constant vertical ocean density gradient following Dinniman et al. (2007):

$$\begin{array}{}\text{(4)}& {\displaystyle }{\displaystyle }{P}_{\text{interface}}=-g\left({\mathit{\rho }}_{o\mathrm{1}}+\mathrm{0.5}{\displaystyle \frac{\mathrm{d}{\mathit{\rho }}_o}{\mathrm{d}z}}{z}_{\mathrm{d}\left[O\right]}\right)z_{\mathrm{d}\left[O\right]},\end{array}$$

where ρ_o is ocean water density, ρ_o1 is ocean water density of the top ocean layer, and the vertical ocean water density gradient, $\frac{\mathrm{d}{\mathit{\rho }}_o}{\mathrm{d}z}$, is set to $\mathrm{8.3927}\times {\mathrm{10}}^{-\mathrm{4}}\mathrm{kg}{\mathrm{m}}^{-\mathrm{4}}$.

2.8 Grounding line evolution

Grounding line movement in FISOC requires that both ice and ocean components support it. Numerical convergence issues place constraints in terms of mesh resolution for representing grounding line movement in ice sheet models (Vieli and Payne, 2005; Pattyn et al., 2006; Gladstone et al., 2010a, b, 2017; Cornford et al., 2013). While FISOC allows ice draft to be passed to the ocean component (Sect. 2.5), FISOC does not impose the ice component grounding line position on the ocean component. Instead, the ocean component uses the evolving cavity geometry to evolve the grounding line.

A recent ice–ocean coupling study (Goldberg et al., 2018) used a “thin film” approach to allow grounding line movement. A thin passive water layer is allowed to exist under the grounded ice, and an activation criterion is imposed to allow the layer to inflate to represent grounding line retreat. The current study takes a conceptually similar approach, modifying the existing wetting and drying schemes independently in both ROMS (Warner et al., 2013) and FVCOM. “Dry” cells are used for the passive water column under grounded ice, and “wet” cells are used for the active water column under floating ice or the open ocean. The wet–dry mask is two dimensional, so while it is conventional to talk about dry or wet cells, this actually refers to dry or wet columns. The grounding line evolves in the two horizontal dimensions and is represented in the ocean component as the vertical surface between dry and wet columns.

The original criterion in both ROMS and FVCOM for a cell to remain dry is given by

$$\begin{array}{}\text{(5)}& {\displaystyle }{\displaystyle }\mathit{\zeta }-z_b<D_{\text{crit}},\end{array}$$

where z_b is the bottom boundary depth (bathymetry, or bedrock depth, positive upward), and D_crit is a critical water column thickness for wet or dry activation. D_crit is a parameter to be set by the user (typical values lie between 1 to 20 m). Thus, cells with a water column thickness less than D_crit are designated as dry. Flux of water into dry cells is allowed, but flux of water out of dry cells is prevented.

The FVCOM criterion for an element to be dry has been modified for the presence of a marine ice sheet and ice shelf system as follows:

$$\begin{array}{}\text{(6)}& {\displaystyle }{\displaystyle }-z_b+z_{\mathrm{d}\left[O\right]}<D_{\text{crit}}.\end{array}$$

This is a purely geometric criterion based entirely on the geometry determined by the ice component. The ROMS criterion for a cell to be dry has been modified for the presence of a marine ice sheet and ice shelf system as follows:

$$\begin{array}{}\text{(7)}& {\displaystyle }{\displaystyle }\mathit{\zeta }-z_b-(z_{s\left[O\right]}-z_{\mathrm{d}\left[O\right]}-D_{\text{crit}}){\displaystyle \frac{{\mathit{\rho }}_{\mathrm{i}}}{{\mathit{\rho }}_{\text{or}}}}\le \mathrm{0},\end{array}$$

where z_s[O] is the ocean representation of ice sheet and ice shelf upper surface height. z_s[O] is needed in this equation because the floatation assumption cannot be made for grounded ice. This equation essentially compares ζ against the height above buoyancy of the grounded ice. In other words, if the dynamic variations in ocean pressure are sufficient to overcome the higher ice pressure due to the positive height above buoyancy, the cell can become ungrounded. The conceptual difference between the FVCOM and ROMS wetting criteria is that ROMS allows dynamic ocean pressure variations to make minor grounding line adjustments relative to the grounding line determined by the ice geometry, whereas FVCOM uses just the ice geometry to determine grounding line position.

FISOC allows the ice component to pass any geometry variables to the ocean, such as ice draft, ice thickness, upper surface elevation or the rates of change of any of these variables. In the event that geometry variables other than z_d are passed to the ocean, the same processing method is used as for z_d, as described in Sect. 2.5. In the current study, $\frac{\mathrm{d}{z}_{\mathrm{d}}}{\mathrm{d}t}$ is passed to the ocean component, and in one case both $\frac{\mathrm{d}{z}_{\mathrm{d}}}{\mathrm{d}t}$ and $\frac{\mathrm{d}{z}_s}{\mathrm{d}t}$ are passed (details in Sect. 3). When $\frac{\mathrm{d}{z}_s}{\mathrm{d}t}$ is passed, $\frac{\mathrm{d}{z}_s}{\mathrm{d}t}$ is processed the same way as $\frac{\mathrm{d}{z}_{\mathrm{d}}}{\mathrm{d}t}$. If the grounding line problem is solved and if z_d is processed for passing to the ocean using the corrected rate method, Eq. (2) is modified to account for the dry water column thickness, which is initialized to D_crit. The correction term changes from $f_{\text{cav}}\left(z_{\mathrm{d}[I,t]}-z_{\mathrm{d}[O,t]}\right)$ to $f_{\text{cav}}\left(max(z_{\mathrm{d}[I,t]},z_b+D_{\text{crit}})-z_{\mathrm{d}[O,t]}\right)$.

There are no connectivity restrictions on wetting and drying for either of the ocean components in the current study. This means that it is possible for individual cells or regions containing multiple cells that are upstream of the grounding line to become wet (i.e. to unground). This occurs on small spatial and temporal scales in ROMS (individual cells a short distance upstream of the grounding line sometimes become temporarily wet) but not at all in FVCOM (likely due to choice of wetting criterion).

3 Verification experiment design

Simulations are carried out on idealized domains as a proof of concept to demonstrate the coupling rather than to address scientific questions. Verification experiment 1 (VE1) aims to assess whether the coupled system conserves mass. Verification experiment 2 (VE2) aims to assess whether the ocean and ice representations of grounding line evolution are consistent.

3.1 Verification experiment 1: floating adjustment

Verification experiment 1 (VE1) is a simple experiment in which a linearly sloping ice shelf is allowed to adjust toward steady state. The experiment is not run long enough to attain steady state but is long enough to demonstrate the evolution of the coupled system; see Table 2 for run length and a summary of other model choices and parameter values used in VE1.

All ice and ocean vertical side boundaries are closed: there is no flow in or out of the domain. There is mass exchange between the ice and ocean (and therefore also heat exchange). The coupling centres on the evolution of ice geometry: the ocean component passes an ice shelf basal melt rate to the ice component, and the ice component passes a rate of change of ice draft to the ocean component.

We expect adjustment toward a uniform-thickness ice shelf to occur by the following two mechanisms.

1. Ice dynamics. The gravitational driving force will tend to cause flow from thicker to thinner regions.

2. Melt and freeze. The greater pressures at greater depth should result in higher melt rates, with the potential for refreezing under thinner regions.

3.1.1 Domain size and meshes

The domain is 30 km across the expected direction of ice flow (y direction) by 100 km along the flow (x direction) for simulation VE1_ER. However, ocean component FVCOM (used in VE1_EF) uses a semi-structured (in principle unstructured but in practice exhibiting some structure) mesh with the dimensions 31 km by 99 km. This results from an auto-generated mesh method using a uniform resolution of 2 km for its triangular elements. FISOC does not in general require that ice and ocean component domains precisely overlap. Indeed the region of overlap is allowed to be small relative to the domains (for example, an Antarctic ice stream interacts with the ocean only in its floating shelf, and the majority of the catchment may be grounded with no possibility to interact with the ocean for the duration of an intended simulation). However, given that we aim to address mass conservation in the coupled system, we choose to require precise domain match between ice and ocean components for the current study. Therefore, for simulations presented in the current study, the ice component has a slightly different domain when coupled to ROMS as compared to when coupled to FVCOM. For VE1_EF the ice component runs on an almost identical mesh to the ocean component. The only difference is at two diametrically opposite corners, where FVCOM prefers to maintain element shape but Elmer/Ice prefers to maintain a strictly rectangular domain (in order to facilitate imposition of consistent boundary conditions at the corners of the domain). These mesh differences are visually summarized in Fig. 3.

3.1.2 Ice component setup

The initial geometry is of an ice shelf at floatation (i.e. hydrostatic equilibrium). The initial ice draft is given in metres by

$$\begin{array}{}\text{(8)}& {\displaystyle }{\displaystyle }{z}_{\mathrm{d}}=-\mathrm{450}+\mathrm{400}\left({\displaystyle \frac{x}{\mathrm{100}\mathrm{000}}}\right),\end{array}$$

where x is distance in metres along the domain. The initial geometry does not vary across the ice flow (y direction). Ice and ocean water densities used in the ice component are ${\mathit{\rho }}_{\mathrm{i}}=\mathrm{910}\mathrm{kg}{\mathrm{m}}^{-\mathrm{3}}$ and ${\mathit{\rho }}_{\text{or}}=\mathrm{1027}\mathrm{kg}{\mathrm{m}}^{-\mathrm{3}}$, respectively. These densities, along with the floatation assumption, determine the ice upper surface.

The pressure acting on the underside of the ice shelf is given by Eq. (3).

Temperature in the ice component is constant through space and time at −5 ^∘C.

VE1 includes ice flow and geometry evolution solving the Stokes equations directly. Glen's power law rheology with n=3 is implemented (Glen, 1952; Gagliardini et al., 2013).

Zero accumulation is prescribed at the upper ice surface. The melt rate from the ocean component is applied at the lower surface. Flow through the vertical side boundaries is not allowed.

Elmer/Ice-specific details. The Stokes equations are solved within Elmer/Ice (Gagliardini et al., 2013). A 2-D horizontal mesh of triangles with an approximate element size of 1 km (VE1_ER) or 2 km (VE1_EF) is extruded in the vertical to give 11 equally spaced terrain-following layers with the bulk element shape being triangular prisms.

3.1.3 Ocean component setup

The ocean bathymetry is set to 500 m throughout the domain. The wet–dry scheme (Sect. 2.8) is not used in this experiment, as the whole domain is ice shelf cavity with no grounded ice. Boundaries are closed and rotation is disabled. Ocean potential temperature is initialized at −1.85 ^∘C, and salinity is initialized at 34.6 on the practical salinity scale. Ice–ocean thermodynamics are captured by means of the three-equation parameterization (Sect. 2.6).

The ocean conditions are chosen to represent a cold cavity ice shelf, such as the Amery Ice Shelf. In this configuration, both basal melting and refreezing can occur.

ROMS-specific details. The horizontal resolution is a constant 1 km. There are 11 vertical layers, with a sigmoidal terrain-following distribution configured to provide increased resolution near the top and bottom surfaces. The ROMS baroclinic (slow) time step size is 200 s, and there are 30 barotropic (fast) time steps for every slow time step. Interior mixing is parameterized with the K-Profile Parameterization (Large et al., 1994). Background vertical mixing coefficients for tracers and momentum are set to constant values of $\mathrm{5.0}\times {\mathrm{10}}^{-\mathrm{5}}$ and $\mathrm{1.0}\times {\mathrm{10}}^{-\mathrm{3}}{\mathrm{m}}^{\mathrm{2}}{\mathrm{s}}^{-\mathrm{1}}$, respectively, while horizontal viscosity and diffusivity are set to 6.0 and 1.0 m² s⁻¹, respectively.

FVCOM-specific details. The horizontal grid resolution is 2 km (defined by the distance between adjacent nodes within a uniform triangular grid), and there are 11 uniformly spaced vertical terrain-following layers. Interior vertical mixing is parameterized using the Mellor and Yamada level 2.5 (Mellor and Yamada, 1982) turbulent closure model (vertical Prandtl number =0.1) together with a constant background viscosity and diffusivity of ${\mathrm{10}}^{-\mathrm{6}}{\mathrm{m}}^{\mathrm{2}}{\mathrm{s}}^{-\mathrm{1}}$. An eddy closure parameterization (Smagorinsky, 1963) is used for the horizontal mixing of momentum (viscosity) and tracers (diffusivity), with both the scaling factor and the Prandtl number being 0.1. Both the barotropic time step and the baroclinic time step sizes are 20 s.

3.1.4 Coupling

The coupling interval is 10 d, the same as the ice component time step size. The cavity update method is rate (Sect. 2.5). For VE1_ER, the regridding method is bilinear using nearest STOD extrapolation for ocean cells that lie outside the ice domain due to grid stagger. For VE1_EF, nearest STOD regridding is used, which results in a one-to-one mapping between ice and ocean nodes due to the meshes being nearly identical (Sect. 3.1.1). There is no grounding line in this experiment.

3.2 Verification experiment 2: grounding line evolution

Verification experiment 2 (VE2) is a modified version of VE1 but with part of the region grounded and a net ice flow through the domain allowed. The setup is identical to VE1 except where stated otherwise in this section. This experiment aims to combine design simplicity with an evolving grounding line rather than to represent a system directly analogous to a real-world example.

3.2.1 Ice component setup

The VE2 initial geometry is given by

$$\begin{array}{}\text{(9)}& {\displaystyle }& {\displaystyle }{z}_{\mathrm{b}}=-\mathrm{20}-\mathrm{980}\left({\displaystyle \frac{x}{\mathrm{100}\mathrm{000}}}\right),\text{(10)}& {\displaystyle }& {\displaystyle }H={\displaystyle \frac{{\mathit{\rho }}_{\text{or}}}{{\mathit{\rho }}_{\mathrm{i}}}}\left(\mathrm{470}-\mathrm{400}\left({\displaystyle \frac{x}{\mathrm{100}\mathrm{000}}}\right)\right),\end{array}$$

where z_b is bedrock elevation relative to sea level and H is ice thickness. Then z_d and z_s are calculated based on floatation and the same densities as in VE1.

The depth-dependent inflow (x=0) and outflow (x=100 km for VE2_ER; x=99 km for VE2_EF) boundary conditions for the ice component are given by

$$\begin{array}{}\text{(11)}& {\displaystyle }& {\displaystyle }{P}_{\text{inflow}}\left(z\right)={\mathit{\rho }}_{\mathrm{i}}g(z_s-z),\text{(12)}& {\displaystyle }& {\displaystyle }{P}_{\text{outflow}}\left(z\right)={\mathit{\rho }}_{o}gz,\end{array}$$

where P_inflow and P_outflow are pressures prescribed at the inflow and outflow boundaries, respectively, and z is height relative to sea level (positive up). Zero normal velocity and free-slip tangential velocity conditions are imposed at the side walls given by y=30 km and either y=0 (for VE2_ER) $y=-\mathrm{1}$ km (for VE2_EF).

The grounding line is allowed to evolve, solving a contact problem (Gagliardini et al., 2013). The pressure acting on the underside of ungrounded ice is given by Eq. (3).

A sliding relation with a simple effective pressure dependency is used under the grounded ice (Budd et al., 1979, 1984; Gladstone et al., 2017),

$$\begin{array}{}\text{(13)}& {\displaystyle }{\displaystyle }{\mathit{\tau }}_b=-C{u}_b^m{z}_{*},\end{array}$$

where τ_b is basal shear stress, u_b is basal ice velocity, z_* is the height above buoyancy (related to effective pressure at the bed, N, by $N={\mathit{\rho }}_{\mathrm{i}}g{z}_{*}$), m is a constant exponent (set to $m=\frac{\mathrm{1}}{\mathrm{3}}$) and C is a constant sliding coefficient (set to $C={\mathrm{10}}^{-\mathrm{4}}\mathrm{MPa}{\mathrm{m}}^{-\frac{\mathrm{4}}{\mathrm{3}}}{\mathrm{a}}^{\frac{\mathrm{1}}{\mathrm{3}}}$).

Height above buoyancy is calculated by

$$\begin{array}{}\text{(14)}& {\displaystyle }{\displaystyle }{z}_{*}=\left\{\begin{array}{ll}H,& \text{if }{z}_{\mathrm{d}}\ge \mathrm{0}\\ H-z_{\mathrm{d}}\frac{{\mathit{\rho }}_{\text{or}}}{{\mathit{\rho }}_{\mathrm{i}}}.& \text{if }{z}_{\mathrm{d}}<\mathrm{0}.\end{array}\right.\end{array}$$

This is equivalent to assuming a sub-glacial hydrology system that is both fully connected to and in pressure balance with the ocean.

3.2.2 Ocean component setup

Ocean bathymetry matches the bedrock prescribed in the ice component (Eq. 9). The wet–dry scheme (Sect. 2.8) is used in this experiment, with a critical water column thickness of D_crit=5 m. Ocean potential temperature is initialized at −1.9 ^∘C, and salinity is initialized at 34.6 on the practical salinity scale.

ROMS-specific details. The ROMS setup is identical to verification experiment 1, except that the baroclinic (slow) time step size is 100 s, with 30 barotropic (fast) time steps for every slow time step.

FVCOM-specific details. The FVCOM model setup is identical to that of verification experiment 1.

3.2.3 Coupling

The cavity update method for VE2_EF is rate (Sect. 2.5). For VE2_ER it is corrected rate, with a correction factor of f_cav=0.01. With the 10 d coupling interval, this equates to a full correction timescale of approximately 3 years. Other coupling details are as in VE1.

Figure 4 Coupled system state after the first (a) and last (b) coupling intervals from the experiment VE1_ER (Table 2). The ice shelf is shown in grey, with basal melt rate computed by the ocean shown in colour on the underside of the ice shelf. Ocean streamlines are shown beneath the ice shelf, with the greyscale indicating magnitude of simulated ocean velocity. The vertical coordinate is given in metres; the horizontal coordinates are given in kilometres. This was a 100-year simulation.

Figure 5 Simulated mass evolution over time for the ocean component (dashed lines), the ice component (dash-dotted lines), and the total across both components (solid lines) from experiments VE1_ER (black) and VE1_EF (red).

Figure 6 Profiles through the centre line for experiment VE2 after the first ice component time step (a) and after 25 years (b). Ice flow speed is shown (flow direction is right to left). Ocean temperature (solid colour) and salinity (contours) are shown after 25 years (these are uniform at the start of the run, hence the solid colour for the ocean in the upper plot). Vertical exaggeration is 50 times. The gap between ocean and ice shelf is half an ocean grid cell and is a plotting artefact (the upper extent of the plotted region for the ocean is the uppermost rho point, which is half an ocean grid cell below the top of the ocean domain).

4 Verification experiment results

4.1 VE1: floating adjustment

Figure 4 summarizes the coupled system state at the start and end of simulation VE1_ER (see also Table 2 for a summary of the experiments). After the first coupling interval (10 d), the ocean component demonstrates a vigorous overturning circulation and high melt rates, especially in the deeper part of the domain. After the last coupling interval (100 years) the combination of melting and ice flow has caused a redistribution of the ice shelf, with an overall reduction in the along-domain gradients. The melt rates and overturning circulation are much weaker than at the start.

The ocean circulation throughout the simulation is predominantly a buoyancy-driven overturning along the domain, with very little cross-domain flow. The peak ocean flow speeds are always located at the top of the ocean domain directly under the ice shelf, where a fast, shallow buoyancy-driven flow from deeper to shallower ice draft is balanced by a much deeper return flow.

Figure 5 shows the evolution over time of the total mass of both ice and ocean components and the total coupled system from experiments VE1_ER and VE1_EF. Note that both ocean models employ the Boussinesq approximation and that the mass in Fig. 5 is calculated as volume multiplied by the reference ocean density from Table 2. Relatively rapid mass transfer from the ice to the ocean occurs during the first few years as the relatively warm ocean water transfers its energy to the ice. After this initial period of net melting, the ocean water temperature is close to freezing point and a long-term freezing trend can be seen that is stronger and more sustained in the ROMS ocean component than FVCOM. In a physically realistic coupled system, the ice and ocean would come into thermodynamic equilibrium and the spatial net mass transfer would approach zero.

The net mass change of the coupled system is more than an order of magnitude smaller than the mass change of the individual components for both experiments VE1_ER and VE1_EF. The current study does not use conservative regridding (Sect. 2.3), and therefore machine precision conservation is not expected. There are additional potential sources of error. The lag of the ocean component behind the ice component (Sect. 2.1.1) will cause a similar lag in total mass evolution. Use of the corrected rate cavity option (Sect. 2.5) prioritizes geometry consistency between components above mass conservation. The aim of analysing mass conservation in the current study is to ensure that the cumulative impact of these potential error sources is small compared to the signal. This has been achieved, and it will be possible to quantify and minimize or eliminate all sources of error in future studies using conservative regridding methods.

Figure 7 Ocean horizontal velocities in the upper layer (black arrows) and basal melt rate (red indicates melting; blue indicates refreezing) after 25 years of simulation VE2_ER (a) and VE2_EF (b). Outputs on the FVCOM mesh were regridded onto a 1 km regular grid. Both FVCOM and ROMS outputs were subsampled at 2 km resolution for this plot.

Figure 8 (a) A comparison of grounded area in the ocean component (total area of dry cells) against grounded area in the ice component (total area of grounded elements). (b) The difference between ocean and ice grounded area. These are from simulation VE2_ER. The green lines are drawn such that their distance apart is equivalent to the area of one row of ocean grid cells.

Figure 9 (a) A comparison of grounded area in the ocean component (total area of dry elements) against grounded area in the ice component (total area of grounded elements). (b) The difference between ocean and ice grounded area. These are from simulation VE2_EF. The green lines are drawn such that their distance apart is equivalent to the area of one row of ocean elements.

4.2 VE2: grounding line evolution

This is a partially grounded experiment in which the ice component boundaries are not closed, a through-domain flow of ice is allowed and the grounding line is allowed to evolve in the coupled system (described in Sect. 3.2). While the initial slope of the lower surface of the ice shelf is the same in both VE1 and VE2, the open inflow and outflow boundaries in the ice component and the relatively shallower ice in the grounded region both lead to a shelf that is much shallower in slope for VE2 than for VE1 for most of the simulation period. Figure 6 illustrates the shape of the ice sheet or ice shelf at the start of the simulation and after 25 years (from simulation VE2_ER). Note that the ice outflow boundary is more active than the inflow, with the flux into the domain through the inflow boundary remaining small and positive throughout the simulation. The ice draft is deepest in the middle of the domain, at around 30 km downstream (in terms of ice flow direction) from the grounding line. The ice draft impacts on circulation and melt, with the strong overturning of VE1_ER not present here. Melting occurs under the deepest ice, with refreezing elsewhere (Fig. 7).

Comparing the coupled simulation VE2_ER to the ice-only simulation (not shown) where the only difference is that the ice component features zero basal melt, it might be expected that the coupled simulation would exhibit a significantly thinner ice shelf due to melting. However, the ice dynamics partially compensate for this in terms of the ice geometry: the melt-induced thinning leads to acceleration in the ice, and the thickness difference is smaller than expected. However, this should not be interpreted as a stabilizing feedback response of ice dynamics to ocean-induced melting, as the increased ice flow would tend to drain the grounded ice more quickly, potentially triggering marine ice sheet instability (Schoof, 2007). Instead this effect may tend to partially mask an ocean-induced ice sheet destabilization if the observational focus is on ice shelf geometry.

As described in Sect. 2.8, the ice and ocean component each evolve the grounding line on their own time step and on their own grid or mesh. There is potential for discrepancy between ice and ocean grounded area due to method of cavity evolution (Sect. 2.5), regridding errors, the inherent differences between grids or meshes, and the methods used to determine grounding line position. While ice geometry is a key determinant of grounding line position, the ice component also tests for a contact force (Gagliardini et al., 2013) and the ocean component ROMS tests height above buoyancy against the free surface variable ζ (Sect. 2.8). Here we look at consistency of grounded area between components.

The evolution of grounded area in both ice and ocean components is shown in Fig. 8 for simulation VE2_ER. While the ice component employs an unstructured mesh of triangular elements (on the lower surface of the 3-D ice body), the ocean component employs a regular grid of square cells. The ocean component appears to exhibit a step-like reduction over time of grounded area. This is due to the row-by-row manifestation of grounding line retreat in the ocean component due to the alignment of grid rows with the linear down-sloping geometry. Grounding line retreat starts at the lateral edges of a row (ungrounding near the sidewall boundary), and the “wetting” of dry cells propagates toward the centre of the row. This step-like behaviour (with the spacing of the green lines in Fig. 8 indicating the total area of a row of cells) explains the main difference between ice and ocean grounded area. The evolution of grounded area is shown in Fig. 9 for simulation VE2_EF. Behaviour is similar to VE2_ER.

The initial rapid reduction in grounded area is due to the initial geometry. A region immediately upstream of the grounding line is initially very lightly grounded, and this region quickly becomes floating. The ocean component lags the ice component in this ungrounding, as can be seen in the first part of the difference plot in Figs. 8 and 9. This lag is in part due to the rate and corrected rate cavity update methods, in which the ocean component uses the most recent two ice component outputs to calculate a rate of change of geometry. This inevitably causes the ocean component to lag by approximately one coupling interval. The discrepancy may also be in part due to the fact that the region in question is close to floatation, and thus the threshold for dry cells to become wet is highly sensitive to ζ, at least for the ROMS implementation. In both experiments, the ice–ocean grounded area discrepancy has a tendency to reduce over time.

The computational time spent in both the ice and ocean components was measured for simulation VE2_ER. The ice component is more expensive than the ocean component during the first coupling interval but is significantly cheaper thereafter. Total time spent in the ice component over the 46-year simulation is approximately one-third that spent in the ocean component. The computational time spent within the central coupling code (calling routines and regridding) was negligible compared to time spent in ice and ocean components. This is with a 10 d coupling interval. If fully synchronous coupling is approached (i.e. if the coupling interval approaches the ocean time step size), the ice component will become much more expensive and it is possible the central coupling code may become significant. We do not anticipate fully synchronous ice–ocean coupling to become practical in the near future, at least not if the ice component directly solves the Stokes equations without simplifying assumptions, as is the case in the current study. The fully synchronous coupling of Goldberg et al. (2018) and Snow et al. (2017) is achieved by using the “shallow shelf approximation” for the ice component and running both components on the same grid.

5 Conclusions

We have presented a flexible coupling framework for ice sheet, ice shelf and ocean models that allows the user to choose between different ice and ocean components. We have demonstrated the functioning of this framework in simple test cases, both with and without a moving grounding line. We have demonstrated conservation of mass and consistency of grounding line evolution using semi-synchronous coupling.

FISOC provides runtime variable exchange on the underside of ice shelves. Providing such variable exchange at vertical ice cliffs, which are more common in Greenland than in Antarctica, will require minor developments to the coupling code, but the ocean components currently coupled through FISOC may need more significant developments in order to represent the buoyant plumes rising up ice cliffs.

Our coupled modelling framework is suitable for studying Antarctic ice sheet, ice shelf and ocean interactions at scales ranging from investigations of ice shelf channels (features with a spatial scale of typically a few km) up to whole Southern Ocean–Antarctic ice sheet coupled evolution. We are currently setting up simulations across this range of scales to address key processes surrounding Antarctic Ice Sheet stability and sea level contribution.