Hostname: page-component-7b7567c8c-2rsb6 Total loading time: 0 Render date: 2026-11-12T09:09:49.603Z Has data issue: false hasContentIssue false
  • English
  • Français

Comparison of numerical schemes for the solution of the advective age equation in ice sheets

Published online by Cambridge University Press:  14 September 2017

Ralf Greve,
Yongqi Wang and
Bernd Mügge
Ralf Greve
Affiliation:
Fachbereich Mechanik, Technische Universität Darmstadt, Hochschulstrasse 1, D-64289 Darmstadt, Germany E-mail: greve@mechanik.tu-darmstadt.de
Yongqi Wang
Affiliation:
Fachbereich Mechanik, Technische Universität Darmstadt, Hochschulstrasse 1, D-64289 Darmstadt, Germany E-mail: greve@mechanik.tu-darmstadt.de
Bernd Mügge
Affiliation:
Fachbereich Mechanik, Technische Universität Darmstadt, Hochschulstrasse 1, D-64289 Darmstadt, Germany E-mail: greve@mechanik.tu-darmstadt.de
Rights & Permissions [Opens in a new window]

Abstract

A one-dimensional model problem for computation of the age field in ice sheets, which is of great importance for dating deep ice cores, is considered. the corresponding partial differential equation (PDE) is of purely advective (hyperbolic) type, which is notoriously difficult to solve numerically. By integrating the PDE over a space–time element in the sense of a finite-volume approach, a general difference equation is constructed from which a hierarchy of solution schemes can be derived. Iteration rules are given explicitly for central differences, first-, second- and third-order (QUICK) upstreaming as well as modified TVD Lax–Friedrichs schemes (TVDLFs). the performance of these schemes in terms of convergence and accuracy is discussed. Second-order upstreaming, themodifiedTVDLF scheme with Minmod slope limiter and, with limitations of the accuracy directly at the base, first-order upstreaming prove to be the most suitable for numerical age computations in ice-sheet models.

Information

Type
Research Article
Information
Annals of Glaciology , Volume 35 , 2002, pp. 487 - 494
Copyright
Copyright © the Author(s) [year] 2002

1. Introduction

In Greenland and Antarctica, several deep ice-core projects have been carried out during recent decades (Camp Century, Dye 3, GRIP and GISP2 in Greenland; Vostok and Byrd Station in Antarctica) or are currently under way (NorthGRIP in Greenland; Dronning Maud Land, Dome C and Dome Fuji in Antarctica). These ice cores reveal a wealth of direct and indirect information on palaeoclimatic conditions such as atmospheric composition, surface temperature and snowfall rate, going back some hundreds of thousands of years, and have therefore greatly improved our knowledge of Earth’s climate and its variability during the Pleistocene and Holocene (e.g. Reference DansgaardDansgaard and others, 1993; Reference PetitPetit and others, 1999; Reference JohnsenJohnsen and others, 2001).

In order to obtain time series of these climatic parameters, proper age–depth relations for the ice and the boreholes are required. In the upper parts, these can be established rather precisely by counting annual signals of isotopic composition and impurities downward. However, beyond ages of a few tens of kyr these stratigraphic techniques fail due to layer thinning and diffusion, so flow models are needed to compute age–depth profiles in the lower parts of the boreholes (cf. Reference JohnsenJohnsen and others, 2001, and references therein).

