**Development and technical paper**
28 Nov 2020

**Development and technical paper** | 28 Nov 2020

# Implementation of the RCIP scheme and its performance for 1-D age computations in ice-sheet models

Fuyuki Saito Takashi Obase and Ayako Abe-Ouchi

^{1},

^{2},

^{2,1}

**Fuyuki Saito et al.**Fuyuki Saito Takashi Obase and Ayako Abe-Ouchi

^{1},

^{2},

^{2,1}

^{1}Japan Agency for Marine-Earth Science and Technology (JAMSTEC), Yokohama, Japan^{2}Atmosphere Ocean Research Institute, University of Tokyo, Kashiwa, Japan

^{1}Japan Agency for Marine-Earth Science and Technology (JAMSTEC), Yokohama, Japan^{2}Atmosphere Ocean Research Institute, University of Tokyo, Kashiwa, Japan

**Correspondence**: Fuyuki Saito (saitofuyuki@jamstec.go.jp)

**Correspondence**: Fuyuki Saito (saitofuyuki@jamstec.go.jp)

Received: 20 Feb 2020 – Discussion started: 04 Mar 2020 – Revised: 30 Sep 2020 – Accepted: 18 Oct 2020 – Published: 28 Nov 2020

Ice-sheet age computations are formulated using an Eulerian advection equation, and there are many schemes that can be used to solve them numerically. Typically, these differ in numerical characteristics such as stability, accuracy, and diffusivity. Furthermore, although various methods have been presented for ice-sheet age computations, the constrained interpolation profile method and its variants have not been examined in this context. The present study introduces one of its variants, a rational function-based constrained interpolation profile (RCIP) scheme, to one-dimensional ice age computation and demonstrates its performance levels via comparisons with those obtained from first- and second-order upwind schemes. Our results show that the RCIP scheme preserves the pattern of input surface mass balance histories in terms of the vertical profile of internal annual layer thickness better than the other schemes.

- Article
(5212 KB) -
Supplement
(15974 KB) - BibTeX
- EndNote

Core samples extracted from ice sheets can provide an archive of past climate history data, and a major issue for researchers attempting to utilize ice-core properties is defining the age of ice along the depth of the ice sheet. This process is often called dating. Dating with numerical ice-flow models is an important approach because it allows researchers to estimate age profiles before actual drilling of ice cores. For example, in Fischer et al. (2013), the authors present an application of ice-flow models to evaluate potential “oldest-ice” study areas.

Various methods for use in ice-sheet model dating have been adopted and compared. Mügge et al. (1999) compared particle tracking (Lagrangian) and Eulerian schemes under simulated steady-state three-dimensional (3-D) velocity fields of the Antarctic ice sheet. That study concluded that the Eulerian scheme works well, except for the bottom part, which encounters problems due to numerical diffusion. In Rybak and Huybrechts (2003), the authors also compared the Lagrangian and Eulerian schemes for simulated Antarctic ice sheets under various schematic steady-state conditions and analytical solutions, as well as under different 3-D velocity fields. Similarly, they concluded that the Lagrangian method produced less error than an Euler approach, although the difference was small over a large part of the domain. Greve et al. (2002) compared several Eulerian schemes such as central difference schemes, first- and second-order upwind schemes, Quadratic Upstream Interpolation for Convective Kinematics (QUICK), and total variation diminished (TVD) Lax–Wendroff (LW) schemes. From comparisons of the one-dimensional (1-D) steady-state age profiles produced by these schemes, they concluded that the second-order upwind and TVD-LW schemes performed well for typical ice-sheet age profiles. Comparisons among semi-Lagrangian schemes have also been performed. Introduction of a semi-Lagrangian trace scheme to ice-sheet modeling was initiated by Clarke and Marshall (2002). They simulate the temporal and spatial variations of water isotopes in the Greenland ice sheet over the past 30 000 years. In Tarasov and Peltier (2003), the authors compared various interpolation schemes in order to compute upwind departure points in a semi-Lagrangian tracer model in terms of preservation of input signal phases and amplitudes, while Lhomme et al. (2005); Clarke et al. (2005) developed a new interpolation method that can be used in a semi-Lagrangian scheme and discussed computed ice-core age–depth relationships for the Greenland Ice Core Project (GRIP) ice core.

To date, various methods have been presented and demonstrated for use in ice-sheet age computations. However, there is still a variety of numerical schemes that have not been examined within this context. These include the constrained interpolation profile (CIP) method (e.g., Yabe et al., 2001) and its variants. Accordingly, the present study introduces a CIP method variant named the rational function-based constrained interpolation profile (RCIP) method (Xiao et al., 1996) for use in 1-D ice age computations and demonstrates the performance of the scheme.

## 1.1 Brief introduction of the RCIP scheme

This section describes a standard algorithm of the CIP scheme family that is used to solve a 1-D advection equation with a non-advection term as follows:

where $f=f(x,t)$ is a free variable to solve,
$u=u(x,t)$ is a velocity field,
$h=h(x,t)$ is an arbitrarily (non-advection) field,
and *t* and *x* are temporal and spatial coordinates, respectively.

As introduced in the previous section, there are three major approaches to solving an advection equation: Eulerian, Lagrangian, and semi-Lagrangian. The CIP scheme family corresponds to a semi-Lagrangian method variation. The basics of the semi-Lagrangian approach, within the context of its comparison with the Lagrangian and Eulerian approaches, have already been presented in a number of past studies. For example, Staniforth and Côté (1990) presented a review of these methods and described the implementation and application of a semi-Lagrangian method in detail. Although a full description of the semi-Lagrangian is not repeated in this paper, its basic principles will be described later in this section.

In CIP schemes, Eq. (1) is solved by performing a time-splitting algorithm (e.g., Yabe and Takei, 1988) in two phases as follows:

Appendix A presents a note on the time-splitting technique.

The primary characteristic of this CIP scheme is the introduction of
an additional equation to solve the spatial derivatives of *f* at the same time.
Differentiation of Eq. (1) provides the equation for
$g(x,t)=\frac{\partial f}{\partial x}$:

Equation (4) is an advection formula that is similar to Eq. (1) with the non-advection function $\widehat{h}(x,t)$ on the right-hand side, which is solved using a time-splitting procedure similar to those used in Eqs. (2) and (3):

The algorithm used to solve the advection phases (Eqs. 2 and 5), which is a core characteristic of the CIP scheme family, is described here, after which the algorithm used to solve the non-advection phases (Eqs. 3 and 6) will be discussed.

In semi-Lagrangian approaches, a particle at $(x,t+\mathrm{\Delta}t)$
originates from the position of the upstream departure point *x*_{dep} such that

where

Figure 1 shows a schematic illustration of
semi-Lagrangian scheme.
The particle at *x*_{j} at time *t*_{m}+Δ*t* originates from a
particle at *x*_{dep}, which is not necessarily on a discretized grid
point *x*_{j}.
Therefore, the free variable *f*(*x*) must be interpolated using the
value on the grid points (represented by color shading in the
figure).

The CIP method constructs an interpolation function *F*_{j}(*x*) for the *f*(*x*)
between two adjacent grid points *x*_{j} and
its upwind point *x*_{j+1} when
*u*_{j}<0 in order to assess the value at the departure point.
Introducing $\langle \mathit{\xi}\rangle ={x}_{\mathrm{dep}}-{x}_{j}$ as the distance to the original
point allows the time evolution of *f*(*x*_{j}) (which is the
original free variable to solve)
and *g*(*x*_{j}) (which is the spatial derivative of *f*
at the grid points *x*_{j})
to be solved
as simple advection equations:

where ${G}_{j}\left(x\right)=\frac{\partial {F}_{j}}{\partial x}$. Note that
computation of distance to the
departure point will be described in
Sect. 1.2. The piecewise interpolation
function *F*_{j}(*x*) for ${x}_{j}\le x\le {x}_{j+\mathrm{1}}$ is defined to be
constrained by the continuity condition at *x*_{j} and *x*_{j+1}
as

A cubic polynomial is chosen in the original CIP scheme, as

where $X=x-{x}_{j}$. The four coefficients (*C*_{0}, *C*_{1}, *C*_{2},
and *C*_{3}) in Eq. (11) are determined to satisfy the
constraints (Eq. 10).
The RCIP scheme framework is occasionally extended to
introduce a rational function (Xiao et al., 1996) such as

The interpolation function is switched from the cubic (Eq. 11) to the rational (Eq. 12) if ${g}_{j}\le {S}_{j}\le {g}_{j+\mathrm{1}}$ or if ${g}_{j}\ge {S}_{j}\ge {g}_{j+\mathrm{1}}$, where

Additionally, the four coefficients (*C*_{0}, *C*_{1}, *C*_{2}, and *D*_{1}) are determined
in order to satisfy the same constraints.
The two interpolation functions (Eqs. 11 and 12) are integrated by introducing
a switching parameter *α*:

The five coefficients used to satisfy the constraints are computed as

The switching parameter $\mathit{\alpha}\in [\mathrm{0},\mathrm{1}]$ is chosen as 1 when it
is necessary to use
rational interpolation. In other cases, 0 is selected.
If *{**f*_{j}*}* and *{**g*_{j}*}* at time *t* are known, the new states
$\mathit{\left\{}{f}_{j}^{*}\mathit{\right\}}$ and $\mathit{\left\{}{g}_{j}^{*}\mathit{\right\}}$ are predicted by shifting by distance along the
characteristics (Eq. 9)
to the departure point 〈*ξ*〉,
as follows:

The solutions above are those of the advection phases
(Eqs. 2 and 5).
The time evolutions of *f* and *g* in the non-advection phases are again
calculated according to Eqs. (3) and (6),
typically by using a forwarding scheme, starting from the solution of the
advection phase $\mathit{\left\{}{f}_{j}^{*}\mathit{\right\}}$ and $\mathit{\left\{}{g}_{j}^{*}\mathit{\right\}}$, as an intermediate solution:

As discussed in Xiao et al. (1996), the formulation of the RCIP scheme possesses attractive properties, such as convexity and monotone preservation, as well as phase speed.

## 1.2 Upstream departure point

The interpolation method used for the field variables, which characterize each scheme, is one of the most important topics in semi-Lagrangian schemes. Another major topic common to the semi-Lagrangian schemes is the method used to compute the departure point.

Equation (8) gives the distance to the departure point:

A simple and primitive way to integrate Eq. (22) is to use the local velocity even if the velocity is a function of time and space (e.g., Toda et al., 2009), such that

Figure 1 shows a trajectory
and departure point
under a constant velocity (dotted line and circle) assumption.
As can be seen, the computed departure point can be
different from a general non-uniform velocity situation.
Another way is to apply the “midpoint rule”, where
both spatial and temporal mean velocities between the target and departure
points replace *u*_{j} in Eq. (23),
which is generally computed in an iterative fashion (Tarasov and Peltier, 2003).
In the present paper, a third approach is adopted.
First, a steady and linear velocity field
between the target and the upstream adjacent points,
*x*_{j} and *x*_{j+1},
is assumed such that

where *u*^{′} is a constant spatial gradient of the velocity.
In order to solve the time evolution of the velocity of a particle at
(*t*_{m},*x*_{j}), Eq. (24) is differentiated by time
*t*:

which is solved as

Introducing Eq. (26), Eq. (22) is integrated as

Based on the above, it can be interpreted that the distance to the departure point is
that of constant velocity case
(Eq. 23 or 28),
multiplied by the bracket term in Eq. (27) as a
correction factor. Here, it should be noted that the correction factor reaches 1
toward the limit of ${u}^{\prime}\to \mathrm{0}$, which definitely corresponds
to the constant velocity case.
The velocity gradient *u*^{′} already appears in
the advection equation of the *g* term (Eq. 4), which is
reused in the departure point computation.

## 2.1 Governing equation

The computation used to determine the age of the ice, i.e., the elapsed time since the ice
deposit, is performed with the pure advection equation^{1}:

where 𝒜 is the age and *t* is time, which is the
Lagrangian approach.
Eq. (29) is then reformulated into the Eulerian equation
for a 1-D problem,

where $\mathcal{A}=\mathcal{A}(z,t)$ and $w=w(z,t)$ are the age and vertical velocity
fields, respectively, and *z* is the vertical coordinate.
Some models introduce
an artificial diffusion term in order to achieve stable
integration (e.g., Mügge et al., 1999). However, the pure
advection form is kept throughout the present paper.
Following most large-scale numerical ice-sheet
models (Greve and Hutter, 1995), the vertical coordinate *z* is scaled with
the local thickness.
Introducing the scaled coordinate *ζ* as

Eq. (30) is reformulated as follows:

where *τ*≡*t* is the corresponding time coordinate in this
system, *b*=*b*(*t*) is the bedrock topography,
and *H*=*H*(*t*) is the ice thickness.
The new velocity term $\mathit{\omega}=\mathit{\omega}(\mathit{\zeta},\mathit{\tau})$ in
*τ*,*ζ* system is computed as

where derivatives of *ζ* are computed as

Since the ice thickness *H*,
which actually reflects the changes in the boundary conditions,
may not be constant throughout the time period,
*H*=*H*(*t*) is prescribed independently of the boundary
conditions in this paper.
The surface mass balance term *M*_{s}
(mass input into the domain), surface evolution, and the vertical
velocity at the surface *z*=*h*(*t*) are related as

which is derived from the kinematic boundary conditions based on the assumption of a flat surface. The spatial derivative of 𝒜 used in the RCIP scheme is derived as follows:

where ${\mathcal{A}}^{\prime}=\frac{\partial \mathcal{A}}{\partial \mathit{\zeta}}$.

In order to solve the time evolution of age and its gradient
(Eqs. 32 and 37), the initial and boundary
conditions are required.
At the free surface *z*=*h*(*t*) (or *ζ*=1), a Dirichlet-type boundary
condition,

holds when the surface mass balance is positive (i.e., *M*_{s}>0).
In contrast, when the surface balance is negative,
the boundary condition is not necessary, because the departure point of
the free surface is inside the ice.
A special treatment is required for the zero mass balance at the
surface, *M*_{s}=0. In this case, the velocity term in *τ*,*ζ* system, *ω* becomes
0, so Eq. (32) is simplified as

which, again, requires no boundary condition for age.
The boundary conditions at the bottom *ζ*=0 simply mirror those at
the surface.

The age derivative, 𝒜^{′}, also satisfies the boundary condition at the free
surface as

when *M*_{s}>0. Conditions similar to age hold for the age
derivative when *M*_{s}<0 and *M*_{s}=0.

In the present study, equivalent but different coefficient representations (Eqs. 15–19) are adopted for the RCIP method implementation, which is described in Appendix B.

## 2.2 Discretization

The spatial discretization of Eqs. (32) and (37) can be either uniform or non-uniform. In the present paper, both types of discretization are examined. Since uniform discretization is a special case of non-uniform discretization, the latter can be described effectively without a loss of generality.

One way to introduce a non-uniform discretization is to apply a non-smooth grid (Shashkov, 1995), which prescribes irregular discretization of the coordinates:

and

Another way to introduce a non-uniform discretization is to apply a smooth grid (Shashkov, 1995), which uses a smooth function to transform the coordinate system. One more coordinate transformation is then performed for a non-uniform smooth-grid system as follows:

where *T* and *Z* are the time and vertical coordinates in the new system. A smooth transformation of *Z*=*Z*(*ζ*) or its inverse *ζ*=*ζ*(*Z*) is prescribed where necessary.
Similarly, a new velocity term $W=W(Z,T)$ in *T*,*Z* system is computed as

Equations (43) and (44), which are the target
equations to solve,
are simply replacements for Eqs. (32) and (37), respectively.
The velocity term $W=W(T,Z)$ is prescribed (as will be explained later). The terms 𝒜, *W*, and 1 on the right-hand side correspond to *f*, *u*, and *h*, respectively, in the RCIP scheme framework (Eq. 1).
Although it is possible to introduce further non-uniform discretization on
the *Z* coordinate, in the present paper, only a uniform
discretization is examined on the smooth-grid discretization:

Actually, the discretization of the *ζ* coordinate corresponds
to the special case of non-uniform smooth discretization with
*Z*≡*ζ*. Therefore, for both uniform and non-uniform discretization, the scheme will be described hereafter using
the *Z* coordinate instead of the *ζ* coordinate.

## 2.3 Comparing other schemes with RCIP schemes

In the present paper, two numerical schemes, the first- and second-order upwind schemes, are examined in comparison with the RCIP schemes. While there are other numerical schemes suitable for such comparisons, including Lagrangian, other semi-Lagrangian, or even higher-order upwind schemes, these have already been reported in past studies (Mügge et al., 1999; Greve et al., 2002; Rybak and Huybrechts, 2003; Clarke and Marshall, 2002; Tarasov and Peltier, 2003; Clarke et al., 2005). Furthermore, since our study focuses on a demonstration of RCIP schemes in relation to the topic of ice dating, a wide range of comparisons is beyond the scope of this paper.

The “first-order” upwind scheme in the present paper evaluates the advection term using the velocity at staggered grid points as follows:

when ${W}_{k+\mathrm{1}/\mathrm{2}}<\mathrm{0}$ and ${W}_{k-\mathrm{1}/\mathrm{2}}<\mathrm{0}$. The velocity at staggered grid points is computed by linear interpolation of the two adjacent velocities at normal grid points. Equation (47) corresponds to numerical integration with the midpoint rule if a Dirichlet-type boundary condition is applied on the upper surface (Eq. 38) and the velocity is kept negative throughout. It is especially notable that, for ice dating at summits, positive (upward) vertical velocity is rarely considered. Therefore, the midpoint rule formulation mentioned above is sufficient for application. On the other hand, a different approach is generally required for a grid point where two velocities at staggered adjacent grid points have opposite signs. In this paper, the velocity term is simply replaced by that at the normal grid point:

where *W*_{k}<0, and ${W}_{k+\mathrm{1}/\mathrm{2}}$ and ${W}_{k-\mathrm{1}/\mathrm{2}}$ have opposite signs.

For the “second-order” upwind scheme, the derivative of the age term is replaced by the second-order upwind difference formulation as

where

for the *W*_{k}<0 case. The age derivative at the surface is required (${{\mathcal{A}}^{\prime}}_{k+\mathrm{1}}$ in Eq. 51),
which is provided as a boundary condition (Eq. 40).
For higher-order numerical schemes, the introduction of a slope
limiter is a standard method for suppressing the development of oscillations near a discontinuity and/or steep gradients (details are described in
Greve et al., 2002). Although it is possible to apply such slope limiters in irregular grids (Murman et al., 2005), an easier approach was adopted instead.
Specifically,
the formulation is switched back to the first-order scheme when
${{A}^{\prime}}_{\mathrm{I}}>\mathrm{0}>{{A}^{\prime}}_{\mathrm{II}}$ or ${{A}^{\prime}}_{\mathrm{I}}<\mathrm{0}<{{A}^{\prime}}_{\mathrm{II}}$.
Although this method may be insufficient to stabilize the solution
near a strong discontinuity, the implementation of more sophisticated slope limiters is beyond the scope of the present paper.

