Ocean Modelling
The objectives of WP1 are to determine the changes in surface forcing required to switch the Antarctic continental shelf seas from cold to warm state (TP1), and to quantify the resulting changes in ocean-induced ice-shelf melting. The work package combines numerical ocean simulations of different complexity and hydrographic observations to (a) decipher the environmental setting at the Antarctic Slope Front necessary to allow on-shore flow of open-ocean waters, (b) determine the changes in external surface forcing required to switch the Antarctic continental shelf seas from cold to warm state – or vice versa, (c) quantify the tipping-point induced changes in ice-shelf basal melting (values transferred to WP2-3), and (d) check the observations for already existing indicators alluding to a tipping point in the Antarctic marginal seas (dissemination in WP4).
FESOM: The Finite Element Sea ice-Ocean Model is a global primitive-equation, hydrostatic ocean model developed at AWI. FESOM comprises a dynamic-thermodynamic sea ice component and the thermodynamic interaction at ice shelf bases. The former is a two-layer model, with the upper layer accounting for snow and the lower layer representing ice. Velocities for the advection of ice and snow volume are derived from an elastic-viscous-plastic solver based on a the three-equation approach with a velocity-dependent parameterization. The model accounts for the thermodynamic coupling, i.e., for the heat loss from the ocean to the ice, and a glacial melt water flux is derived and applied to the salinity in the layer adjacent to the ice shelf base. Despite this coupling, which results in mass loss and gain at the base, the ice shelf cavity geometry remains constant in time.
The governing equations of FESOM are solved on a horizontally unstructured mesh all at the same location, which is analogous to the Arakawa A-grid. The horizontal mesh consists of triangular elements at the surface while the 3D domain is represented by tetrahedral elements. In the vertical, nodes are aligned below the surface nodes. In addition, a hybrid grid is used with 25 terrain-following (sigma) layers down to a pre-defined isobath (e.g., 2500 m) around the Antarctic continent, while 36 z-layers resolve the ocean north of the isobath (Fig. 6). A finite-element grid allows for local refinement of the horizontal resolution, which is particularly useful in high latitudes, where the internal Rossby radius of deformation becomes very small (O[3 km]). Examples for this refinement are shown in Fig. 7 – (a) low resolution from 4 km in the southern and western Weddell Sea to 73 km in the northeastern corner of the Weddell gyre, (b) high resolution from around 15 km over the whole Weddell gyre, including its eastern boundary, to a mean resolution of 3.5 km on the southern Weddell Sea continental shelf including the Filchner-Ronne Ice Shelf cavity. Although is hat been pointed out that a model resolution O(1 km) is needed to fully resolve mesoscale processes, it has to be noted that the (global) mesh with a total of 11 million nodes is at the limit of computational capacity at the time of writing. Therefore, TiPACCs also includes regional models (see below), confined to areas of interest and forced at the open boundaries with the output of the global model, which allow for the resolution needed to consider the relevant processes in the Antarctic marginal seas.
Website: https://fesom.de
NEMO: The Nucleus for European Modelling of the Ocean (NEMO) is a state-of-the-art modeling framework providing a flexible tool for studying the ocean and its interactions with the other components of the earth climate system over a wide range of space and time scales. It is developed and maintained by a pan-European consortium. The core physical models include ocean (OPA, Océan PArallélisé) and sea ice (LIM, Louvain-la-Neuve Sea Ice Model, currently being merged into SI3, Sea Ice modeling Integrated Initiative), along with optional biogeochemistry, adaptive mesh refinement, and data assimilation components.
The NEMO community ocean model is a primitive equation dynamic/thermodynamic model adaptable to both regional and global ocean circulation problems down to kilometric scales (Madec, 2016). In the horizontal direction, the model uses a curvilinear orthogonal grid, while in the vertical direction, coordinate choices are full or partial step z– or z*, s, or hybrid z/s. The distribution of variables is a three-dimensional Arakawa C-type grid, staggered in all dimensions such that velocity components are carried at intermediate grid points surrounding the grid points that carry the scalar properties.
The sea-ice component of NEMO is a state-of-the-art, finite-difference implementation, solved on curvilinear orthogonal coordinates with variables arranged on a C-grid. It combines the conservation of momentum for an elastic-viscous-plastic continuum, energy and salt-conserving halo-thermodynamics using one layer of snow and five of ice, and an explicit representation of sub-grid scale ice thickness variations using a multi-category distribution.
The NEMO system includes standard global configurations with differing horizontal resolution (“ORCA” grids), from which regional domains can readily be created. All configurations are based on an isotropic Mercator grid, which provides improved horizontal resolution at high latitudes. Recent work undertaken at UKRI-BAS has added the capability of simulating the circulation beneath ice shelves to the NEMO ocean model. This capability has been used in several global and regional configurations47,48. Another important capability for the Antarctic seas is a Lagrangian iceberg module.
Website: https://www.nemo-ocean.eu/
MITgcm: The MIT General Circulation Model (MITgcm) is a numerical model designed for the study of atmosphere, ocean and climate. It solves the primitive equations using finite volume techniques that support the treatment of irregular geometries using horizontal orthogonal curvilinear grids and vertical partial cells. Its non-hydrostatic formulation enables it to simulate fluid phenomena over a wide range of scales from planetary down to sub-kilometre. Forward and tangent linear counterparts of MITgcm are supported using an “automatic adjoint compiler”. These can be used in parameter sensitivity and data assimilation studies. By employing fluid isomorphisms, one hydrodynamical kernel can be used to simulate flow in both the atmosphere and ocean.
The code includes specialised packages that optionally deal with biogeochemical processes, ocean interior and boundary layer processes, including a range of sub-grid-scale mixing parameterisations, atmospheric processes, sea-ice, ice shelves, coupled simulations and state estimation. The MITgcm sea ice model (MITgcm/sim) is based on a variant of the viscous-plastic (VP) dynamic-thermodynamic sea ice model with many important aspects of the original code modified and improved. The code has been rewritten for an Arakawa C-grid, includes three different options for solving the nonlinear momentum equations, different options for the ice-ocean stress and a three-layer thermodynamic model. The “shelfice” provides an additional pressure exerted by the floating ice shelf on the top of the ocean and a thermodynamic model for basal melting of the ice.
All three ocean models have the capability to consider the influence of tides on on-shore flow and cavity circulation with consequences for the ice shelf basal mass balance.
Website: http://mitgcm.org/
Ice Modelling
The objective of WP2 is to determine the stability regime of the grounding lines of the Antarctic Ice Sheet and the existence of tipping points with respect to ice-shelf melt (TP2). The work package involves, primarily, the use of numerical ice-flow models.
Using three state-of-the-art large-scale ice-flow models (PISM, Úa, Elmer/Ice) we will determine if perturbations in applied ice-shelf melt gives rise to an unstable and irreversible grounding-line retreat. The existence of ice-sheet tipping points will be investigated using several different approaches, as outlined above in the section Definition of tipping points and how they will be assessed. Firstly, as exemplified in Fig. 4a-c, transitions to unstable and self-sustained retreat will be identified by calculating rates of grounding-lines migrations for individual ice-streams following step-wise perturbations in applied ice-shelf melt. Secondly, the stability of grounding lines will be determined by calculating the functional relationship between ice flux and ice thickness. Thirdly, for those numerical models using a fully-implicit forward time integration, the stability will be diagnosed directly from the eigenvalues of the forward time operator of the numerical models. We will furthermore assess if the tipping points were crossed in the past, during the Last Interglacial, and will quantify the related impact of Antarctic Ice Sheet retreat on global mean sea level.
PISM: The Parallel Ice Sheet Model PISM, developed at University of Alaska (UAF) and the Potsdam Institute of Climate Impact Research (PIK) is an open-source ice dynamics model, with the code available at http://www.pism-docs.org. The model is thermomechanically coupled and uses a hybrid shallow approximation of ice flow. PISM ice velocities therefore transition smoothly from vertical shear-dominated, slow-moving regions in the interior to the faster flow in ice streams and shelves. The efficient superposition allows for resolutions of up to few kilometres for whole Antarctica simulations.
Temperatures within PISM are determined based on an energy-conserving enthalpy scheme including a thin sub-glacial water layer and a layer of the thermal bedrock. PISM further includes a subglacial hydrology model that allows for spatially and temporally variable shear stresses. The subglacial hydrology model transports water under mass conversation to the grounding line or land-terminating ice margins.
Both the grounding line and the calving front are simulated at subgrid scale in PISM and evolve according to the physical boundary conditions. Grounding-line motion is reversible and shown to be consistent with full-Stokes simulations for higher resolutions. The improved representation of grounding line motion was confirmed within the framework of the MISMIP-3D model inter-comparison project. PISM has participated in a number of model intercomparison projects, including the MISMIP+ project, ABU-MIP as well as LAR-MIP. Furthermore, it takes part in the Ice Sheet Model Intercomparison Project ISMIP6.
For the coupling of continental ice and ocean, PISM includes the Potsdam Ice shelf Cavity mOdel PICO. PICO simulates sub-shelf melt rates consistent with the vertical overturning circulation underneath ice shelves. It can adapt to the advance and retreat of grounding lines and is able to reproduce the wide range of melt rates currently observed for Antarctic ice shelves. PICO can therefore serve as a coupler between ice-sheet and ocean models.
Due to adaptive time-stepping and the employment of the SIA and SSA, PISM is computationally very efficient and capable of simulating large ensembles on multi-millennial time-scales. Previous simulations of the Antarctic Ice Sheet have shown rapid transitions and threshold behaviour, in idealistic setups and for Antarctic drainage basins which makes the model well-suited for the project proposed here.
Website: https://pism-docs.org/wiki/doku.php
Elmer/Ice: The open source Finite Element software for ice-sheets, glaciers and ice flow modelling Elmer/Ice is the glaciological add-on package to Elmer, which is a multi-physics finite-element suite mainly developed by numerical experts at CSC in Finland. Elmer/Ice was initially developed to study processes at a local scale but it has progressively proved its ability to investigate large-scale ice flow evolution. Elmer/Ice has been providing reference results for the model inter-comparison projects ISMIP-HOM, MISMIP and MISMIP3 and actively participated to MISMIP+. Since 2004, more than 100 publications were based on simulations computed with Elmer/Ice. The worldwide recognition of Elmer/Ice is due to its physics completeness, its ability to be initialized to the actual ice-sheet state and simulate coastal changes (grounding line, calving). Based on the finite element method, unstructured meshes can be produced to better capture ice flow in regions of interest (e.g. grounding line area, fast flowing ice stream). For the current project, both SSA and full-Stokes ice flow simulations will be performed. The comprehensive list of Elmer/Ice capabilities in terms of physics and numerical methods is detailed on a dedicated website (http://elmerice.elmerfem.org/).
Úa is a finite-element ice-flow model developed at the University of Northumbria, Newcastle, UK, by Hilmar Gudmundsson. The model is based on a vertically integrated formulation of the momentum equations and can be used to simulate the flow of large ice sheets such as the Antarctic and the Greenland Ice Sheets, ice caps and mountain glaciers (see Fig. 8 for an example). The ice-flow equations are solved on an unstructured mesh consisting of linear, quadratic or cubic triangular elements. Various meshing options are available, including automated mesh refinement and coarsening. When simulating the flow of marine ice sheets these meshing options allow, for example, the areas around grounding lines to be automatically highly resolved as the grounding-lines migrate through the computational domain. Elements can be activated and deactivated. This enables the computational domain itself to change in the course of a run, for example when simulating the growth and decay of a large group of mountain glaciers. Inversion for model parameters is done using the adjoint method. Both Bayesian and Tikhonov formulations are supported. Ice-thickness positivity is enforced using the active set-method. The model has been used extensively to study impacts of ice-shelf buttressing on grounding-line dynamics and has recently been coupled with MITgcm.
Website: https://github.com/GHilmarG/UaSource
Coupled Ocean-Ice Modelling
The objective of WP3 is to determine the impacts that a switch in the ocean state from cold to warm (TP1) has on the stability regime of the grounding lines (TP2) and the resulting implications for global sea level. The work package involves, primarily, the use of coupled numerical ice-flow and ocean circulation models.
Using three coupled ice-ocean models (PISM+FESOM, Úa+MITgcm, Elmer/Ice+NEMO) we will determine potential interations between the ocean and the ice-sheet tipping points, TP1 and TP2. In particular, we will detetermine if TP1 gives rise to TP2, i.e. if ice-shelf melt rates correponding to warm oceans give rise to an unstable and self-sustained grounding-line retreats. Furthermore, we will determine how changes in ice-shelf cavities in turn affect the ocean circulation, and if both TP1 and TP2 arise in fully coupled simulations.
FESOM-PISM: The coupling of ice flow model PISM and ocean model FESOM follows an ‘asynchronous’ approach. For each coupling time-step of typically a few months to one year, PISM uses the model output at the ice margins provided by FESOM as boundary conditions to compute the ice flow in sheet and shelves.
For the experiments proposed here, sub-shelf melt rates will generally be provided to PISM directly by FESOM. In regions where the model grids and cavity geometries differ significantly in the ice and ocean model, however, we will follow a complimentary approach, using the ice-shelf cavity model PICO as a coupler. FESOM will then provide PICO – which is readily available as a module in PISM – with the necessary ocean temperatures and salinities outside of the ice shelf cavities, while sub-shelf melt rates are computed by PICO, consistent with the vertical overturning circulation for each ice-shelf cavity. Following this approach, PISM+PICO has been shown to reproduce the wide range of sub-shelf melt rates for present-day Antarctica and can easily adapt to changes in the grounding line and ice-front positions due to the interaction with the surrounding ocean. From the ice draft and grounding line location at the end of each coupled PISM simulation time-step, a new cavity geometry is derived and an updated FESOM mesh generated. FESOM’s prognostic variables are then projected to the new mesh and FESOM is run forward for another coupling time-step, and the cycle repeats.
NEMO-Elmer/Ice: Elmer/Ice and NEMO have recently been coupled for the MISOMIP inter-comparison project. The two components are coupled continuously with a typical coupling period of a few months. At every asynchronous coupling step, NEMO provides melt rates averaged over the coupling period to Elmer/Ice, then Elmer/Ice runs over the coupling period and provides the updated ice-shelf geometry to NEMO. Ocean temperature, salinity and velocity are preserved through coupling steps. To handle changes in water column thickness while avoiding the generation of spurious barotropic waves (tsunamis), we nonetheless impose the conservation of the local barotropic transport at every coupling step. Initialization of the coupled simulations is done by running NEMO with a static ice shelf and imposing the resulting melt rates in the initialization of Elmer/Ice. In this H2020 project, we will move towards more realistic configurations at the regional or global scales, building upon recently developed NEMO configurations that include ice-shelf cavities.
MITgcm-Úa: The coupling between Úa and MITgcm follows an ‘asynchronous’ approach. This means that Úa and MITgcm are stepped forward continuously, and information is exchanged at fixed intervals (coupling time-step), which is typically, but not necessarily, larger than the individual ice and ocean time steps. At each coupling time step, updated meltrates are delivered from the ocean to the ice model, and the new shape of the iceshelf cavity is transferred from the ice to the ocean model. To accommodate the differences in model grids and step-wise changes in cavity shape between coupling time-steps, Úa nodes are chosen to include all MITgcm nodes. Therefore no interpolation is required when transferring data from Úa to MITgcm. Melt rates for Úa nodes are obtained through linear interpolation of the MITgcm grid. The thermohaline history of the ocean and the momentum budget are retained at each coupling time-step. Newly opened ocean cells are filled with tracer values from the nearest flooded cell and initial velocities are set to zero. At each coupling time-step, the ice-shelf induced pressure loading onto the ocean is recalculated to avoid the formation of unphysical waves (‘tsunamis’).
The coupling time-step is a model parameter and needs to be chosen at the start of a coupled simulation. Typical values are on the order of several months to years. Note that the coupling approach used here reflects recent developments and is different from the previous ‘discontinuous’ approach where ocean history was not retained and the ocean circulation was spun up from initial conditions at each coupling time-step.