Applications of the three-dimensional dynamic/thermodynamic ice-sheet model SICOPOLIS to compute ages for the Greenland Summit (GRIP/GISP2) region were discussed by Reference GreveGreve (1997) and Greve and others (Reference Greve, Weis and Hutter1998, Reference Greve, Mügge, Baral, Albrecht, Savvin, Hutter, Wang and Beer1999), and for Dronning Maud Land (where a deep ice core is planned within the European Project for Ice Coring in Antarctica (EPICA)) by Reference Calov, Savvin, Greve, Hansen and HutterCalov and others (1998) and Reference Savvin, Greve, Calov, Mügge and HutterSavvin and others (2000). In these studies, the present state of the Greenland and Antarctic ice sheets was obtained via palaeoclimatic simulations over two glacial/interglacial cycles, and the three-dimensional age field was computed by solving the advective age evolution equation along with the thickness, temperature and flow-velocity equations. A major limitation in those calculations was the need for some artificial diffusion in the vertical direction in order to achieve numerical stability, which requires an unphysical boundary condition for the age at the ice base. Reference Calov, Savvin, Greve, Hansen and HutterCalov and others (1998) show that this leads to a basal boundary layer of ~15% of the ice thickness in which the computed ages are affected by this arbitrary boundary condition and therefore erroneous.

In this study, we consider a simple one-dimensional steady-state approximation for the age computation. This simplification allows an analytical solution for the age– depth relation against which numerical methods can be checked. After showing a general method to construct numerical solution schemes, we discuss a variety of different schemes in terms of their accuracy and convergence properties.The objective is to find out which scheme is most suitable for applications in full three-dimensional time-dependent simulations in order to provide age fields as precisely as possible over the whole depth of ice cores.

2. Age Equation In Ice Sheets

2.1. General equation

Let us consider an ice sheet in Cartesian coordinates x, y (which span the horizontal plane), z (vertical) and time t. the equation which describes the evolution of the age field A(x, y, z, t) in the accumulation zone of the ice sheet is then

(1)

where d/dt is the material time derivative which follows the motion of the ice particles with velocity v = (vx, vy, vz), and z = h(x, y, t) is the free surface of the ice sheet where the ice particles settle as snowflakes. In an Eulerian frame,

(2)

where ∂/∂t is the local time derivative. Assuming incompressibility, div v =0, Equation (2) can be written in conservative form as

(3)

This equation is of hyperbolic type with a constant source term and a homogeneous Dirichlet boundary condition.

2.2. One-dimensional steady-state approximation

We now assume steady-state conditions and neglect horizontal advection in Equation (2), which is a coarse approximation of the flow conditions in the vicinity of many ice-core positions in Greenland and Antarctica situated on or close to ice domes or ridges (GRIP, GISP2, Dome C, DML05, etc.). Equation (2) then simplifies to

(4)

or equivalently in conservative form,

(5)

With the scales

(6)

where H is the local ice thickness and as is the net accumulation rate, dimensionless quantities are introduced as

(7)

If we choose z = 0 for the local ice base, then for the free surface z = h = H, or and the dimensionless forms of Equations (4) and (5) become

(8)

In the following, all quantities are taken dimensionless, and tildes are omitted for simplicity of notation.

For the vertical velocity, a Dansgaard–Johnsen type distribution (Reference Dansgaard and Johnsen.Dansgaard andJohnsen, 1969) is assumed, which consists of a constant vertical strain rate ∂vz/∂z from the free surface z = 1 down to a position z = z*, and a linearly decreasing vertical strain rate ∂vz/∂z ∝ z below. With vz(1) = –1 (vertical velocity balances accumulation at the free surface), (small offset at the base in order to avoid an infinite age) and continuity of vz and ∂vz/∂z at z = z*, this yields

(9)

where

(10)

Note that all constants c1–c4 are positive provided and In the following, we will set

(11)

The resulting profile is plotted in Figure 1.

Fig. 1

Analytical profiles of the vertical velocity vz (Equation (9)) and the age A (Equation (12)).

With the velocity distribution (9), an analytic solution of the age Equation (8) 1 can readily be obtained. It reads

(12)

This solution is plotted in Figure 1. At the bottom, A = 20.76, which for the GRIP ice core (H = 3028m, as = 0.23m ice equiv. a–1; Reference Dahl-Jensen, Johnsen, Hammer, Clausen, Jouzel and PeltierDahl-Jensen and others, 1993) corresponds to 273.2 kyr. That age is in reasonable agreement with current estimates of about 250 kyr (Reference DansgaardDansgaard and others, 1993; Reference GreveGreve, 1997).

