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

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.

1 Introduction

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:

$$\begin{array}{}\text{(1)}& {\displaystyle \frac{\partial f}{\partial t}}+u(x,t){\displaystyle \frac{\partial f}{\partial x}}=h(x,t),\end{array}$$

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:

$$\begin{array}{}\text{(2)}& {\displaystyle }& {\displaystyle \frac{\partial f}{\partial t}}+u(x,t){\displaystyle \frac{\partial f}{\partial x}}=\mathrm{0},& {\displaystyle }& {\displaystyle }\text{the advection phase}\text{(3)}& {\displaystyle }& {\displaystyle \frac{\partial f}{\partial t}}=h(x,t)& {\displaystyle }& {\displaystyle }\text{the non-advection phase}.\end{array}$$

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}$:

$$\begin{array}{}\text{(4)}& {\displaystyle \frac{\partial g}{\partial t}}+u(x,t){\displaystyle \frac{\partial g}{\partial x}}=\widehat{h}(x,t)={\displaystyle \frac{\partial h}{\partial x}}-g{\displaystyle \frac{\partial u}{\partial x}}.\end{array}$$

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):

$$\begin{array}{}\text{(5)}& {\displaystyle }& {\displaystyle \frac{\partial g}{\partial t}}+u(x,t){\displaystyle \frac{\partial g}{\partial x}}=\mathrm{0},& {\displaystyle }& {\displaystyle }\text{advection phase}\text{(6)}& {\displaystyle }& {\displaystyle \frac{\partial g}{\partial t}}=\widehat{h}(x,t)& {\displaystyle }& {\displaystyle }\text{non-advection phase}.\end{array}$$

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

$$\begin{array}{}\text{(7)}& f(x,t+\mathrm{\Delta }t)=f(x_{\mathrm{dep}},t),\end{array}$$

where

$$\begin{array}{}\text{(8)}& x_{\mathrm{dep}}=x+\underset{t+\mathrm{\Delta }t}{\overset{t}{\int }}\mathrm{d}tu(x,t).\end{array}$$

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).

Figure 1 Schematic illustration of advection and semi-Lagrangian scheme. The new state computation for a target point x_j from time t_m in the case of u_j<0 is presented. The colors symbolically express the value of field variables. The boxes correspond to model grid points. The solid arrow is the trajectory of one particle and the solid circle is the departure point. In a semi-Lagrangian scheme, the distance to the departure point, 〈ξ〉, is computed using an assumed trajectory. Interpolating the state at departure point x_dep, the value is advected to the arrival point $(x_j,t_m+\mathrm{\Delta }t)$. The dotted arrow and circle correspond to the trajectory and departure point in a different case while assuming a constant velocity, which may lead to a different state.

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:

