Semi-Lagrangian transport of oxygen isotopes in polythermal ice sheets : implementation and first results

Modelling the evolution of the Earth system on long timescales requires the knowledge and understanding of driving mechanisms, such as the hydrological cycle. This is dominant in all components of the Earth’s system, such as atmosphere, ocean, land surfaces/vegetation and the cryosphere. Observations and measurements of stable water isotopes in climate archives can help to decipher and reconstruct climate change and its regional variations. For the cryosphere, theδ18O cycle in the current generation of Earth system models is missing and an efficient and accurate tracer transport scheme is required. We describe ISOPOLIS 1.0, a modular semi-Lagrangian transport scheme of second-order accuracy, which is coupled to the polythermal and thermomechanical ice sheet model SICOPOLIS (version 2.9). Model skill is demonstrated by experiments with a simplified ice sheet geometry and by comparisons of simulated ice cores with data from Greenland (GRIP) and Antarctica (Vostok). The presented method is a valuable tool to investigate the transport of any kind of passive tracer inside the cold ice part of a polythermal ice sheet and is an important step to model the wholeδ18O cycle.


Introduction
Oxygen isotopes are an important proxy for the reconstruction of temperatures of the past.Air temperature is related to stable isotopic composition of precipitation as indicated by observations (e.g.Dansgaard, 1964;Gat, 1996;Jouzel et al., 1997;Gornitz, 2009).When assuming this relationship remains the same in the past, it is possible to reconstruct past temperatures from ice or sediment cores.To model the whole isotopic cycle in an Earth system model (ESM), components for the atmosphere, ocean and cryosphere are required.The biggest components of today's cryosphere are the two huge ice sheets of Greenland and Antarctica which together contain about 99 % of the Earth's ice volume and about 80 % of global fresh water.Therein and in smaller glaciers ice appears in two different states; as "cold" ice with temperatures below the pressure-melting point and "temperate" or "warm" ice where temperatures are at the pressure melting point.Ice masses neither pure "cold" or "temperate" are called polythermal ice.Temperate ice may have liquid water in it and, therefore, may be considered as at least a two-component fluid.In ice sheets this temperate ice may exit as a thin layer near the base.Although it is much less in volume than cold ice, temperate ice has significant consequences on the ice dynamics (Lliboutry and Duval, 1985) and therefore on isotope and tracer transport in general.
Models for oxygen isotopes in the atmosphere, ocean and land-biosphere have existed for some time (e.g.Jouzel et al., 1987;Hoffmann et al., 1998;Sturm et al., 2010) but the cryospheric part was first successfully modelled by Clarke and Marshall (2002) for Greenland.Until then, thermomechanical ice sheet models were mainly unsuccessful or had large limitations in their dimension, spatial extent, temporal coverage or steady-state assumptions.Problems were mainly associated with the Eulerian scheme which is commonly used to solve for advection processes in numerical ice sheet models.The main problems are discontinuities of the advected property which is usually tackled by an artificial diffusion term.This leads to an inaccurate solution near the base, and reliable solutions for tracer dispersion are only optioned for the upper half of the ice sheet (Rybak and Huybrechts, 2003).Clarke and Marshall (2002) and the follow-up papers (Clarke et al., 2005;Lhomme et al., 2005a, b), therefore, use a semi-Lagrangian scheme to track passive tracers.Tarasov and Peltier (2003) also uses a semi-Lagrangian scheme to model the transport of δ 18 O directly while Clarke and Marshall derived the δ 18 O value indirectly via ice age and ice provenance.In a semi-Lagrangian scheme a regular grid is used and particles are usually tracked back to its origin from one time step before.This backtracking is of importance for the overall performance and accuracy.A review by Staniforth and Côté (1991) recommends at least second-order accuracy for the backtracking scheme.Clarke and Marshall (2002), Clarke et al. (2005) and Lhomme et al. (2005a, b) use a firstorder backtracking scheme in their ice sheet model, considering cold ice only.
In this study we are using the second-order backtracking scheme of de Almeida et al. (2009) together with a polythermal ice sheet model to simulate the δ 18 O distribution in ice sheets.The whole model consists of three components: a polythermal ice sheet model, a semi-Lagrangian transport module and a post-processing tool which are described in detail in Sect. 2. As polythermal ice sheet model we use SICOPOLIS (Version 2.9.) which is based on Greve (1997a, b).SICOTRACE and SICOSTRAT are two new components designed for semi-Lagrangian transport and the reconstruction of the stratigraphy.While SICOTRACE calculates the transport variables, SICOSTRAT can be used to generate cross sections along grid points and ice cores at arbitrary locations with information such as δ 18 O concentration, place of origin and age of the ice.The performance of the model is demonstrated with a simplified geometry of the EISMINT (Payne et al., 2000) inter-comparison project and later applied to Greenland and Antarctica (Sect.3).A summary of the paper and general conclusions are given in Sect. 4.

