the Creative Commons Attribution 4.0 License.
the Creative Commons Attribution 4.0 License.
Implementation of the RCIP scheme and its performance for 1D age computations in icesheet models
Takashi Obase
Ayako AbeOuchi
Icesheet 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 icesheet 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 functionbased constrained interpolation profile (RCIP) scheme, to onedimensional ice age computation and demonstrates its performance levels via comparisons with those obtained from first and secondorder 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)  Fulltext XML

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 icecore properties is defining the age of ice along the depth of the ice sheet. This process is often called dating. Dating with numerical iceflow 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 iceflow models to evaluate potential “oldestice” study areas.
Various methods for use in icesheet model dating have been adopted and compared. Mügge et al. (1999) compared particle tracking (Lagrangian) and Eulerian schemes under simulated steadystate threedimensional (3D) 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 steadystate conditions and analytical solutions, as well as under different 3D 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 secondorder upwind schemes, Quadratic Upstream Interpolation for Convective Kinematics (QUICK), and total variation diminished (TVD) Lax–Wendroff (LW) schemes. From comparisons of the onedimensional (1D) steadystate age profiles produced by these schemes, they concluded that the secondorder upwind and TVDLW schemes performed well for typical icesheet age profiles. Comparisons among semiLagrangian schemes have also been performed. Introduction of a semiLagrangian trace scheme to icesheet 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 semiLagrangian 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 semiLagrangian scheme and discussed computed icecore age–depth relationships for the Greenland Ice Core Project (GRIP) ice core.
To date, various methods have been presented and demonstrated for use in icesheet 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 functionbased constrained interpolation profile (RCIP) method (Xiao et al., 1996) for use in 1D 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 1D advection equation with a nonadvection 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 (nonadvection) 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 semiLagrangian. The CIP scheme family corresponds to a semiLagrangian method variation. The basics of the semiLagrangian 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 semiLagrangian method in detail. Although a full description of the semiLagrangian 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 timesplitting algorithm (e.g., Yabe and Takei, 1988) in two phases as follows:
Appendix A presents a note on the timesplitting 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 nonadvection function $\widehat{h}(x,t)$ on the righthand side, which is solved using a timesplitting 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 nonadvection phases (Eqs. 3 and 6) will be discussed.
In semiLagrangian 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 semiLagrangian 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 nonadvection 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 semiLagrangian schemes. Another major topic common to the semiLagrangian 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 nonuniform 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 1D 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 largescale numerical icesheet 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 Dirichlettype 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 nonuniform. In the present paper, both types of discretization are examined. Since uniform discretization is a special case of nonuniform discretization, the latter can be described effectively without a loss of generality.
One way to introduce a nonuniform discretization is to apply a nonsmooth grid (Shashkov, 1995), which prescribes irregular discretization of the coordinates:
and
Another way to introduce a nonuniform 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 nonuniform smoothgrid 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 righthand side correspond to f, u, and h, respectively, in the RCIP scheme framework (Eq. 1). Although it is possible to introduce further nonuniform discretization on the Z coordinate, in the present paper, only a uniform discretization is examined on the smoothgrid discretization:
Actually, the discretization of the ζ coordinate corresponds to the special case of nonuniform smooth discretization with Z≡ζ. Therefore, for both uniform and nonuniform 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 secondorder upwind schemes, are examined in comparison with the RCIP schemes. While there are other numerical schemes suitable for such comparisons, including Lagrangian, other semiLagrangian, or even higherorder 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 “firstorder” 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 Dirichlettype 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 “secondorder” upwind scheme, the derivative of the age term is replaced by the secondorder 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 higherorder 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 firstorder 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 1D 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 steadystate 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 timerelated terms as 1∕10.
All the computations in our present study were performed on a personal computer (PC) equipped with an Intel Xeon E52609 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 secondorder upwind scheme (UP2), and the firstorder upwind scheme with midpoint rule (UP1). Additionally, a firstorder scheme without a midpoint rule (UP1n) is sometimes used. Multiple 1Dcolumn 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 UP1, UP2, 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 UP2 or UP1 run.
3.2 A verification experiment using uniform velocity
Before performing an experiment under a typical icesheet 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 UP1 scheme deviates the most from the solution. Specifically, it deviates 1 year from around twothirds 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 UP1 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 UP2 is even better at 7.5 to 154 years (0.04 % to 0.77 %).
3.3 Experiment with steady nonuniform vertical velocity
Hereafter, nonuniform 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 UP1 scheme shows the secondbest 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 UP1, 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 onedirection 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 secondorder accuracy. A “true” firstorder 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 firstorder accuracy. Figures 4b and 5b also contain results obtained using such a scheme, marked as UP1n. However, as expected, when such a normalgrid velocity is introduced for the advection equation, the results have less accuracy than those of the second order upwind (UP2). Furthermore, as shown in Greve et al. (2002), the improved performance of the UP1 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 UP2 scheme shows intermediate errors between RCIP and UP1 at the bottom.
3.4 Nonsteady surface mass balance experiments
This section presents the results of experiments conducted with nonsteady velocity profiles, which were performed with the prescribed surface mass balance time series. First, a very simple squarewave 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 lowvalue 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 squarewave 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 UP2 and UP1 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 UP2 and UP1 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 UP2 and UP1 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 UP2 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 UP2 scheme, which may exaggerate the oscillation of age gradients more than the simple firstorder Taylor expansion. For the UP1 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 cosinewave formulation of the surface mass balance (Fig. 6), which is relatively more continuous than the squarewave version. Similar performance levels were obtained by the UP2 and RCIP+corr (RCIP) methods for the small amplitude case (Fig. 9c). Instability also arises at lowtohigh transitions when the ratio of hightolow accumulation is larger (archived in the Supplement).
Figure 10 shows the results obtained by squarewave 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 UP2 and UP1 schemes. Dissipation at the discontinuity blur the squarewave shapes, particularly at the deeper part, but the phases are still maintained better by the RCIP+corr (RCIP) scheme than by the UP2 scheme.
3.5 Nonsteady 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 nonsteady 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 steadystate 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 squarewave or cosinewave 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 squarewave 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 nonsteady 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 nonpositive 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 icecore experiments, where the interpretation of icecore 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 steadystate icesheet 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 squarewave 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 squarewave shape in Fig. 10). The results obtained by UP2 show loss of variation at similar or shallower depths, and those by UP1 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, higherfrequency 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 highfrequency input (Fig. 15c). The squarewave 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 squarewave 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 higherresolution 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 largescale 3D icesheet 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 onefourth the resolution of the reference experiment. The time step for the lowerresolution experiments was set to 200 years.
When compared to the higherresolution cases, the annual layer thickness patterns seem to be less preserved. The squarewave pattern in the results of UP1 has already disappeared at around 1400 m depth, while those in the other schemes are almost the same, although less than the 129level 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 129level 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 highfrequency 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 Nonuniform 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 nonuniform 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 squarewave 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 nonuniform discretization types are adopted in this section. One is a nonsmooth grid and the other is a smooth grid (introduced in Sect. 2.2). Various methods can be applied for nonsmooth 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 nonuniform 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 Zcoordinate system. Figure 14a shows the uniform grid spacing of Z, which corresponds to the nonuniform ζ 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 uniformspacing (Z: 513e), smoothgrid (Z: 513p), and nonsmoothgrid (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 UP1 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 UP2 were preserved at slightly deeper depths than those of UP1. 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 lowerresolution experiment were obtained by the other methods. The UP1 and UP2 methods deviate from the RCIP method by around 1 and 10 kyr at most, while RCIP+corr deviates by around 100 years. Using nonuniform discretization, preservation of the input shape is further extended to the deeper part (Fig. 18). As shown in Fig. 14, the nonsmooth nonuniform 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 nonuniform 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 lowerresolution experiment, Fig. 19 shows that the RCIP+corr (RCIP) scheme performs relatively better in terms of the phases than those with the UP2 and UP1 schemes.
The present study demonstrates a method for performing 1D age computations of ice sheets under constant velocity, variable velocity responding to transient changes in surface mass balance, and/or changes in icesheet thickness. Herein, comparisons of the vertical profiles of computed ages, as well as annual layer thicknesses, were examined among the RCIP schemes (semiLagrangian) and upwind schemes (Eulerian). Although the experiments in the present study were limited to 1D 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 cosinewavetype 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 secondorder upwind scheme with a more suitable filter will change the characteristics. The introduction of slope limiters on general nonuniform discretization for higherorder 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 secondorder 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 firstorder 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 firstorder 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 firstorder 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 secondorder upwind scheme. Furthermore, as discussed above, the upstream correction significantly improves the RCIP solution, which suggests that it is important to consider the nonconstant velocity between the arrival and departure points. Since the midpoint rule formulation in the firstorder 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 firstorder approximation. In any case, it is clear that some or all of these points contribute to the higher performance of the firstorder 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 firstorder 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 secondorder scheme preserves the history better than the firstorder 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 secondorder scheme are still slightly better than those for the firstorder scheme throughout most of this study's experiments.
Greve et al. (2002) presented “practical suggestions” for numerical dating schemes: the secondorder, the total variation diminishing Lax–Friedrichs (TVDLF) scheme with the minimum modulus (minmod) filter, and even the firstorder 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 timerelated 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., millennialscale 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 squarewave shape pattern is well preserved, at least around the 1700 m with P_{T}=5 (a) even though the higherfrequency 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 1D age computations, implementation of the RCIP scheme for 3D computation of the age field is also a suitable subject for future discussions.
Extension to 3D would require the consideration of complex 3D 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 3D 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 nonzero 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 semiLagrangian schemes using higherorder 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 backtracing method. Huybrechts et al. (2007) suggested a very effective backtracing method, which can be sufficient by itself for icecore dating. The small but primary advantage of the RCIP method over the powerful backtracing 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 highprecision 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 icesheet modeling. The evolutions of icesheet 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 multidimensional conservative formulation of CIP schemes has already been proposed (Yabe et al., 2002). Accordingly, the implementation of the scheme to 3D age and temperature fields in numerical icesheet models has already been set as the next target of our development.
A timesplitting 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 nonadvection 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 semiLagrangian 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 nonadvection 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 floatingpoint 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 floatingpoint 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 roundingup 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 nonzero 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
IcIES2/JP version 0, which is a subset package of the
Ice sheet model for Integrated Earth system Studies II (IcIES2).
IcIES2 is available at
https://github.com/saitofuyuki/icies2 (last access: 26 September 2020, Saito and AbeOuchi, 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_gmd2020
in the git repository.
The supplement related to this article is available online at: https://doi.org/10.5194/gmd1358752020supplement.
FS developed the icesheet 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/S02773791(01)001111, 2002. a, b, c, d
Clarke, G. K. C., Lhomme, N., and Marshall, S. J.: Tracer transport in the Greenland ice sheet: threedimensional 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., DahlJensen, 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 ”OldestIce” ice core, Clim. Past, 9, 2489–2505, https://doi.org/10.5194/cp924892013, 2013. a, b
Goelzer, H., Nowicki, S., Payne, A., Larour, E., Seroussi, H., Lipscomb, W. H., Gregory, J., AbeOuchi, 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 sealevel contribution of the Greenland ice sheet: a multimodel ensemble study of ISMIP6, The Cryosphere, 14, 3071–3096, https://doi.org/10.5194/tc1430712020, 2020. a
Greve, R. and Hutter, K.: Polythermal threedimensional 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 nonclimatic 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/cp35772007, 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 ColdRegion 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 NAS05007, NAS, 2005. a, b
Parrenin, F., Dreyfus, G., Durand, G., Fujita, S., Gagliardini, O., Gillet, F., Jouzel, J., Kawamura, K., Lhomme, N., MassonDelmotte, V., Ritz, C., Schwander, J., Shoji, H., Uemura, R., Watanabe, O., and Yoshida, N.: 1Dice flow modelling at EPICA Dome C and Dome Fuji, East Antarctica, Clim. Past, 3, 243–259, https://doi.org/10.5194/cp32432007, 2007. a, b, c
Rybak, O. and Huybrechts, P.: A comparison of Eulerian and Lagrangian methods for dating in numerical icesheet models, Ann. Glaciol., 37, 150–158, https://doi.org/10.3189/172756403781815393, 2003. a, b, c, d
Saito, F. and AbeOuchi, A.: Ice sheet model for Integrated Earth system Studies II (IcIES2), GitHub, available at: https://github.com/saitofuyuki/icies2, last access: 26 September 2020. a
Saito, F., Obase, T., and AbeOuchi, A.: Resources relating to Saito et al. (2020, GMD) – IcIES2/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., AbeOuchi, A., Agosta, C., Albrecht, T., AsayDavis, X., Barthel, A., Calov, R., Cullather, R., Dumas, C., GaltonFenzi, 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 multimodel ensemble of the Antarctic ice sheet evolution over the 21st century, The Cryosphere, 14, 3033–3070, https://doi.org/10.5194/tc1430332020, 2020. a
Shashkov, M.: Conservative FiniteDifference 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.: SemiLagrangian Integration Schemes for Atmospheric Models: A Review, Mon. Weather Rev., 119, 2206–2223, https://doi.org/10.1175/15200493(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.: Multidimensional conservative semiLagrangian 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/00104655(95)001247, 1996. a, b, c
Yabe, T. and Takei, E.: A New HigherOrder 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 semiLagrangian solver for solid, liquid and gas, J. Comput. Appl. Math., 149, 267–277, https://doi.org/10.1016/S03770427(02)005356, 2002. a, b
Some models adopt 0 for the righthand 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