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 semisynchronous 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-iceshelf 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.

nisms 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 global sea level. In order to better quantify thesemechanisms using computer models, feedbacks need to be incorporated into the modelling system. To achieve this ocean and ice dynamic models must be coupled, allowing run time 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) 10 provides a modular approach to online 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 multiple methods for coupling the sub ice shelf cavity geometry. Thermodynamic coupling may also be activated. We present idealised simulations using FISOC with a Stokes flow ice dynamic model coupled to a regional ocean model. We demonstrate the modularity 15 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 idealised coupled ice -ocean system, both with and without grounding line movement.

Introduction
The Antarctic and Greenland ice sheets have the potential to provide the greatest contributions to global sea level rise on cen-20 tury timescales (Church et al., 2013;Moore et al., 2013), with the greatest uncertainty in projections being due to the 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 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 evi-25 dence (Christianson et al., 2016) and model studies (Favier et al., 2014;Gladstone et al., 2012). While MISI is fundamenatally 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 (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 ice sheet model 30 based studies use relatively simple parameterisations for calculating basal melt rates beneath ice shelves, a recent comparison has highlighted limitations of this approach (Favier et al., 2019). In particular, melt parameterisations as a function of depth or thermal driving do not impose conservation of heat in the system, and none of the parameterisations 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 Coupling projects take different approaches to handling the different timescales of ice and ocean processes. An ice sheet flowline model coupled to a five box ocean model allows large ensemble simulations to be carried out, but is limited in terms of 40 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 offline coupling with partial restarts for the ocean model (Seroussi et al., 2017;De Rydt and Gudmundsson, 2016;Thoma et al., 2015).

45
Here we present a new flexible Framework for Ice Sheet -Ocean Coupling (FISOC, Section 2). FISOC allows runtime coupling between ice and ocean components with a user choice of synchronicity options. Adopting Earth System Modelling Framework terminology (ESMF; Section 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 3D ocean models to an ice dynamic model and present idealised simulations demonstrating mass conservation and consistent 50 grounding line behaviour (Section 3). FISOC is also currently being used to contribute to ISOMIP+ and MISOMIP1.
FISOC is an open source coupling framework built using the Earth System Modelling Framework (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 55 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 (Section 2.1) facilitates swapping between different ice components or between different ocean components according to purpose (Section 2.2).

60
-Access to ESMF tools allows multiple regridding and interpolation options, including between regular grids and unstructured meshes (Section 2.3).
-Multiple options for handling differing ice and ocean time scales include fully synchronous coupling, passing rates of change, time averaging of variables (Sections 2.4 and 2.5).
-Flexible run-time control over the exchange of variables allows specific coupling modes to be (de)activated as required, e.g. geometric coupling, thermodynamic coupling.
-Grounding line movement (Section 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 (Section 2.5).
-Flexibility for parallelisation options. Currently sequential coupling is implemented, but any combination of sequential and concurrent parallelisation is possible with minimal coding effort.

70
-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 75 repository (see "code availability" at the end of this paper).

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 80 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 minimise the additional complexity.  -Call the component's initialise, run and finalise routines as required.
-Convert the component's grid or mesh to ESMF format, using ESMF data structures. FISOC adopts a fail-fast approach. Errors are generally considered to be fatal, in which case FISOC will log error information and finalise 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 run-time configuration file, and these can be checked for consistency in the model- There are a few cases where ESMF errors can be handled at run time. Details can be found in the FISOC manual.

Components
FISOC is designed to facilitate swapping between different ocean or ice components. Currently two different ocean components 115 and one ice component are available through FISOC. Table 1 summarises components currently coupled into FISOC. In some cases a non-standard build of the component is required for FISOC compatibility, and these are described in the FISOC manual, which can be obtained through the FISOC repository (Section 2.1).
The ice component Elmer/Ice (Gagliardini et al., 2013) is a powerful, flexible, state-of-the-art ice dynamic model. (2005)) is a 3D terrain-following sigma 120 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 nonhydrostatic options, 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). As stated above, FISOC provides coupling on a horizontal plane onto which the lower surface of an ice shelf can be projected.

The Regional Ocean Modelling System (ROMS; Shchepetkin and McWilliams
It is this plane on which ice and ocean properties are exchanged through the FISOC framework. 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 3D ice-ocean interface geometries are challenging not only for FISOC but also for the current generation of ice sheet and ocean models.

130
FISOC has access to all the run-time 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.
When using FISOC to couple Elmer/Ice to ROMS, the ROMS staggered grid, including ghost cells, extends beyond the 135 domain of the Elmer/Ice mesh. ESMF regridding methods also 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 with nearest STOD to extrapolate to points outside the domain.

Coupling timescales
The timescales for sub-shelf cavity circulation behaviour are in general much shorter than the timescales for ice flow and 140 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 being much larger. We define relevant terminology for coupling timescales: In principal, FISOC supports all three synchronicity options, though fully synchronous coupling is not practical to achieve when solving the Stokes equations for the ice. The experiments carried out for this paper use semi synchronous coupling with cavity geometry evolution as described in Section 2.5. Corrected rate. The same as above, except that a drift correction is applied to ensure ice and ocean representations of cavity 185 geometry do not diverge.
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 next coupling time-step.
Linear interpolation. The ocean representation of the ice draft is given by temporal linear interpolation between the two 190 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 initialised 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 195 FISOC initialisation. This ensures consistency.
Handling cavity evolution is a little more complicated in the case of an evolving grounding line, as discussed in Section 2.8 below.

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 200 adopt the three-equation formulation for thermodynamic exchange (Hellmer and Olbers, 1989;Holland and Jenkins, 1999).
This parameterisation 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. 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 205 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 the coupling time step (including the effect of melting/freezing, as well as any ice dynamical response), which is used to update the ocean model ice shelf cavity shape (Section 2.5).
For some applications, conductive heat fluxes into the ice shelf due to vertical temperature gradients in the ice at the ice-210 ocean interface are required by the three-equation parameterisation to calculate the flux balance at the ice ocean interface.
While ice-ocean thermodynamic parameterisations 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.

Interface pressure 215
Aside from the geometry evolution, an ocean boundary condition for pressure at the ice-ocean interface, P interf ace , 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 (and for dry cells, see Section 2.8) than ocean models can typically handle (Goldberg et al., 2018). In the current study, the ocean component uses ice draft (i.e. the depth of the ice-ocean interface) to estimate a floation pressure. ROMS assumes a 220 constant reference ocean density: where g is acceleration due to gravity, ρ or is a reference ocean density and D [O] is the ocean representation of ice draft. For the current study, all simulations with ROMS use ρ or = 1027 kg m −3 . FVCOM assumes a constant vertical ocean density gradient following Dinniman et al. (2007): where ρ o is ocean water density, ρ o1 is ocean water density of the top ocean layer and the vertical ocean water density gradient, dρo dz , is set to 8.3927 × 10 −4 kgm −4 .

Grounding line evolution
Grounding line movement in FISOC requires that both ice and ocean components support it. Numerical convergence issues A recent ice-ocean coupling study (Goldberg et al., 2018) used a "thin film" approach to allow grounding line movement.

235
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)  The original criterion in both ROMS and FVCOM for a cell to remain dry is given by: where η is the free surface variable, z b is the bottom boundary depth (bathymetry, aka bedrock depth), and D crit is a critical 245 water column thickness (or depth) for wet/dry activation. D crit is a parameter to be set by the user (typical values lie between 1 to 20m). Both η and z b are defined relative to sea level. Thus cells with a water column thickness less than D crit are designated 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/shelf system as follows: The ROMS criterion for a cell to be dry has been modified for the presence of a marine ice sheet/shelf system as follows: where S [O] is the ocean representation of ice sheet/shelf upper surface height and D [O] is the ocean representation of ice draft.
FISOC allows the ice component to pass any geometry variables to the ocean, such as ice draft, ice thickness, upper surface 255 elevation, or rates of change of any of these. In the event that geometry variables other than D are passed to the ocean, the same processing method is used as for D, as described in Section 2.8.
In the current study dD dt is passed to the ocean component, and in one case both dD dt and dS dt are passed (details in Section 3). When dS dt is passed, dS dt is processed the same way as dD dt . If the grounding line problem is solved, and if D is processed for passing to the ocean using the Corrected rate method, 260 equation 2 is modified to account for the dry water column thickness, which is initialised to D crit . The correction term changes from f cav D [I,t] − D [O,t] to f cav max(D [I,t] , zb + D crit ) − D [O,t] .

Verification experiment design
Simulations are carried out on idealised 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 265 experiment 2 (VE2) aims to assess whether the ocean and ice representations of grounding line evolution are consistent.

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 enough to demonstrate evolution of the coupled system. See Table 2 for run length and a summary of other model choices and parameter values used in VE1.

270
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 is purely geometric in that the ocean component Table 2. Model choices and input parameters used in verification experiment 1 (VE1, Section 3.1) and verification experiment 2 (VE2, Section 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 principal unstructured, but in practice structure can be seen (See Figure 2 middle and  Cavity correction factor, fcav n/a n/a 0.01 n/a Minimum water column Dcrit n/a n/a 5m 5m Ocean We expect adjustment toward a uniform thickness ice shelf to occur by two mechanisms: 275 1. Ice dynamics. The gravitational driving force will tend to cause flow from thicker to thinner regions.
2. Melt/freeze. The greater pressures at greater depth should result in higher melt rates, with the potential for refreezing under thinner regions.

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 280 simulation VE1_ER. However, ocean component FVCOM (used in VE1_EF) uses a semi-structured (in principal structured but in practice exhibiting some structure) mesh with dimensions 31 km by 99 km. 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 285 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 the ice component for VE1_EF 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 summarised in Figure 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 m by The pressure, P , acting on the underside of the ice shelf is given by where z is height relative to sea level (positive upward) and g is gravitational acceleration (set to 9.81 m s −2 ).
Temperature in the ice component is constant through space and time at −5°C. 300 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 net 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 2D horizontal 305 mesh of triangles with an approximate element size of 1km (VE1_ER) or 2km (VE1_EF) is extruded in the vertical to give 11 equally spaced terrain-following layers with the bulk element shape being triangular prisms.

Ocean component setup
The ocean bathymetry is set to 500 m throughout the domain. The wet/dry scheme (Section 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. Temperature 310 is initialised at −1.85°C and salinity at 34.6 psu. Ice-ocean thermodynamics are captured by means of the three-equation parameterisation (Section 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- 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 320 Number = 0.1) together with a constant background viscosity and diffusivity of 10 −6 m 2 s −1 . An eddy closure parameterisation (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.

Coupling
The coupling interval is 10 days, the same as the ice component time-step. Cavity update method is Rate (Section 2.5). For

325
VE1_ER the regridding method is bilinear with 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 (Section 3.1.1). There is no grounding line in this experiment.

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 330 the domain allowed. The setup is identical to VE1 except where stated otherwise in this section.

Ice component setup
The VE2 initial geometry is given by

335
where z b is bedrock elevation relative to sea level and H is ice thickness. Then D and S are calculated based on floatation and the same densities as in VE1.
No restrictions to ice flow are imposed at the upstream and down stream boundaries given by x = 0 and either x = 100 km (for VE2_ER) or x = 99 km (for VE2_EF). 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 = −1 km (for VE2_EF).

340
The grounding line is allowed to evolve solving a contact problem (Gagliardini et al., 2013).
A sliding relation with a simple effective pressure dependency is used under the grounded ice (Budd et al., 1979(Budd et al., , 1984Gladstone et al., 2017), 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, Height above buoyancy is calculated by: This is equivalent to assuming a sub-glacial hydrology system fully connected to, and in pressure balance with, the ocean. FVCOM specific details. The FVCOM model setup is identical to that of Verification Experiment 1.

Coupling
The cavity update method for VE2_EF is Rate (Section 2.5). For VE2_ER it is Corrected Rate with a correction factor of f cav = 0.01. With the 10 day coupling interval this equates to a full correction timescale of approximately 3 years. Other 360 coupling details are as in VE1.

VE1: Floating adjustment
This is a floating only experiment in which ice and ocean components adjust together towards an equilibrium state (described in Section 3.1). We do not run long enough to achieve equilibrium, but instead investigate the conservation of mass in the 365 coupled system as it evolves.   Figure 4 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. 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, 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 380 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.

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 385 allowed, and the grounding line is allowed to evolve in the coupled system (described in Section 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 5 illustrates the shape of the ice sheet/shelf after 25 years from simulation VE2_ER. The ice draft is deepest in the middle of the domain, at around 30km downstream (in terms of ice flow 390 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 (Figure 6).
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 395 thinning leads to acceleration in the ice and the thickness difference is smaller than expected. However, this should not be interpreted as a stabilising 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 destabilisation if the observational focus is on ice shelf geometry.
As described in Section 2.8 the ice and ocean component each evolve the grounding line on their own time-step and on their 400 own grid or mesh. There is potential for discrepancy between ice and ocean grounded area due to method of cavity evolution (Section¸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 η (Section 2.8). Here we look at consistency of grounded area between components.  behaviour (with the spacing of the green lines in Figure 7 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 Figure 8 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 415 line is initially very lightly grounded, and this region quickly becomes floating. The ocean component lags the ice component in this un-grounding, as can be seen in the first part of the difference plot in Figures 7 and 8. 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, thus the 420 threshold for dry cells to become wet is highly sensitive to η, at least for the ROMS implementation.
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 425 time spent in ice and ocean components. This is with a 10 day coupling interval. If fully synchronous coupling is approached (i.e. if the coupling time-step approaches the ocean time-step) 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); Snow et al. (2017) 430 is achieved by using the "shallow shelf approximation" for the ice component and running both components on the same grid.

Conclusions
We have presented a flexible coupling framework for ice sheet/shelf and ocean models which 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 435 using semi-synchronous coupling.
FISOC provides run time 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.