3. Numerical Schemes

3.1. General considerations

We now turn to the numerical solution of the one-dimensional age Equation (8) 2. In order to investigate schemes that can be adopted for non-steady-state conditions, we do not attempt to solve this equation directly, but iterate

(13)

with the initial condition A(t = 0) = 0 until steady state is reached. By introducing the flux F and source Q as

(14)

Equation (13) can be written in compact fashion,

(15)

Let us introduce a discretization of space and time,

(16)

where Δz and Δt are the grid spacing and the time-step, respectively, and define averaged ages fluxes and sources

(17)

Then integration of Equation (15) over the space–time rectangle in the sense of a finite-volume approach yields the difference equation

(18)

which is still exact.

The bottom gridpoint k = 0 needs to be treated separately. Here,

(19)

and by integration of Equation (15) over [z0, z1/2]×[tn, tn+1],

(20)

Numerical solution schemes can be constructed by assigning approximate values to the averaged quantities in Equations (18) and (20). We will only consider explicit schemes with

(21)

where denote an optional dissipative limiter. Thus, the general iteration rule is

(22)

We apply a numerical grid for which the ages A and sources Q are defined on gridpoints (nodes) k, and the velocities vz on cell boundaries k ± 1/2 (staggered grid, Fig. 2). In addition, the basal velocity (at k = 0) is assumed to be known. Therefore, in Equation (22) the ages do not conform to the grid, and the several numerical schemes derived from Equation (22) differ only in the method of computing the ages from neighbouring ages and in the selection of dissipative limiters

Fig. 2

Staggered numerical grid for the one-dimensional age computation.

3.2. Errors

We assess the error of several solution techniques by comparison with the analytic solution (Equation (12)). Near-basal ice is of particular interest, so we consider the relative error of the basal age,

(23)

where A0 is the steady-state basal age of Equation (22) andis the exact analytic value given by Equation (12). We also introduce the relative error for an arbitrary position in the column,

(24)

Note that, in contrast to does not contain the sign information because of the logarithmic plotting to be applied.

3.3. Central differences

The simplest calculation of

(25)

leads to the forward-time central-space scheme (FTCS)

(26)

which is of first order in time and of second order in space. However, the resulting iteration is always unstable (Reference Morton and Mayers.Morton and Mayers, 1994). Stability can be achieved by extending the time integration over two time-steps (``leap-frog’’),

(27)

which is known as the central-time central-space scheme (CTCS, of second order in time and space). In order to restrict the growth of spurious numerical modes due to decoupling of even and odd iteration steps n, every 101st time-step the CTCS scheme (Equation (27)) is replaced by the FTCS scheme (Equation (26)). This number was found to be sufficient to avoid decoupling. the calculation is also initialized with the FTCS scheme.

3.4. First-order upstream

The CTCS scheme tends to produce numerical oscillations due to the unphysical central discretization of the advection term. At any point in space, the flow of information is from the upstream direction only. to correct for this, non-central differences pointing in the upstream direction can be used. Schemes of that kind are called upstream (upwind) schemes, 0 and the simplest one is constructed by choosing

(28)

For our problem, (Equation (9); Fig. 1), so that

(29)

This is the first- order upstream scheme (UP1), which is of first order in time and space.

By subtracting Equations (29)2 and (26)2, one obtains the difference

which corresponds to an additional diffusion term

(30)

(numerical diffusivity λnum) on the righthand side of the age equation (15), which is discretized by central differences. the numerical diffusion efficiently dampens numerical oscillations, but a disadvantage is that regions with steep gradients within a smooth solution are smeared out over the spatial domain.

3.5. Second-order upstream

An improvement can be made in the numerical diffusion by replacing Equation (28) with a linear upstream interpolation,

(31)

With and central interpolation for the second-order upstream scheme (UP2)