Model description
The diffusion of δ 18 O over multi-annual periods is considered to be negligible in cold ice (e.g.Jean-Baptiste et al., 1998).Therefore, it can be considered as a passive tracer in cold ice which is not altered chemically or physically on its way through the ice sheet and does not influence the flow.
When ice crosses the cold-temperate surface (CTS) and enters the temperate regime the passive trace assumption is no longer true since diffusion is possible in the liquid water which my be present there.Therefore it would be possible to introduce a marker which indicates whether the ice crossed the CTS at some time in its flow history.This would possibly aid ice core interpretation but will not be implemented at this stage in the model.
The advection of such a passive tracer in Eulerian description is where the transported property itself and its gradient must be given.This is not the case when discontinuities are present, e.g. when two ice-flow branches are merging or when the surface mass balance or conditions near the bedrock are changing (Clarke et al., 2005).
In contrast, a Lagrangian description is not influenced by such discontinuities as individual particles are tracked, although the Lagrangian description encounters other problems such as error accumulation along the trajectory because of the required velocity interpolation (Rybak and Huybrechts, 2003).Due to the divergence of ice flow there are areas with very low particle density which lead to very low information density at the same time.This necessitates having a large number of modelled particles whereby new ones are constantly introduced at the ice surface and old ones are removed at the ice sheet's base (Lhomme, 2004).These and other practicalities linked to the irregular grid render Lagrangian schemes unsuitable in ice sheet modelling.
To overcome the drawbacks of the Eulerian and Lagrangian scheme, the semi-Lagrangian scheme tries to combine the best of both, namely the regular grid of the Eulerian and the better stability of the Lagrangian scheme.

Temperature reconstruction with water isotopes
Today's δ 18 O ratio in snow at any given location is controlled by the local temperature, the amount of precipitation, the distance from the coast and the altitude (e.g.Gornitz, 2009).Analyses indicate a strong correlation between δ 18 O and surface temperature, especially for temperatures below 20 • C in mid-to-high latitudes and with the amount of precipitation at low latitudes.
The measured correlation between δ 18 O and temperature in polar regions is strongest for annual means (e.g.Gornitz, 2009).Observational data indicate that the relationship of temperature and δ 18 O are not spatially uniform and that it can be described as a linear function of the surface temperature T s (Eq.2) or a function of the present δ 18 O value, local surface elevation S and of temperature change T over the ice sheet (Cuffey, 2000) (Eq.3): (2) (3) Clarke and Marshall (2002) formulate an indirect approach by transporting so-called "depositional provenance labels" t, x, y and by maintaining a "depositional archive" with information about surface temperatures, ice topography and mass balance (see Fig. 1).Hence, it is possible to reconstruct the δ 18 O value with a transfer function in the form of Eqs.

General framework
(2) or (3) and get the additional benefits of information of age and transport.
If the history of the surface temperature is stored and the time and place of origin is known, the δ 18 O can be calculated by inverting Eqs.(2) or (3).The history of the surface temperature is a climate variable which can either be an external input from an atmospheric model or prescribed as is the case here.
The other information of the origin (x, y) and time (t) of deposition is transported with the ice flow and is not altered on its way, i.e. it is a passive tracer.In order to calculate a passive tracer, the velocity field needs to be known, which is calculated by the ice sheet model SICOPOLIS in our study.Tarasov and Peltier (2003) modelled δ 18 O directly and use a very high vertical resolution of up to 4096 layers on a subgrid near ice core locations.Their approach is not feasible for our purpose because we need to know the δ 18 O value in the whole ice sheet in order to close the hydrological cycle in the ESM framework.
The new transport model consists of three separate modules: SICOPOLIS which is the polythermal ice sheet model, SICOTRACE which is the tracer transport model and SICOSTRAT which reconstructs the stratigraphy.Figure 1 illustrates the framework with in-and outputs of each component.