$$\begin{array}{}\text{(9)}& \left\{\begin{array}{rl}{\displaystyle }f(x_j,t+\mathrm{\Delta }t)& {\displaystyle }=f(x_j+\langle \mathit{\xi }\rangle ,t)=F_j(x_j+\langle \mathit{\xi }\rangle ),\\ {\displaystyle }g(x_j,t+\mathrm{\Delta }t)& {\displaystyle }=g(x_j+\langle \mathit{\xi }\rangle ,t)=G_j(x_j+\langle \mathit{\xi }\rangle ),\end{array}\right.\end{array}$$

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

$$\begin{array}{}\text{(10)}& \left\{\begin{array}{rlrl}{\displaystyle }{F}_j\left(x_j\right)& {\displaystyle }=f\left(x_j\right),& {\displaystyle }{F}_j\left(x_{j+\mathrm{1}}\right)& {\displaystyle }=f\left(x_{j+\mathrm{1}}\right),\\ {\displaystyle }{G}_j\left(x_j\right)& {\displaystyle }=g\left(x_j\right),& {\displaystyle }{G}_j\left(x_{j+\mathrm{1}}\right)& {\displaystyle }=g\left(x_{j+\mathrm{1}}\right).\end{array}\right.\end{array}$$

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

$$\begin{array}{}\text{(11)}& F_j\left(X\right)=C_{\mathrm{0}}+C_{\mathrm{1}}X+C_{\mathrm{2}}{X}^{\mathrm{2}}+C_{\mathrm{3}}{X}^{\mathrm{3}},\end{array}$$

where $X=x-x_j$. The four coefficients (C₀, C₁, C₂, and C₃) 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

$$\begin{array}{}\text{(12)}& F_j\left(X\right)={\displaystyle \frac{C_{\mathrm{0}}+C_{\mathrm{1}}X+C_{\mathrm{2}}{X}^{\mathrm{2}}}{\mathrm{1}+D_{\mathrm{1}}X}}.\end{array}$$

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

$$\begin{array}{}\text{(13)}& S_j={\displaystyle \frac{f_{j+\mathrm{1}}-f_j}{\mathrm{\Delta }{x}_{j+\frac{\mathrm{1}}{\mathrm{2}}}}},\mathrm{\Delta }{x}_{j+\frac{\mathrm{1}}{\mathrm{2}}}=x_{j+\mathrm{1}}-x_j.\end{array}$$

Additionally, the four coefficients (C₀, C₁, C₂, and D₁) are determined in order to satisfy the same constraints. The two interpolation functions (Eqs. 11 and 12) are integrated by introducing a switching parameter α:

$$\begin{array}{}\text{(14)}& F_j\left(X\right)={\displaystyle \frac{C_{\mathrm{0}}+C_{\mathrm{1}}X+C_{\mathrm{2}}{X}^{\mathrm{2}}+C_{\mathrm{3}}{X}^{\mathrm{3}}}{\mathrm{1}+\mathit{\alpha }{D}_{\mathrm{1}}X}}.\end{array}$$

The five coefficients used to satisfy the constraints are computed as

$$\begin{array}{}\text{(15)}& {\displaystyle }{D}_{\mathrm{1}}& {\displaystyle }={\displaystyle \frac{\mathrm{1}}{\mathrm{\Delta }{x}_{j+\frac{\mathrm{1}}{\mathrm{2}}}}}\left[\left|{\displaystyle \frac{S_j-g_j}{g_{j+\mathrm{1}}-S_j}}\right|-\mathrm{1}\right],\text{(16)}& {\displaystyle }{C}_{\mathrm{3}}& {\displaystyle }={\displaystyle \frac{g_j-S_j+(g_{j+\mathrm{1}}-S_j)(\mathrm{1}+\mathit{\alpha }{D}_{\mathrm{1}}\mathrm{\Delta }{x}_{j+\frac{\mathrm{1}}{\mathrm{2}}})}{\mathrm{\Delta }{x_{j+\frac{\mathrm{1}}{\mathrm{2}}}}^{\mathrm{2}}}},\text{(17)}& {\displaystyle }{C}_{\mathrm{2}}& {\displaystyle }=S_j\mathit{\alpha }{D}_{\mathrm{1}}+{\displaystyle \frac{S_j-g_j}{\mathrm{\Delta }{x}_{j+\frac{\mathrm{1}}{\mathrm{2}}}}}-C_{\mathrm{3}}\mathrm{\Delta }{x}_{j+\frac{\mathrm{1}}{\mathrm{2}}},\text{(18)}& {\displaystyle }{C}_{\mathrm{1}}& {\displaystyle }=g_j+f_j\mathit{\alpha }{D}_{\mathrm{1}},\text{(19)}& {\displaystyle }{C}_{\mathrm{0}}& {\displaystyle }=f_j.\end{array}$$

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:

$$\begin{array}{}\text{(20)}& \left\{\begin{array}{rl}{\displaystyle }{f}_j^{*}& {\displaystyle }=F_j(\langle \mathit{\xi }\rangle )={\displaystyle \frac{C_{\mathrm{0}}+C_{\mathrm{1}}\langle \mathit{\xi }\rangle +C_{\mathrm{2}}\langle \mathit{\xi }{\rangle }^{\mathrm{2}}+C_{\mathrm{3}}\langle \mathit{\xi }{\rangle }^{\mathrm{3}}}{\mathrm{1}+\mathit{\alpha }{D}_{\mathrm{1}}\langle \mathit{\xi }\rangle }},\\ {\displaystyle }{g}_j^{*}& {\displaystyle }=G_j(\langle \mathit{\xi }\rangle )={\displaystyle \frac{C_{\mathrm{1}}+\mathrm{2}{C}_{\mathrm{2}}\langle \mathit{\xi }\rangle +\mathrm{3}{C}_{\mathrm{3}}\langle \mathit{\xi }{\rangle }^{\mathrm{2}}}{\mathrm{1}+\mathit{\alpha }{D}_{\mathrm{1}}\langle \mathit{\xi }\rangle }}\\ {\displaystyle }& {\displaystyle }-{\displaystyle \frac{\mathit{\alpha }{D}_{\mathrm{1}}}{\mathrm{1}+\mathit{\alpha }{D}_{\mathrm{1}}\langle \mathit{\xi }\rangle }}{f}_j^{*}.\end{array}\right.\end{array}$$

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:

$$\begin{array}{}\text{(21)}& \left\{\begin{array}{rl}{\displaystyle \frac{f_j(t+\mathrm{\Delta }t)-f_j^{*}}{\mathrm{\Delta }t}}& {\displaystyle }=h_j\\ {\displaystyle \frac{g_j(t+\mathrm{\Delta }t)-g_j^{*}}{\mathrm{\Delta }t}}& {\displaystyle }={\widehat{h}}_j.\end{array}\right.\end{array}$$

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:

$$\begin{array}{}\text{(22)}& \langle \mathit{\xi }\rangle =-\underset{t}{\overset{t+\mathrm{\Delta }t}{\int }}\mathrm{d}tu(x,t).\end{array}$$

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

$$\begin{array}{}\text{(23)}& \langle \mathit{\xi }\rangle =-u_j\mathrm{\Delta }t.\end{array}$$

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

$$\begin{array}{}\text{(24)}& u\left(x\right)=u\left(x_j\right)+(x-x_j)u^{\prime }\text{for}{x}_j\le x\le x_{j+\mathrm{1}},\end{array}$$

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:

$$\begin{array}{}\text{(25)}& {\displaystyle \frac{\mathrm{d}u}{\mathrm{d}t}}={\displaystyle \frac{\mathrm{d}x}{\mathrm{d}t}}{u}^{\prime }=u^{\prime }u,\end{array}$$

which is solved as

$$\begin{array}{}\text{(26)}& u\left(t\right)=u\left(t_m\right)\exp \left[u^{\prime }(t-t_m)\right].\end{array}$$

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

$$\begin{array}{ll}\text{(27)}& {\displaystyle }\langle \mathit{\xi }\rangle =-\underset{t_m}{\overset{t_m+\mathrm{\Delta }t}{\int }}\mathrm{d}tu\left(t\right)& {\displaystyle }=-u\left(t_m\right)\mathrm{\Delta }t\left[{\displaystyle \frac{\exp \left(u^{\prime }\mathrm{\Delta }t\right)-\mathrm{1}}{u^{\prime }\mathrm{\Delta }t}}\right],{\displaystyle }& {\displaystyle }(\text{when}{u}^{\prime }\ne \mathrm{0})\\ \text{(28)}& {\displaystyle }\langle \mathit{\xi }\rangle =-u\left(t_m\right)\mathrm{\Delta }t,& {\displaystyle }(\text{when}{u}^{\prime }=\mathrm{0}).\end{array}$$

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 Model description

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¹:

$$\begin{array}{}\text{(29)}& {\displaystyle \frac{\mathrm{d}\mathcal{A}}{\mathrm{d}t}}=\mathrm{1},\end{array}$$

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,

$$\begin{array}{}\text{(30)}& {\displaystyle \frac{\partial \mathcal{A}}{\partial t}}+w(z,t){\displaystyle \frac{\partial \mathcal{A}}{\partial z}}=\mathrm{1},\end{array}$$

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

$$\begin{array}{}\text{(31)}& \mathit{\zeta }={\displaystyle \frac{z-b}{H}},\end{array}$$

Eq. (30) is reformulated as follows:

$$\begin{array}{}\text{(32)}& {\displaystyle \frac{\partial \mathcal{A}}{\partial \mathit{\tau }}}+\mathit{\omega }{\displaystyle \frac{\partial \mathcal{A}}{\partial \mathit{\zeta }}}=\mathrm{1},\end{array}$$

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

$$\begin{array}{}\text{(33)}& \mathit{\omega }=w{\displaystyle \frac{\partial \mathit{\zeta }}{\partial z}}+{\displaystyle \frac{\partial \mathit{\zeta }}{\partial t}},\end{array}$$

where derivatives of ζ are computed as

$$\begin{array}{}\text{(34)}& {\displaystyle \frac{\partial \mathit{\zeta }}{\partial z}}& {\displaystyle }={\displaystyle \frac{\mathrm{1}}{H}},\text{(35)}& {\displaystyle \frac{\partial \mathit{\zeta }}{\partial t}}& {\displaystyle }=-{\displaystyle \frac{\mathrm{1}}{H}}\left[{\displaystyle \frac{\partial b}{\partial t}}+\mathit{\zeta }{\displaystyle \frac{\partial H}{\partial t}}\right].\end{array}$$

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

$$\begin{array}{}\text{(36)}& w(z=h(t\left)\right)={\displaystyle \frac{\partial h}{\partial t}}+M_{\mathrm{s}}\left(t\right),\end{array}$$

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:

$$\begin{array}{}\text{(37)}& {\displaystyle \frac{\partial {\mathcal{A}}^{\prime }}{\partial \mathit{\tau }}}+\mathit{\omega }{\displaystyle \frac{\partial {\mathcal{A}}^{\prime }}{\partial \mathit{\zeta }}}=-{\displaystyle \frac{\partial \mathit{\omega }}{\partial \mathit{\zeta }}}{\mathcal{A}}^{\prime },\end{array}$$

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,

$$\begin{array}{}\text{(38)}& \mathcal{A}(\mathit{\zeta }=\mathrm{1})=\mathrm{0},\end{array}$$

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

$$\begin{array}{}\text{(39)}& {\displaystyle \frac{\partial \mathcal{A}}{\partial \mathit{\tau }}}=\mathrm{1},\end{array}$$

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

$$\begin{array}{}\text{(40)}& {\mathcal{A}}^{\prime }(\mathit{\zeta }=\mathrm{1})=-{\displaystyle \frac{\mathrm{1}}{M_{\mathrm{s}}}},\end{array}$$

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:

$$\begin{array}{}\text{(41)}& \mathrm{0}\equiv {\mathit{\zeta }}_{\mathrm{0}}<{\mathit{\zeta }}_{\mathrm{1}}<\mathrm{\dots }<{\mathit{\zeta }}_{N_k-\mathrm{1}}\equiv \mathrm{1},\end{array}$$

and

$$\begin{array}{}\text{(42)}& \mathrm{\Delta }{\mathit{\zeta }}_{k+\mathrm{1}/\mathrm{2}}={\mathit{\zeta }}_{k+\mathrm{1}}-{\mathit{\zeta }}_k\text{for}k=\mathrm{0},\mathrm{\dots },N_k-\mathrm{2}.\end{array}$$

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:

$$\begin{array}{}\text{(43)}& {\displaystyle }& {\displaystyle \frac{\partial \mathcal{A}}{\partial T}}+W{\displaystyle \frac{\partial \mathcal{A}}{\partial Z}}=\mathrm{1},\text{(44)}& {\displaystyle }& {\displaystyle \frac{\partial {\mathcal{A}}^{\prime }}{\partial T}}+W{\displaystyle \frac{\partial {\mathcal{A}}^{\prime }}{\partial Z}}=-{\displaystyle \frac{\partial W}{\partial Z}}{\mathcal{A}}^{\prime },\end{array}$$

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

$$\begin{array}{}\text{(45)}& W=\mathit{\omega }{\displaystyle \frac{\partial Z}{\partial \mathit{\zeta }}}.\end{array}$$

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:

$$\begin{array}{}\text{(46)}& Z_k={\displaystyle \frac{k}{N_k-\mathrm{1}}},\text{for}k=\mathrm{0},\mathrm{\dots },N_k-\mathrm{1}.\end{array}$$

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:

$$\begin{array}{}\text{(47)}& {\left.W{\displaystyle \frac{\partial \mathcal{A}}{\partial Z}}\right|}_{Z=Z_k}\simeq W_{k+\mathrm{1}/\mathrm{2}}{\displaystyle \frac{{\mathcal{A}}_{k+\mathrm{1}}-{\mathcal{A}}_k}{\mathrm{\Delta }{Z}_{k+\frac{\mathrm{1}}{\mathrm{2}}}}}=W_{k+\mathrm{1}/\mathrm{2}}{A^{\prime }}_{\mathrm{I}}\left(Z_{k+\frac{\mathrm{1}}{\mathrm{2}}}\right),\end{array}$$

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:

$$\begin{array}{}\text{(48)}& {\left.W{\displaystyle \frac{\partial \mathcal{A}}{\partial Z}}\right|}_{Z=Z_k}\simeq W_k{A^{\prime }}_{\mathrm{I}}\left(Z_{k+\frac{\mathrm{1}}{\mathrm{2}}}\right),\end{array}$$

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

$$\begin{array}{}\text{(49)}& {\left.W{\displaystyle \frac{\partial \mathcal{A}}{\partial Z}}\right|}_{Z=Z_k}\simeq W_k{A^{\prime }}_{\mathrm{II}}\left(Z_k\right),\end{array}$$

where

$$\begin{array}{ll}\text{(50)}& {\displaystyle }& {\displaystyle }{A^{\prime }}_{\mathrm{II}}\left(Z_k\right){\displaystyle }& {\displaystyle }={\displaystyle \frac{\left(\mathrm{2}\mathrm{\Delta }{Z}_{k+\frac{\mathrm{1}}{\mathrm{2}}}+\mathrm{\Delta }{Z}_{k+\frac{\mathrm{3}}{\mathrm{2}}}\right){A^{\prime }}_{\mathrm{I}}\left(Z_{k+\frac{\mathrm{1}}{\mathrm{2}}}\right)-\mathrm{\Delta }{Z}_{k+\frac{\mathrm{1}}{\mathrm{2}}}{A^{\prime }}_{\mathrm{I}}\left(Z_{k+\frac{\mathrm{3}}{\mathrm{2}}}\right)}{\mathrm{\Delta }{Z}_{k+\frac{\mathrm{1}}{\mathrm{2}}}+\mathrm{\Delta }{Z}_{k+\frac{\mathrm{3}}{\mathrm{2}}}}}\\ {\displaystyle }& {\displaystyle }\text{for}k<N_k-\mathrm{2},\\ \text{(51)}& {\displaystyle }& {\displaystyle }{A^{\prime }}_{\mathrm{II}}\left(Z_k\right)=\mathrm{2}{A^{\prime }}_{\mathrm{I}}\left(Z_k\right)-{{\mathcal{A}}^{\prime }}_{k+\mathrm{1}},{\displaystyle }& {\displaystyle }\text{for}k=N_k-\mathrm{2},\end{array}$$

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 Experiment and results

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:

$$\begin{array}{}\text{(52)}& w\left(\mathit{\zeta }\right)=-\left[\left(M_{\mathrm{s}}+M_{\mathrm{b}}-{\displaystyle \frac{\partial H}{\partial t}}\right)\stackrel{\mathrm{̃}}{w}\left(\mathit{\zeta }\right)-M_{\mathrm{b}}\right],\end{array}$$

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{̃}}{w}\left(\mathit{\zeta }\right)$ is the normalized velocity profile. Assuming no basal sliding, $\stackrel{\mathrm{̃}}{w}\left(\mathit{\zeta }\right)$ can be approximated by

$$\begin{array}{}\text{(53)}& \stackrel{\mathrm{̃}}{w}\left(\mathit{\zeta }\right)=\mathrm{1}-{\displaystyle \frac{p+\mathrm{2}}{p+\mathrm{1}}}(\mathrm{1}-\mathit{\zeta })+{\displaystyle \frac{\mathrm{1}}{p+\mathrm{1}}}(\mathrm{1}-\mathit{\zeta }{)}^{p+\mathrm{2}},\end{array}$$

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{̃}}{w}$, as

$$\begin{array}{}\text{(54)}& {\displaystyle \frac{\partial \stackrel{\mathrm{̃}}{w}}{\partial \mathit{\zeta }}}={\displaystyle \frac{p+\mathrm{2}}{p+\mathrm{1}}}\left[\mathrm{1}-(\mathrm{1}-\mathit{\zeta }{)}^{p+\mathrm{1}}\right].\end{array}$$

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⁻¹ for the surface and basal mass balance, respectively, has the same normalized shape as that under 30 and 0 cm yr⁻¹, 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 Experimental results obtained using a uniform velocity of w=15 cm yr⁻¹. Snapshots of the computed vertical age profiles obtained by RCIP with correction, RCIP, second-, and first-order upwind schemes at t=5, 10, 15, 20 kyr are shown. Since the “correction factor” of the departure points is 1 (Eq. 28), the results of RCIP with correction are identical to those of RCIP. The results of the second-order upwind scheme are close to those of the RCIP, which are barely visible at this scale. The solution is also shown as a benchmark (thick gray line). Symbols are plotted for every eight vertical levels.

Figure 2 shows the computed age profile under the uniform velocity of −15 cm yr⁻¹ 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

$$\begin{array}{}\text{(55)}& \mathcal{A}(z,t)=min\left[t,\underset{h}{\overset{z}{\int }}\mathrm{d}{z}^{\prime }{\displaystyle \frac{\mathrm{1}}{w}}\right];\end{array}$$

thus, the exact solution for an uniform velocity is

$$\begin{array}{}\text{(56)}& \mathcal{A}(z,t)=min(t,-z/w_c),\end{array}$$

where $w_c=-\mathrm{15}$ cm yr⁻¹. 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 %).

Figure 3 Experimental results obtained using a uniform velocity of w=15 cm yr⁻¹. Snapshots of computed vertical age profiles obtained by RCIP with correction, RCIP, second-, and first-order upwind schemes at t=20 kyr relative to the exact solution are shown. The results of different time steps of 100, 50, and 25 years are shown for each scheme. The results of RCIP with correction are identical to those of RCIP. Age differences are shown on a logarithmic scale, except for those near 0, which are shown on a linear scale. Symbols are plotted for every eight vertical levels.

Figure 4 Experimental results obtained under steady vertical velocity profiles with H=3000 m, M_s=3 cm yr⁻¹, and M_b=0 mm yr⁻¹. (a) Computed vertical age profiles by RCIP with departure point correction (RCIP+corr), the RCIP scheme, the second-order upwind scheme (UP-2), the first-order upwind scheme (UP-1), and the first-order upwind scheme without the midpoint rule (UP-1n) at t=2000 kyr, (b) those relative to the benchmark profile obtained by numerical integration, and (c) a zoomed-in portion of the bottom part of panel (a) showing that the differences among the experiments are on the order of 10 kyr. Age differences are shown on a logarithmic scale, except for those near 0, which are shown on a linear scale. Uniform grid spacing of 129 levels is adopted in this simulation.

Figure 5 Same as Fig. 4a and b except for the experiment with M_b=3 mm yr⁻¹.

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⁻¹, 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⁻¹), 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⁻¹. The other two parameters are fixed. Surface mass balance is set as 3 cm yr⁻¹ 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.