## 3.1 Experimental design

Following some modeling studies on the dating of deep drilling sites that used simplified 1-D vertical ice flow models (e.g., Parrenin et al., 2007), the present study adopts an analytical vertical velocity profile under the assumption that there are no horizontal variations in the bedrock elevation, surface, and basal mass balances:

where *M*_{s} and *M*_{b} are the surface and basal mass balance (positive is input), respectively, *H* is the ice thickness, and $\stackrel{\mathrm{\u0303}}{w}\left(\mathit{\zeta}\right)$ is the normalized
velocity profile. Assuming no basal sliding, $\stackrel{\mathrm{\u0303}}{w}\left(\mathit{\zeta}\right)$ can be approximated by

where *p* is a parameter for the profile (Parrenin et al., 2007).
Under Glen's flow law with a steady-state isotropic ice condition,
*p* is equal to the flow law exponent *n* (typically *n*=3).
In addition, the RCIP scheme requires the derivative of the velocity,
which is computed using the derivative of $\stackrel{\mathrm{\u0303}}{w}$, as

In addition to the vertical velocity, the time evolutions of the surface and basal mass balances and the ice thickness are required for the age computations. These will be presented in each of the following sections.

The initial conditions for the 𝒜 and 𝒜^{′} fields are set to 0 for all our experiments. In these cases, the age derivative 𝒜^{′} is kept 0 under the level at which the age reaches the integration time. Starting from the 0 field, time integration is computed for 2000 kyr for most of our experiments.

It is worth mentioning that formulations like
Eq. (52), which is a function of normalized depth,
make it possible to interpret example results
for different configurations when appropriate spatial/temporal
dimension scaling is used. In this case, the spatial and temporal characteristic scales can be defined, for example, by the ice thickness and the surface mass balance. This means that the age solution under the configuration
of 3 and 0 cm yr^{−1} for the surface and basal
mass balance, respectively, has the same normalized shape as that
under 30 and 0 cm yr^{−1},
by scaling all the time-related terms as 1∕10.

All the computations in our present study were performed on a personal computer (PC) equipped with an Intel Xeon E5-2609 central processing unit (CPU) and compiled with GNU Fortran. Each surface/basal mass balance, ice thickness, and vertical resolution configuration is repeated using four numerical schemes: the RCIP with departure correction (RCIP+corr), the RCIP without correction (RCIP), the second-order upwind scheme (UP-2), and the first-order upwind scheme with midpoint rule (UP-1). Additionally, a first-order scheme without a midpoint rule (UP-1n) is sometimes used. Multiple 1-D-column experiments with different boundary configurations using one numerical scheme are examined simultaneously in one run. For example, the mean computational costs for one run (with 28 different configurations) in the case of 129 levels over 200 kyr are 30, 28, 32, and 34 s, using UP-1, UP-2, RCIP, and RCIP+corr, respectively. Those in the case of 513 levels are 338, 296, 364, and 392 s, respectively. Details differ among the configurations, and it takes 30 % to 40 % more time to perform a RCIP+corr run than to perform a UP-2 or UP-1 run.

## 3.2 A verification experiment using uniform velocity

Before performing an experiment under a typical ice-sheet configuration,
verification of the numerical model used in the present study is
presented under further simplified conditions, namely the constant velocity
case.
This is easily performed using Eq. (52), in which the parenthesis
term equals 0, in other words, by keeping *H* constant and
setting ${M}_{\mathrm{s}}\equiv -{M}_{\mathrm{b}}$ for
arbitrarily *p*.

Figure 2 shows the computed age profile under
the uniform
velocity of −15 cm yr^{−1} and *H*=3000 m.
Uniform grid spacing of 129 levels is adopted,
which corresponds to *Z*≡*ζ* and $\mathrm{\Delta}Z=\mathrm{\Delta}\mathit{\zeta}=\mathrm{1}/{\mathrm{2}}^{\mathrm{7}}$
(i.e., Δ*z*=23.4375 m) using the smooth grid.
The time step is set as 100 years,
which corresponds to the Courant–Friedrichs–Lewy (CFL) condition ∼0.64. The vertical age profile is formulated as

thus, the exact solution for an uniform velocity is

where ${w}_{c}=-\mathrm{15}$ cm yr^{−1}.
For completion purposes, the results of the RCIP scheme are plotted in the figure, which is (by definition) identical to those of RCIP+corr scheme. For the steady state, a linear age profile from 0 years at the surface
and 20 kyr at the bottom is expected (corresponding to the
thick gray line in Fig. 2),
which is obtained by all the methods after integration of
around 27 kyr (not shown).
In contrast, the transient states are different among the results of
the four schemes examined. Figure 3 shows the computed age profile relative to the exact solution, with three different time steps (100, 50, and 25 years) for each scheme.
The results of RCIP+corr (and thus RCIP) are shown to be less sensitive to the time step than the upwind schemes, which reflects the fact that both the interpolation and the departure point calculation are successful. At 20 kyr, a linear age profile should be obtained, but
all four results show ages that are younger than the exact solution,
due to numerical diffusion. Additionally, while all of the schemes show relatively good performance for the upper part, the result obtained by the UP-1 scheme deviates the most from the solution.
Specifically, it deviates
1 year from around two-thirds of the total depth and reaches almost
1 kyr at the bottom, which is already visible in
Fig. 2. In contrast, the other results
deviate from the solution only near the bottom $\sim \mathrm{9}/\mathrm{10}$ of the
total depth and reach ∼100 years, or even less, at the bottom.
The error at the bottom of UP-1 is 759 to 902 years (3.8 % to
4.5 %),
while that of RCIP is 76 to 98 years (0.38 % to 0.49 %),
and the best of UP-2 is even better at 7.5 to 154 years (0.04 % to
0.77 %).

## 3.3 Experiment with steady non-uniform vertical velocity

Hereafter, non-uniform velocity experiments are performed using $p=n=\mathrm{3}$ in Eq. (53). First, simple cases with constant surface/basal mass balances, as well as thicknesses that correspond to steady vertical velocity profiles, are shown. Since Eq. (55) cannot be solved using Eq. (52), profiles created by numerical integration (the Runge–Kutta scheme) are used as “benchmark” solutions in this and the following section.

The ice thickness and the accumulation rate chosen in this and the following sections are ∼ 3000 m and ∼ 3 cm yr^{−1}, respectively, which correspond to typical quantities for the East Antarctic Plateau during the glacial period (e.g., Parrenin et al., 2007; Fischer et al., 2013). On the other hand,
there are other cases with similar thicknesses and 10 times higher accumulation (∼ 30 cm yr^{−1}), typically in
Greenland or West Antarctica (e.g., Clarke and Marshall, 2002).
As mentioned in Sect. 3.1,
with proper scaling,
the results in these sections can be interpreted as results
with much higher accumulation rates.
This will be examined in the discussion (Sect. 4).

Two sets of basal melting are presented: no basal melting
and 3 mm yr^{−1}.
The other two parameters are fixed. Surface mass balance is set as
3 cm yr^{−1} and thickness is set as *H*=3000 m.
We use a uniform grid spacing of 129 levels (Δ*z*=23.4375 m), and the time step is set as 100 years,
which is the same configuration used in the previous section.

Figures 4a and 5a
show computed age profiles at
*t*=2000 kyr for all the schemes along with the benchmark age
profiles. Very few differences can be seen among the profiles over most
parts of the figures under this scale.
Deviations from the benchmark are shown in
Figs. 4b and 5b.
The results of each scheme show larger errors
near the bottom than near the upper part.
Some results show sudden increases in the error at certain depths,
which correspond to the depths around where the age should
reach the time of integration.

The RCIP+corr scheme shows the best result for all depths. The UP-1 scheme shows the second-best result, which is even better than the RCIP scheme around the depth of 2600 m. However, it also shows the largest errors among all the schemes examined at deeper depths. The good performance of UP-1, in spite of its smallest spatial accuracy, which is attributed to the cancellation of errors due to discretization and numerical diffusion, has already been presented in Greve et al. (2002). The midpoint rule formulation (Eq. 47) also plays a role in the increased accuracy. Due to simple situations, such as the one-direction advection and the constant upper boundary conditions, the age profile computation can be formulated as a vertical integration from top to bottom. This means that the midpoint rule integration actually has second-order accuracy. A “true” first-order upwind scheme can be applied by using Eq. (48) over the whole domain. In this case, vertical integration from top to bottom corresponds to an Euler integration, which has first-order accuracy. Figures 4b and 5b also contain results obtained using such a scheme, marked as UP-1n. However, as expected, when such a normal-grid velocity is introduced for the advection equation, the results have less accuracy than those of the second order upwind (UP-2). Furthermore, as shown in Greve et al. (2002), the improved performance of the UP-1 scheme is limited to the upper part, and the errors become larger as the depth increases. The results of the RCIP scheme show relatively larger errors than the other methods, except for the top and the bottom parts, which highlights the importance of accurate departure point calculations. The result of the UP-2 scheme shows intermediate errors between RCIP and UP-1 at the bottom.