Polythermal ice sheet model SICOPOLIS
SICOPOLIS (Greve, 1997a) is a three-dimensional polythermal ice sheet model and is based on the shallow ice approximation.The model has been applied previously for Greenland, Antarctica, the polar ice caps of Mars and the entire Northern Hemisphere.As input it uses the surface mass balance, mean air temperature, eustatic sea level and geothermal heat flux (e.g.Greve, 1997b).As output, temperature, water Hence, there are two separate sigma-transformed regular grids above each other, and they share their horizontal coordinate x → ξ and y → η (not shown here).The grid spacing ζ c and ζ t for cold and temperate are usually different and also the sigma-transformation is different for both domains (see Greve, 1997a).CTS denotes the cold-temperate surface.content in temperate regions, ice extent and thickness, ice velocities, isostatic displacement and the temperature of the lithosphere are calculated.In this study, the velocities in the cold ice v c and temperate ice v t , and the respective thickness H c and H t are of major importance.The two domains are separated by the CTS at z = z m (see Fig. 2).As an input for the transport component SICOTRACE the three-dimensional velocity field in the cold and temperate region as well as the ice extent and thickness is stored for every time step (Fig. 1) in a netCDF file (Rew and Davis, 1990).

The numerical grid
SICOPOLIS uses three grids with different sigmatransformations on top of each other (see Fig. 2).A grid for the lithosphere is overlaid by a grid for temperate ice below cold ice.In general, values in the cold ice domain are denoted with subscript c and in the temperate region with subscript t.
For cold ice the σ -transformation is where a is a stretching parameter which is 5 herein (Greve, 1997a).
The transformation in temperate ice with bedrock elevation b (see Fig. 2) is With this sigma transformation the depositional provenance labels are ξ d , η d , τ d .
The two domains for cold and temperate ice use an Arakawa C-Grid (Arakawa, 1997) with the velocity components at intersections between grid point centres.In this grid the velocity components in cold and temperate ice, beside gradients, are placed in between grid points.All the other values, like temperature or water content for example are defined on the grid centre, itself.
The two different grids in the polythermal ice sheet model and the fact that the ice extent and the thicknesses of cold and temperate ice changes with time makes the semi-Lagrangian transport more difficult for glacier models than for atmospheric or ocean models.This will be discussed in more detail in the next section.

Semi-Lagrangian transport module SICOTRACE
SICOTRACE stands for SImulation COde of TRACErs and is a separate program which reads SICOPOLIS outputs and calculates the transport of the three provenance variables.In this section the semi-Lagrangian transport with backward trajectories is described as well as the problems and remedies associated with the specifics of ice sheet models, especially for polythermal ice.

Semi-Lagrangian scheme
The basic idea of a semi-Lagrangian scheme is illustrated in Fig. 3.A transported property is constant over time and is not changing along its trajectory: If the place of origin at time t − t is known, then the value of (t) at the current time step is also known.A standard semi-Lagrangian scheme uses a regular spaced grid, and for each grid point the back trajectory is calculated to get the position of the particle at one time step before.This position is usually not on the grid.Therefore, an interpolation is required to get the value (t− t) for which we used a trilinear interpolation (Press et al., 1996).
In order to reach high accuracy, a review of semi-Lagrangian schemes for atmospheric models by Staniforth and Côté (1991) recommends the use of a backtracking method with at least second-order accuracy.However, in the original work by Clarke and Marshall (2002) and subsequent papers (Clarke et al., 2005;Lhomme et al., 2005a) a firstorder approximation was used.Much research has been done in developing accurate, cheap and robust backtracking methods, (e.g.McGregor, 1993;Nair et al., 2003;Purser and Leslie, 1994;Staniforth et al., 2003;Hortal, 2002).For our application we opt for the scheme of de Almeida et al. ( 2009) because it is of second-order accuracy and its robustness for small and large Courant number in the presence of weak and strong flow curvature makes it well suited for a transient and three-dimensional ice sheet model.