(32)

is produced. This scheme has first-order accuracy in time and second-order accuracy in space.

3.6. QUICK

Another upstream scheme, the Quadratic Upstream Interpolation for Convective Kinematics (QUICK; Reference LeonardLeonard, 1979), is defined by

(33)

As in UP2, with and central interpolation forthe corresponding iteration rule is

(34)

This scheme is of first order in time and of third order in space.

3.7. TVD Lax–Friedrichs

The concept of Total Variation Diminishing (TVD) schemes was introduced by Reference HartenHarten (1983). the main idea of TVD methods is to combine the advantages of accurate high-order schemes and dissipative first-order schemes like UP1. This is achieved by introducing a dissipative limiter which ensures that no spurious oscillations develop near a discontinuity or a zone with steep gradients, and high-order accuracy is retained elsewhere. for instance, the solution can be second- or third-order accurate in the smooth parts, whereas the scheme possesses only first-order accuracy at local extrema or discontinuities. In this way only negligible numerical diffusion is introduced, and advection problems even with discontinuities like travelling shock waves can be well modelled. for more details on TVD and related schemes see, for example, Reference YeeYee (1989), Nessyahu andTadmor (1990) and Reference Jiang, Levy, Lin, Osher and TadmorJiang and others (1998).

TVD schemes make use of a piecewise linear reconstruction of quantities within gridcells. for our problem this means that for a gridcell in place of the single age value as was used in the above schemes (this corresponds to a piecewise constant step reconstruction), a linear reconstruction

(35)

is applied. for the slope limiter

(36)

which yields a weighted average of the left-sided and right-sided gradients, and at the boundaries the one-sided gradients

(37)

are used. the function ϕ(θ) in Equation (36) must satisfy certain conditions in order to yield second-order-accurate cell reconstructions and satisfy the TVD property (Reference SwebySweby, 1984). We consider three limiters. the Superbee limiter

(38)

produces the largest slopes and the smallest numerical diffusion, while the Minmod limiter

(39)

produces the smallest slopes and the largest numerical diffusion. the Woodward limiter

(40)

lies between those extremes. With all three functions, the slope limiter (Equation (36)) is symmetric with respect to the one-sided gradients and and vanishes for or sgn n sgn

A typical TVD method is the TVD Lax–Friedrichs scheme (TVDLF), defined by

(41)

where

(42)

are the values at the cell boundary k+1/2 resulting from reconstruction (35) for the adjacent left (index k) and right (index k + 1) cell, respectively (with either Superbee (TVDLF/S), Minmod (TVDLF/M) orWoodward (TVDLF/ W)), and

(43)

A difficulty with the TVDLF scheme is that the dissipative limiter(41)2, which is chosen in analogy to the usual Lax–Friedrichs scheme, results in rather large numerical diffusion. An alternative we therefore also consider is the modified TVD Lax–Friedrichs scheme (MTVDLF/S,M,W) in which the Courant number is used as a multiplier1 (Reference Tόth and OdstrčilTόth and Odstrčil, 1996), so that

(44)

This leads to the scheme

(45)

which possesses first-order accuracy in time and up to second-order accuracy in space.

4. Discussion

In order to compute the steady-state solution of the age equation (15), we iterate from t = 0 until tf = 1000, which is, considering the maximum age of the analytical solution A = 20.76 at the bottom, easily sufficient to reach steady state. We test the solution schemes using a variety of domain grid densities, kmax =20,40,60,80,100, which correspond to grid spacings Δz =1/20,1/40,1/60,1/80 and 100respectively. the time-step is Δt =0.5Δz, Friedrichs–Lewy condition so that with the Courant– (Reference Morton and Mayers.Morton and Mayers, 1994) is fulfilled.