$$\begin{array}{}\text{(57)}& \left\{\begin{array}{rlrl}{\displaystyle }{M}_{\mathrm{s}}\left(t\right)& {\displaystyle }=a_{\mathrm{H}},& {\displaystyle }& {\displaystyle }(\mathrm{0}\le \text{mod}(t,P_{\mathrm{T}})<P_{\mathrm{H}})\\ {\displaystyle }{M}_{\mathrm{s}}\left(t\right)& {\displaystyle }=a_{\mathrm{L}},& {\displaystyle }& {\displaystyle }(P_{\mathrm{H}}\le \text{mod}(t,P_{\mathrm{T}})<P_{\mathrm{T}}=P_{\mathrm{H}}+P_{\mathrm{L}}),\end{array}\right.\end{array}$$

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⁻¹. The other patterns examined in this paper are provided as the Supplement to this paper.

Figure 6 Schematic figure showing the time evolution of the surface mass balance adopted in these experiments. Only the first two cycles are plotted.

Figure 7 Results of transient experiments with square-wave surface mass balances of a_H, $a_{\mathrm{L}}=\mathrm{3},\mathrm{1.5}$ cm yr⁻¹, P_T=100 kyr, $P_{\mathrm{H}}:P_{\mathrm{L}}=\mathrm{1}:\mathrm{1}=\mathrm{50}:\mathrm{50}$ kyr, square-wave, M_b=0, and constant H=3000 m. (a) Vertical profiles of the computed age and (c) annual layer thickness at 1000 kyr using RCIP+corr, RCIP (overlapped on RCIP+corr), UP-2, and UP-1 are shown. The last 500 years are clipped from panel (c), where the age gradient (inverse of λ) is close to 0, reflecting the initial condition. (b) The computed age differences at the same depth relative to the result of the RCIP+corr case are shown on a logarithmic scale, except for those near 0, which are shown on a linear scale. For reference purposes, the gray lines indicate benchmark solutions for the constant surface mass balance cases of a_H, a_L, and a_M. Uniform grid spacing of 129 levels is adopted in this simulation.

Figure 8 Same as Fig. 7 but for the results of transient experiments using a square-wave surface mass balance of a_L=3 and 0.75 cm yr⁻¹.

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⁻¹ and 0.75 cm yr⁻¹. The higher surface mass balances and the basal are set as a_H=3 cm yr⁻¹ 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.

Figure 9 Same as Fig. 7 but for the results of transient experiments conducted with a cosine-wave surface mass balance of $a_{\mathrm{H}},a_{\mathrm{L}}=\mathrm{3}$ and 1.5 cm yr⁻¹.

Figure 10 Results of transient experiments with square-wave surface mass balances of (a) $a_{\mathrm{H}},a_{\mathrm{L}}=\mathrm{3}$ and 1.5 cm yr⁻¹ and (b) 3 and 0.75 cm yr⁻¹ with high- and low-value phase durations (Eq. 57) set as $P_{\mathrm{H}}:P_{\mathrm{L}}=\mathrm{1}:\mathrm{1}=\mathrm{50}:\mathrm{50}$ kyr. The basal mass balances are set as 0. The computed annual layer thickness of λ against the computing age is shown (RCIP overlaps with RCIP+corr). The gray lines indicate benchmark solutions of the constant surface mass balance cases of a_H and a_L, while the black lines $a_{\mathrm{M}}=(a_{\mathrm{H}}+a_{\mathrm{L}})/\mathrm{2}$ are provided as references. Uniform grid spacing of 129 levels is adopted in this simulation.

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

$$\begin{array}{}\text{(58)}& {\displaystyle \frac{\partial \mathit{\lambda }}{\partial t}}={\displaystyle \frac{\partial w}{\partial z}}\mathit{\lambda },\end{array}$$

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

$$\begin{array}{}\text{(59)}& {\displaystyle \frac{\partial w}{\partial z}}=-M_{\mathrm{s}}{\displaystyle \frac{\partial \stackrel{\mathrm{̃}}{w}}{\partial \mathit{\zeta }}}.\end{array}$$

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

$$\begin{array}{}\text{(60)}& {\displaystyle \frac{\partial \mathit{\lambda }}{\partial z}}={\displaystyle \frac{\mathit{\lambda }}{\stackrel{\mathrm{̃}}{w}}}{\displaystyle \frac{\partial \stackrel{\mathrm{̃}}{w}}{\partial \mathit{\zeta }}},\end{array}$$

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.

Figure 11 Prescribed time evolution patterns of ice thickness adopted in the non-steady thickness experiment. The thickness evolution is computed using an e-folding time of 10 kyr against square-wave and cosine-wave formulation of the surface mass balance with $a_{\mathrm{H}},a_{\mathrm{L}}=\mathrm{3}$ and 1.5 cm yr⁻¹ (Fig. 6) provided as an example.

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:

$$\begin{array}{}\text{(61)}& {\displaystyle \frac{\partial H}{\partial t}}=-{\displaystyle \frac{\mathrm{1}}{{\mathit{\tau }}_{\mathrm{H}}}}\left\{H-H_{\mathrm{ref}}\left(M_{\mathrm{s}}\right)\right\},\end{array}$$

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

$$\begin{array}{}\text{(62)}& H_{\mathrm{ref}}\left(M_{\mathrm{s}}\right)=H(t=\mathrm{0}){\left[{\displaystyle \frac{M_{\mathrm{s}}\left(t\right)}{M_{\mathrm{s}}(t=\mathrm{0})}}\right]}^{\mathrm{1}/(\mathrm{2}n+\mathrm{2})}.\end{array}$$

For cases where $H(t=\mathrm{0})=\mathrm{3000}$ m, $a_{\mathrm{H}},a_{\mathrm{L}}=\mathrm{3},\mathrm{1.5}$ cm yr⁻¹ 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⁻¹ or 3,0.75 cm yr⁻¹; τ_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⁻¹. 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.

Figure 12 Same as Fig. 7 but for the results of non-steady thickness experiments conducted with the response timescale set as τ_H=10 kyr. The reference thickness values H for the benchmark profiles are explained in the text.

Figure 13 Same as Fig. 7 but for the results of transient experiments conducted using a square-wave surface mass balance with a_L=0 cm yr⁻¹.

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⁻¹ 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⁻¹ 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⁻¹, 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.

Figure 14 Vertical discretization adopted in the present study: (a) Z vs. ζ (b) Δζ vs. ζ. Five patterns are shown: uniform grid spacing of 129 levels (129e), that of 513 levels (513e), a smooth non-uniform discretization (513p; see text), a non-smooth, non-uniform discretization (477o; see text), and uniform grid spacing of 33 levels (33e). Symbols are plotted for every 16 vertical levels. The four gray lines in panel (b) correspond to the layer thickness necessary to resolve 10, 5, 2, and 1 kyr differences under the condition of M_s=1.5 cm yr⁻¹, M_b=0, and H=3000 m. The two horizontal magenta lines correspond to the depth needed to reach 1000 and 2000 kyr under the same conditions.

Figure 15 Results of transient experiments conducted with the square-wave surface mass balances of $a_{\mathrm{H}},a_{\mathrm{L}}=\mathrm{3}$ and 1.5 cm yr⁻¹ and $P_{\mathrm{H}}:P_{\mathrm{L}}=\mathrm{1}:\mathrm{1}$ (Eq. 57), and the total duration as (a) P_T=50 kyr, (b) P_T=20 kyr, and (c) P_T=10 kyr. The basal mass balance is set as 0 for all the experiments in the figure. The vertical profiles of the annual layer thickness λ at 1000 kyr using RCIP+corr, RCIP (which overlaps with RCIP+corr), UP-2, and UP-1 are shown. The gray and black lines indicate benchmark solutions of constant surface mass balance cases with a_H, a_L, and $a_{\mathrm{M}}=(a_{\mathrm{H}}+a_{\mathrm{L}})/\mathrm{2}$ given as references. Uniform grid spacing of 129 levels is adopted in this simulation. The results covering the depth from 1000 to 2600 m are shown.

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⁻¹, 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⁻¹, and that they are interpreted as T_res=5, 2.5, 1, and 0.5 kyr, respectively, while M_s=0.75 cm yr⁻¹ 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 16 Same as Fig. 15, except for the vertical resolution, which is shown as a uniform grid spacing of 513 levels. The results at 2000 kyr covering the depth from 1000 to 2600 m are shown.

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∕10P_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.

Figure 17 Same as Fig. 7, except for the experiment with the lower vertical resolution of 33 uniform grid spacing. Symbols are plotted for every two vertical levels.

Figure 18 Same as Fig. 7c, except for the vertical resolution as (a) uniform grid spacing of 513 levels, (b) smooth non-uniform grid spacing of 513 levels, and (c) non-smooth non-uniform grid spacing of 477 levels. The results at 2000 kyr covering the depths from 2300 to 2900 m are shown.

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⁻¹. 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⁻¹ 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⁻¹ and/or P_T=50 kyr, which are archived in the Supplement).