Numerical grid for the semi-Lagrangian transport
In general it is possible to perform semi-Lagrangian transport on the same grid as used by the Eulerian ice sheet model.In a polythermal model, however, two layers with σ -transformed coordinates are stacked on each other.Therefore, a new grid consisting of the same horizontal grid as in SICOPOLIS but with a σ -transformation including both cold and temperate ice (see Fig. 4) has been formulated: The time step for the semi-Lagrangian transport remains the same as in SICOPOLIS.In general, it could be longer than for the Eulerian ice sheet model but for now it is set to be the same as in SICOPOLIS.

Trajectory calculation
The backtracking method is important for the overall performance of a standard semi-Lagrangian transport scheme.Here, we use the second-order accuracy scheme of de Almeida et al. ( 2009) which allows the velocity and acceleration vectors of particles to vary between t and t − t. Figure 5 shows the scheme in one dimension.The scheme is based upon the general idea of multistage methods, where each integration step requires the estimate of the dependent variable at several intermediate times, similar to Runge-Kutta (de Almeida et al., 2009).
Compared to first-order approximations the second-order scheme is known to be more accurate and exhibits better conservation properties (e.g.McGregor, 1993;Staniforth and Pudykiewicz, 1985).A discussion of the differences in our application is given below (Sect.3.1.2).
Point A is the desired ending point on the grid point ξ(i).The red line EA is the exact trajectory and the dashed black line E A is the approximate one.We know the property at all grid points ξ i at one (t − t) and two time steps (t − 2 t) beforehand and we wish to obtain values at the same grid points at time t.
In this method the trajectory is split into two steps.In the first step, starting from grid point ξ(i) at time t the particle is displaced backwards for a time interval t with the velocity calculated at an intermediate position at point B for time t − t.In the second step, starting from the particle position C at time t − t calculated in the previous step, the particle is displaced backwards for another t, with the velocity calculated for the intermediate position D, also at time t − t.The intermediate points B and D are obtained by considering displacements of the particle for a time interval t/2.
The velocities at point B and D need to be interpolated.For the interpolation of the velocity field we use a trilinear scheme which is sufficiently accurate (Behrens, 1996).The positions of the points themselves rely on the velocities and, hence, the whole process needs to be repeated until a convergence criterion of = 0.001 m is met.
In three dimensions Eq. ( 6) is approximated by Eq. ( 8) by applying the described backtracking for all three co-ordinates (ξ , η, ζ s ): where (ξ , η , ζ s , t −2 t) is usually not on the grid and is interpolated.In three dimensions the first step is and the second step is Equations ( 9) and ( 10) are valid for cold and temperate ice with different velocities v c , v t and sigma transformations σ c , σ t .There are some difficulties associated with semi-Lagrangian transport and sigma coordinates in polythermal ice sheet models which are not encountered in models of other compartments of the climate system.These are discussed in the next section.

Polythermal ice sheet and sigma transformations
The sigma transformations σ c , σ t are dependent on the ice thickness of cold and temperate ice, the position of the coldtemperate surface and the bedrock elevation (see Eqs. 4c and 5c).These quantities are all dynamic which makes it necessary to calculate them at times t − t 2 and t − 3 t 2 in order to compute ζ * and ζ * * (see Eqs. 9a-c and 10a-c): Geosci.Model Dev., 7, 1395-1408, 2014 1401 Therefore the fields of H c , H t , z m , b are calculated at times t − t 2 and t − 3 t 2 and then bilinearly interpolated at location ξ i − a * (m) , η j − b * (m) .This is an additional complication encountered in ice sheet models whereas in for example ocean models ζ * would simply be σ −1 (z = z k − c * (m) , t − t 2 ).In addition it is possible that the tracked particle is crossing the CTS and therefore it is necessary to check at each iteration in which domain the particle is.This is done by checking whether the z coordinate is below or above the z m at the current position ξ and η.

Boundary conditions at the ice-bedrock interface and the ice surface
Basal melting and isostatic adjustment are included in the velocity fields.Hence, no special treatment is required here.
On the other hand, the boundary condition at the ice surface is expressed by the surface mass balance.A positive mass balance means that new ice is forming on top of the ice sheet while a negative mass balance means loss of ice.In the case of a positive mass balance the three provenance variables are the current values at that grid point, whereas t d is the current time and x d , y d are the coordinates at the grid point.Since the applied backtracking scheme uses two time levels, it can arise that point C is at a position which was ice free at the last time step.In such a case the backtracking stops and the current value of provenance variable is interpolated from the field (t − t).