Table 1 shows the relative errors of the basal age, of the computed age profiles for the numerical schemes we wish to test. for all schemes and all grid spacings, stable integration is achieved, and all schemes show clear convergence toward the analytical solution with decreasing grid spacing (increasing number of gridpoints). Surprisingly, for the smallest number of gridpoints (kmax = 20), the UP1 scheme, which has the smallest spatial accuracy (first order) and the largest numerical diffusion, gives by far the best result with less than 4%, whereas for the other schemes is greater than 20%. Apparently, the discretization error and the error due to numerical diffusion cancel each other out to a large extent. However, this balance disappears for larger numbers of gridpoints, so that for kmax ≥ 60 the result of UP1 is worst, and even for kmax = 100 the error is larger than for kmax = 20.

Table 1.

Relative errors rb of computed basal ages (see Equation (23))

The spatially second-order schemes UP2 and MTVDLF/ M give the best results throughout (except at kmax = 20, as noted above), and further, the results of the two schemes are identical. Closer inspection shows that this surprising finding does not hold in general, but is due to our particular problem (see Appendix). the errors rb of both schemes are very small, ≤1% for kmax ≥ 60. Even for the third-order scheme QUICK, rb is distinctly greater than for UP2 and MTVDLF/M at all resolutions. As for the TVD schemes, application of the less diffusive limiters (MTVDLF/S, MTVDLF/W) reduces the accuracy, so that some numerical diffusion is evidently desirable.

The CTCS scheme is particularly bad for small numbers of gridpoints, kmax = 20, 40. This is because it tends to produce some oscillations over the whole spatial domain, which may even lead to unstable integrations. These oscillations become evident in Figure 3, which shows the relative errors of the computed age profiles as functions of z for kmax = 20, 100 (schemes MTVDLF/S,W not considered). An attempt may be made to minimize this problem by adjusting the frequency of FTCS steps (section 3.3), but nevertheless the CTCS scheme cannot be recommended due to this susceptibility to instabilities.

Fig. 3

Relative errors r of computed ages (see Equation (24)) for (a)> kmax = 20 and (b) kmax = 100.

A comparison of the relative errors of UP1, UP2, MTVDLF/M and QUICK over the whole column is instructive (see Fig. 3). for UP1 the error is smallest close to the surface and increases essentially monotonically toward the base. By contrast, for UP2, MTVDLF/M and QUICK the errors are somewhat larger close to the surface, reach a minimum at approximately half the depth and increase again further down. This is probably because UP2, MTVDLF/M and QUICK use central interpolation like CTCS for the age see Equations (32)2 and (34)2), so that the near-surface error resembles that of CTCS. the errors of UP1 and QUICK are smaller than those of UP2 and MTVDLF/M for (small age gradients), whereas for z < z* (large age gradients) UP2 and MTVDLF/M perform better. Nevertheless, except for the very basal gridpoint at z = 0 (see above discussion of the basal error) the simplest UP1 scheme yields results comparable to those of UP2 and MTVDLF/M for high and low resolutions.

At the transition point where the vertical velocity is discontinuous in the second and the age in the third z derivative, the errors of UP1, UP2, MTVDLF/M and QUICK show some jumps, naturally more pronounced for high resolution. Further, some oscillations appear for high resolution close to the base. Generally, among those schemes the most artificial structure is introduced to the age profile by QUICK due to its high spatial order.

It is worth noting that in other similar models, different schemes are favoured. Reference Wang, Straughan, Greve, Ehrentraut and WangWang and Hutter (2001) investigate the one-dimensional problem of a discontinuity (Heaviside step) of some field quantity which travels with constant velocity c in direction x in an infinitely extended system governed by the advective wave equation Apart from the missing source term, the constant velocity and the missing boundaries, this model equation is equivalent to our age equation (15). However, the moving discontinuity poses a challenge to numerical schemes that is not present for our continuous age profile. Reference Wang, Straughan, Greve, Ehrentraut and WangWang and Hutter (2001) show that numerical diffusion in the UP1 scheme largely smears out the discontinuity, and CTCS, UP2 and QUICK produce some numerical oscillations in the vicinity of the discontinuity. Best reproduction is achieved by the MTVDLF schemes; however, in contrast to our problem the less diffusive limiters (MTVDLF/S, MTVDLF/W) are favourable because they best preserve the shape of the Heaviside step.