Figure 19 Same as Fig. 10 except for the vertical resolution as (a) uniform grid spacing of 513 levels, (b) smooth non-uniform grid spacing of 513 levels, and (c) non-smooth non-uniform grid spacing of 477 levels. The results at 2000 kyr are shown.

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⁻¹¹. 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

$$\begin{array}{}\text{(63)}& \mathit{\zeta }={\displaystyle \frac{Z+\mathit{\gamma }{Z}^{\mathit{\psi }}}{\mathrm{1}+\mathit{\gamma }}},\end{array}$$

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:

$$\begin{array}{}\text{(64)}& \mathit{\zeta }\equiv {\displaystyle \frac{\left(Z+\mathrm{4}{Z}^{\mathrm{14}}\right)}{\mathrm{5}}},Z_k={\displaystyle \frac{k}{N_k-\mathrm{1}}}\text{for}k=\mathrm{0},\mathrm{\dots },N_k-\mathrm{1}.\end{array}$$

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⁻¹, 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⁻¹ 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.

4 Discussion and conclusion

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⁻¹) 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⁰ to 10² mm, respectively (remember that the units of annual layer thickness λ are substantially mm yr⁻¹). 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⁰ to 10² mm. In summary, the results in the present paper can be interpreted as cases of $\sim \mathrm{30}\mathrm{cm}{\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.

Appendix A: Notes on time splitting

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

$$\begin{array}{}\text{(A1)}& f^{*}(x,t)=f(x,t)-\int \mathrm{d}th(x,t)\simeq f(x,t)-\int \mathrm{d}t{h}^{*}\left(t\right).\end{array}$$

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

$$\begin{array}{}\text{(A2)}& {\displaystyle \frac{\partial f^{*}}{\partial t}}+u(x,t){\displaystyle \frac{\partial f^{*}}{\partial x}}=\mathrm{0},\end{array}$$

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

$$\begin{array}{}\text{(A3)}& f(x,t+\mathrm{\Delta }t)=f^{*}(x,t+\mathrm{\Delta }t)+\underset{t}{\overset{t+\mathrm{\Delta }t}{\int }}\mathrm{d}t{h}^{*}\left(t\right),\end{array}$$

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^{*}$.

Appendix B: Implementation of the RCIP method in the present paper

“Machine epsilon” is defined as the smallest ϵ in a computer such that $\mathrm{1}+\mathit{ϵ}>\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:

$$\begin{array}{}\text{(B1)}& F_j\left(X\right)=f_j+{\displaystyle \frac{g_{j}X+C_{\mathrm{2}}{X}^{\mathrm{2}}+C_{\mathrm{3}}{X}^{\mathrm{3}}}{\mathrm{1}+\mathit{\alpha }{D}_{\mathrm{1}}X}}=f_j+\mathit{\delta }f,\end{array}$$

where C₁ and C₀ 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

$$\begin{array}{}\text{(B2)}& F_j\left(X\right)=f_j+{\displaystyle \frac{{\widehat{C}}_{\mathrm{1}}X+{\widehat{C}}_{\mathrm{2}}{X}^{\mathrm{2}}+{\widehat{C}}_{\mathrm{3}}{X}^{\mathrm{3}}}{{\widehat{D}}_{\mathrm{0}}+\mathit{\alpha }{\widehat{D}}_{\mathrm{1}}X}},\end{array}$$

where new series of constants are

$$\begin{array}{ll}\text{(B3)}& {\displaystyle }{\widehat{D}}_{\mathrm{0}}& {\displaystyle }=\left|g_{j+\mathrm{1}}-S_j\right|,\text{(B4)}& {\displaystyle }{\widehat{D}}_{\mathrm{1}}& {\displaystyle }=D_{\mathrm{1}}\left|g_{j+\mathrm{1}}-S_j\right|={\displaystyle \frac{\left|S_j-g_j\right|-\left|g_{j+\mathrm{1}}-S_j\right|}{\mathrm{\Delta }{x}_{j+\frac{\mathrm{1}}{\mathrm{2}}}}},\text{(B5)}& {\displaystyle }{\widehat{C}}_{\mathrm{1}}& {\displaystyle }=g_j\left|g_{j+\mathrm{1}}-S_j\right|,\text{(B6)}& {\displaystyle }{\widehat{C}}_{\mathrm{2}}& {\displaystyle }=C_{\mathrm{2}}\left|g_{j+\mathrm{1}}-S_j\right|{\displaystyle }& {\displaystyle }=S_j\mathit{\alpha }{\widehat{D}}_{\mathrm{1}}+{\displaystyle \frac{(S_j-g_j)\left|g_{j+\mathrm{1}}-S_j\right|}{\mathrm{\Delta }{x}_{j+\frac{\mathrm{1}}{\mathrm{2}}}}}-{\widehat{C}}_{\mathrm{3}},\\ \text{(B7)}& {\displaystyle }{\widehat{C}}_{\mathrm{3}}& {\displaystyle }=C_{\mathrm{3}}\left|g_{j+\mathrm{1}}-S_j\right|{\displaystyle }& {\displaystyle }={\displaystyle \frac{\left|g_{j+\mathrm{1}}-S_j\right|}{\mathrm{\Delta }{x}_{j+\frac{\mathrm{1}}{\mathrm{2}}}}}\left[g_j-S_j+(g_{j+\mathrm{1}}-S_j)+\mathit{\alpha }{\widehat{D}}_{\mathrm{1}}\mathrm{\Delta }{x}_{j+\frac{\mathrm{1}}{\mathrm{2}}}\right],\end{array}$$

respectively. When $g_{j+\mathrm{1}}-S_j=\mathrm{0}$ (and α=1), the coefficients lead to

$$\begin{array}{}\text{(B8)}& {\displaystyle }{\widehat{C}}_{\mathrm{3}}& {\displaystyle }=\mathrm{0},{\widehat{C}}_{\mathrm{1}}=\mathrm{0},{\widehat{D}}_{\mathrm{0}}=\mathrm{0},\text{(B9)}& {\displaystyle }{\widehat{C}}_{\mathrm{2}}& {\displaystyle }=S_j{\displaystyle \frac{\left|S_j-g_j\right|}{\mathrm{\Delta }{x}_{j+\frac{\mathrm{1}}{\mathrm{2}}}}},\text{(B10)}& {\displaystyle }{\widehat{D}}_{\mathrm{1}}& {\displaystyle }={\displaystyle \frac{\left|S_j-g_j\right|}{\mathrm{\Delta }{x}_{j+\frac{\mathrm{1}}{\mathrm{2}}}}},\end{array}$$

respectively, and using this combination, F_j(X) is formulated as

$$\begin{array}{}\text{(B11)}& F_j\left(X\right)=f_j+{\displaystyle \frac{{\widehat{C}}_{\mathrm{2}}{X}^{\mathrm{2}}}{{\widehat{D}}_{\mathrm{1}}X}}=f_j+S_{j}X,\end{array}$$

which means a linear profile is adopted, regardless of g_j.