Deriving the stratigraphy with SICOSTRAT
SICOSTRAT (Simulation COde for STRATigraphy) generates the stratigraphy of δ 18 O everywhere within the ice sheet by inversion of Eqs.
(2) or (3) and with the provenance labels as well as the values stored in the depositional archive.
To calculate the δ 18 O value with the simple relationship in the form of Eq. ( 2), where δ 18 O only depends on the surface temperature, the routine is as follows: 1. Convert z to ζ s , then use trilinear interpolation to get ξ d , η d and t d .
2. T s (x, y, t) can be interpolated bilinear in space and linear in time to calculate the surface temperature the particle had when it was deposited on the ice sheet surface.
3. Calculate δ 18 O with the transfer function.
This procedure can be used to generate the isotopic stratigraphy of the whole ice sheet which can be written as a netCDF and later be used by an Earth system model.SICOSTRAT uses generic mapping tools (GMTs) to generate plots of ice cores at arbitrary locations and cross sections along grid lines.This makes it possible to validate the model against ice core data and to get a general overview of the isotopic composition.

Results
In a first study, the model is applied to the EISMINT intercomparison project and later to the ice sheets of Greenland and Antarctica.

EISMINT
The EISMINT experiment phase two (Payne et al., 2000) is a simplified geometry experiment with regular boundary conditions to compare thermomechanical ice-sheet models.All boundary conditions are symmetrical and time independent.Twelve experiments were defined and the experiment K is used here to test the transport model SICOTRACE and the post-processing tool SICOSTRAT.In experiment K the bedrock consists of a regular array of 500 m high mounds and with zero ice initial condition.
This EISMINT setup uses a 25 km×25 km horizontal grid in the model domain of the size 1500 km × 1500 km and for the semi-Lagrangian grid we use 100 vertical layers for ζ s .The time step in the SICOPOLIS simulation was 200 yr.The accumulation/ablation rate is a function of geographical position, which changes its sign in a given distance from the summit (for details see Payne et al., 2000).Figure 6 shows the steady-state ice surface at time 300 kyr.In addition, two ice-core locations C1 and C2 as well as one cross section are marked.The core C1 is located at the ice divide (x = 750 km, y = 750 km) and the results are shown in Fig. 11.Core C2 is located at the border of the region with positive mass balance at x = y = 1000 km, and the results of the isotope modelling are given in Fig. 12.For the cross section the two depositional coordinates, the depositional age and the resulting δ 18 O distributions are shown in Figs.7-10.

Discussion for EISMINT
The simulated mass transport in the x direction at the cross section is small compared to the one in the y direction, which is expected due to the general ice flow from the centre to the outer regions for the given geometry.This can also be seen in Fig. 8 for the depositional provenance of y, changing from central deposition to radial origin.
In Fig. 8 the whole ice body is yellow which means that there is little transport perpendicular to the cross section.Inside the ice body the contour lines indicate the depositional coordinate x = 750 km, which is a combination of topographical effect due to the mounts and numerical variations in the order of less than 10 km.
At the area with positive mass balance the oldest ice can be found close to the ice sheet base (Fig. 9).This is the ice deposited during the initial phase of the build-up process which remains located there for the whole time.It can also be seen in Fig. 8 where the ice has nearly the same depositional coordinate as the actual coordinate.Since the surface temperature is constant over time, the δ 18 O ratio is only influenced by the depositional coordinates and not by the age.This can be seen in the cross section Fig. 10 and the ice cores Figs.11d and 12d.Especially in the ice cores, the δ 18 O value is a mere combination of the depositional coordinates.This is due to the symmetry of the surface temperature with its lowest value at the centre of the ice sheet.The stepped behaviour of the profile at the upper levels are due to the step size of the vertical ζ s coordinate used for semi-Lagrangian transport and are also influenced by the used time step.This behaviour disappears with shorter time step but for our experiments we chose a longer one in order to save computational time.