To conclude, our tests show clearly that second-order upstreaming (UP2), modified TVD Lax–Friedrichs (MTVDLF/M) solution techniques and, with limitations of the accuracy directly at the base for higher resolutions, even first-order upstreaming (UP1) yield the best results for typical age profiles in ice sheets (here approximated by a one-dimensional model problem). UP1 and UP2 are less expensive computationally, but spatial discontinuities of the age field (which are not included in our model problem) will either be smeared out (UP1) or will be accompanied by artificial oscillations (UP2). MTVDLF/M has the potential to reproduce such jumps more accurately at the cost of more computing time.

Generalization of the methods presented in section 3 to three-dimensional time-dependent situations described by Equation (3) is straightforward and will be carried out in the near future within the ice-sheet model SICOPOLIS. the applicability to such problems was successfully demonstrated by Reference Wang, Straughan, Greve, Ehrentraut and WangWang (2001) in the context of wind-induced circulations in lakes, where three-dimensional advection occurs in the evolution equations for the flow velocity and the temperature field.

Acknowledgements

The very thorough and constructive reviews of C. L. Hulbe and an anonymous referee are gratefully acknowledged. We also thank K. Hutter for reading and correcting the original manuscript. the scientific editor was T.H. Jacka.

Appendix

Equivalence Between MTVDLF/M AND UP2

The equivalence between the two schemes MTVDLF/M and UP2 for our particular problem can be proved in the following manner. for the whole spatial domain 0 ≤ z ≤ 1, the analytical steady-state solution of the age equation (15) has the properties

(A1)

(see Fig. 1). Provided the numerical solution at time tn is sufficiently close to the analytical steady state so that these relations also hold for the numerical solution, then, due to (A1)3,

(A2)

and with (A1)2,

(A3)

so that the slope limiter (Equations (36) and (37)) computed with Minmod is

(A4)

Thus, from Equation (42),

(A5)

When this result and condition (A1)1 in the form |vz| = –vz are inserted into the MTVDLF/M scheme (Equation (45)), it reduces to the UP2 scheme (Equation (32)) for all grid-points k = 0 . . . kmax.

The equivalence of MTVDLF/Mand UP2 can be similarly shown for the cases

(A6)

(A7)

(A8)

For case (A6) the proof is analogous to (A1), whereas for (A7) and (A8),

(A9)

so that the slope limiter is

(A10)

which is (except for k = 0) the upstream direction for positive velocities vz.

References