## 3.4 Non-steady surface mass balance experiments

This section presents the results of experiments conducted with non-steady velocity profiles, which were performed with the prescribed surface mass balance time series. First, a very simple square-wave formulation is adopted for the time evolution of the surface mass balance.

where *a*_{H} and *a*_{L} are the prescribed high and low surface mass balance terms, *P*_{H} and *P*_{L} are the durations with high- and low-value phases, and *P*_{T} is the duration of one cycle. Figure 6 shows the
time evolution of a normalized surface mass balance with *P*_{T}=100 kyr cycles
and a phase pattern of ${P}_{\mathrm{H}},{P}_{\mathrm{L}}=\mathrm{1}:\mathrm{1}$ as an example.
Several experimental configuration combinations are examined, including
${P}_{\mathrm{H}},{P}_{\mathrm{L}}=\mathrm{1}:\mathrm{1}$, 7:1, or 1:7;
and *M*_{b}=0, 0.3, or 3 mm yr^{−1}.
The other patterns examined in this paper are provided as the
Supplement to this paper.

Figures 7 and 8
show computed age profiles at
*t*=1000 kyr under the square-wave surface mass balance,
where the lower surface mass balances are set as
*a*_{L}=1.5 cm yr^{−1} and 0.75 cm yr^{−1}.
The higher surface mass balances and the basal are set as
*a*_{H}=3 cm yr^{−1} and *M*_{b}=0, respectively.
For reference purposes, the benchmark solutions with constant surface mass balances of *a*_{H} and *a*_{L} are shown with gray lines.
The black line is the benchmark solution with the constant
surface mass balances of the mean, ${a}_{\mathrm{M}}=({a}_{\mathrm{H}}+{a}_{\mathrm{L}})/\mathrm{2}$.
As shown in these figures, the computed age profiles are close to the
benchmark solution with *a*_{M}, particularly at the bottom.
For the upper part, the computed age profiles
are along the benchmark solutions for the constant surface mass balance cases of *a*_{L}.

Since there are few visible differences among the computed ages, the computed age profiles relative to the one produced by the RCIP+corr scheme (Fig. 7b) are shown. The figures show comparable relative performance levels in spite of the different input surface mass balance histories. The age profiles produced by the RCIP scheme deviate systematically from RCIP+corr by less than 1 kyr throughout the depth range, which reflects the differences in computing the departure points. The other two schemes deviate around 10 kyr at most. The age difference oscillations seen in the UP-2 and UP-1 schemes are visible near the age corresponding to the time when switching was conducted between the high and low surface mass balances. (Δ𝒜 vs. 𝒜 plots are presented in the Supplement). These oscillations reflect the characteristics of the UP-2 and UP-1 schemes at the discontinuities.

Figures 7c and 8c
show the computed annual layer thickness, *λ*, against the depth.
In the present paper, the annual layer thickness is defined as the
inverse of the age gradient. For the RCIP+corr and RCIP methods, the computed field of the age derivative itself (𝒜^{′} in
Eq. 44) is used with the coordinate transformation.
On the other hand, for the UP-2 and UP-1 methods,
the diagnosed field is used (${{A}^{\prime}}_{\mathrm{I}}$ or ${{A}^{\prime}}_{\mathrm{II}}$
in Eqs. 47, 50, and 51,
respectively).
An infinite or a very large annual layer thickness may be present near the
bottom, due to the zero gradient of age as a consequence of the initial
experiment conditions.
In this study, the last 500 years are clipped from the figures.

The annual layer thickness has the following relationship in terms of the thinning rate:

while assuming that layers remain horizontal (Cuffey and Paterson, 2010). When the basal mass balance is 0 and the thickness is constant, the vertical velocity gradient can be formulated from Eq. (52) as

Finally, after some derivation, the vertical gradient of annual layer thickness can be formulated as

which is a function of *λ* and the normalized vertical velocity shape.
This experiment was conducted with zero basal mass balance, constant thickness,
and the same normalized velocity. Therefore, the vertical profile of
the annual layer thickness should go back and forth on the two lines
produced by those computed using the constant surface mass balances.

In terms of computed annual layer thickness profiles, the RCIP+corr and RCIP (which overlap with RCIP+corr in the figure)
methods show particularly good performance over the upper part,
as shown in Figs. 7c and 8c.
Dissipation at the discontinuity becomes larger towards the bottom,
but the solution of RCIP+corr (RCIP) is somewhat more stable on the two benchmark lines than the other schemes.
Overshooting at the discontinuity is shown for the solution by the UP-2 scheme, which becomes larger as the difference between the high and low surface mass balances increases.
In the present study, this is considered to be a consequence of
an inadequate slope filter. In addition, the annual layer thickness is diagnosed with Eq. (50) for the UP-2 scheme, which may exaggerate the oscillation of age gradients more than the simple first-order Taylor expansion. For the UP-1 scheme results,
the annual layer thickness diffuses with the depth level and approaches
the constant accumulation case of its mean. In deeper areas, the annual layer thickness is found in the vicinity (above or below) of the mean *a*_{M} benchmark profile in all of the numerical schemes.

The same exercises were performed using a different shape for the time evolution of the surface mass balance. Figure 9 shows the results for an experiment conducted using the cosine-wave formulation of the surface mass balance (Fig. 6), which is relatively more continuous than the square-wave version. Similar performance levels were obtained by the UP-2 and RCIP+corr (RCIP) methods for the small amplitude case (Fig. 9c). Instability also arises at low-to-high transitions when the ratio of high-to-low accumulation is larger (archived in the Supplement).

Figure 10 shows the results obtained by
square-wave forcing in terms of computed annual layer thickness,
*λ*, against the computed age for all the schemes, obtained by the
relative duration of the ${P}_{\mathrm{H}}:{P}_{\mathrm{L}}=\mathrm{1}:\mathrm{1}$ case (similar figures obtained by other experimental configurations are archived in the Supplement). Since the periodicity of the input cycle is 100 kyr in this experiment, the annual layer thickness profiles should show the same periodicity. The obtained results show relatively good performance for the RCIP+corr (RCIP) scheme in terms of
the phases when compared with the UP-2 and UP-1 schemes.
Dissipation at the discontinuity blur the square-wave shapes, particularly
at the deeper part, but the phases are still maintained better by the
RCIP+corr (RCIP) scheme than by the UP-2 scheme.

## 3.5 Non-steady thickness experiments

The time evolution of the surface mass balance often involves the evolution of ice thickness as a response. In this section, age computation performance levels under non-steady mass balance and ice thickness conditions are presented. In the present paper, the time evolution of thickness is computed as follows:

where *H*_{ref}(*M*_{s}) is the reference thickness as a function that
depends solely on the surface mass balance and *τ*_{H} is the response
thickness timescale.
Under ideal conditions, the steady-state ice thickness at the summit is
proportional to the $\mathrm{1}/(\mathrm{2}n+\mathrm{2})$ power of the surface mass balance,
where *n* is Glen's flow law exponent (Cuffey and Paterson, 2010).
Following this relationship, the reference thickness is formulated as

For cases where $H(t=\mathrm{0})=\mathrm{3000}$ m,
${a}_{\mathrm{H}},{a}_{\mathrm{L}}=\mathrm{3},\mathrm{1.5}$ cm yr^{−1} and *P*_{T}=100 kyr, the
evolution of *H* over the first two cycles can be computed as shown in
Fig. 11 using Eqs. (61) and (62). The lower thickness limit in this case is
$\mathrm{3000}\times (\mathrm{1.5}/\mathrm{3.0}{)}^{\mathrm{1}/\mathrm{8}}\sim \mathrm{2751.01}$ m.

Several experimental configuration combinations are examined.
These include square-wave or cosine-wave forcing;
${a}_{\mathrm{H}},{a}_{\mathrm{L}}=\mathrm{3},\mathrm{1.5}$ cm yr^{−1} or 3,0.75 cm yr^{−1};
*τ*_{H}=3 kyr or 10 kyr;
${P}_{\mathrm{H}},{P}_{\mathrm{L}}=\mathrm{1}:\mathrm{1}$, 7:1, or 1:7;
*M*_{b}=0, 0.3, or 3 mm yr^{−1}.
Figure 12 shows the result
of experiments conducted with response timescales of 10 kyr
for the 100 kyr cycle square-wave case provided as an example. The gray lines in the figures are the benchmark solution with constant
surface mass balances of *a*_{H} and *a*_{L} and their corresponding reference thickness of *H*_{ref}.
The black line is computed using the mean surface mass balance
${a}_{\mathrm{M}}=({a}_{\mathrm{H}}+{a}_{\mathrm{L}})/\mathrm{2}$, and the mean thickness over the last cycle of its evolution.

A comparison with the fixed thickness experiments (Fig. 7 vs. Fig. 12, or other combinations archived in the Supplement) shows no significant differences. The preservation of discontinuity in the annual layer thickness is similar to that seen in the non-steady thickness case. The differences in computed age, as well as the performance levels of the phases in the annual layer thickness, are qualitatively the same. In addition, all of the combinations examined in this paper show corresponding results that are qualitatively similar to those obtained in the fixed thickness experiments.

## 3.6 Occasional non-positive surface mass balance experiments