Comparison of backtracking schemes
Figure 13 illustrates the difference between the first-order backtracking method as used in the papers based on Clarke and Marshall (2002) and the backtracking scheme from de Almeida et al. (2009).The figure is based on the EIS-MINT experiment with flat topography (experiment A in Payne et al., 2000).We chose the flat topography because the effects of the different backtracking schemes are more obvious with flat ground than with the mound topography.The apparent feature in the figure is that the differences  are close to zero in the accumulation zone, which is within a radial distance of 450 km from the summit (x = y = 750 km).However in the ablation zone the deviations between the backtracking schemes are substantial.In general the first order scheme calculates older ice in the ablation zone, as can be seen in additional plots in the supplement.While in the previous studies the δ 18 O value in the ablation zone was of no particular interest, the values close to the margin are crucial in order to close the hydrological cycle.
The backtracking scheme from de Almeida et al. ( 2009) is numerically more expensive than the first-order backtracking, mainly because of the required iterations and interpolations of the de Almeida et al. scheme.
These deviations between the backtracking methods in the ablation zone are likely associated with the greater velocity gradients near the margin.In addition, during ice sheet buildup the velocities vary more in the ablation zone, which is better handled by the two-level time scheme with second-order accuracy.Studies of atmospheric models by Staniforth and Pudykiewicz (1985) and McGregor (1993) found that firstorder schemes are inaccurate for large Courant numbers and exhibit poor conservation properties and that a first-order scheme with straight lines and velocities taken at the end point produces an error of 4 % each time step for trajectories in a solid-body rotation problem.

Greenland and Antarctica
In order to apply the isotope transport model to real geometries we present some simulation results of the Greenland and Antarctic ice sheets.Both simulations start at 422 kyr before the present (pre-industrial at 1950) with no initial ice and relaxed bedrock (in respect to glacio-isostatic adjustment).The horizontal grid is 20 km × 20 km for Greenland and 40 km × 40 km for Antarctica.This leads to 83 × 141 × 101 grid points in the semi-Lagrangian grid for Greenland and .
141 × 121 × 101 points for Antarctica, respectively.We use 11 levels in the bedrock and temperate ice domain and 81 in the cold ice domain.
A glacial index gi(t) is used to vary the air temperature and precipitation distribution by interpolating between the present and the last glacial maximum (LGM) conditions.This index is defined gi = 1 for conditions at the LGM and gi = 0 for present conditions (Forsstrom et al., 2003) and is based on data derived from the δ 18 O GRIP ice core record from Greenland and from the δD Vostok ice core record from Antarctica (Dansgaard et al., 1993;Johnsen et al., 1995;Petit et al., 1999).For the first 100 kyr BP the GRIP record is used and prior to that the Vostok ice core.This is necessary since the GRIP record is believed to be corrupted due to iceflow irregularities (Greve, 2005).In total, the glacial index reaches 422 kyr back in time.
With these settings the total computational time including the semi-Lagrangian transport and writing of the depositional archive on one core of a 2.8 GHz dual-core AMD Opteron was 141.3 h for Antarctica and 214.7 h for Greenland.
As an example of a simulated Greenland ice core the GRIP core is shown in Fig. 15a which is located in the central region of the Greenland Ice Sheet (see Fig. 14a for its location).For Antarctica the Vostok Station has been chosen, which is located in central East Antarctic Ice Sheet (see Fig. 14b).The modelled δ 18 O depth profiles are compared to observational data by Johnsen et al. (1997) for the GRIP ice core and by Petit et al. (1999) for the Vostok ice core.Since there is a difference between the modelled and the observed ice thickness we defined level 0 to be the modelled height of the real ice surface.
Figure 16a shows two cross sections through the (a) Greenland and (b) Antarctic ice sheets and the vertical δ 18 O distribution close to the chosen ice-core locations.The Antarctic cross section reveals a broader range of δ 18 O variation with more depleted values than the Greenland cross section.Hence, the colour bars are chosen differently.