Calov, R., Savvin, A., Greve, R., Hansen, I. and Hutter, K.. 1998. Simulation of the Antarctic ice sheet with a three-dimensional polythermal ice-sheet model, in support of the EPICA project. Ann. Glaciol., 27, 201–206.Google Scholar
Dahl-Jensen, D., Johnsen, S.J., Hammer, C.U., Clausen, H. B. and Jouzel, J.. 1993. Past accumulation rates derived from observed annual layers in the GRIP ice core from Summit, central Greenland. In Peltier, W.R., ed. Ice in the climate system. Berlin, etc., Springer-Verlag, 517–532. (NATO ASI Series I: Global Environmental Change 12.)Google Scholar
Dansgaard, W. and Johnsen., S.J. 1969. A flow model and a time scale for the ice core from Camp Century, Greenland. J. Glaciol., 8(53), 215–223.Google Scholar
Dansgaard, W. and 10 others. 1993. Evidence for general instability of past climate from a 250-kyr ice-core record. Nature, 364(6434), 218–220.CrossRefGoogle Scholar
Greve, R. 1997. Large-scale ice-sheet modelling as a means of dating deep ice cores in Greenland. J. Glaciol., 43(144), 307–310. (Erratum: 43(145), p. 597–600.)CrossRefGoogle Scholar
Greve, R.,Weis, M. and Hutter, K.. 1998. Palaeoclimatic evolution and present conditions of the Greenland ice sheet in the vicinity of Summit: an approach by large-scale modelling. Palaeoclimates, 2(2–3),133–161.Google Scholar
Greve, R., Mügge, B., Baral, D., Albrecht, O. and Savvin, A.. 1999. Nested high-resolution modelling of the Greenland Summit region. In Hutter, K.,Wang, Y. and Beer, H.. eds. Advances in cold-region thermal engineering and sciences: technological, environmental, and climatological impact. Berlin, etc., Springer-Verlag, 285–306. (Lecture Notes in Physics 533.)Google Scholar
Harten, A. 1983. High resolution schemes for hyperbolic conservation laws. J. Comput. Phys., 49(3), 357–393.Google Scholar
Jiang, G.-S., Levy, D., Lin, C.-T., Osher, S. and Tadmor, E..1998. High-resolution non-oscillatory central schemes with non-staggered grids for hyperbolic conservation laws. SIAMJ. Numer. Anal., 35(6), 2147–2168.Google Scholar
Johnsen, S. J. and 8 others. 2001. Oxygen isotope and palaeo-temperature records from six Greenland ice-core stations: Camp Century, Dye-3, GRIP, GISP2, Renland and NorthGRIP. J. Quat. Sci., 16(4), 299–307.Google Scholar
Leonard, B. P. 1979. A stable and accurate convective modelling procedure based on quadratic upstream interpolation. Comput. Methods Appl. Mech. Eng., 19(1), 59–98.Google Scholar
Morton, K.W. and Mayers., D.F. 1994. Numerical solution of partial differential equations. Cambridge, Cambridge University Press.Google Scholar
Petit, J.-R. and 18 others. 1999. Climate and atmospheric history of the past 420,000 years from the Vostok ice core, Antarctica. Nature, 399(6735), 429–436.Google Scholar
Savvin, A., Greve, R., Calov, R., Mügge, B. and Hutter, K.. 2000. Simulation of the Antarctic ice sheet with a three-dimensional polythermal ice-sheet model, in support of the EPICA project. II. Nested high-resolution treatment of Dronning Maud Land, Antarctica. Ann. Glaciol., 30,69–75.Google Scholar
Sweby, P.K. 1984. High resolution schemes using flux limiters for hyperbolic conservation laws. SIAMJ. Numer. Anal., 21(5), 995–1101.Google Scholar
Tόth, G. and Odstrčil, D.. 1996. Comparison of some flux corrected transport and total variation diminishing numerical schemes for hydrodynamic and magnetohydrodynamic problems. J. Comput. Phys., 128(1), 82–100.Google Scholar
Wang, Y. 2001. Comparing different numerical treatments of advection terms for wind-induced circulations in Lake Constance. In Straughan, B., Greve, R., Ehrentraut, H. and Wang, Y.. eds. Continuum mechanics and applications in geophysics and the environment. Berlin, etc., Springer-Verlag, 368–393.Google Scholar
Wang, Y. and Hutter, K.. 2001. Comparisons of several numerical methods with respect to convectively-dominated problems. Int. J. Num. Methods Fluids, 37(6),721–745.Google Scholar
Yee, H. C. 1989. A class of high-resolution explicit and implicit shock-capturing methods. U.S. Nat. Aeron. SpaceAdmin. Tech.Mem. 101088.Google Scholar
Figure 0

Fig. 1Analytical profiles of the vertical velocity vz (Equation (9)) and the age A (Equation (12)).

Figure 1

Fig. 2Staggered numerical grid for the one-dimensional age computation.

Figure 2

Table 1.Relative errors rb of computed basal ages (see Equation (23))

Figure 3

Fig. 3Relative errors r of computed ages (see Equation (24)) for (a)> kmax = 20 and (b) kmax = 100.