So far, the surface mass balance values adopted in our experiments
have been positive (corresponding to the accumulation zone).
This limitation is sufficient for the usual topics relating to deep
ice-core experiments, where the interpretation of ice-core data may become too complex. However, in order to provide a complete demonstration of the
performance levels of numerical age computations for more general cases, it is worthwhile to examine other cases. Although it may be considered pointless to examine steady negative mass balance cases because they simply mirror the steady
positive cases presented above,
the surface mass balance level adopted in this section is examined with
zero or negative *a*_{L} in Eq. (57) and
Fig. 6.
One encountered difficulty is computing vertical velocity when the surface mass
balance is negative.
Strictly speaking, it is possible to apply negative *M*_{s} to Eqs. (52) and (53), but the validity of such a formulation may be questionable because it is based on an idealized steady-state ice-sheet solution under positive surface mass balance conditions (e.g., Rybak and Huybrechts, 2003). However, for the sake of simplicity,
the vertical velocity profiles in the current study are
prescribed using the same set of equations for both positive and negative
surface mass balances. This is considered to be sufficient,
particularly for evaluations of the numerical performance levels of different schemes.

The results of the transient experiments that were conducted under a square-wave surface mass balance of *a*_{L}=0 cm yr^{−1}
are presented in Fig. 13,
while the other configuration is the same as in Figs. 7 and 8.
The results obtained under a configuration with ${a}_{\mathrm{L}}=-\mathrm{1.5}$ cm yr^{−1} are archived in the Supplement.
For both experiments, the thickness is fixed as 3000 m,
the mass balances are *a*_{H}=3.0 cm yr^{−1},
*M*_{b}=0, and the phases are
${P}_{\mathrm{H}},{P}_{\mathrm{L}}=\mathrm{50}$ and 50 kyr.

Several experimental configuration combinations are examined.
In a comparison involving the experimental results of the positive mass
balance cases examined in this paper, qualitatively similar results are presented. As the prescribed surface mass balance at the lower *a*_{L} becomes smaller, errors in the annual layer thickness become clearer around the middle depth.
The *λ* for the RCIP+corr (RCIP) scheme at the
1600 m
depth and below does not extend to the reference gray line
(Fig. 13c).
This is due to the lack of sufficient vertical resolution when capturing the
variation.
However, the results are still better than those obtained from the other schemes.

## 3.7 Resolution

Annual layer thickness becomes smaller with depth, which reflects the vertical velocity profile. Therefore, differences in age between two neighboring levels become larger with increasing depth. At a certain depth, the grid spacing becomes insufficient to hold the variation of the input age cycles, which means that the preservation of the input variation is lost below that depth. Typically, in the experiments shown above, 100 kyr cycle properties of input surface mass balance are maintained at around 300 to 500 kyr with the RCIP+corr (RCIP) scheme, and the computed age becomes smoother before that age (e.g., the square-wave shape in Fig. 10). The results obtained by UP-2 show loss of variation at similar or shallower depths, and those by UP-1 do so at even shallower depths, which results from the numerical diffusion of the schemes.

In the same manner as computing an approximate depth–age solution under
constant surface/basal mass balance and constant thickness (e.g., the gray
benchmark lines in Fig. 7a),
the inverse age–depth solution can be also computed.
Using this solution, the vertical profiles of layers that are
sufficient to hold a constant
age difference (*T*_{res}) can be obtained.
Figure 14b shows four gray lines, which correspond
to Δ*ζ* sufficient to hold the 10, 5, 2, and 1 kyr
differences when the experiment configuration is constant at *M*_{s}=1.5 cm yr^{−1},
*M*_{b}=0, and *H*=3000 m.
It is worth mentioning that for other *M*_{s} cases with the same
*M*_{b} and *H* constant, the four reference *ζ*–Δ*ζ*
relationships correspond to those with *T*_{res} divided by the factor of
*M*_{s} for *M*_{s}=3 cm yr^{−1}, and that they are interpreted as *T*_{res}=5, 2.5,
1, and 0.5 kyr, respectively, while
*M*_{s}=0.75 cm yr^{−1} are interpreted as *T*_{res}=20, 10
4, and 2 kyr, respectively.
Therefore, the Δ*ζ* limit estimated by using the lower surface mass
balance of the experiment should be considered.
For example, if the grid size at a certain *ζ* is larger than
the *T*_{res}=2 kyr line, characteristics with higher frequencies than
*T*_{res} cannot be sampled.
The vertical line marked as Z: 129e in Fig. 14b
corresponds to uniform grid spacing of 129 levels adopted in the
experiments conducted thus far. As the figure shows, this discretization
can hold 2, 5, and 10 kyr differences by *ζ*∼0.82, 0.44, and
0.29, respectively, and the 1 kyr is not resolved.

Figure 15 is the same as
Fig. 7c, except for the results using the different *P*_{T} of 50, 20, and 10 kyr.
As shown in the figure, higher-frequency properties disappear at
shallower depths. The results of RCIP (RCIP+corr) keep the
oscillation relatively stable, but the computed annual layer thicknesses are
not on the lines of the constant mass balance cases (gray lines), even at shallow
depths, for high-frequency input (Fig. 15c).
The square-wave shape pattern seems to be well preserved, at least around
the 1700 m depth (*ζ*∼0.44) in
Fig. 15a, and around the 600 m depth
(*ζ*∼0.8) in the case of Fig. 15b (which is
beyond the range of the figure but presented in the Supplement).
Therefore, by comparison with Fig. 14, it can be roughly
estimated that *T*_{res}=5 kyr and
*T*_{res}=2 kyr or longer are necessary to resolve the
*P*_{T}=50 kyr and *P*_{T}=20 kyr square-wave shapes,
respectively, which correspond to 1∕10*P*_{T}.

Here, the same series of experiments is repeated using a higher resolution
and a uniform grid spacing of 513 levels, which is 4 times
the resolution of the previous experiments.
The vertical line marked as Z: 519e in Fig. 14b
corresponds to this grid spacing. As the figure shows, this discretization can hold 5, 2, and 1 kyr differences by *ζ*∼0.19, 0.33, and 0.51, respectively, corresponding to 2430, 2010, and 1470 m in depth, respectively. The time step for higher-resolution experiments
conducted hereafter is set as 25 years.

Figure 16 is the same as
Fig. 15, except for the vertical grid spacing
adopted and with zooming shown near the bottom part.
The patterns seem to be well preserved by around the depths corresponding to
the *ζ* above, but the computed ages produced by the RCIP+corr (RCIP) schemes are not on the line *M*_{s}=*a*_{L} below the corresponding depths.

The number of vertical layers presented above exceeds 100,
which is substantially more than those used in many operational large-scale
3-D ice-sheet models. The typical number of layers is 30 or even less
(e.g., Goelzer et al., 2020; Seroussi et al., 2020).
Therefore, since it would be helpful to evaluate performance levels
at such a lower resolution for a broader range of applications,
a series of experiments was performed using a lower resolution.
Figure 17 is the same as
Fig. 7,
except that the results are provided using a uniform grid spacing of 33 levels (i.e., Δ*z*=93.75 m), which is one-fourth
the resolution of the reference experiment.
The time step for the lower-resolution experiments was set to
200 years.

When compared to the higher-resolution cases, the annual layer thickness patterns seem to be less preserved. The square-wave pattern in the results of UP-1 has already disappeared at around 1400 m depth, while those in the other schemes are almost the same, although less than the 129-level cases.

Figure 14b also contains a line corresponding
to a uniform grid spacing of 33 levels marked
as Z: 33e. As the figure shows, this discretization can hold
10 kyr differences by *ζ*∼0.67, corresponding to
990 m in depth. Similarly, it holds 20 kyr differences for
*ζ*∼0.43 (1710 m, not shown in the figure). Figure 17c shows that the result
patterns obtained for the RCIP+corr (RCIP) schemes are well preserved by around the 1400 m depth, which is between the two vertical levels of these estimations. Thus, like the 129-level cases, it can be roughly estimated that 1∕10 (or slightly more)
duration of the input cycle is necessary
to resolve by one grid.

A comparison between Figs. 17b and 7b shows that
differences in computed age among the schemes are within a comparable range
(*≲* 10 kyr) and thus are neither exaggerated nor
converged by reducing this resolution.
The high-frequency oscillation seen near the bottom in
Fig. 17b
is not in Fig. 17b,
which reflects the fact that
even the RCIP+corr scheme cannot
preserve the input shape near the bottom with the lower resolution.
Nevertheless, it should be emphasized that there are still systematic 1 to 10 kyr biases are left by the upwind schemes.

## 3.8 Non-uniform discretization

So far, all of the experiments were performed with uniform
discretization of either 129 or 513 levels.
For most cases in the present study,
it is reasonable to adopt non-uniform discretization,
which means large spacing toward the top and small spacing toward the bottom.
Since it was previously estimated that at least 1∕10
duration of the input cycle is necessary to resolve one grid,
the discretization can be optimized according to the Δ*ζ*
profile computed with the minimum surface mass balance of the experiment.
Here, for example, the target experiments are set as the square-wave surface
mass balance with *P*_{T}=100 kyr, ${P}_{\mathrm{H}}:{P}_{\mathrm{L}}=\mathrm{1}:\mathrm{1}$,
*M*_{b}=0, *H*=3000 m constant, and
${a}_{\mathrm{H}},{a}_{\mathrm{L}}=\mathrm{3},\mathrm{1.5}$ cm yr^{−1}.
This is the same configuration seen in Fig. 7.
For this configuration, the combination of ${T}_{\text{res}}={P}_{\mathrm{T}}/\mathrm{20}=\mathrm{5}$ kyr
and *M*_{s}=1.5 cm yr^{−1} is adopted in order to compute the reference profile in the same way as Fig. 14.
This number, which is half the number of the estimates used in
the discussion above, was chosen for safety and to facilitate additional experiments with other configurations (e.g., *a*_{L}=0.75 cm yr^{−1} and/or
*P*_{T}=50 kyr, which are archived in the Supplement).