Discussion for Greenland and Antarctica
The comparison between the simulated cores and observational data shows in general a good agreement of the isotope records.Taking into consideration that the core data stems from a single ice core with high vertical and horizontal resolution (on the cm scale) and the simulated core is based on 20 km × 20 km and 40 km × 40 km model simulations for Greenland and Antarctica, respectively, the overall coincidence is satisfying.Looking into more detail, however, the observational data shows more high-frequency variability and a shift in the main signals for the Antarctic record (Fig. 15).There are several reasons which need to be discussed when comparing the records.
Firstly, the ice dynamic model SICOPOLIS was taken as a given tool and no effort was made to tune the model for the present day.This was beyond the scope of the study and we used the standard setting here.As a second argument, the isotope boundary condition comes into play which is taken from transfer functions and is, therefore, not a correct local function.While this works properly for Greenland, a mismatch for the Vostok location of about 8 per mille can be found between observations and model results.Third, the time step for the Greenland simulation is five years and for Antarctica is ten years while the measurements resolve the seasonal scale.In addition, the glacial index of the surface temperature forcing is also smoothed to 100 yr so the small-scale variations cannot be resolved.
On the other hand, the overall variability of δ 18 O is comparable with the measurements, while the absolute position of the spikes is also influenced by the difference of the ice thickness.In the cross section for the Greenland Ice Sheet the sequences of values between −40 per mille (green) and −35 per mille (yellow) indicates the different glacial-interglacial cycles of the past, also seen in the lower section of the GRIP ice core.The profile is just a cross section near the GRIP  location but SICOSTRAT generates a complete netCDF output of the whole three-dimensional isotope field.
For the Vostok ice core (Fig. 15b) the modelled values show the same variations as the measurements but the values are generally higher and the features are up to 500 m deeper than in the measurements.Here, the ice dynamics need to be adjusted and more care has to be given to the ice dynamic model performance.The very low accumulation rates in central East Antarctica are challenging to the ice dynamics and lead to too low signals in the simulation (e.g. at 2000 m depth).In addition, the 40 km grid does not resolve details in the bedrock topography, which has a strong influence on the stratigraphy of the Vostok ice core (Parrenin et al., 2004).

Conclusions
In this paper, an oxygen isotope transport model for polythermal ice sheets has been presented, which makes it possible to study the oxygen isotope ratio inside polythermal ice sheets where the shallow ice approximation is valid.The model was applied to the EISMINT inter-comparison phase 2 project and applied to the Greenland and Antarctic ice sheets.As an example, one simulated ice core of each ice sheet has been compared to measured ice core data.The model has been developed to tackle two goals -firstly, to use local comparison of vertical profiles at drill core sites to validate and improve the ice sheet model, and secondly, to close the oxygen isotopic budget in the hydrological cycle of a fully coupled Earth system model approach, when coupling the ice model to an atmosphere-ocean-land surface general circulation model (e.g.Werner et al., 2011;Xu et al., 2012).
In our study the backtracking scheme from de Almeida et al. ( 2009) with second-order accuracy is used which requires the three-dimensional velocity fields to be interpolated on intermediate steps.The sigma coordinate formulation of the ice sheet model make it necessary to evaluate and interpolate ice topography variables in space and time.This adds to the computational costs and overall complexity of the approach.
If the focus of a study with semi-Lagrangian transport is on ice core locations in the ice sheet interior a first-order backtracking may be sufficiently accurate, with the additional benefit of lower computational lost.Otherwise, if the values of the transported property near the margin are of interest, a second-order backtracking scheme such as the one of de Almeida et al. ( 2009) may be better suited for the task.We will address the comparison between different semi-Lagrangian schemes as well as Lagrangian and Eulerian methods focusing on values near the margin in a subsequent paper.
The indirect semi-Lagrangian approach with provenance transport and three different programs has some advantages and disadvantages.
The advantages are that transfer function between surface temperature and δ 18 O can be changed without the need of re-running the whole simulation again.This is important for long-term simulation since SICOPOLIS and SICO-TRACE are computationally expensive and, therefore, it allows to experiment with different transfer functions.The transfer function for δ 18 O could also include local changes of surface elevation and changes in mean surface temperature (Cuffey, 2000;Langebroek et al., 2010).A downside of the three different programs is that SICOTRACE needs the three-dimensional velocity field in the cold and temperate ice domain and information about the ice topography for each time step.For high resolution and simulations on palaeo-timescales this leads to a high amount of data which could be of the order of terabytes.Since the information of the velocity field and the whole history of the ice sheet evolution is in general of no interest, the data can be deleted after SICOTRACE has performed the transport calculation.For paleo-runs, a shell script subsequently runs SICOPOLIS and SICOTRACE to generate a depositional archive and provenance archive for the whole time period for SICOSTRAT.This leads to a further complication of file handling and an additional layer of complexity in the whole workflow.
The semi-Lagrangian transport can be used for other species, such as deuterium, but if the species involve some feedback with the ice dynamics, then the approach of three different programs is not feasible because of all the overhead with file transfer and initialisation of programs only to run for one time step.
This chosen approach is in any case an important step in the direction of fully coupled Earth System Models for investigating the climate system and comparing model output and in situ measurements.It is also the basis for studies involving the transport of passive tracers.

