1. Introduction
Interpretation of palaeoclimatic records in ice sheets requires accurate dating. Conventional stratigraphic techniques fail below the depth where annual signals in chemical and physical properties become indistinguishable due to layer thinning and diffusion (e.g. Reference Reeh, Oeschger and LangwayReeh, 1989; Reference JohnsenJohnsen and others, 2001). Here models present the only recourse for computing age–depth profiles. General knowledge of the age distribution in large ice sheets is also required to incorporate physical characteristics of the ice in flow laws which depend on the age of the ice or its strain history. An example is the hardness difference between Holocene and Wisconsin ice found in the Greenland ice sheet, which has to be taken into account to model thinning rates in the basal layers correctly (Reference HuybrechtsHuybrechts, 1994).
The dating problem can be solved analytically only for a limited number of cases where the horizontal advection of ice can be neglected or prescribed in terms of a simple model. Other limitations include the difficulty of dealing with temporal variations in accumulation rate, flow pattern, ice thickness, temperature and ice rheology (Reference Reeh, Oeschger and LangwayReeh, 1989). Elaborate three-dimensional ice-sheet models are needed to deal with these parameters. In such flow models, the age distribution is derived from the numerical solution of a purely advective equation. There are essentially two ways to solve such an equation. The Eulerian method considers age as a transportable quantity by solving the advective equation in a frame fixed in space. Due to the nature of the equation, straightforward finite-difference schemes of the type central in space are known to be unconditionally unstable, requiring the use of diffusive schemes. First attempts to model the age distribution in large ice-sheet models added an artificial diffusion term to the numerical solution to damp the instabilities (Reference HuybrechtsHuybrechts, 1994; Reference GreveGreve, 1997; Reference Greve, Mugge, Baral, Albrecht, Savvin, Hutter, Wang and BeerGreve and others, 1999). However, this method is not very accurate in the basal layers, and reliable solutions can only be obtained for the upper half of an ice sheet. In recent work, Reference Greve, Wang and MuggeGreve and others (2002) have shown that for the one-dimensional case a second-order upwinding scheme for the vertical advection term is stable and yields an accuracy of 499% for a grid in excess of 60 points.
The alternative way to deal with the age equation is to consider a frame following the motion, which is a numerically diffusion-free scheme. In the Lagrangian approach, the path of a particle is tracked through the numerical grid as the ice sheet evolves through time. Its application in a three-dimensional numerical ice-sheet model was discussed by Reference Clarke and S. J.Clarke and Marshall (2002). They use the term“marker” instead of “tracer”, and each particle has its own“birth certificate”, since it can carry information on any conservative characteristic (e.g. time of deposition, place of deposition, water isotopes). Moreover, the Lagrangian approach makes it possible to reverse time and to carry out backward tracing. The Lagrangian approach discussed in this paper is distinct from the semi-Lagrangian technique as widely used in, for instance, meteorological applications. In a semi-Lagrangian approach, the solution is provided on a regular grid, but this is achieved by using a different set of particles at each time-step, the set of particles being chosen such that they arrive exactly at the points of a Cartesian mesh at the end of the time-step (Reference Bates and McdonaldBates and McDonald, 1982; Reference Staniforth and CoteStaniforth and Cote, 1991). The main advantage of a semi-Lagrangian numerical scheme over an Eulerian one is the greatly improved stability properties, because it makes it possible to overcome the Courant–Friedrichs– Levy (CFL) restriction on time-steps (Reference Machenhauer, Olk, Lin, Laprise and RitchieMachenhauer and Olk, 1997). It is therefore often used for problems with large Courant numbers, where the product of velocity and time-step is significant compared to grid size. In ice sheets, however, Courant numbers are much smaller, and such large time-steps are not needed as they would imply an unacceptable loss of accuracy. In addition, the semi-Lagrangian method does not allow reconstruction of the trajectory of an individual particle.
The objective of this paper is to assess the merits and shortcomings of both the Eulerian and Lagrangian methods and to perform a detailed quantitative comparison of their (relative) accuracy. In section 2 we give a general outline of both numerical approaches. In section 3 we consider a series of simple schematic two-dimensional ice-sheet models in the vertical plane and describe the results of a comparison of the date distribution under the different model set-ups. The ice divide is a special case, for which a perfect analytical solution is available to compare with the numerical results. We also test a linear and a cubic-spline interpolation method for the Lagrangian scheme. Problems with the application of both interpolation methods are discussed with respect to the horizontal velocity field. In section 4 we apply the Eulerian and Lagrangian methods in a three-dimensional Antarctic ice-sheet model. For this purpose, we consider an experiment under present-day climatic conditions. Results of the comparison of the two approaches, and the main conclusions follow in section 5. In this paper, we limit ourselves to stationary problems, as the emphasis is on investigating basic differences between the two approaches for the simplest cases. Time-dependent applications will be discussed in future work.
2. Numerical Dating Methods
Hereafter we shall understand under the quantity A the absolute moment in time when a particle was deposited at the surface. This quantity is negative in the past. Since all results discussed in the paper are for integrations ending at the present day, these dates of deposition can also be interpreted as negative ages BP (before present), in accordance with the more usual definitions given in the literature.
2.1. Eulerian approach
We consider a Cartesian coordinate system with x, y the horizontal plane, z the vertical, and t time. The fundamental equation states that date of deposition A is a quantity that is conserved with time:
where d/dt is the full (total) time derivative. In an Eulerian frame, this equation can be rewritten with a partial (local) time derivative as:
where vx, vy, vz are components of the three-dimensional velocity field. Boundary conditions are A(x, y, z = H + h, t) = t at the surface, with H ice thickness and h bed elevation. For an ice sheet frozen to the bottom, the date of deposition is undefined at the base (or physically equal to the time the first permanent ice cover settled at that particular location). In the case of basal melting, a flux boundary condition can be set at the base as follows:
where B is the basal melting rate (m a–1), positive in the case of basal melting. In all runs discussed in this paper, we use second-order upwinding for all the advection terms. Since we only consider cases with basal melting, implying strict downward motion perpendicular to the ice layers, the basal boundary condition on A actually stands on its own and is never invoked to compute the date of deposition elsewhere. This prevents an arbitrary choice of basal ice age from polluting the solution at another location in the ice sheet. The case with negative B (basal freeze-on) cannot be handled with Equation (3) and is excluded from the analysis in this paper.
2.2. Lagrangian approach
In a Lagrangian frame which follows the particles, Equation (1) can be discretized in a straightforward way without any artificial diffusion. The approach consists of integration of a vector evolution equation:
where X(t) is the position vector of the ice particle and V(x, y, z, t) is the velocity field in the ice sheet. The initial condition for Equation (4) states that the trajectory of an ice particle deposited at the surface starts at time t0 = A at the free surface H + h:
A trajectory ends when a particle crosses an arbitrary conventional bottom (we use 99% of relative depth) or leaves the ice sheet through its lateral margins.
The general problem in the Lagrangian approach consists of proper numerical interpolation of the velocity field. If piecewise linear interpolation is used, the velocity field is computed as a weighted sum of the field values in the surrounding gridpoints (four in two-dimensional space, eight in three-dimensional space). Cubic-spline interpolation requires knowledge of the velocity values in all or a subset of the model domain (Reference Press, Teukolsky, Vetterling and FlanneryPress and others, 1992).
3. Schematic Tests in a Two-Dimensional Vertical Plane
First we assume a series of schematic experiments representing ideal cases to compare the Eulerian and Lagrangian approaches. The first of these experiments is for a flat bedrock under a steady state (standard model). Further schematic experiments test the effects of basal sliding, basal melting and an undulating bedrock.
3.1. Design of the standard experiment
We consider a two-dimensional ice sheet on a flat bedrock at z = 0 for which the surface profile can be given as:
where
This is the well-known Nye–Vialov solution (Reference PatersonPaterson, 1994), in which H is ice thickness (m), equated here to surface elevation, H0 is the height of the ice sheet at the divide (m), L is the half-width of the profile (106 m), A is the rate factor of the flow law (10–16 Pa–3 a–1), n is the flow-law exponent (taken to be 3), p is ice density (910 kg m–3), g is acceleration of gravity (9.81 m s–2) and M is the mass balance (m a–1). With a typical value of M = 0.1 ma–1, corresponding to a large polar ice sheet, we obtain a value for H0 of 3598.4 m (see Table 1).
Table 1. Parameter values for the various schematic tests in a two-dimensional plane
The horizontal velocity field as a function of depth can be derived along the flowline from Equations (6) and (7) to yield the following expression:
For the vertical velocity vz a simple linear expression is used. This expression is not consistent with vx as provided by Equation (8) but is used for simplicity because it allows for a straightforward analytical comparison at the divide:
Equations (6–9) entirely define a well-circumscribed problem for the purpose of comparing the dating methods. In practice, we use a dimensionless vertical coordinate scaled to ice thickness, and transform Equations (8) and (9) to this new system. Both the Lagrangian and Eulerian calculations take place directly in the transformed system.
The domain of the ice sheet given by Equations (6) and (7) is discretized by 101 6 101 gridpoints, equally spaced in dimensionless space. The experiments run for 500 kyr with a time-step of 1 year. Since the ice sheet is symmetrical with respect to the ice divide, we consider only the righthand domain, which consists of horizontal gridpoints 51–101. The Lagrangian method uses cubic-spline interpolation for the position integration. The results are shown in Figure 1. They are only valid for the upper 99% of ice thickness and between the divide and x = 980 km, because both boundaries are singular for one or both of the methods.
Fig. 1. Results for the standard model run with flat bedrock, no basal melting and no basal sliding (experiment A). The streamlines are a direct output of the Lagrangian method and coincide with particle paths for a stationary state. The lower panels show the calculated dates of deposition from the Lagrangian and Eulerian methods, as well as their relative difference. Results are shown in the square that arises from the dimensionless vertical coordinate scaled to ice thickness.
In steady state, particle trajectories are equivalent to streamlines. These are shown in the top panel of Figure 1. Only 6 tracers out of 51 reach a relative depth of 99 %, while the others finish at the margin. This immediately demonstrates a basic problem with the Lagrangian method: in the deeper layers of the sheet, it yields less information than the Eulerian method, leading to an inherent loss of detail with depth. In our comparison, we solve this problem by launching 45 additional tracers between gridpoints 51 and 56, giving a density of 51 markers at 99 % depth, comparable to the Eulerian method. For clarity, however, these are not shown on the streamline plot.
On the scale of the plots, the two dating methods are almost indistinguishable. We therefore define a relative difference on date of deposition (“Eulerian error”) as follows:
where AE is Eulerian date and AL Lagrangian date. Negative relative differences indicate that the Lagrangian ice is older (i.e. the date of deposition is more negative) than the Eulerian ice for the same position. The relative difference between the two methods is very small (< 0.1 %) for most of the ice sheet, except at the model boundaries. The probable reasons for these differences are discussed in section 3.4 further below.
3.2. Comparison with an analytical solution at the divide
For the simple stationary case of no bottom melting, no horizontal advection, and constant vertical strain rate, a simple logarithmic expression for the date-of-deposition profile with depth can be derived, corresponding to our standard experiment at the divide:
where G is the surface accumulation rate, here equal to the mass balance M. This analytical derivation is due to Reference NyeNye (1963) and Reference HaefeliHaefeli (1963), and is known as the Nye–Haefeli time-scale.
We compare the analytical result at the divide with the numerical solution of the standard experiment A. In this case, the relative error of the numerical method is defined as:
where Anum is the date of deposition computed by either numerical method (Eulerian or Lagrangian), and Aan is given by Equation (11). The result is displayed in Figure 2. Clearly, the Lagrangian method produces less error, as is evident from the logarithmic scale required to show the relative errors. All along the vertical, the Lagrangian method slightly underestimates the date of deposition (i.e. produces an older date) but with a relative error always below 0.1% and often as good as 0.001 %. This error depends on the time-step and the grid size and is due to the accumulation of errors during the subsequent position interpolations of the ice particle travelling downwards. Relative to the date, this error becomes smaller the smaller the vertical velocity becomes nearer to the bottom. The error on the Eulerian calculations is also acceptably small, but has a different profile with depth. Except for the upper two vertical gridpoints, where second-to order upwinding cannot be used and a two-point upwinding scheme is only accurate to first order, the relative error is always < 0.1 % down to a relative depth of approximately 93%. For most of this part, Enum is negative, that is, the Eulerian date is more recent than the analytical solution, but Enum changes sign at 8 4% depth, when it becomes positive. Near to the bottom, the Eulerian error increases rather abruptly to a substantial level of several per cent. This originates from an inability of the numerical scheme to produce a precise result close to the lower boundary. There, the accuracy of the finite-difference scheme changes from second to first order, and gradients also become very large.
Fig. 2. The error on date of deposition for both methods relative the analytical solution at the ice divide. The upper logarithmic axis is for the full line (Lagrangian), and the lower linear axis refers to the stippled line (Eulerian).
Notes: (U) is the vertically averaged horizontal velocity corresponding to experiment A. The last column is for a parameter that gives the amplitude of bedrock undulations in Equation (13).
3.3. Design of additional schematic tests
To test more challenging conditions more like real world-conditions, we modified the standard model set-up to investigate the effects of basal sliding, basal melting and an undulating bedrock. The modifications to the standard model are listed in Table 1. I n all of the additional experiments, the requirement was that the vertically integrated ice flux along the flowline was identical to that in the standard model. Also the surface profile was retained from the standard experiment. Although these additional tests violate basic principles underlying the Nye–Vialov solution, this is unproblematic as long as identical velocity and ice-thickness fields are used to compare the two dating methods. The basic design of the additional four experiments can be derived from the first column of streamline plots in Figure 3.
Fig. 3. Results for the additional tests performed for a schematic two-dimensional stationary ice sheet showing streamlines, Lagrangian date and the relative difference between Lagrangian and Eulerian dates of deposition. The bold letters in the upper right corner refer to the specific experiment: B, basal melting; C, basal sliding; D, undulating bedrock; E, undulating bedrock and basal sliding (cf. Table 1 and the text for more details).
To meet the requirement of flux conservation, the basal melting experiment (experiment B) has to increase the surface accumulation by the same amount as the basal melting rate of B = 0.02 m a– 1 . The major effect of this modification is a non-zero vertical velocity at the base, resulting in more recent ice in the basal layers as compared to the standard experiment. By contrast, the introduction of basal sliding (experiment C) causes a flattening of the trajectories to resemble hyperbolas such that only the tracer starting at the divide reaches a depth of 99%, whereas the other tracers leave at the margin at x = 980 km. In this case, a 50–50 partitioning was put forward between horizontal movement due to deformation and due to basal sliding. To investigate the effect of a non-flat bed, an undulating bedrock profile (experiment D) was obtained as follows:
where c is a constant between 0 and 1 that gives the amplitude of the waves as a fraction of ice thickness H(x) in the standard experiment. To preserve the flux, the horizontal velocities are now scaled to the difference between the surface profile of the standard experiment and the revised bed elevation. Particles now follow an undulating pattern, with higher speeds over thinner ice and lower speeds over thicker ice, compared to the standard experiment. Finally, experiment E combines both an undulating bedrock and basal sliding.
3.4. Diversity in ice-age patterns produced by the two dating methods
From close inspection of Figures 1and 3 a number of general inferences can be drawn. The overriding conclusion is that the difference between the Lagrangian and Eulerian methods is very small. In general, the relative difference does not exceed 0.5% for most of the ice sheet. Larger relative differences only occur at the top and basal boundaries, and near to the lateral margin. This pattern is qualitatively similar for all schematic experiments. The reason for these systematic differences originates from the different methods and different sources of numerical errors. We cannot assess which method is closer to perfection, since no exact solution exists except at the divide. Nevertheless, the relative comparisons permit some conclusions.
It is reasonable to assume that the difference between the two methods near to the base and at the top is due to imperfections in the Eulerian scheme, since that is what the divide comparison brought to light. However, near the lateral margin (righthand margin of Fig. 4) the difference could also be attributed to the Lagrangian scheme due to the cumulative nature of its errors. In the Lagrangian approach, errors are proportional to travel time and to the velocity magnitude. Therefore the largest error will be at the end of the trajectory of a particle at the margin, where both quantities are largest.
Fig. 4. (a) Relative error of interpolation (%) of the linear and the cubic-spline methods for the horizontal velocity in the experiment. (b) Cumulative error on date of deposition using linear interpolation in the Lagrangian scheme at a position from the divide.
3.5. Problems of interpolation in the Lagrangian approach
The problem of applying the Lagrangian approach and comparing it with the Eulerian method is fundamentally a problem of interpolation. Since the velocity field is prescribed on a regular grid, and a particle’s position needs to be obtained in any point during a model run, one has to interpolate the velocity field between gridpoints to obtain spatial increments of a particle’s position. Further, since the position of a particle is computed step by step, these errors are cumulative. The second problem arises from the comparison itself and is due to uneven spatial coverage of tracers. The only solution to the latter problem is a higher abundance of initial tracers. A preliminary analysis of the flow is therefore necessary to decide at which location, and at which frequency in time-dependent applications, additional tracers need to be launched.
Concerning the first problem, it is clear that bi-cubicspline interpolation is better than piecewise bilinear interpolation. This is demonstrated in Figure 4. In Figure 4a we compute the modulus of the relative error of both methods of interpolation for the horizontal velocity with respect to the analytical result as given by Equation (8), and average the result over the vertical coordinate. That provides a good measure of interpolation errors. It can be seen that spline interpolation of the horizontal velocity is three to as much as seven orders of magnitude more accurate than linear interpolation. In our schematic set-up, the interpolation error on velocity almost exclusively originates from the horizontal velocity component because the vertical velocity is linear, and therefore both interpolation methods give similar results for the vertical to within computer accuracy. Based on the results shown in Figure 4a, we can safely consider the spline interpolation results as the reference state to assess the accuracy of the linear interpolation method.
We find that linear interpolation generally leads to underestimation of the velocity and therefore to smaller displacements and a calculated date of deposition that is earlier as compared to the spline method. This fact is illustrated in Figure 4b, where the cumulative error on date is shown for the transect at the margin at x = 980 km. Still, the absolute value of accumulated error on date of deposition of approximately 1200 years in the lower layers corresponds to a relative error of 51.5%. Similar results were obtained at other locations closer to the divide, giving a typical relative error on date by linear interpolation of about 1% at 90% depth.
The difference between the two methods of date estimation is therefore small, and linear interpolation can be expected to yield acceptable results. This is the more so in view of the large central processing unit (CPU) time required for spline interpolation whose magnitude makes it hardly applicable in a three-dimensional ice-sheet model. The main reason is that spline interpolation is not a local interpolation technique but requires calculation of two partial derivatives at the tracer’s position at every time-step (Reference Press, Teukolsky, Vetterling and FlanneryPress and others, 1992). In three-dimensional applications, the additional CPU cost is at least two orders of magnitude larger than for linear interpolation.
4. Application in a Three-Dimensional Antarctic Ice-Sheet Model
4.1. Experimental set-up
To assess the applicability of the two dating methods under real conditions, we implemented the age-calculation routines in a three-dimensional thermomechanical model of the Antarctic ice sheet (Reference HuybrechtsHuybrechts and de Wolde, 1999; Reference Huybrechts and de WoldeHuybrechts, 2002). This model is time-dependent, calculates the coupled velocity and temperature fields, and has components describing basal sliding, basal melting, the surface climate, isostatic adjustment, and the flow in a coupled ice shelf. The horizontal resolution is 20 km and there are 31 layers in the vertical, gradually decreasing in thickness towards the base. The model employs a dimensionless vertical coordinate scaled to ice thickness with grid spacings of 5% of ice thickness at the surface and 1% at the base.
The model was run for stationary present-day climatic conditions. The only constraint compared to recent applications of the model is that the grounding line was fixed to its current position. Ice thickness is, however, still a free variable allowed to establish an equilibrium with the climate forcing and the flow law. Consequently, the vertical velocity is obtained from the continuity equation and retains full consistency with the horizontal velocity components and the surface mass balance, which is positive over the entire model domain. Basal melting and basal sliding occur in the model whenever basal temperature reaches the pressure-melting point. The Antarctic experiment therefore combines all of the features discussed in the schematic experiments above, with the notable exception that vertical velocity is no longer linear.
For the Eulerian ice-date computations, the same second-order upwinding scheme was used as in the schematic experiments, and, for the Lagrangian method, piecewise linear interpolation was applied for tracking tracer position in between gridpoints. The calculations took place directly within the transformed coordinate system, ensuring that vertical movement is always directed downwards over the entire model domain, thus avoiding basal boundary-condition problems in the upwinding advection scheme. For the Lagrangian computations, tracers were launched only once at the surface at the initial time for every grounded gridpoint. The total number of tracers therefore equalled approximately 30500 out of a total number of 281×281=78961 horizontal gridpoints. The total integration time was 500 kyr for both methods, with a time-step of 5 years. This is a long enough period to ensure that a stationary solution is established for the advection equation and that all tracers that have remained within the grounded ice sheet are well below 90% depth. For the Eulerian calculations, the initial condition on date of deposition consisted of the analytical result given by Equation (11). To facilitate the experiment, the Lagrangian run was actually run in offline mode with the velocity output from the Antarctic model for the present time.
4.2. Results
Figure 5 shows an example of the resulting distributions of the date of deposition at an intermediate depth of 60%. These results should not be confused with the“real”date distribution in the Antarctic ice sheet since ice thickness and accumulation rate have not been allowed to vary over the full length of the integration period in response to climatic changes (though the accumulation rate has changed in accordance with surface elevation changes relative to the observations). Nevertheless, the pattern is expected to be qualitatively similar to reality, with old ice concentrated in the interior of East Antarctica roughly inversely proportional to the surface accumulation rate, together with much younger ice at comparable relative depth in West Antarctica and at the margin. Also clear are the effects of flow divergence/convergence and of the flow chanelling into the outlet glaciers that drain most of the ice sheet at the margin. In general, the two methods give very similar results, with absolute date differences of maximally ±10 kyr, but well below ±2 kyr over most of the area, with the Lagrangian method giving some what older ice at the margin and somewhat younger ice in the interior (Fig. 5).
Fig. 5. Date of deposition computed in the Antarctic ice sheet for present-day stationary climatic conditions by both methods, as well as their absolute difference. This example is for the layer at 60% depth.
However, the largest differences occur in zones of large gradients and can therefore be explained by diffusive effects in the Eulerian scheme and, probably more importantly, the irregular distribution of tracers that results for deeper layers in the Lagrangian method. This is demonstrated in more detail in Figure 6, marking the positions where tracers have crossed a given depth layer during the course of the integration. At a relative depth of 60%, only half of the original tracers remain in the ice sheet, and the other half have migrated into the surrounding ice shelves or have been washed out in the ocean. At 90% depth, less than a quarter of the tracers remain, and they are distributed unevenly in space. It should, however, be noted that this problem will be less severe when tracers are launched at regular time intervals in time-depend-ent applications with changing climatic input and ice-sheet geometry.
Fig. 6. Distribution of tracer crossings through the 60% and 90% depth layers in the Antarctic experiment, showing how tracers are gradually concentrated and washed out towards the margin in deeper layers. In this experiment, tracers were launched only once for every ice-sheet gridpoint on the 2816281 numerical mesh.
The inverse problem, in which the tracers are used to compute the relative depth of an isochrone surface, does not depend on the rhythm by which tracers are launched at the surface. Figure 7 shows the distribution of the final positions of tracers 100 kyr after they were launched at the surface. In our steady-state experiment, these are the points one would have to use to reconstruct the depth of the 100 kyr isochrone surface. Almost all of the remaining tracers concentrate within interior East Antarctica, with a very sparse distribution at the margin and in West Antarctica, and virtually none in the Antarctic Peninsula area. In this case, the isochrone surface of 100 kyr is marked by only 310 9 tracers out of the 30548 initially inserted at the surface.
Fig. 7. Distribution of tracers at the depth of the isochrone surface of 100 kyr in the Antarctic model, demonstrating the problem of tracer density for interpolating age horizons on a regular numerical grid.
5. Concluding Remarks
A comparison of Lagrangian and Eulerian approaches for ice dating in numerical ice-sheet models has enabled us to establish the advantages and disadvantages of the two methods. From experiments with a schematic ice sheet and an application in a three-dimensional thermomechanical model of the Antarctic ice sheet, our main findings can be summarized as follows:
1. The traditional Eulerian method, in which ice age is estimated by solving an advection equation with a second-order upwinding finite-difference scheme, generally produces larger errors than a Lagrangian method on the same grid, as evident from the divide solution where an exact analytical result exists. Relative differences between the two methods are, however, small (51%) for the largest part of the ice sheet, except near the base, where errors increase exponentially in the Eulerian method.
2. Errors in the Lagrangian computations essentially find their origin in interpolation procedures between grid-points. They accumulate along the trajectory of an ice particle and are proportional to the age of a particle. These errors are fundamentally different from those in the Eulerian approach which arise from numerical diffusive effects. The accuracy problem arising from using linear interpolation of tracer distributions onto a regular numerical grid can be alleviated by a higher frequency of tracer launches.
3. The application of a piecewise linear interpolation algorithm in the Lagrangian approach requires less CPU time than the numerical solution of the advection equation. Interpolation by cubic splines is more precise but is hardly applicable in large-scale ice-sheet models with present-day computers because of the very large CPU times required. In time-dependent applications, the number of data needing to be stored can increase dramatically, but limits on storage are easier to overcome than limitations on physical time.
4. A clear advantage of the Eulerian approach is its ease of implementation and computation. It provides a continuous solution on the same grid as the one on which the ice-velocity components are provided in a numerical model, with only little data storage. The Lagrangian approach inherently leads to loss of information with depth due to the divergence of ice flow. High basal sliding and comparatively low vertical velocity may cause shade effects in the deeper layers where a particle cannot penetrate or their final distribution is very sparse. Preliminary analysis of bottom layers can help to select the optimal amount of tracers that need to be launched from the surface.
5. A disadvantage of the Eulerian approach is that it cannot deal with hiatus in the layer sequence due to surface or bottom melting. The practical application of the Eulerian method is also limited to smooth and continuous fields such as ice age or date of deposition. The Lagrangian computation can in principle provide the distribution of any desired conservative property that is carried with the flow, such as chemical composition and place and elevation of deposition of the ice. Those are the properties that need to be studied in detail in connection with the drilling of deep ice cores.
Acknowledgements
This work is a contribution to the European Project for Ice Coring in Antarctica (Epica III) funded by the European Union under the FP5 Environment and Climate Programme, contract No. EVK2-CT-2000-00077.We also thank N. Lhomme and M. Siegert for their thoughtful comments which markedly improved the quality of the paper.