Two non-uniform discretization types are adopted in this section.
One is a non-smooth grid and the other is a smooth grid (introduced in
Sect. 2.2). Various methods can be applied for non-smooth discretization. A very simple method, which is adopted in this study, calls for starting an initial spacing from the top and then keeping the same grid spacing as long as it is smaller than the reference profile. When the spacing exceeds the reference, it is halved
from the coordinates and maintained at that size until it exceeds the
reference again. It is necessary to limit the minimum grid spacing in order to avoid an infinite number of discretizations.
The line marked as Z: 477o in Fig. 14
is a computed profile following the method mentioned above, which runs from
$\mathrm{\Delta}\mathit{\zeta}={\mathrm{2}}^{-\mathrm{6}}$ to 2^{−11}. It contains 477 levels from
*ζ*=0 to *ζ*=1. The vertical coordinate system *Z* (model
logical coordinate; see Sect. 2.2) is identical
to *ζ*, and the series of *ζ*_{k} is shown in
Fig. 14a.

For non-uniform smooth discretization, a transformation function that
follows the reference profile is necessary between *ζ* and *Z*. Since there is no fixed method for choosing the formulation of *Z*, the following three constraints are adopted
for this paper:
(i) $\mathit{\zeta}(Z=\mathrm{0})=\mathrm{0}$,
(ii) $\mathit{\zeta}(Z=\mathrm{1})=\mathrm{1}$,
and (iii) $\frac{\mathrm{d}\mathit{\zeta}}{\mathrm{d}Z}>\mathrm{0}$ for $\mathrm{0}\le Z\le \mathrm{1}$.
A simple formulation to satisfy these constraints is

where the two parameters *γ* (a weight) and *ψ* (a power) control the
shape of transformation.
The linear term *Z* in Eq. (63) is
needed to avoid infinite $\frac{\partial Z}{\partial \mathit{\zeta}}$ at $\mathit{\zeta}=Z=\mathrm{0}$, which is
used for in Eq. (45).
After some trial and error with changing *γ* and *ψ* chosen from
integer numbers,
we found the following formulation can be used to maintain a grid spacing
that is less than or equal to the reference profile until approaching the bottom:

The line marked as Z: 513p in
Fig. 14 is the profile obtained by
Eq. (64).
Here, the same number of levels (513) is adopted to discretize under
the *Z*-coordinate system.
Figure 14a
shows the uniform grid spacing of *Z*, which
corresponds to the non-uniform *ζ* grid spacing achieved by this method.
Additionally, Fig. 14 marks the two vertical coordinates as references in order to reach 1000 and 2000 kyr, respectively, under the constant condition *M*_{s}=1.5 cm yr^{−1}, *M*_{b}=0, and *H*=3000 m.

Figures 18
and 19 show
the results obtained using the uniform-spacing (Z: 513e), smooth-grid (Z: 513p), and non-smooth-grid (Z: 477o) discretization methods,
in terms of *λ* vs. depth, and *λ* vs. age.
A comparison with Fig. 7 shows that the latter
preserves the input shape deeper than the former.
As shown in Fig. 14, the uniform discretization
case (marked as Z: 513e) is expected to fail to resolve 10 kyr
at age 1000 kyr, which is presented in
Fig. 19.
The results of UP-1 preserved the input shape deeper than the lower
resolution, which are almost only half of those achieved with the RCIP+corr (RCIP) methods. The results of UP-2 were preserved at slightly deeper depths than those of UP-1.
However, their phases are shown to be shifted from those of the RCIP+corr (RCIP)
methods, particularly at the deeper part.
For differences in the computed age from the RCIP method (Fig. 19a and b),
quantitatively, the same performance levels as the lower-resolution experiment were obtained by the other methods. The UP-1 and UP-2 methods deviate from the RCIP method
by around 1 and 10 kyr at most,
while RCIP+corr deviates by around 100 years.
Using non-uniform discretization, preservation of the input shape is
further extended to the deeper part (Fig. 18).
As shown in Fig. 14, the non-smooth non-uniform
discretization case (marked as Z: 477o) crosses the *T*_{res}=10 kyr
line at *ζ*∼0.07 (i.e., 2790 m depth)
for the *a*_{L}=1.5 cm yr^{−1} experiment, which is observed in
Fig. 19.
The smooth non-uniform discretization case (marked as Z: 513p)
crosses the *T*_{res}=10 kyr line slightly below that depth, *ζ*∼0.06 (2820 m depth), which is
observed in Fig. 19 again.
In addition, similar to the lower-resolution experiment,
Fig. 19 shows that the RCIP+corr (RCIP) scheme performs relatively better in terms of the phases than those with the UP-2 and UP-1 schemes.

The present study demonstrates a method for performing 1-D age computations of ice sheets under constant velocity, variable velocity responding to transient changes in surface mass balance, and/or changes in ice-sheet thickness. Herein, comparisons of the vertical profiles of computed ages, as well as annual layer thicknesses, were examined among the RCIP schemes (semi-Lagrangian) and upwind schemes (Eulerian). Although the experiments in the present study were limited to 1-D computations under summits, we believe the characteristics of the RCIP schemes have been presented sufficiently to allow evaluations of their performance levels.

Overall, the RCIP schemes show the best performance levels among the schemes examined in the present study. In particular, the computed vertical profiles of the annual layer thicknesses produced by RCIP schemes follow the expected depth profiles more reasonably than the other methods. This advantage reflects the design of the RCIP scheme, which explicitly computes the evolution of the age derivative, i.e., the inverse of annual layer thickness, using an advection equation that is similar to the one used to compute the age itself. Using the other schemes, the computed vertical profiles of annual layer thickness either show more smoothing at shallower depths than that were found with the RCIP scheme or the development of oscillation at steep changes in the input surface mass balance. Such oscillation development is shown even when the input is a smooth cosine-wave-type pattern and the amplitude is large. Since the slope filter adopted in this study is extremely simple, it is possible that the results obtained by the use of a second-order upwind scheme with a more suitable filter will change the characteristics. The introduction of slope limiters on general non-uniform discretization for higher-order upwind schemes is possible (Murman et al., 2005), but the conditions used for switching between a cubic polynomial and a rational form (Eqs. 11 and 12) in the RCIP scheme may be simpler and easier to implement. Under some configurations, oscillation development is not shown by the second-order upwind scheme. However, the phases of annual layer thickness against the age are shifted from those expected from the initial inputs, which again demonstrates the advantages of the RCIP scheme.

We examined two methods of computing the departure points in our RCIP scheme experiments. Under a constant velocity case, the results obtained by the simpler method show even less accurate solutions than the first-order upwind scheme, while the other “correction” method shows the best performance. The computed age difference between the two RCIP methods is 1000 years at most for all the configurations examined in the present study, including the vertical resolution. As a result, the simpler method still performs well if the expected accuracy of the application is less than that period. Under an evolving surface mass balance, the solution of the upwind scheme deviation is around 10 kyr, which is slightly larger.

As has already been discussed in previous studies (Greve et al., 2002), the first-order upwind scheme shows somewhat better performance than other schemes in some experiments. Greve et al. (2002) attributes this result to the cancellation of errors between discretization and numerical diffusion. In addition, from comparisons between results obtained via the first-order upwind scheme, with and without the midpoint rule (Fig. 4), we find that the midpoint rule does provide an advantage because the results obtained without the rule are worse by 1 order of magnitude than those obtained via the second-order upwind scheme. Furthermore, as discussed above, the upstream correction significantly improves the RCIP solution, which suggests that it is important to consider the non-constant velocity between the arrival and departure points. Since the midpoint rule formulation in the first-order scheme, in principle, corresponds to this upstream correction, they are consistent. The shape of the normalized vertical velocity profile also may play a role in the relative performance levels. For example, the upper part is more linear than the bottom part, which may increase the accuracy of the first-order approximation. In any case, it is clear that some or all of these points contribute to the higher performance of the first-order scheme, except for the bottom part.

As long as the annual layer thickness is not a concern, we feel that the classical upwind schemes are acceptable choices for use when dating. Note that using a first-order upwind scheme causes the structural details of the surface mass balance history to disappear very rapidly, but average features will compute quite well, except for near the bottom. The second-order scheme preserves the history better than the first-order scheme, but without an effective slope limiter, strange oscillations can appear in the results, as we have demonstrated in the present paper. However, in spite of these oscillations in the annual layer thickness, the results achieved by the second-order scheme are still slightly better than those for the first-order scheme throughout most of this study's experiments.

Greve et al. (2002) presented “practical suggestions” for numerical dating schemes: the second-order, the total variation diminishing Lax–Friedrichs (TVDLF) scheme with the minimum modulus (min-mod) filter, and even the first-order upwind schemes. In line with those recommendations, we would like to add the following additional practical suggestions. If good performance is required from the annual layer thickness computation, we strongly recommend the application of RCIP. We also strongly recommend the application of RCIP if age computations near the bottom are required to be within the error range of, e.g., 10 kyr. In other cases, the classical upwind schemes are acceptable choices.