Figure 1 .
Figure 1.Flow chart of the three programs included in ISOPOLIS: SICOPOLIS -the polythermal ice sheet model.SICOTRACE the tracing program which uses the output from SICOPOLIS and calculates the semi-Lagrangian tracer transport method.SICOSTRAT is the plotting routine for the construct of the stratigraphy.

Figure 2 .
Figure 2. Terrain-following sigma transformation in the polythermal ice sheet.If temperate ice is present it is always below the cold ice.Hence, there are two separate sigma-transformed regular grids above each other, and they share their horizontal coordinate x → ξ and y → η (not shown here).The grid spacing ζ c and ζ t for cold and temperate are usually different and also the sigma-transformation is different for both domains (seeGreve, 1997a).CTS denotes the cold-temperate surface.

Figure 3 .
Figure 3. Basic principle of the semi-Lagrangian method explained in the cold ice grid.If the value of (t − t) at the departure point (open circle) is known, the value (t) at the desired grid point (i, j, k c ) can be calculated, where i and j are the indices of the horizontal grid and k c the one for the vertical sigma coordinates in cold ice.

Figure 4 .
Figure 4. Semi-Lagrangian sigma-transformed grid.The original cold and temperate grid on the left are both covered by the semi-Lagrangian grid.

Figure 5 .
Figure 5. Schematic for the two-step three time level scheme for one dimension (ξ ).The red solid curve is the actual trajectory of a particle and the dashed line is the approximation.a * is the displacement to the intermediate position where the velocity is evaluated and a * * is the displacement of the second step.

Figure 6 .
Figure 6.Steady-state surface elevation of the EISMINT experiment K at simulation time 300 kyr.The two crosses mark the core locations C1 and C2, and the dashed red line indicates the cross section.This cross section cuts through the ice centre dome where core C1 is located.

Figure 7 .
Figure 7. Cross section of the depositional provenance of x at time 300 kyr with the mound topography and the dashed contour of the temperate ice thickness H t .The colour bar indicates the origin of the ice where blue for example means that the ice is coming from the x coordinate of 0 km.

Figure 8 .
Figure 8. Cross section of the depositional provenance of y at time 300 kyr.The same colour coding as in Fig. 7 is applied here but for the y coordinate.

Figure 9 .Fig. 10 .Figure 10 .
Figure 9. Age of the ice derived by calculating t final −t d where t final is the final time of the simulation.The older ice is located near the base and the majority of the ice is younger than 30 kyr.

Figure 13 .
Figure 13.Difference of calculated age between first-order and de Almeida et al. (2009) backtracking with the EIS-MINT experiment A with flat topography.The colours indicate the vertical sum over all layers in the semi-Lagrangian grid of the absolute value of the differences at each grid cell: ksmax ks=1 first order(ks)−de Almeida et al.(ks) first order(ks)

Figure 14 .
Figure 14.Simulated present-day surface topography with 200 m contours of elevation in km a.s.l.Major ice core locations are marked with crosses and the location of the cross sections are marked with a red dashed line.(a) Greenland and (b) Antarctica, without the ice shelves because SICOPOLIS 2.9 does not include ice shelf dynamics.

Figure 15 .
Figure 15.Comparison of modelled δ 18 O values in black and observed values in blue.The modelled data in the GRIP core (a) have an offset of 10 per mille to the right in order to make the comparison more easy to read.The vertical axis is the ice core depth and level 0 is defined as the ice surface in both, measured and modelled ice core.For the GRIP core (a) the modelled depth is 2945 compared to observed 3029 m with reliable δ 18 O data until 2983.2m, and for the Vostok core (b) the modelled depth is 3549 m compared to 3623 m in reality and data are available until 3310 m.

Figure 16 .
Figure 16.Two cross sections of the (a) Greenland and (b) Antarctic ice sheets showing the δ 18 O distribution close to the GRIP and Vostok ice cores.