The ice thickness and accumulation rate values used in the present
paper correspond to typical values found on the East Antarctic Plateau,
and the values used for the cycles in surface mass balance are 10 to
100 kyr. Providing appropriate scaling is used,
all of the results can be interpreted in the same way as those with different
configurations.
In order to simplify the situation, the set of
calculations will herein adopt a configuration
with a different magnitude and surface mass balance cycle
while keeping the same thickness and zero mass balance.
Figure 7 is taken as an example.
In this experiment, we will use a figure with the same shape
while replacing all the time-related variables to 1∕10 –
10 times higher accumulation rates (${a}_{\mathrm{H}},{a}_{\mathrm{L}}=\mathrm{30}$ and
15 cm yr^{−1}) with 1∕10 cycles, *P*_{T}=10 kyr,
${P}_{\mathrm{H}},{P}_{\mathrm{L}}=\mathrm{5}$ and 5 kyr.
The range of the horizontal axis needs to be adjusted from (a) 𝒜 from 0 to 100 kyr, (b) Δ𝒜 from −1 to
1 kyr, and (c) *λ* from 10^{0} to 10^{2} mm,
respectively (remember that the units of annual layer
thickness *λ* are substantially mm yr^{−1}).
Similarly, Fig. 15 can be interpreted as
the results of (a) *P*_{T}=5 kyr, (b) *P*_{T}=2 kyr and (c) *P*_{T}=1 kyr, respectively, providing that the horizontal axis is adjusted *λ* from 10^{0} to 10^{2} mm. In summary, the results in the present paper can be interpreted as cases of $\sim \mathrm{30}\phantom{\rule{0.125em}{0ex}}\mathrm{cm}\phantom{\rule{0.125em}{0ex}}{\mathrm{yr}}^{-\mathrm{1}}$ surface mass balance with cycles of 1 to 10 kyr,
i.e., millennial-scale climate oscillations on a typical Greenland
site. Under the scaled configuration,
it can be interpreted from Fig. 7b
that the age profiles produced by the RCIP scheme deviate
from RCIP+corr by less than 100 years throughout the depth range, which reflects the differences in computing the departure points.
The other two schemes deviate by around 1 kyr at most.
In addition, from Fig. 15a,
it can be seen that the square-wave shape pattern is well preserved, at least
around the 1700 m with *P*_{T}=5 (a) even though the
higher-frequency properties disappear at shallower depths with
*P*_{T}=1 kyr case (c). Moreover, from examining a series of experiments with 10 times higher accumulation under 1∕10 shorter cycles, we confirmed the same normalized shape of the results (not shown).

Although the focus of the present study is limited to 1-D age computations, implementation of the RCIP scheme for 3-D computation of the age field is also a suitable subject for future discussions.

Extension to 3-D would require the consideration of complex 3-D flow fields and typically much lower horizontal ice age gradients. In addition, the negative mass balance experiment demonstrated in the present study is too simple to be compatible with the 3-D situation. One important characteristic of the CIP scheme family is that the spatial gradient of the field variable (age in this case) is not a diagnostic (passive) value but is instead a prognostic field. Yabe et al. (2002) argued that even in an extreme case where values of the three adjacent grid points are zero, one wave still can exist, and thus non-zero spatial gradients can be held at these grid points. Therefore, it is speculated that the accuracy of the RCIP approach is not worse than that of other semi-Lagrangian schemes using higher-order interpolation techniques over the field variables, which have been discussed in past studies (Clarke and Marshall, 2002; Clarke et al., 2005). As described in the present study, RCIP is an effective scheme for preserving the flux information at deposition (annual layer thickness in the case of dating). However, detection of “points of origin” requires another technique, e.g., the back-tracing method. Huybrechts et al. (2007) suggested a very effective back-tracing method, which can be sufficient by itself for ice-core dating. The small but primary advantage of the RCIP method over the powerful back-tracing method is that it is a forward scheme. This means that it is not necessary to record all the past velocity field data during the simulation. Therefore, we consider the combination of the high-precision forward scheme and the powerful backward scheme to be a good choice when the objective is to obtain rough and detailed pictures of ice age fields.

Furthermore, it is expected that the RCIP scheme will be applicable to other advection problems in ice-sheet modeling. The evolutions of ice-sheet thickness and temperature are formulated using transport or advection equations, which are also good candidates for extending the discussion of this study. For such cases, researchers may be interested in mass or energy conservation in the field. Actually, a multi-dimensional conservative formulation of CIP schemes has already been proposed (Yabe et al., 2002). Accordingly, the implementation of the scheme to 3-D age and temperature fields in numerical ice-sheet models has already been set as the next target of our development.

A time-splitting technique (Eqs. 2 and 3) is at the core of the CIP schemes
and is somewhat difficult to understand at first glance.
We will attempt to clarify matters with the following
simple explanation.
If it is assumed that the non-advection term *h*(*x*,*t*) satisfies at least
locally $h(x,t)={h}^{*}\left(t\right)$, i.e., it is not dependent on *x*,
a new variable ${f}^{*}(x,t)$ can be introduced such that

Then, by introducing Eq. (A1) into the original advection
equation (Eq. 1), a pure advection form of *f*^{*} can be
obtained:

which is the same form as the advection phase equation (Eq. 2).
Using a semi-Lagrangian algorithm, solving Eq. (A2)
for *f*^{*} at time *t*+Δ*t* requires ${f}^{*}(x,t)$, which is
identical to *f*(*x*,*t*) by cancellation of the integral term of
Eq. (A1). Therefore, Eq. (A2) is
solved by the identical procedure used for Eq. (2).
After solving ${f}^{*}(t+\mathrm{\Delta}t)$, *f*(*t*+Δ*t*) can be computed using
Eq. (A1), such that

which is the same as the non-advection
equation (Eq. 3) and the solution (Eq. 21)
where *f*(*x*) is integrated with the initial condition ${f}_{j}^{*}$.

“Machine epsilon” is defined as the smallest *ϵ* in a computer
such that $\mathrm{1}+\mathit{\u03f5}>\mathrm{1}$ under floating-point arithmetic.
Similarly, an arbitrarily number *f* has the corresponding smallest number
(hereafter *ε*_{f}) which satisfies $f+{\mathit{\epsilon}}_{f}>f$.
In very rare cases, the authors observed that the age at the upwind grid point becomes close to the value at the target grid point, which differs by
*ε*_{f} (i.e., ${f}_{j+\mathrm{1}}={f}_{j}+{\mathit{\epsilon}}_{f}$).
Since no representative value exists between *f*_{j+1} and
*f*_{j} under floating-point arithmetic, the upwind value is either
*f*_{j} or *f*_{j+1}. Sometimes, *f*_{j+1} corresponds to
a value at a grid point that is too far away to be transported.
If there is an accumulation of errors of this type, the computed age
may show unexpected oscillations.

Although rounding up very small differences may be a possible solution for such cases, a different approach was adopted in the present study. After some trials, the authors finally adopted the following procedure for avoiding such oscillations, which (to the degree they used it) worked better than the rounding-up procedure. In the numerical model of the present paper, Eq. (14) is transformed as follows:

where *C*_{1} and *C*_{0} are substituted using
Eqs. (18) and (19).
The second term *δ**f* can be computed as the difference
between *f*_{j} and *f*_{j+1}. When *δ**f* is non-zero but
sufficiently small, i.e., less than *ε*_{f},
the value *f*_{j}+*δ**f* is maintained as *f*_{j}.
After simple reformulation, *F*_{j}(*X*) in the model code is finally
formulated as

where new series of constants are

respectively. When ${g}_{j+\mathrm{1}}-{S}_{j}=\mathrm{0}$ (and *α*=1), the coefficients lead to

respectively, and using this combination, *F*_{j}(*X*) is formulated as

which means a linear profile is adopted, regardless of *g*_{j}.

All the numerical experiments in the present paper are performed with
IcIES-2/JP version 0, which is a subset package of the
Ice sheet model for Integrated Earth system Studies II (IcIES-2).
IcIES-2 is available at
https://github.com/saitofuyuki/icies2 (last access: 26 September 2020, Saito and Abe-Ouchi, 2020) under Apache license
version 2.0.
The exact version of the full code, including the scripts used to run
the model for all the simulations in this paper, is archived on Zenodo
(https://doi.org/10.5281/zenodo.4034557, Saito et al., 2020).
Also it is tagged as `archive_gmd-2020`

in the git repository.

The supplement related to this article is available online at: https://doi.org/10.5194/gmd-13-5875-2020-supplement.

FS developed the ice-sheet model and then implemented the RCIP and other dating schemes in the model. FS performed numerical experiments designed by all the authors. The manuscript was written by FS with contributions from TO and AAO.

The authors declare that they have no conflicts of interest.

We would like to thank Shawn Marshall and an anonymous referee for their valuable comments, which have substantially improved our manuscript. This study was supported by the Japan Society for the Promotion of Science (JSPS) KAKENHI under grant nos. 17K05664, 17H06323, and 17H06104.

This research has been supported by the Japan Society for the Promotion of Science (JSPS) KAKENHI (grant nos. 17K05664, 17H06323, and 17H06104).

This paper was edited by Philippe Huybrechts and reviewed by Shawn Marshall and one anonymous referee.

Clarke, G. K. and Marshall, S. J.: Isotopic balance of the Greenland Ice Sheet: modelled concentrations of water isotopes from 30,000 BP to present, Quatern. Sci. Rev., 21, 419–430, https://doi.org/10.1016/S0277-3791(01)00111-1, 2002. a, b, c, d

Clarke, G. K. C., Lhomme, N., and Marshall, S. J.: Tracer transport in the Greenland ice sheet: three-dimensional isotopic stratigraphy, Quatern. Sci. Rev., 24, 155–171, https://doi.org/10.1016/j.quascirev.2004.08.021, 2005. a, b, c

Cuffey, K. M. and Paterson, W. S. B.: The Physics of Glaciers, Academic Press, 4th edn., 2010. a, b

Fischer, H., Severinghaus, J., Brook, E., Wolff, E., Albert, M., Alemany, O., Arthern, R., Bentley, C., Blankenship, D., Chappellaz, J., Creyts, T., Dahl-Jensen, D., Dinn, M., Frezzotti, M., Fujita, S., Gallee, H., Hindmarsh, R., Hudspeth, D., Jugie, G., Kawamura, K., Lipenkov, V., Miller, H., Mulvaney, R., Parrenin, F., Pattyn, F., Ritz, C., Schwander, J., Steinhage, D., van Ommen, T., and Wilhelms, F.: Where to find 1.5 million yr old ice for the IPICS ”Oldest-Ice” ice core, Clim. Past, 9, 2489–2505, https://doi.org/10.5194/cp-9-2489-2013, 2013. a, b

Goelzer, H., Nowicki, S., Payne, A., Larour, E., Seroussi, H., Lipscomb, W. H., Gregory, J., Abe-Ouchi, A., Shepherd, A., Simon, E., Agosta, C., Alexander, P., Aschwanden, A., Barthel, A., Calov, R., Chambers, C., Choi, Y., Cuzzone, J., Dumas, C., Edwards, T., Felikson, D., Fettweis, X., Golledge, N. R., Greve, R., Humbert, A., Huybrechts, P., Le clec'h, S., Lee, V., Leguy, G., Little, C., Lowry, D. P., Morlighem, M., Nias, I., Quiquet, A., Rückamp, M., Schlegel, N.-J., Slater, D. A., Smith, R. S., Straneo, F., Tarasov, L., van de Wal, R., and van den Broeke, M.: The future sea-level contribution of the Greenland ice sheet: a multi-model ensemble study of ISMIP6, The Cryosphere, 14, 3071–3096, https://doi.org/10.5194/tc-14-3071-2020, 2020. a

Greve, R. and Hutter, K.: Polythermal three-dimensional modelling of the Greenland ice sheet with varied geothermal heat flux, Ann. Glaciol., 21, 8–12, 1995. a

Greve, R., Wang, Y., and Mügge, B.: Comparison of numerical schemes for the solution of the advective age equation in ice sheets, Ann. Glaciol., 35, 487–494, https://doi.org/10.3189/172756402781817112, 2002. a, b, c, d, e, f, g, h

Huybrechts, P., Rybak, O., Pattyn, F., Ruth, U., and Steinhage, D.: Ice thinning, upstream advection, and non-climatic biases for the upper 89 % of the EDML ice core from a nested model of the Antarctic ice sheet, Clim. Past, 3, 577–589, https://doi.org/10.5194/cp-3-577-2007, 2007. a

Lhomme, N., Clarke, G. K. C., and Marshall, S. J.: Tracer transport in the Greenland Ice Sheet: constraints on ice cores and glacial history, Quatern. Sci. Rev., 24, 173–194, https://doi.org/10.1016/j.quascirev.2004.08.020, 2005. a

Mügge, B., Savvin, A., Calov, R., and Greve, R.: Numerical Age Computation of the Antarctic Ice Sheet for Dating Deep Ice Cores, in: Advances in Cold-Region Thermal Engineering and Sciences, edited by: Hutter, Y. W. K. and Beer, H., Springer, 307–318, 1999. a, b, c

Murman, S., Berger, M., and Aftosmis, M.: Analysis of slope limiters on irregular grids, in: Technical Report NAS-05-007, NAS, 2005. a, b

Parrenin, F., Dreyfus, G., Durand, G., Fujita, S., Gagliardini, O., Gillet, F., Jouzel, J., Kawamura, K., Lhomme, N., Masson-Delmotte, V., Ritz, C., Schwander, J., Shoji, H., Uemura, R., Watanabe, O., and Yoshida, N.: 1-D-ice flow modelling at EPICA Dome C and Dome Fuji, East Antarctica, Clim. Past, 3, 243–259, https://doi.org/10.5194/cp-3-243-2007, 2007. a, b, c

Rybak, O. and Huybrechts, P.: A comparison of Eulerian and Lagrangian methods for dating in numerical ice-sheet models, Ann. Glaciol., 37, 150–158, https://doi.org/10.3189/172756403781815393, 2003. a, b, c, d

Saito, F. and Abe-Ouchi, A.: Ice sheet model for Integrated Earth system Studies II (IcIES-2), GitHub, available at: https://github.com/saitofuyuki/icies2, last access: 26 September 2020. a

Saito, F., Obase, T., and Abe-Ouchi, A.: Resources relating to Saito et al. (2020, GMD) – IcIES-2/JP version 0 and results (Version 7879532a), Zenodo, https://doi.org/10.5281/zenodo.4034557, 2020. a

Seroussi, H., Nowicki, S., Payne, A. J., Goelzer, H., Lipscomb, W. H., Abe-Ouchi, A., Agosta, C., Albrecht, T., Asay-Davis, X., Barthel, A., Calov, R., Cullather, R., Dumas, C., Galton-Fenzi, B. K., Gladstone, R., Golledge, N. R., Gregory, J. M., Greve, R., Hattermann, T., Hoffman, M. J., Humbert, A., Huybrechts, P., Jourdain, N. C., Kleiner, T., Larour, E., Leguy, G. R., Lowry, D. P., Little, C. M., Morlighem, M., Pattyn, F., Pelle, T., Price, S. F., Quiquet, A., Reese, R., Schlegel, N.-J., Shepherd, A., Simon, E., Smith, R. S., Straneo, F., Sun, S., Trusel, L. D., Van Breedam, J., van de Wal, R. S. W., Winkelmann, R., Zhao, C., Zhang, T., and Zwinger, T.: ISMIP6 Antarctica: a multi-model ensemble of the Antarctic ice sheet evolution over the 21st century, The Cryosphere, 14, 3033–3070, https://doi.org/10.5194/tc-14-3033-2020, 2020. a

Shashkov, M.: Conservative Finite-Difference Methods on General Grids, CRC press, Boca Raton, 1st edn., https://doi.org/10.1201/9781315140209, 1995. a, b

Staniforth, A. and Côté, J.: Semi-Lagrangian Integration Schemes for Atmospheric Models: A Review, Mon. Weather Rev., 119, 2206–2223, https://doi.org/10.1175/1520-0493(1991)119<2206:SLISFA>2.0.CO;2, 1990. a

Tarasov, L. and Peltier, W. R.: Greenland glacial history, borehole constraints, and Eemian extent, J. Geophys. Res., 108, 2143, https://doi.org/10.1029/2001JB001731, 2003. a, b, c

Toda, K., Ogata, Y., and Yabe, T.: Multi-dimensional conservative semi-Lagrangian method of characteristics {CIP} for the shallow water equations, J. Comput. Phys., 228, 4917–4944, https://doi.org/10.1016/j.jcp.2009.04.003, 2009. a

Xiao, F., Yabe, T., and Ito, T.: Constructing oscillation preventing scheme for advection equation by rational function, Comp. Phys. Comm., 93, 1–12, https://doi.org/10.1016/0010-4655(95)00124-7, 1996. a, b, c

Yabe, T. and Takei, E.: A New Higher-Order Godunov Method for General Hyperbolic Equations, J. Phys. Soc. Jpn., 57, 2598–2601, https://doi.org/10.1143/JPSJ.57.2598, 1988. a

Yabe, T., Xiao, F., and Utsumi, T.: The Constrained Interpolation Profile Method for Multiphase Analysis, J. Comput. Phys., 169, 556–593, https://doi.org/10.1006/jcph.2000.6625, 2001. a

Yabe, T., Ogata, Y., Takizawa, K., Kawai, T., Segawa, A., and Sakurai, K.: The next generation CIP as a conservative semi-Lagrangian solver for solid, liquid and gas, J. Comput. Appl. Math., 149, 267–277, https://doi.org/10.1016/S0377-0427(02)00535-6, 2002. a, b

^{1}

Some models adopt 0 for the right-hand side
(e.g., Rybak and Huybrechts, 2003) simply because they use a different
age definition. For such cases, redefining 𝒜 as 𝒜−*t*
results in an equation that is identical to Eq. (29).

- Abstract
- Introduction
- Model description
- Experiment and results
- Discussion and conclusion
- Appendix A: Notes on time splitting
- Appendix B: Implementation of the RCIP method in the present paper
- Code availability
- Author contributions
- Competing interests
- Acknowledgements
- Financial support
- Review statement
- References
- Supplement

- Abstract
- Introduction
- Model description
- Experiment and results
- Discussion and conclusion
- Appendix A: Notes on time splitting
- Appendix B: Implementation of the RCIP method in the present paper
- Code availability
- Author contributions
- Competing interests
- Acknowledgements
- Financial support
- Review statement
- References
- Supplement